Magnetization processes of zigzag states on the honeycomb lattice:
Identifying spin models for RuCl and NaIrO
Abstract
We study the fieldinduced magnetization processes of extended HeisenbergKitaev models on the honeycomb lattice, taking into account offdiagonal and longerrange exchange interactions, using a combination of MonteCarlo simulations, classical energy minimization, and spinwave theory. We consider a number of different parameter sets, previously proposed to describe the magnetic behavior of RuCl and NaIrO with their antiferromagnetic zigzag ground states. By classifying these parameter sets, we reveal the existence of three distinct mechanisms to stabilize zigzag states, which differ in the sign of the nearestneighbor Kitaev interaction, the role of longerrange interactions, and the magnitude of the offdiagonal interaction. While experimentally hardly distinguishable at zero field, we find that the three different scenarios lead to significantly different magnetization processes in applied magnetic fields. In particular, we show that a sizable offdiagonal interaction naturally explains the strongly anisotropic field responses observed in RuCl without the need for a strong anisotropy in the effective tensor. Moreover, for a generic field direction, it leads to a highfield state with a finite transversal magnetization, which should be observable in RuCl.
I Introduction
Stronglycorrelated electron systems in the regime of strong spinorbit coupling have recently attracted a lot of attention, luring with the prospect of material realizations of nontrivial phases with possibly exotic excitations. Candidate systems that are presently intensely studied, both experimentally and theoretically, are based on transition metal ions with partially filled or shells, such as RuCl and IrO ().trebst2017 In these systems, the transitionmetal ions are caged in chlorine RuCl and oxygen IrO octahedra, respectively. These cages share an edge with each of their three neighboring octahedra, such that the central ions form weakly coupled layers of honeycomb lattices. The combined effect of level splitting due to the octahedral crystal field and spinorbit coupling generates moments. Additional Coulomb repulsion drives the system into a Mott insulating state already above room temperature, suggesting a localmoment description as a lowenergy model.plumb2014 ; choi2012 In the simplest, nondistorted cubic geometry, the edgesharing octrahedra allow two 90° RuClRu and IrOIr exchange paths, respectively. This leads to a partial destructive interference of the symmetric Heisenberg interaction in favor of a dominant bonddirectional Isingtype exchange of the form of the celebrated Kitaev compass model,jackeli2009 ; chaloupka2010 fueling the hope to realize Kitaev’s spin liquid phase of Majorana fermions.kitaev2006
Experimentally, RuCl and NaIrO show an ordering transition towards an antiferromagnetic state of collinear “zigzag” type,choi2012 ; sears2015 ; johnson2015 while LiIrO exhibits an incommensurate spiral magnetic ground state. williams2016 ; footnote1 In all cases, the Néel temperature is small compared to the CurieWeiss temperature, indicating substantial frustration. A suitable perturbation that allows to stabilize the spinliquid state experimentally is therefore much soughtafter.yadav2016 Interestingly, in RuCl a number of very recent measurements indeed appear to be consistent with a, possibly continuous,wolter2017 transition from the ordered phase towards a quantum paramagnetic phase, driven by magnetic fieldleahy2017 ; baek2017 ; sears2017 ; zheng2017 ; hentrich2017 ; banerjee2017 or external pressure.wang2017 ; cui2017
At zero field and ambient pressure, magnetic order must be caused by interactions beyond the nearestneighbor Kitaev model,chaloupka2010 and it has been pointed out that longerrangesingh2012 as well as offdiagonal interactionsrau2014 may play a role. An additional complication is that the materials appear to realize a lowtemperature monoclinic crystal structure instead of the idealized cubic one,choi2012 ; johnson2015 however, the actual size and role of the trigonal distortion is disputed.yadav2016 ; park2016 ; tjeng2017 Furthermore, when subject to an external magnetic field, RuCl, for instance, unveils a huge anisotropy between inplane and outofplane magnetic responses. johnson2015 ; majumder2015 ; kubota2015 In the spin models with only Kitaev and Heisenberg interactions, such large anisotropy can only arise from largely anisotropic factors.kubota2015 This, however, would require a substantial trigonal distortion, not easily consistent with crystallographic measurements.cao2016 ; park2016 ; tjeng2017 Evidently, a more natural explanation for the magnetic anisotropy would be an effective spin model which displays a strong intrinsic magnetic anisotropy, present even if the tensor were entirely isotropic. To the best of our knowledge, such has so far not been demonstrated.
To date, the debate about the most appropriate spin model for RuCl and NaIrO is by no means settled. There are basically (at least) three scenarios that were suggested to stabilize the experimentally observed zigzag order in these systems:
 Scenario 1

If the most dominant interactions are between nearest neighbors only, a zigzagordered state occurs if the Kitaev interaction is antiferromagnetic, .chaloupka2013
 Scenario 2

Zigzag order can be stabilized also for ferromagnetic Kitaev coupling, , if a sizable thirdneighbor antiferromagnetic Heisenberg interaction, , is present.fouet2001 ; winter2016
 Scenario 3
For RuCl, for instance, inelastic neutronscattering data have been argued to be consistent with antiferromagnetic Kitaev coupling (Scenario 1).banerjee2016a ; banerjee2016b ; gohlke2017 In contrast, abinitio studies typically find a ferromagnetic and have emphasized the significance of longerrange couplings (Scenario 2).winter2016 ; kim2016 ; yadav2016 Finally, some recent works, based on abinitio calculations or fits to neutronscattering data, suggest by contrast a nearestneighbor  model with ferromagnetic and strong positive as a minimal description of RuCl (Scenario 3).ran2017 ; wang2016 (Notably, we will demonstrate below that the plain  model is insufficient to stabilize zigzag order and needs to be supplemented with, e.g., a finite ferromagnetic Heisenberg coupling .) A similar  model with and strong , in this case supplemented with small Heisenberg interactions and (stabilizing zigzag order) has recently also been proposed.winter2017 Experimentally, it is difficult to differentiate between these scenarios at zero field. A possible indicator is the moment direction in the ordered state;sizyuk2016 ; chaloupka2016 however, this observable is difficult to access and requires the use of polarized neutronscao2016 or magnetic xray scattering data.chun2015 Clearly, further theoretical predictions that distinguish between these different mechanisms are necessary to resolve the debate.
In this paper, we demonstrate that the response of the system to an external field differs substantially for the different scenarios of stabilizing the zigzag state, and we argue that this can be used to narrow down the experimentally relevant parameter range of the models. Recently, we have mapped out the magnetization processes in the nearestneighbor HeisenbergKitaev model.janssen2016 Here, we contrast these findings to new results that are obtained within the extended HeisenbergKitaev models including longerrange and offdiagonal interactions. Given the huge parameter space, we mainly restrict our attention to those parameter regimes which lead to zigzag order in zero field, but some of the analytical considerations are more general, as will be noted below. Our methodology is concerned with the semiclassical limit of large spin , but we also compute quantum corrections to the highfield magnetization, which we evaluate for .
We devote a specific discussion to the effect of large in the presence of a magnetic field. We argue that this naturally leads to a strongly anisotropic field response, and we also discuss under which circumstances the term induces a finite transversal magnetization even in the highfield phase.
The remainder of the paper is organized as follows: In Sec. II, we introduce the general honeycomblattice spin Hamiltonian together with the different parameter sets considered in this paper. Sec. III quickly summarizes the different methods used to tackle the problem. Secs. IV, V, and VI discuss the influence of a magnetic field on the zigzag states stabilized by the three different mechanisms, corresponding to different parameter regimes of the general model. In all cases, we discuss degeneracies and the expected lowfield behavior using analytical arguments. We also show extensive numerical results for the magnetization processes for various concrete parameter sets. Sec. VII contains a specific discussion of the effects of a large offdiagonal interaction. A summary of our results, together with a discussion of their experimental relevance and suggestions for future experiments, closes the main part of the paper. In the appendix, we demonstrate that a nearestneighbor model with ferromagnetic and strong stabilizes zigzag order only when supplemented with a finite .
Ii Models
ii.1 Generic Hamiltonian
Out of the large variety of spin Hamiltonians proposed for RuCl and NaIrO, most can be written as extensions of the HeisenbergKitaev model originally proposed in Ref. chaloupka2010, . Consequently, we consider models on the honeycomb lattice of the general form
(1) 
Here, correspond to first, second, and thirdneighbor Heisenberg couplings, similarly are Kitaev couplings, and is a symmetric offdiagonal coupling. and denote first and secondneighbor bonds, respectively, with . On bonds , with cyclic permutation for and bonds. The thirdneighbor bonds , and neglect possible trigonal distortions. are along opposite points of the same hexagon. We assume cubic symmetry, corresponding to a perfect honeycomb structure in
We are interested in the effects of a uniform external magnetic field, which couples linearly to the spins via
(2) 
where the constants have been absorbed into the field , with the (possibly anisotropic) effective tensor and the Bohr magneton. The total Hamiltonian then is .
ii.2 Conventions for field directions
Due to the broken spin rotation symmetry, and as will be demonstrated explicitly below, the response of the system to an external field crucially depends on the field direction . For convenience, we will give all field directions in the cubic spin basis and label them in the form , i.e.,
(3) 
In RuCl, the cubic axes , , and point along nearestneighbor RuCl bonds, see Fig. 1(b). Consequently, the two inplane field directions along the crystallographic and axes correspond to and , respectively. The direction perpendicular to the honeycomb layer, which is sometimes referred to as axis,johnson2015 is labelled as the direction in our conventions. In the same way, denotes an intermediate direction which lies in the plane and is tilted 55° away from the axis towards the axis.
ii.3 Parameter sets
Set  Material  [meV]  [meV]  [meV]  [meV]  [meV]  [meV]  Method  Ref.  Year 

1  RuCl  —  —  —  —  fit to neutron scattering  banerjee2016a, , banerjee2016b,  2016  
1’  NaIrO  —  —  —  —  fit to susceptibility & neutron scattering  chaloupka2013,  2013  
1+  RuCl  —  —  —  DFT + expansion  kim2015,  2015  
2  NaIrO  —  DFT + exact diagonalization  winter2016,  2016  
2+  NaIrO  DFT + expansion, direction of moments  sizyuk2016, , sizyuk2014,  2016  
(2+)’  NaIrO  —  MRCI, fit to  katukuri2014,  2014  
(2+)”  RuCl  —  MRCI, fit to magnetization  yadav2016,  2016  
2/3  RuCl  —  DFT + exact diagonalization  winter2016,  2016  
3  RuCl  —  —  —  —  fit to neutron scattering  ran2017,  2017  
3’  RuCl  —  —  —  —  DFT + expansion  wang2016,  2016  
3”  RuCl  —  —  —  DFT + expansion  kim2016,  2016  
3+  RuCl  —  —  fit to neutron scattering  winter2017,  2017 
In Table 1, we list popular parameter sets for the couplings , , and , that were suggested either on the basis of abinitio calculations or by fitting the predictions of different simplified versions of the above model to experimental data. The parameter sets can roughly be divided into three groups, corresponding to the three scenarios listed in the introduction:

Dominant antiferromagnetic , supplemented by ferromagnetic and small longerranged interactions,

Dominant ferromagnetic together with large antiferromagnetic ,

Strong in conjunction with ferromagnetic and only small Heisenberg couplings.
The classical energy contributions to a zerofield zigzag state within three different minimal models, representative for the three scenarios, are illustrated in Fig. 1. The observed zigzag states in RuCl (Refs. sears2015, ; johnson2015, ) and NaIrO (Ref. choi2012, ) are in principle compatible with all three groups. We reiterate, though, that the plain  model suggested in Refs. ran2017, ; wang2016, needs to be supplemented with, e.g., a finite ferromagnetic Heisenberg coupling in order to stabilize zigzag order, see Sec. VI and the appendix for details. Parameter Sets 3 and 3’ therefore do not lead to a zigzag ground state, at least on the classical level.
In the absence of an external magnetic field, it appears difficult to experimentally rule out any of the various (mutually incompatible) parameter sets. In what follows, we argue that the behavior in an external field qualitatively differs for the three scenarios. In particular, we propose the magnetic field response as an experimental indicator of the sign of and demonstrate that the value of can be obtained (i) by measuring the angle dependence of the magnetic susceptibility or critical field strength, or (ii) by applying a strong magnetic field provided that the finite transversal magnetization is accessible in the experimental setup.
Iii Methods
iii.1 MonteCarlo simulations
To obtain concrete results for the magnetization processes of different microscopic models, we have employed classical MonteCarlo (MC) simulations, combining singlesite and paralleltempering updates in order to equilibrate the spin configurations at low ; for details see Ref. janssen2016, . Here, we have performed field scans on lattices with up to sites at temperatures down to and have computed the static spin structure factor , where is the lattice vector at site . The lowest configurations have been further cooled down to obtain the zerotemperature magnetization in the respective classical ground state.
iii.2 Variational approach
The MC simulations are augmented with a semianalytical variational calculation in which we have assumed a foursite magnetic unit cell with two pairs of pairwise parallel spins according to the three proposed , , and zigzag patterns (Fig. 1). The classical energy has been minimized with respect to the four independent spherical angles and of the two pairs of spins for each zigzag pattern. We have crosschecked the variational results with our MC data points. For all cases studied the ground state is indeed either a (uniformly or nonuniformly) canted zigzag state or a homogeneous highfield state, and we find perfect agreement within a relative uncertainty of typically (except for a few data points that are very close to firstorder transitions). This way, we obtain the full classical magnetization curve as function of continuous for various field directions including (corresponding to the axis in RuCl), ( axis), ( axis), and (intermediate direction in the plane).
iii.3 Spinwave theory
As another crosscheck, and in order to obtain an estimate on the influence of quantum fluctuations in the respective systems, we employ spinwave theory in the homogeneous highfield phase. For large , the magnon spectrum is gapped. Upon lowering the field strength, the gap decreases. If the transition towards the canted zigzag phase at is continuous, it can be understood as condensation of magnons, and the critical field strength at which the gap closes can be computed analytically.janssen2016 In linear spinwave theory, it precisely agrees with the classical continuous transition point, and we have crosschecked our spinwave results against the MC simulations and the variational technique. In the cases of continuous transitions, we have checked that the instability wave vector, at which the magnon gap closes, agrees with the ordering wave vector of the canted zigzag state, as it should be as long as there are no intermediate phases with other ordering wave vectors.janssen2016 Furthermore, we have computed the leadingorder correction to the classical magnetization in the highfield phase for and for , the sizes of which allow us to assess the validity of our classical approximation. Details on the spinwave calculations are given in the Supplemental Materials to Refs. janssen2016, ; wolter2017, .
Iv Scenario 1: Zigzag order from antiferromagnetic and ferromagnetic
In this and the following two sections, we discuss and contrast the different stabilization mechanisms of zigzag order corresponding to the three scenarios, and map out the consequential different magnetization processes. We remind the reader that zigzag order is characterized by ferromagnetic chains along zigzag lines and antiferromagnetic bonds between neighboring zigzag chains, see Fig. 1(a) for a case with antiferromagnetic bonds (dubbed “zigzag”). Symmetryequivalent zigzag states (“zigzag” and “zigzag”) are obtained by a rotation about one site in real space in conjunction with a rotation about the axis in spin space;footnote4 in addition there is a trivial global spinflip (Ising) symmetry. At zero external field, any zigzag state is therefore at least sixfold degenerate.
Let us start with the case of the nearestneighbor HeisenbergKitaev model with , as originally proposed in Ref. chaloupka2013, . In this scenario, the zigzag chains are stabilized by a ferromagnetic Heisenberg interaction, while the antiferromagnetic ordering between neighboring zigzag chains is induced by an antiferromagnetic Kitaev interaction. In the situation with antiferromagnetically ordered bonds and ferromagnetic chains along and bonds [“zigzag state”, Fig. 1(a)], the Kitaev coupling is satisfied only when the spins are aligned along the axis in spin space, . This is illustrated in Figs. 1(c) and (d). Equivalently, for the  ()zigzag states, the spins point along (). In Scenario 1, there are therefore six possible zigzag ground states at zero field. We denote them as “cubicaxis zigzag.”
In general, when a collinear antiferromagnet is brought into a small external magnetic field, the fielddependent part of the energy is minimized when the staggered order is perpendicular to . In a Heisenberg system with spin rotation invariance, such alignment is always possible already at infinitesimal field strength. For increasing field strength, all spins then cant homogeneously, i.e., with a common “canting angle” towards the magnetic field axis, until the system reaches saturation at some critical field strength . In the classical limit the magnetization curve is linear below .
This situation drastically changes when the spin symmetry is broken. In the present Scenario 1, with the possible zerofield spin alignment along the cubic axes only, it is impossible for the spins to align perpendicular to a small external field, unless the magnetic field is orthogonal to one of the three cubic axes , , and . At finite field strength, the spins therefore generically cant inhomogeneously towards the magnetic field with two different canting angles, and the classical magnetization curve becomes nonlinear. Furthermore, the nonuniformly canted zigzag states then compete with other states that allow for a uniform canting mechanism. This happens, for example, for a field in the outofplane direction. As we have demonstrated recently, the magnetization process for undergoes a sequence of metamagnetic transitions, and eventually (if ) exhibits a continuous transition from a vortexlike state with uniform canting angle towards the polarized state. On the other hand, for , there is a direct firstorder transition from the nonuniformly canted zigzag state towards the polarized state.janssen2016
By contrast, for an external field orthogonal to one of the cubic axes, e.g., , the spins can align perpendicular to an infinitesimal . The spins therefore cant towards with a uniform canting angle and a linear classical magnetization curve, until eventually a continuous transition towards the polarized phase occurs.
While the nature of the magnetization process itself is highly anisotropic in this scenario, both the magnitude of the magnetic susceptibility at small fields and the critical field strength (at which the transition to the highfield state occurs) only slightly depend on the field direction. Classically, one finds, for instance,
(4) 
as long as the highfield transition is continuous. Note that is even larger for the inplane direction than for the outofplane direction [in contrast to the anisotropy in the field response observed in RuCl (Refs. johnson2015, ; majumder2015, )].
Our numerical simulations for the prototype Parameter Set 1 are fully consistent with these expectations, see Fig. 2. For [Fig. 2(b)] and [Fig. 2(c)] the spins cant homogeneously, i.e., with a uniform canting angle, towards the magnetic field axis, resulting in a linear classical magnetization curve in the intermediatefield regime. At the critical field strength , there is a continuous transition towards the polarized highfield phase with uniform . For [Fig. 2(a)] and [Fig. 2(d)], for which no spin alignment perpendicular to at small field is possible, the canting becomes inhomogenous and the magnetization curve becomes nonlinear with a positive curvature in the intermediatefield regime. For this particular Parameter Set 1 with , and these field directions, there are no metamagnetic transitions, and, instead, there is a direct firstorder transition towards the polarized phase. Magnetization curves for other parameter sets with , showing various novel largeunitcell and vortexlike phases at intermediate field strengths are displayed in Ref. janssen2016, . In the classical limit, , the highfield phase is fully polarized, , while quantum fluctuations reduce the magnetization below full saturation, , for .
V Scenario 2: Zigzag order from ferromagnetic and antiferromagnetic
The second possible mechanism, suggested to stabilize the zigzag state in both RuCl and NaIrO, is characterized by a ferromagnetic Kitaev coupling in conjunction with a thirdneighbor antiferromagnetic Heisenberg interaction between opposite sites of the same hexagon. In this case, the ferromagnetic zigzag chains (along, e.g., and bonds) are stabilized by the Kitaev coupling, while the antiferromagnetic ordering between neighboring zigzag chains is stabilized by . This is illustrated in Figs. 1(e) and (f). Classically, the spins in the zigzag state can therefore point anywhere in the plane, and analogously for the  and zigzag states. Consequently, the system has an accidental continuous classical groundstate degeneracy of . We denote these states as “cubicplane zigzag.” In the real systems, the degeneracy is partially lifted by three generically competing effects: (i) Quantum fluctuations drive an orderfromdisorder mechanism which favors a spin alignment along the cubicaxes direction, i.e., or for the zigzag state and analogously for the  and zigzag state.baskaran2008 ; sizyuk2016 ; rousochatzakis2017b (ii) A small offdiagonal interaction () lifts the degeneracy, e.g., in the zigzag state in favor of a state in which ().chaloupka2015 ; sizyuk2016 (iii) A small external magnetic field favors a state in which the spins are aligned perpendicular to .
For conceptual clarity, let us first discuss the magnetization process within an idealized classical  model (using, e.g., Parameter Set 2). For any field direction, the plane perpendicular to intersects all three cubic planes , , and in (at least) one axis. The zigzag spins can therefore always align perpendicular to a small external field. They cant uniformly towards the magnetic field and exhibit in the classical limit a direct continuous transition towards the polarized state in the minimal  model at
(5) 
independent of . Furthermore, there is a residual groundstate degeneracy at finite field for all field directions. This is in perfect agreement with our numerical results for Parameter Set 2, Fig. 3. Note that the presence of further Heisenberg or Kitaev interactions, such as , , and modify Eq. (5), but do not spoil the isotropy of the magnetization process in the classical limit and when .
Let us now include the effects of perturbations away from this idealized limit. Both quantum fluctuations and offdiagonal interactions lift the continuous groundstate degeneracy in favor of one or two preferred axes, which, as discussed above, differ among the two cases. At zero field, the spins align along the axis that corresponds to the stronger one of the two perturbations.sizyuk2016 For a generic field direction that is not perpendicular to this axis, the zigzag spins will then cant nonuniformly towards the field axis for small fields , where is the characteristic energy scale of the orderfromdisorder effects. If both orderfromdisorder effects and offdiagonal interactions are weak (e.g., for large or large , and small ), there will be a transition (or crossover) at some finite field strength towards the canted version of the cubicplane zigzag state with uniform canting angle. We confirm this expectation numerically for the Parameter Set 2+ in Fig. 4. At zero field, the spins align along the , , or direction, with . Upon switching on a small magnetic field in, e.g., the or direction, they cannot align perpendicular to , and therefore cant nonuniformly towards it, leading to a nonlinear magnetization curve. At a critical field strength , however, Figs. 4(a) and (b) show a transition towards an intermediate phase with a uniform canting and a linear magnetization curve. For this parameter set, there is no true transition for the two inplane directions displayed in Figs. 4(c) and (d), although the latter case can be understood as a (very weak) crossover.
Hence, the zigzag state that is stabilized by a ferromagnetic and antiferromagnetic , in the presence of only weak , shows a metamagnetic transition at intermediate field strength for generic (but not all) field directions. The transition towards the polarized state, by contrast, occurs at a field strength that is, up to small corrections due to quantum fluctuations and perturbations from the offdiagonal interaction, independent of the field axis. For stronger perturbations as, for instance, in the Parameter Set 2/3 with sizable , the “wouldbe” transition towards the canted cubicplane zigzag state may be preempted by a firstorder transition towards the polarized state, see Fig. 5. Consequently, the transition towards the highfield phase at is discontinuous for [Fig. 5(a)] and (b), but remains continuous for (c) and (d).
Vi Scenario 3: Zigzag order from positive and ferromagnetic ,
Finally, we consider the zigzag state that is stabilized within a nearestneighbor model with strong positive offdiagonal interaction . We emphasize that the classical ground state of the plain  model with and is a multi state with six major and three minor Bragg peaks in the first Brillouin zone. The zigzag ground state can, however, be recovered for nearestneighbor interactions only, when the model is supplemented with a ferromagnetic Heisenberg interaction , provided that the ratio is large enough. This is demonstrated in detail in the appendix. The minimal description of Scenario 3 is therefore given by a  model. In this model, the spins in, e.g., the zigzag state are aligned for zero external field along an direction determined by , with
(6) 
and thus . Note that the value of the Heisenberg interaction does not affect the direction of the ordered moments at zero field. The same applies to further isotropic interactions (, ), such that Eq. (6) is valid for all parameter sets in Table 1 with , , and . For small , we have , such that . For large , on the other hand, one finds , and hence . In this latter case, the effective moments in RuCl and NaIrO would point towards the centers of opposite faces of the surrounding chlorine or oxygen octahedra. This is illustrated in Figs. 1(g) and (h). To stay within a sufficiently simple naming scheme, we therefore denote these states, in a slight abuse of notation for finite , as “facecenter zigzag.” In the  and zigzag states, the spins are analogously aligned along an and direction. The zerofield ground state in this scenario is therefore sixfold degenerate.
Due to the alignment of the zigzag spins at zero field, the magnetization process crucially depends on the direction of the external field. If the magnetic field is perpendicular to the zerofield orderedmoment direction, a simple canting mechanism with uniform canting angle and, in the classical limit, a linear magnetization curve can occur. This happens, for example, for the inplane direction . At a critical field strength with (in the simplest  model)
(7) 
there is a direct and continuous transition towards the polarized state. There is another linearly independent direction of the form , perpendicular to the zerofield moments in the zigzag state. A uniform canting and a continuous highfield transition are expected in this field direction as well. Analogous perpendicular directions and exist for the  and zigzag states. All these directions, however, depend on the ratio via the function .
For other field directions, however, the zigzag spins cannot align perpendicular to the field axis for small . Consequently, they cant with two different angles towards the magnetic field. We demonstrate numerically in Figs. 6(a), (b), and (d) that the magnetization curve is nonlinear and the transition to the highfield state for these field directions is typically first order. By contrast, the magnetization curve is linear and the highfield transition is continuous for [Figs. 6(c)]. Fig. 6 also shows that the critical field strength , at which the transition to the highfield phase occurs, now critically depends on the field axis, with the largest among the four shown here for the outofplane direction [Fig. 6(a)]. This is in sharp contrast to the situation with small shown in Figs. 2–4. Furthermore, the highfield phase for shown in Fig. 6(b) for is characterized by a nonsaturated classical magnetization in field direction, corresponding to a finite transversal magnetization perpendicular to . We elaborate on these key aspects of the models with strong in the following section.
Vii Consequences of strong
The presence of a sizable offdiagonal interaction as suggested within Scenario 3 has two characteristic consequences, which are in sharp contrast to the effects of the Heisenberg and Kitaev interactions and which allow to test this scenario in future experiments.
vii.1 Anisotropic magnetic response
First, we note that identifying the sign of as “ferromagnetic” or “antiferromagnetic” is inappropriate. The nature of the offdiagonal interaction for fixed sign of in fact depends on the direction of the ordered moments. Let us assume for simplicity a uniform spin state in an arbitrary direction, as it would occur for a very strong magnetic field. The argument, however, can be made in a similar way also for components of the spins in the direction of a small external field. We find it convenient to choose the new orthonormal basis in which the coordinate axes are aligned along the crystallographic , , and axes in Fig. 1(b), i.e., , , and . The uniform state can then be written as
(8) 
where and are the spherical angles in the basis, such that corresponds to the axes and corresponds to the plane. In this basis, the fieldindependent part of the total classical energy becomes particularly simple,
(9) 
Notably, the energy becomes independent of the azimuthal angle , and only the offdiagonal interaction term depends on the polar angle . For any inplane direction, this anisotropic part becomes , while for the outofplane direction . In general, the dependent contribution is negative for and positive otherwise (assuming that ). Within the honeycomb plane, a positive therefore favors a ferromagnetic spin alignment over an antiferromagnetic one. If the spins are forced to lie perpendicular to the honeycomb plane, on the other hand, a positive favors an antiferromagnetic spin pattern over a ferromagnetic alignment.
Consequently, the response of a system with a large offdiagonal interaction to an external magnetic field crucially depends on the field axis: For an outofplane direction, the offdiagonal interaction has an antiferromagnetic character. Increasing lifts the energy of the polarized state, leads to an increase of the critical field strength at which the transition to the highfield state occurs, and makes the system less susceptible to the external field. For an inplane direction, by contrast, a positive behaves as if the interaction were ferromagnetic. Increasing reduces the critical field strength, and makes the system more susceptible to the external field. We illustrate this behavior for the lowtemperature susceptibility at small inplane and outofplane fields as well as the critical field strength in Fig. 7. Here, we have used Parameter Set 3+ as an example, and have solved the finitefield problem in the classical limit by employing the methods described in Sec. III. Without any factor anisotropy, the ratio between inplane and outofplane susceptibility for this parameter set would be , and the ratio of critical field strengths would be . Between different inplane directions, on the other hand, the anisotropy is much smaller, e.g., for the ratio of the critical field strengths. Note that a negative has precisely the opposite effect: It leads to a ferromagnetic character of the offdiagonal interaction for an outofplane magnetic field and an antiferromagnet character for an inplane field . Hence, a positive offdiagonal interaction naturally explains the strong magnetic anisotropy observed experimentally in RuCl.footnote6 This is further discussed in Sec. VIII.
vii.2 Highfield transversal magnetization
Second, we investigate the nature of the highfield phase in the presence of a finite . To this end, it is instructive to recall the symmetries of our (idealized cubic) system.jiang2011 Due to the bonddependent interactions, reflecting strong spinorbit coupling, the continuous spin rotation symmetry known from isotropic Heisenberg models is replaced by a residual discrete symmetry of rotation around the axis. This corresponds to the cyclic permutation in spin space together with a rotation about one site in real space. There is also a discrete rotation symmetry about the inplane axis with angle , which corresponds to in spin space and a reflection in real space. Furthermore, there is a reflection symmetry through the plane perpendicular to the axis, corresponding to the spin transformation and a reflection in real space. The composition maps a vector in spin space onto its inverse, , which corresponds to the global spinflip (Ising) symmetry mentioned earlier. However, if the offdiagonal interaction vanishes, , the model has an enhanced symmetry which is not protected by the lattice point group: It is given by the individual inversion of components of , i.e., , , or .
An external magnetic field breaks part of the symmetry. For instance, a field in the direction breaks and , but leaves invariant. If , inverting the and components of individually remains a symmetry, but not for . At high fields beyond a certain critical field strength , one expects a ground state that is adiabatically connected the fully polarized state at , without spontaneous symmetry breaking. Consequently, without the offdiagonal interaction, the spin expectation value in the highfield state is locked to the field axis,
(10) 
As a side note, we recall that even in this case the state at is not fully polarized, , since quantum fluctuations reduce the magnetization as a consequence of the broken spin symmetry, possibly even substantially so, see the red/light curves in Figs. 2–6, as well as Ref. janssen2016, .
In the presence of a finite offdiagonal interaction , however, there is no symmetry which forbids a highfield state that is not parallel to , such as, e.g., with for . If this indeed happened, it would mean that the magnetization would exhibit a transversal component which is perpendicular to the external field for all finite , and only in the strict limit would the total magnetization become parallel to the external field . In order to examine for which parameters this indeed happens in our model, and if so, for which field directions, we again consider a uniform classical highfield state with the fieldindependent part of the total energy given in Eq. (VII.1). We first note that becomes entirely independent of and when . This is consistent with our symmetryreasoned expectation formulated above, stating that is then parallel to the magnetic field for all field directions. With the offdiagonal interaction present, however, whether or not and are collinear depends on the magneticfield direction and the ratio . In any case, since is independent of , will always lie in the plane spanned by and . Let us denote the polar angle of in the basis by . The angle between and is then given by . To see whether has a finite expectation value in the highfield state, we have to minimize the total energy , where
(11) 
This yields a conditional equation for the energyminimizing as function of , reading
(12) 
and which can be solved exactly for any given . From this we find for or and large enough . For , by contrast, is finite for all . Consequently, we find a finite transversal magnetization in the highfield state for all field directions which are neither parallel nor perpendicular to the axis. In the classical limit, the relative size of the transversal magnetization is given by
(13) 
with determined through Eq. (12). These results are consistent with the numerical data presented in Figs. 2–6: The uniform highfield phase is characterized by a finite transversal magnetization (indicated by a nonsaturated classical magnetization in field direction) if and only if and for an external field that is neither parallel nor perpendicular to the honeycomb plane. We have explicitly checked that our numerical results agree with Eq. (13) also quantitatively. For instance, at in the direction (i.e., ), we find and for the Parameter Sets 2/3 and 3+, respectively.
Viii Discussion
viii.1 Susceptibility anisotropy
RuCl (Ref. kubota2015, ) and, to a somewhat lesser extent, NaIrO (Ref. singh2010, ) show a distinguished anisotropy in the magnetic response. The inplane magnetic susceptibility in RuCl at low fields and temperatures is about four times larger than its outofplane counterpart .kubota2015 Similarly, at a field strength of, e.g., 15 T, the magnetization is about 5–6 times larger when the external field is applied within the plane as compared to when .johnson2015 Given that the trigonal distortion is expected to be relatively small,cao2016 ; park2016 ; tjeng2017 it appears quite unnatural to attribute such large anisotropy exclusively to the effects of an anisotropic tensor. In fact, abinitio calculations rather point to a fairly isotropic tensor,majumder2015 ; tjeng2017 suggesting an intrinsic mechanism originating from bonddependent spinspin exchange interactions.
In this work, we have demonstrated that a model with a strong offdiagonal interaction naturally leads to a strongly anisotropic field response, without the need to assume a largely anisotropic tensor. For the particular Parameter Set 3+, the ratio of susceptibilities between an inplane field direction and the outofplane direction becomes of the order of (depending on the particular inplane direction chosen), which agrees with the experimental value for RuCl. We therefore suggest that the offdiagonal interaction in RuCl is positive and large, possibly of the order of , as given by Parameter Set 3+. As an aside, we note that the direction of the ordered moments at zero field as given by Eq. (6) for follows as , which corresponds to an axis in the plane that is tilted away from the axis towards the axis. Remarkably, this appears to be in fair agreement with one of the two possible options expected for RuCl.cao2016
The anisotropy within the honeycomb plane is much smaller, but finite, with a characteristic periodicity resulting from the discrete lattice rotation symmetry. We predict that the inplane susceptibility in both RuCl and NaIrO is maximal for a field direction along RuRu bonds and IrIr bonds, respectively (i.e., for being parallel to the crystallographic axis in our conventions), and minimal for the inplane direction that is perpendicular to this axis ( axis in our conventions).
viii.2 Fieldinduced transitions
Recent experiments report a fieldinduced transition in RuCl from a canted ordered phase into a putative quantum disordered phase when a field of around 7 T (10 T) is applied along an inplaneleahy2017 ; sears2017 ; zheng2017 ; hentrich2017 ; banerjee2017 (intermediatebaek2017 ) direction. While at least one, wolter2017 or more, banerjee2017 fieldinduced phase transitions appear to be observable independent of the specific sample, the nature of the adjacent phases has so far not been established beyond doubt.
We have argued that the magnetization process differentiates between the three scenarios for stabilizing a zigzag state. While the transition towards the uniform highfield phase is generically continuous when the field is directed along a RuRu bond (e.g., along the axis), this is different when the field is applied along an inplane direction orthogonal to a RuRu bond (e.g., axis): Here, the highfield transition is continuous if and only if is ferromagnetic and is not too large, with the Parameter Set 2/3 appearing to describe a rough upper bound for the value of the offdiagonal interaction. For larger values of , such as in the Parameter Set 3+, the transition becomes discontinuous. When the actual nature of the Kitaev interaction in RuCl is antiferromagnetic, as suggested within Scenario 1, the transition from the canted zigzag state into the highfield state becomes continuous only if the magnetic field lies within one of the three cubic planes, which in RuCl are the planes perpendicular to the RuCl bonds. Detailed magnetization measurements, using high (pulsed) magnetic fields, for different field directions are clearly called for.
Our results feature phase transitions at intermediate fields only for Parameter Set 2+, there between different canted states, and possibly in Scenario 1, there towards more complicated vortex and multi states.janssen2016 This situation might change upon accounting for (i) quantum effects beyond the semiclassical limit (possibly leading to fieldinduced spinliquid behavior) and (ii) deviations from the idealized cubic structure. Investigations in these directions are left for future work.
Recent neutron diffraction measurements in RuCl in small magnetic field applied along an inplane direction that is orthogonal to a RuRu bond (e.g., along the axis in our notation), find a depopulation of zigzag domains with ordering wavevectors .sears2017 ; banerjee2017 This observation is consistent with our calculations in the models with (Scenario 3) or (Scenario 1): A magnetic field along the inplane direction, which is perpendicular to a bond in our conventions, favors the  and zigzag states to the detriment of the zigzag state.
viii.3 Transversal magnetization at high fields
We have furthermore shown that the offdiagonal interaction leads to a finite transversal magnetization in the uniform highfield phase if the magnetic field axis is neither parallel nor perpendicular to the honeycomb plane. As the offdiagonal interaction already breaks part of the symmetry of the original HeisenbergKitaev model explicitly, this is possible without spontaneous symmetry breaking. For a magnetic field along RuCl bonds, for instance, the transversal magnetization becomes comparatively large and should in principle be observable in RuCl. Its magnitude is a direct measure of the size of , see Eqs. (12) and (13).
Ix Conclusions and outlook
In this paper, we have studied the magneticfield behavior of various extended HeisenbergKitaev models on the honeycomb lattice. The magnetization processes sensitively depend on the signs and relative sizes of the model parameters, and we have assessed the validity of the models in light of the experimental observations in RuCl and NaIrO. We have also predicted a number of further measurable consequences.
With these theoretical predictions, we believe that a detailed experimental characterization of the magnetization processes of the zigzag states should allow to fix the sign of the Kitaev interaction and determine the magnitude of the offdiagonal interaction in RuCl and NaIrO. Once an appropriate effective spin Hamiltonian has been established, a full quantum computation that goes beyond the present semiclassical analysis is indispensable. This should also allow to understand and characterize the fieldinduced phases found experimentally in RuCl.leahy2017 ; baek2017 ; sears2017 ; zheng2017 ; hentrich2017 ; banerjee2017
Acknowledgements.
We thank B. Büchner, P. LampenKelley, R. Moessner, S. Nagler, N. Perkins, S. Rachel, R. Valentí, J. van den Brink, A. Vishwanath, Z. Wang, S. Winter, and A. Wolter for illuminating discussions and collaboration on related work. This research was supported by the DFG through SFB 1143 and GRK 1621. ECA was supported by FAPESP (Brazil) Grant No. 2013/006818 and CNPq (Brazil) Grant No. 302065/20164.Appendix A Zigzag state in nearestneigbor  model
Recently, a minimal nearestneighbor  model with and was suggested to describe the physics of RuCl.ran2017 ; wang2016 In this appendix, we demonstrate that such an oversimplified model—without any Heisenbergtype interactions—is (at least in the classical limit) not sufficient to stabilize a zigzag antiferromagnet when and . A finite ferromagnetic Heisenberg interaction and a large enough ratio are needed to recover the zigzag ground state. This result is in contrast to the exploratory study of Ref. rau2014, , which used a combined LuttingerTisza and single analysis for the classical energy minimization. Here, we employ classical MonteCarlo simulations which are capable to also detect multi and incommensurate phases in a sufficiently unbiased way.
Fig. 8(a) shows the static structure factor from lowtemperature MC simulations for the Parameter Set 3 at zero external field. The ground state is characterized by six major Bragg peaks at , , as well as three minor peaks . We have checked that this remains true also upon increasing the ratio . By contrast, a zigzag antiferromagnet would lead to only three Bragg peaks at the three points of the first Brillouin zone. As a further check, we have explicitly verified that the multi state as found in the MC simulations has indeed a lower energy as the zigzag state when both are cooled down to , as suggested by the low structure factor. Consequently, the classical ground state of the  model with