Phase Diagram and Quantum Order by Disorder in the Kitaev $K_1$-$K_2$ Honeycomb Magnet

We show that the topological Kitaev spin liquid on the honeycomb lattice is extremely fragile against the second-neighbor Kitaev coupling $K_2$, which has recently been shown to be the dominant perturbation away from the nearest-neighbor model in iridate Na$_2$IrO$_3$, and may also play a role in $\alpha$-RuCl$_3$ and Li$_2$IrO$_3$. This coupling naturally explains the zigzag ordering (without introducing unrealistically large longer-range Heisenberg exchange terms) and the special entanglement between real and spin space observed recently in Na$_2$IrO$_3$. Moreover, the minimal $K_1$-$K_2$ model that we present here holds the unique property that the classical and quantum phase diagrams and their respective order-by-disorder mechanisms are qualitatively different due to the fundamentally different symmetries of the classical and quantum counterparts.


I. Introduction
The search for novel quantum states of matter arising from the interplay of strong electronic correlations, spin-orbit coupling (SOC), and crystal field splitting has recently gained strong impetus in the context of 4d and 5d transition metal oxides [1]. The layered iridates of the A 2 IrO 3 (A=Na,Li) family [2][3][4][5][6][7] have been at the center of this search because of the prediction [8,9] that the dominant interactions in these magnets constitute the celebrated Kitaev model on the honeycomb lattice, one of the few exactly solvable models hosting gapped and gapless quantum spin liquids (QSLs) [10]. This aspect together with the realization that the Kitaev spin liquid is stable with respect to moderate Heisenberg-like perturbations [9,11] has triggered a lot of experimental activity on A 2 IrO 3 and, more recently, on the similar α-RuCl 3 compound [12][13][14].
In the layered A 2 IrO 3 magnets, the single-ion ground state configuration of Ir 4+ is an effective pseudospin J eff = 1/2 doublet, where spin and orbital angular momenta are intertwined due to the strong SOC. In the original Kitaev-Heisenberg model proposed by Jackeli and Khaliullin [8], the pseudospins couple via two competing nearest neighbor (NN) interactions: An isotropic antiferromagnetic (AFM) Heisenberg exchange, J 1 , and a highly anisotropic Kitaev interaction, K 1 , which is strong and ferromagnetic, a fact that is also confirmed by ab-initio quantum chemistry calculations by Katukuri et al [15,16]. Nevertheless, neither Na 2 IrO 3 nor Li 2 IrO 3 are found to be in the spin liquid state at low temperatures. Instead, they show, respectively, AFM zigzag and incommensurate long-range magnetic orders, none of which is actually present in the Kitaev-Heisenberg model for FM K 1 coupling.
The most natural way to obtain these magnetic states is by including further neighbor Heisenberg couplings [15][16][17][18], which are non-negligible due to extended nature of the 5dorbitals of Ir 4+ ions [6,19]. In addition, recent calculations by Sizyuk et al [20] based on the ab-initio density-functional data of Foyevtsova et al [21] have shown that, for Na 2 IrO 3 , the next nearest neighbor (NNN) exchange paths must also give rise to an anisotropic, Kitaev-like coupling K 2 , which turns out to be AFM. More importantly, this coupling is the largest interaction after K 1 . It has also been argued [22] that K 2 plays an important role in the stabilization of the IC spiral state in Li 2 IrO 3 and might be deduced from the strong-coupling limit of Hubbard model with topological band structure [23,24].
Recent structural [12] and magnetic [13] studies have shown that the layered honeycomb magnet α-RuCl 3 is another example of a strong SOC Mott insulator, where the Ru 3+ ions are again described by effective J eff = 1/2 doublets. At low T , this magnet exhibits zigzag ordering as in Na 2 IrO 3 . Furthermore, the superexchange derivations [25,26] based on the ab initio tight-binding parameters show that the NNN coupling K 2 is again appreciable, and the signs of both K 1 and K 2 are reversed compared to Na 2 IrO 3 (i.e., K 1 is AFM and K 2 is FM). However, a strong off-diagonal symmetric NN exchange Γ term [15,16,27], which is allowed by symmetry, is also present [25,26], together with a much smaller J 1 cou- The model we consider here is described by the effective spin-1/2 Hamiltonian where ij (respectively ij ) label NN (NNN) spins on the honeycomb lattice, S a j defines the ath cartesian component of the spin operator at site j, and γ ij (λ ij ) define the type of Ising coupling for the bond (ij), see Fig. 1. This model interpolates between two well known limits, the exactly solvable Kitaev spin liquid [10] at K 2 = 0, and the triangular Kitaev model at K 1 = 0 [28][29][30][31][32]. It is easy to see that a finite K 2 ruins the exact solvability of the NN Kitaev model because the flux operators [10] W p = 2 6 S z 1 S x 2 S y 3 S z 4 S x 5 S y 6 (see site-labeling convention in Fig. 5, top left), around hexagons p are no longer conserved.
In the following we parametrize K 1 = cos ψ and K 2 = sin ψ, and take ψ ∈ [0, 2π). It turns out that the physics actually remains the same under a simultaneous sign change of K 1 and K 2 , because this can be gauged away by an operation Here we show two members (where spins point along the z-axis, blue/red circles denote spin up/down) that are related to each other by flipping the spins in every second ladder (shaded strips) of Fig. 1. The Bragg peaks corresponding to S z i S z j correlations are also shown in the extended Brillouin zone (assuming the same magnetic form factor in the two unit cell sublattices). The corresponding Bragg reflections for S x i S x j and S y i S y j are related to S z i S z j by C6v spin-orbit rotations [33].
, which is the product of π-rotations around the y, z, and x axis, respectively, for the B, C, and D sublattices of Fig. 1. This hidden duality is a very common feature in many spin-orbital models [9,34,35] but does not exist when Heisenberg couplings are also present (in contrast to the symmetry H xyz discussed below). Here it reduces our study to the first two quadrants of the unit circle of ψ. Figure 2 shows the quantum phase diagram as found by exact diagonalizations (ED) on finite clusters, see discussion below and numerical data shown in Fig. 3. There are six different regimes as a function of the angle ψ: the two quantum spin liquids (QSLs) regions (which have been enlarged for better visibility) around the exactly solvable Kitaev points (ψ = 0 and π) and four long-range magnetic regions (I-IV), hosting FM, Neel, stripy, as well as the zig-zag phases that are relevant for Na 2 IrO 3 (II) and α-RuCl 3 (IV). Under the duality transformation H yzx , the two QSLs map to each other, I maps to III, and II maps to IV.
Each of the magnetic regions actually hosts twelve degenerate quantum states, some of which are even qualitatively different among themselves, with very distinct Bragg reflections. For example, the region III hosts six FM and six stripy AFM ground states, and IV hosts six Néel and six zigzag AFM ground states. This striking aspect stems from a non-global symmetry, H xyz , which is the product of π-rotations around the x, y, and z axis, respectively, for the B, C, and D sublattices of Fig. 1. The two states shown in each magnetic region FIG. 3. (Color online) (a-b) Exact low-energy spectra (measured from the ground state energy E0) of the 24-site (a) and 32-site (b) clusters, defined, respectively [33], by the spanning vectors (2t1−4t2, 4t1−2t2) and (2t1−4t2, 4t1). A non-linear x-axis is used in order to highlight all regions of interest equally. The states are labeled by momenta k in the first BZ, parity ("e" for even, "o" for odd) under inversion through hexagon centers, and parity under global spin π-rotations around the x-axis ("Sze" for even, "Szo" for odd). The (red) numbers in (a) denote the multiplicity of the lowest five levels in regions I and II, and the ground state degeneracy at ψ = 0 and π. (c) Ground state expectation value Wp of Kitaev's flux operators. (d) Square root of the 'symmetrized' ground state spin structure factor S(Q) (see text), along with the spin length calculated from a self-consistent non-linear spin-wave theory (NLSWT).
of Fig. 2 are related to each other by this symmetry, which for these particular states amounts to flipping the z-component of the spins in every second shaded ladder of Fig. 1. The remaining ten states of the quantum ground state manifold arise by applying the global symmetries of the model: i) the double cover C 6v of C 6v , and ii) the double cover D 2 of the D 2 group of global π rotations in spin space.
Let us now turn to the numerical spectra shown in Fig. 3 (a,b). First, the QSL regions are extremely narrow: They survive in a tiny window of δψ = 0.05π around the exact Kitaev points, which is confirmed by the comparison of ED against large scale pseudofermion functional renormalization group (PFFRG) calculations [36][37][38][39]. So the QSLs are extremely fragile against K 2 .
Second, Fig. 3 (a,b) show very dense spectral features in the QSL regions, reflecting the continuum structure of fractionalized excitations above the Kitaev spin liquid. More specifically, for finite systems the ground state degeneracy at the exact Kitaev points [40] is lifted by K 2 . Still, for small enough |K 2 |, the QSLs must be gapless in the thermodynamic limit, because K 2 respects time reversal symmetry and is therefore not expected [10] to open a gap in the Majorana spectrum [41].
Third, unlike the QSL regions, the low-energy spectrum inside the magnetic regions is very discrete. In addition, most of the low-lying states within the energy window shown in Figs. 3 (a,b) correspond precisely to the twelve quantum ground states discussed above. For finite systems, these states are admixed by a finite tunneling, leading to twelve symmetric eigenstates with quantum numbers corresponding to the decomposition of the symmetry broken states. This decomposition is worked out in detail in [33] and is indeed fully consistent with the ED data. So the lowest twelve states in each magnetic region of Figs. 3 (a,b) will collapse to zero energy in the thermodynamic limit, leaving the true magnon excitations with a large anisotropy gap (modulo finite size corrections), reflecting the anisotropic, Ising-like character of the magnetic model.
Fourth, the magnetic instabilities, which serve as good examples of deconfinement-confinement transitions [42][43][44][45] for the underlying spinons, are of first order, as they are accompanied by finite, abrupt changes [46] in several ground state properties, e.g., in W p , and in the spin-spin correlations. Specifically, at ψ = 0 and π, all fluxes W p have a value of +1 [10]. A finite K 2 admixes sectors of different W p , and so W p drops continuously as we depart from the exact Kitaev's points, until it jumps to very low absolute values when we enter the magnetic phases, see Fig. 3 (c).
Turning to the spin-spin correlations, their abrupt change at the transition can be seen in the behavior of the 'symmetrized' spin structure factor S(Q) shown in Fig. 3 (d), which is defined as where N is the number of sites, Q (α) is the ordering wavevector (see below) of the α-th component of the spins (α = x, y, z), and the extra factor of 2 in this definition accounts for the fact [33] that, for finite systems, there are no correlations between NN ladders like the ones shaded in Fig. 1, due to the non-global symmetry H xyz discussed above. These data show clearly the short-range (long-range) character of spinspin correlations inside (outside) the QSL regions. This aspect can be seen more directly in Fig. 4, which shows the real-space spin-spin correlation profiles S α i S α j , in the three channels α = x, y, z, as calculated in the ground state of the 32-site cluster, inside the first QSL phase and slightly ψ =0.01π ψ =0.028π Real-space spin-spin correlation profiles evaluated at the ground state of the N = 32 cluster, inside the first QSL phase (ψ = 0.01π, left column) and inside the magnetic phase I (ψ = 0.028π, right column). Different rows correspond to the three different channels S α i S α j , α = x, y and z. The reference site i is indicated by the small black open circle. Positive (negative) correlations are shown by filled blue (filled red) circles, whose radius scales with the magnitude of the correlation. The difference between α = z and α = x, y stems from the fact that the 32-site cluster does not have the full point-group symmetry of the infinite lattice, and the momentum point Mz is not equivalent by symmetry to Mx and My, see [33]. outside (magnetic phase I). The results show clearly the ultra short-range nature of the correlations inside the QSL region, and the long-range nature outside.
Finally, the spin-spin correlation profiles demonstrate the special anisotropic character of the correlations, whereby different spin components α are correlated along different directions of the lattice (or, equivalently, different spin components α order at different ordering wavevectors Q (α) , see also Fig. 2), reflecting the locking between spin and orbital degrees of freedom in this model. Similar behavior is found for all other magnetic phases, including the zig-zag phases that are relevant for Na 2 IrO 3 and α-RuCl 3 . Such a signature of directional dependent Kitaev [7]; see also last paragraph of Sec. VII.
In the following we shall probe the physical mechanism of the spin liquid instabilities by taking one step back and examining the classical limit first.

III. Classical limit
For classical spins, the frustration introduced by the K 2 coupling is different from the one of the pure K 1 model studied by Baskaran et al [47]. A straightforward classical minimization in momentum space [33] gives lines of energy minima instead of a whole branch of minima [47], suggesting a sub-extensive ground state manifold structure, in analogy to compass-like models [48] or other special frustrated antiferromagnets [49].
We can construct one class of ground states by satisfying one of the three types of Ising bonds. We can choose for example the horizontal zz-bonds and align the spins along the z-axis with relative orientations dictated by the signs of K 1 and K 2 . The energy of the resulting configuration saturates the lower energy bound [33] and is therefore one of the ground states. We can then generate other ground states by noting that K 1 and K 2 fix the relative signs of the spin projections S z only within the vertical 2-leg ladders of the lattice (shaded strips in Fig. 1), but do not fix the relative orientation between different ladders, because these couple only via xx and yy Ising interactions which drop out at the mean field level. This freedom leads to 2 nlad ground states, where n lad ∝ √ N is the number of vertical ladders. This sub-extensive degeneracy stems from the presence of non-global, sliding operations [48,[50][51][52] of flipping S z → −S z for all spins belonging to one vertical ladder. Similarly, we can saturate the xx or the yy bonds, leading to 2-leg ladders running along the diagonal directions of the lattice. In total, this procedure delivers 3×2 nlad classical ground states.
These states are actually connected in parameter space by valleys formed by other, continuous families of ground states that can be generated by global SO(3) rotations of the discrete states [33]. The degeneracy associated with these valleys is accidental and can therefore be lifted by fluctuations. This is in fact the situation at finite T where thermal fluctuations select one of the three types of discrete ground states, thereby breaking the three-fold symmetry of the model in the combined spin-orbit space. This corresponds to a finite-T nematic phase where spins point along one of the three cubic axes but still sample all of the 2 nlad corresponding states, without any long-range magnetic order. To achieve the latter one needs to spontaneously break all sliding symmetries and this cannot happen at finite T , according to the generalized Elitzur's theorem of Batista and Nussinov [50]. The sliding symmetries can break spontaneously only at T = 0 and in all possible ways, which is reflected in the divergence of the spin structure factor along lines in momentum space.

IV. Quantum spins & Strong-coupling expansion
Turning to quantum spins, the situation is fundamentally different because the sliding symmetries are absent from the beginning: To flip one component of the spin we must combine a π-rotation in spin space and the time reversal operation [53]. The latter, however, involves the complex conjugation which cannot be constrained to act locally only on one ladder. Essentially, this means that the ladders must couple to each other dynamically by virtual quantum-mechanical processes, which in turn opens the possibility for long-range magnetic ordering even at finite T .
The natural way to understand the dynamical coupling between the ladders is to perform a perturbative expansion around one of the three strong coupling limits where the above discrete states become true quantum-mechanical ground states. Consider for example the limit where the xx and yy couplings, denoted by K , are much smaller than the zz couplings, K z 1 and K z 2 . Let us also parametrize K x(y) 1,2 = rK z 1,2 , K z 1 = cos ψ and K z 2 = sin ψ. For r = 0 we have n lad decoupled vertical ladders, and 2 nlad quantum ground states. Degenerate perturbation theory [33] then shows that the degeneracy is first lifted at fourth order in r via three, loop-four virtual processes that involve: (i) only K perturbations, see the top panel of Fig. 5.
The processes (i) give rise to intra-ladder, six-body terms which are nothing else than the flux operators W p . As shown by Kitaev [10], these terms can be mapped to the square lattice Toric code [54] which has a gapped spin liquid ground state. Next, the processes (ii) and (iii) give rise to effective, NNN inter-ladder couplings of the form JS z i S z j , where i and j have the same (ii) or different (iii) sublattice unit cell indices, see top panel of Fig. 5. To fourth-order in r, the corresponding couplings J W (i), J 1 (ii), and J 2 (iii) read Note that J 2 is always AFM and competes with J 1 in the regions I and III of Fig. 2. We also emphasize that there is no S z i S z j coupling when i and j belong to NN ladders. This is actually true to all orders in perturbation theory, because of the above non-global symmetry H xyz , which changes the sign of S z on every second vertical ladder (B and C sites of Fig. 1).
The main panel of Fig. 5 shows the behavior of |J W |/r 4 , 2|J 1 |/r 4 , and J 2 /r 4 as a function of the angle ψ, where the relative factor of 2 between |J 1 | and J 2 accounts for their relative contribution to the total classical energy. Close to the exactly solvable points ψ = 0 and π, the physics is dominated by the flux terms W p which, as mentioned above, lead to the gapped Toric code QSL [10,54]. The gapless QSL at r = 1 is eventually stabilized by off-diagonal processes that necessarily admix states outside the lowest manifold of the r = 0 point [55].
The four magnetic phases I-IV of Fig. 2 are all stabilized by J 1 which, according to Fig. 5, is the dominant coupling in a wide region away from ψ = 0 and π. Note that there are also two windows (shaded in Fig. 5) in the beginning of regions I and III where the two inter-ladder terms compete and 2|J 1 | < J 2 . This opens the possibility for two more states (the ones favored by J 2 ) in these regions. This scenario is however not confirmed by our ED spectra and spin structure factors (especially for the 32-site cluster which is commensurate with both types of competing phases), showing that these phases are eventually preempted by the QSLs and the phases I and III at higher values of r.
We remark here that the 1-loop formulation of PFFRG delivers the J 2 but not the J 1 processes because, in a diagrammatic formulation of Abrikosov fermions, these processes relate to 3-particle vertex contributions, which require a 2-loop formulation. However, for ψ around 0 and π, where J 1 is small, a 1-loop formulation already yields good agreement.

V. Semiclassical picture
The magnetic phases of the model can be captured by a standard semiclassical expansion, but this has to go beyond the non-interacting spin-wave level. Indeed, the zero-point energy of the quadratic theory lifts the accidental continuous degeneracy of the problem (selecting the cubic axes for the global direction in spin space, see Ref. [33]), but fails to lift the discrete 2 nlad degeneracy (the spectrum has lines of zero modes corresponding to the soft classical twists along individual ladders), and does not deliver a finite spin length, in analogy to several frustrated models [31,49,56,57]. The spurious zero modes are gapped out by spin-wave interactions, leading to the expected anisotropy gap and a finite spin length. The latter (obtained here from a self-consistent treatment of the quartic theory; details will be given elsewhere) tracks closely the behavior of the spin length extracted from the ED 'symmetrized' spin structure factor [58] S(Q), see Fig. 3 (d). Furthermore, both methods give values that are very close to the classical value of 1/2 inside the magnetic regions, showing that these phases are very robust. The quartic spin wave expansion is however insensitive to the proximity of the QSLs, most likely due to the first-order character of the transitions.

VI. Triangular Kitaev points
At ψ = ± π 2 the system decomposes into two interpenetrating triangular sublattices, where the K 2 coupling plays the role of a NN Kitaev coupling. This problem has been studied for both classical [28,29] and quantum spins [30][31][32]. The above analysis for the magnetic phases still holds here, the only difference being that the two legs of each ladder decouple, since they belong to different triangular sublattices. The ordering between the legs belonging to the same sublattice stems from the effective coupling J 1 , which is the only one surviving at K 1 = 0. This coupling connects NNN legs only, leading to twelve states in each sublattice and thus 12 2 states in total, instead of 12 for finite K 1 . The accumulation of such extra states at low energies can be clearly seen in Fig. 3(a-b) at ψ = ± π 2 . Note that while the ED spectra are broadly independent of system size, significant differences between the two cluster sizes are apparent near ψ = ±π/2. These differences, e.g. on the ground state multiplicity, can be easily traced back to the different point group symmetry of the two clusters, see detailed explanation in [33].
Finally we would like to point out that the origin of the ordering mechanism at the triangular Kitaev points has also been discussed independently in a recent paper by G. Jackeli and A. Avella [31].

VII. Discussion
Charting out the stability region of the Kitaev spin liquid is an extremely relevant endeavor for the synthesis and characterization of new materials. One of the counterintuitive results of this study is that the frustrating (with respect to longrange magnetic order) NNN coupling K 2 , which has exactly the same anisotropic form and symmetry structure as the K 1 term, destabilizes the Kitaev spin liquid much faster than the non-frustrating isotropic Heisenberg J 1 coupling. This finding gives a very useful hint in the search of realistic materials that exhibit the Kitaev spin liquid physics. In A 2 IrO 3 materials, for example, the role of the size of the central ion (Na in Na 2 IrO 3 , or Li in Li 2 IrO 3 ) in mediating the K 2 coupling (see also below) is a key aspect that can be easily controlled by experimentalists [59, 60].
On a more conceptual note, the physical mechanism underpinning the magnetic long range ordering in the present model is a novel example of order-by-disorder. Unlike many other classical states, here the ordering manifests only for quantum spins and not for classical spins. This striking contrast between classical and quantum spins is even more surprising in the light of the fact that all these phases have a strong classical character with local pseudo-spin lengths that are very close to the maximum classical value of 1/2.
On this issue, we should stress that there is no discrepancy between the very large pseudo-spin length that we report here and the small length of the magnetic moments extracted from magnetic reflections, e.g., in Na 2 IrO 3 [5]. Such an apparent discrepancy can be explained by the value of the g-factor which can be significantly smaller then 2, because the orbital angular momentum is not quenched in strong SOC compounds. For the ideal cubic symmetry, for example, the well-known Landé formula gives g = 2/3, and similar values could be expected for lower symmetry.
Let us now elucidate further our main reasons on why the K 2 coupling must play an important role in Na 2 IrO 3 , and can be relevant in Li 2 IrO 3 and α-RuCl 3 : i) The super-exchange expansion of [20] shows clearly that the NNN Kitaev coupling is the second largest term in Na 2 IrO 3 , with K 2 7-9 meV. All other perturbations are at most 1-2 meV, consistent with the numbers given by the large-scale ab initio quantum chemistry study of [15]. The mechanism behind the large magnitude of K 2 in Na 2 IrO 3 is physically very clear: It originates from the large diffusive Na ions that reside in the middle of the exchange pathways, and the constructive interference of a large number of four pathways [20].
In Li 2 IrO 3 , the K 2 interaction comes from the same mechanism but it is relatively smaller because of the smaller size of Li ions [26]. Still, as discussed in [22], this coupling can be important to explain the current experimental evidence in terms of magnetic susceptibility profile, Curie-Weiss temperature, and the relevant range of couplings.
Finally, in α-RuCl 3 , the analogous super-exchange path is absent, but an appreciable K 2 still arises from the anisotropy of diagonal interactions originated from the interplay between different hopping processes [26]. However, as we already pointed out in the Introduction, the second largest coupling in α-RuCl 3 is the anisotropic exchange Γ [15,27]. According to the study of J. Rau et al. [27], a positive Γ seems to compete with K 2 for positive K 1 [26]. However, the situation is still unclear since the Bragg peaks of the states favored by Γ do not reside at the M points of the BZ found experimentally by J. A. Sears et al. [13], whereas such Bragg peaks are naturally present in the zig-zag phases favored by K 2 , or even by a negative J 1 . So a lot more work is needed to clarify the relative importance of Γ, K 2 , and J 1 in α-RuCl 3 .
ii) The K 2 coupling explains naturally the zig-zag ordering in Na 2 IrO 3 . This phase cannot arise in the original J 1 -K 1 model, because this would require an AFM coupling K 1 , whereas it is widely accepted that K 1 is FM and large in magnitude, see e.g. [16]. Also, the much smaller Γ terms, which are positive, also favor the zig-zag phase and do not compete with K 2 , according to [27].
iii) The K 2 coupling can provide in addition the basis to resolve the long-standing puzzle of the large AFM Curie-Weiss temperature [2,3,6], without incorporating unrealis-tically large values of longer-range Heisenberg couplings J 2 and J 3 .
iii) The recent diffusive x-ray scattering experiments by S. H. Chun et al. [7] have provided direct evidence for the predominant role of anisotropic, bond directional interactions in Na 2 IrO 3 . In conjunction with the above discussion and the results of Fig. 4, the K 2 term then emerges naturally as the number one anisotropic candidate term that can drive the zigzag ordering and the directional dependence of the scattering found in [7].
An aspect that remains to be discussed in the context of Na 2 IrO 3 is the direction of the magnetic moments which, according to the x-ray scattering data of S. H. Chun et al. [7], do not point along the cubic axes but along the face diagonals. As discussed above, the K 2 coupling stabilizes the zig-zag phase but it is unable to lock the direction of the moments at the mean-field level due to an infinite accidental degeneracy. The fact that the locking along the cubic axes in the K 1 -K 2 model eventually proceeds via a quantum order-by-disorder process (see Ref. [33]) renders this result very susceptible to much smaller anisotropic interactions that can pin the direction of the moments already at the mean field level. A very small positive anisotropic Γ term can for example play such a role and can account for the locking along the face diagonals, as can be directly seen by a straightforward minimization of the classical energy. An alternative scenario involves a competing order-by-disorder effect within a more extended model that includes weak longer-range exchange interactions [26].
Acknowledgements. We acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota and the Max Planck Institute for the Physics of Complex Systems, Dresden, where a large part of the numerical computations took place. We are also grateful to R. Moessner, C.  [58] The extra factor of 2 in this definition accounts for the fact that there are no correlations between NN ladders for finite systems, due to the symmetry Hxyz, see also [33]. i

Supplemental material
In this Supplementing material we provide auxiliary information and technical details and derivations. Specifically, Sec. A deals with the Luttinger-Tisza minimization of the classical energy in momentum space (1), and the order-by-disorder process by harmonic spin-waves (2). Sec. B gives details about our finite-size ED study, including the symmetry analysis of the low-energy spectra in regions I and II of the phase diagram (3), and the definition of the 'symmetrized' spin structure factor S(Q). In Sec. C we provide results from the pseudofermion functional renormalization group (PFFRG) approach. Finally, in Sec. D we provide the derivation of the effective Hamiltonian around the strong coupling limit of K x(y) 1,2 = 0.
A. Semiclassical analysis

Lutinger-Tisza minimization
We choose the primitive vectors of the honeycomb lattice as t 1 = ay and t 2 = (− √ 3 2 x+ 1 2 y)a, where a is a lattice constant, see Fig. 1 of the main paper. We also define t 3 = t 1 −t 2 = √ 3 2 x+ 1 2 y. In the following, we label the Bravais lattice vectors as R = nt 1 +mt 2 , where n and m are integers. We also denote the two sites in the unit cell by a sublattice index i = 1-2. The total classical energy of the K 1 -K 2 model reads where N uc = N/2, is the number of unit cells, and the matrices Λ (α) (where α = x, y, z) are given by To find the classical minimum we need to minimize the energy under the strong constraints S 2 R,i = S 2 , ∀(R, i). The Luttinger-Tisza method [1][2][3][4] amounts to relax the strong constraints with the weaker one R,i S 2 R,i = N S 2 , or equivalently k,i S k,i · FIG. 1. The first two Brilouin zones of the honeycomb lattice, along with the special lines in momentum space Q (x) , Q (y) , and Q (z) (respectively Q (x) , Q (y) , and Q (z) ) corresponding to the minima of the classical energy for K2 > 0 (< 0), see text.
ii S −k,i = S 2 . If we can find a minimum under the weak constraint that also satisfies the strong constraints then we have solved the problem. To this end, we minimize the function with respect to {S α −k,i }, which gives a set of three eigenvalue problems for the Λ matrices: If we can satisfy these three relations (plus the strong constraint) with a single eigenvalue λ, then = λS 2 . So the energy minimum corresponds to the minimum over the three eigenvalues λ (α) of the matrices Λ (α) (−k), and over the whole Brillouin zone (BZ). The eigenvalues of these matrices and the corresponding eigenvectros are: For K 2 positive, the minima of λ ± are located on the lines Q (x) = r(G 1 +G 2 )+(l+ 1 2 )G 2 , Q (y) = rG 1 +(l+ 1 2 )G 2 , and Q (z) = rG 2 + (l + 1 2 )G 1 , respectively, where l is any integer and r ∈ (− 1 2 , 1 2 ). On the other hand, for K 2 negative, the minima are located on the lines: Q (x) = r(G 1 +G 2 )+lG 2 , Q (y) = rG 1 +lG 2 , and Q (z) = rG 2 +lG 1 . Both sets of lines are shown in Fig. 1.
Let us now try to build a ground state from the minima of the above eigenvectors for the case K 1,2 > 0, by using the line of minima Q (z) as follows: where we used the relation R = nt 1 + mt 2 and have defined ξ m ≡ 1/2 −1/2 drf (r)e i2πmr , which is the Fourier transform of the envelope function f (r). We still need to satisfy the spin length constraint, which imposes a condition that the inverse Fourier transform of f (r) takes only the values ±1. This freedom corresponds to the sliding symmetries of flipping individual vertical ladders, and leads to 2 nlad degenerate states (where n lad is the number of vertical ladders), as discussed in the main text.
Similarly we can construct another 2 × 2 nlad states by using the lines Q (x) or Q (y) in momentum space, which correspond to decoupled ladders running along the diagonal directions of the lattice. Altogether, we have found the 3 × 2 nlad discrete classical ground states discussed in the main text by using the Luttinger-Tisza minimization method.
Finally, it is easy to see that we can also combine the three types of states into a continuous family of other ground states that include coplanar and non-coplanar states. This family can be parametrized by two angles θ and φ as follows, where i = 1, 2 and S x R,i , S y R,i and S z R,i denote the three type of discrete solutions discussed above.

Harmonic order-by-disorder
As we claimed in the main text, harmonic spin waves lift the accidental continuous degeneracy of the classical ground state manifold and select the discrete 3 × 2 nlad states, whereby spins point along the cubic axes. Here we shall demonstrate this result by considering a one-parameter family of coplanar states obtained by linearly combining two zigzag states and two stripy states with spins pointing along the cubic axes. In the resulting family of states, spins are pointing in some direction on the zx-plane. Figure 2 shows the two zigzag and two stripy phases with spins pointing along the cubic axes. Here "yz-zigzag//x" denotes a zigzag state with FM zig-zag lines running along the yy and zz bonds of the Kitaev Hamiltonian, and the spins point along the x-axis. Similarly, "x-stripy//z" denotes a stripy state with FM ladders formed by the xx bonds of the Kitaev Hamiltonian, and the spins point along the z-axis. Specifically, these states can be written as: 3 , π and M y = π √ 3 , π (see Fig. 1) and R = nt 1 + mt 2 . The one-parameter family of classical ground states are obtained by linear combinations of the above states: where ζ = 1 for the zigzag case and ζ = −1 for the stripy case. The effect of harmonic spin waves can be found by a standard linear spin-wave expansion around the corresponding states for each value of θ. Figs. 2 (e-f) show the zero-point energy correction (per number of unit cells) as a function of the angle θ for a representative point inside region II (ψ = 0.8π, ζ = 1) and another point inside region I (ψ = 0.3π, ζ = −1). The data show clearly that harmonic fluctuations select the states with the spins pointing along the cubic axes (θ = 0, ±π/2, and π).
We have checked that the result is the same for the corresponding order-by-disorder process for the one-parameter family of states obtained by combining two states with the same wavevector, such as the "zx-zigzag // z" and "zx-zigzag // x". B. Technical details about the ED study

The symmetry group of the Hamiltonian
The full symmetry group of the K 1 -K 2 model, for half-integer spins, is T × C 6v × D 2 , which consists of: iv  A1 1 1  1  1  1  1  1  1  1  A2 1 1  1  1  1  1 The translation group T generated by the primitive translation vectors t 1 and t 2 , see Fig. 1 of the main text.
2. The double cover C 6v of the group C 6v ⊂ SO(3) in the combined spin and real space, where the six-fold axis goes through one of hexagon centers. This group is generated by two operations: the six-fold rotation C 6 around [111], whose spin part maps the components (x, y, z) → (y, z, x), and the reflection plane (110) that passes through the zz-bonds of the model, whose spin part maps (x, y, z) → (−y, −x, −z).
3. The double cover D 2 of the point group D 2 ⊂ SO (3), which consists of three π-rotations C 2x , C 2y , and C 2z in spin space. The first maps the spin components (x, y, z) → (x, −y, −z), etc.

Finite clusters
In our ED study we considered two clusters with periodic boundary conditions, one with 24 and another with 32 sites, with spanning vectors (2t 1 −4t 2 , 4t 1 −2t 2 ) and (2t 1 −4t 2 , 4t 1 ), respectively. These clusters are shown in Fig. 3 (a, c). The 24-site cluster has the full point group symmetry of the infinite lattice, i.e. C 6v × D 2 , whereas the 32-site cluster has the lower symmetry C 2v × D 2 , where C 2v contains the reflection planes (110) and (110). Turning to translational symmetry, the allowed momenta for each cluster are shown in Fig. 3(b, d). Both clusters accommodate the three M points of the Brillouin zone (BZ) and are therefore commensurate with all magnetic states of the phase diagram. The difference between the two clusters is that the three M points are degenerate for N = 24 but not for N = 32.
In our ED study we have exploited: i) translations, ii) the C 2 subgroup of full C 6v point group (which is equivalent to the inversion I in real space through the hexagon centers), and iii) the global spin inversion which maps the local S z basis states | ↑ → | ↓ . This operation is described by i σ x i , which is nothing else than the global π-rotation C 2x in spin space, divided by

Mz
My Mx  For N = 24 and 32, the product of all these phase factors give +1.
Consequently, the energy eigenstates are labeled by: i) the momentum k, ii) the parity under C 2 ('e' for even, 'o' for odd), and iii) the parity under S z spin inversion ('Sze' for even, 'Szo' for odd).

Symmetry spectroscopy of classical phases
Here we derive the symmetry decomposition of the twelve magnetic states of region I and II of the phase diagram. As explained in the main paper, the other two regions, III and IV, map to I and II, respectively, by the hidden duality of H yxz followed by a simultaneous change of sign in K 1 and K 2 .
a. Phase I In the following, |str, α β denotes the stripy state with FM ladders running along the direction of the α-bonds, and the spins pointing along β in spin space. The twelve magnetic states of region I of the phase diagram can be split into four groups: S1 = {|str, x z , |str, y x , |str, z y }, S1 = {|str, x −z , |str, y −x , |str, z −y }, Table II shows how these twelve states transform under some of the symmetry operations of the group. Let us first examine the translation group. We have, ∀β: Tt 1 · |str, y β = |str, y −β , Tt 2 · |str, y β = |str, y β , Tt 1 · |str, z β = |str, z β , Tt 2 · |str, z β = |str, z −β . Next, let us examine the parities with respect to the C 2 rotation in real space and the C 2x rotation in spin space. It is easy to see that the first symmetry is not broken by any of the twelve states, while the second is broken when β = y and z. So all twelve vi states are even with respect to C 2 , the β = x are even with respect to C 2x , while β = y and z must decompose into both even and odd parities with respect to C 2x . Altogether: 'Extra' degeneracy at the M points for N = 24. The above quantum numbers for the M points are fully consistent with what we find in the low-energy spectra of Fig. 3 (a) of the main paper. For the symmetric, N = 24 cluster, the three M points are degenerate due to the six-fold symmetry. However we see that the two sets of M points are also degenerate with respect to each other, i.e. we have a six-fold degeneracy. This extra degeneracy comes from the D 2 symmetry in spin space. To see this, let us relabel the spin inversion part of (B1) using the actual IR of the group D 2 (see Table I, right), instead of the parity with respect to C 2x (which contains less information about the state): We see that the two states belonging to a given M point transform differently under D 2 , so the Hamiltonian does not couple the two states. Yet, these states are mapped to each other by one of the reflection planes of C 6v , so they must be degenerate, leading to an overall six-fold degeneracy at the M points.
Degeneracies at the Γ point for N = 24. The little group of the Γ point is the full point group C 6v × D 2 . However, all of the above six states that belong to the Γ point belong to the identity IR of D 2 , so it is enough to decompose them with respect to the C 6v part of the little group. To this end we use the well known formula from group theory [5] which gives the number of times m α that the α-th IR of C 6v appears in the decomposition of the 6 × 6 representation formed by the six states belonging to the Γ point. Here X(g) gives the character of this representation, while χ α (g) is the character of the α-th IR of C 6v , see Table I (left). From Table II it follows that X(g) is finite only for the elements E, E, C 2 , and C 2 , and using the characters of Table I (left) we find that the only finite m α are the following: m A1 = m A2 = 1, m E2 = 2, namely i.e. we expect two singlets and two doublets. All states are found in the low-energy spectra shown in Fig. 3 (a) of the main paper, where the degeneracy of the E 2 levels has been confirmed numerically.

b. Phase II
Here we denote by |zig, αα β the zigzag state with FM lines formed by consecutive α and α type of bonds, and the spins pointing along β in spin space. The twelve magnetic states of region II can be split into four groups: Under T and C 2x in spin space, these states transform in analogous way with the twelve states of region I, see (B1). The difference is that the present states break the C 2 rotation around the hexagon centers, and therefore the decomposition will contain both even and odd parities with respect to C 2 . Specifically, In analogy with region I, for the symmetric 24-site cluster, the six states belonging to the M points are degenerate due to the additional D 2 symmetry, and the six states belonging to the Γ point decompose as in (B4), namely 6Γ → A 1 ⊕ A 2 ⊕ 2E 2 . Again, all states are found in the low-energy spectra shown in Fig. 3 (a) of the main paper. vii c. Special points ψ = ±π/2: Different ground state structure for N = 24 and N = 32 As shown in Figs. 3(a) and (b) of the main text, the ED results are broadly independent of system size but significant differences between the two cluster sizes are apparent for the GS structure near ψ = ±π/2. The reason behind this difference lies in the different point group symmetry of the two clusters. The 24-site cluster has the full point group symmetry of the infinite lattice, whereas the 32-site cluster does not. This is also true for the two triangular sublattices of each cluster at ψ = ±π/2, where they become independent from each other.
Due to the high symmetry, each of the 12-site sublattices of the 24-site cluster have a two-fold degenerate ground state at ψ = ±π/2; let us denote them by |α and |β . On the other hand, the lower symmetry of the 16-site sublattices of the 32-site cluster leads to a single, non-degenerate ground state; let us denote it by |γ . Now, the global ground state structure of the two clusters at ψ = ±π/2 follows simply by taking the tensor product of the ground state manifolds in each sublattice. The 24-site cluster has four ground states: |α sub1 ⊗ |α sub2 , |β sub1 ⊗ |β sub2 , |α sub1 ⊗ |β sub2 , |β sub1 ⊗ |α sub2 . (B6) The first two states belong to the representation Γ.e.Sze, i.e. they have even parity with respect to inversion through the middle of the hexagons (this operation maps one sublattice to the other), and the same is true for the combination 1 √ 2 (|α sub1 ⊗ |β sub2 + |β sub1 ⊗ |α sub2 ). The remaining, antisymmetric combination, 1 √ 2 (|α sub1 ⊗ |β sub2 − |β sub1 ⊗ |α sub2 ), belongs to Γ.o.Sze, i.e. it has odd parity. This is in perfect agreement with the ED data.
For the 32-site cluster on the other hand, there is only one global ground state, namely |γ sub1 ⊗ |γ sub2 , which has even parity, again in agreement with the ED data.
Of course, as we discuss in the main text, in the thermodynamic limit a large number of states (12 2 ) will collapse to the ground state, which is how the corresponding symmetry-broken (classical) states are eventually formed.

'Symmetrized' spin structure factor and spin length
Here we discuss the 'symmetrized' spin structure factor S(Q) and explain the overall normalization factor that we use to extract the spin length. As we discuss in the main text, NN ladders do not couple by the symmetry H xyz , and so the quantum ground state of a finite cluster contains both relative orientations of the two sets of ladders L 1 and L 2 with equal amplitude. As a result, the spin-spin correlations between two spins that belong to L 1 and L 2 are zero for any finite cluster. If we wish to calculate the local spin lengths from the ground state spin-spin correlation data we can calculate the 'symmetrized' spin structure factor for one of the two subsets of ladders only, say L 1 : where N 1 = N/2 is the number of sites inside the sublattice L 1 , and Q (a) is the ordering wavevector corresponding to the spin component α = {x, y, z}. By translation symmetry, where we have chosen a reference site r = 0. The local spin length m is then given by m 2 = 2 N S 1 (Q). By contrast, the corresponding 'symmetrized' spin structure factor of the full lattice S(Q), defined by S(Q) = 1 N 2 α r,r ∈ L1∪L2 S α r S α r e iQ (α) ·(r−r ) , would give in the present case and the corresponding local spin lengths would be off by a factor of √ 2.
viii FIG. 4. Static spin-structure factor χ zz (k) plotted in the extended Brillouin zone (black lines inside the plotted region mark the boundaries of first Brillouin zone) for various values of ψ in the Kitaev spin-liquid phase. Note that χ xx (k) (χ yy (k)) are related to χ zz (k) by clockwise (counterclockwise) 2π/3-rotations in k-space.

C. Pseudofermion functional renormalization group (PFFRG) approach
In addition to ED we studied the K 1 -K 2 honeycomb model using the pseudofermion functional renormalization group (PF-FRG) approach. Rewriting the spin operators in terms of Abrikosov auxiliary fermions, the resulting fermionic model can be efficiently treated using a one loop functional renormalization group procedure. This technique calculates diagrammatic contributions to the spin-spin correlation function in infinite order in the exchange couplings, including terms in different interaction channels: The inclusion of direct particle-hole terms insures the correct treatment of the large spin limit S → ∞ while the crossed particle-hole and particle-particle terms lead to exact results in the large N limit. This allows to study the competition between magnetic order tendencies and quantum fluctuations in an unbiased way. For details we refer to reader to Ref. [6].
The PFFRG method calculates the static spin-structure factor as given by with S α (k, τ ) = 1 √ N i e −ikri e Hτ S z i e −Hτ , where τ denotes the imaginary time and T τ is the corresponding time-ordering operator. Being able to treat large system sizes (calculations for the K 1 -K 2 model are performed for a spin cluster with 265 sites) the PFFRG yields results close to the thermodynamic limit. Fig. 4 shows three representative plots for the momentum resolved spin-structure factor χ zz (k) in the Kitaev spin-liquid phase in the vicinity of ψ = 0. While in the exact Kitaev limit ψ = 0 the PFFRG reproduces the well-known nearest neighbor correlations as indicated by a single harmonics profile of the spin-structure factor, deviations from ψ = 0 lead to longer-range correlations and a more diverse spin-structure factor.