Kitaev Spin Liquid in 3d Transition Metal Compounds

We study the exchange interactions and resulting magnetic phases in the honeycomb cobaltates. For a broad range of trigonal crystal fields acting on Co2+ ions, the low-energy pseudospin-1/2 Hamiltonian is dominated by bond-dependent Ising couplings that constitute the Kitaev model. The non-Kitaev terms nearly vanish at small values of trigonal field \Delta, resulting in spin liquid ground state. Considering Na3Co2SbO6 as an example, we find that this compound is proximate to a Kitaev spin liquid phase, and can be driven into it by slightly reducing \Delta by \sim 20 meV, e.g., via strain or pressure control. We argue that due to the more localized nature of the magnetic electrons in 3d compounds, cobaltates offer the most promising search area for Kitaev model physics.

We study the exchange interactions and resulting magnetic phases in the honeycomb cobaltates. For a broad range of trigonal crystal fields acting on Co 2+ ions, the low-energy pseudospin-1/2 Hamiltonian is dominated by bond-dependent Ising couplings that constitute the Kitaev model. The non-Kitaev terms nearly vanish at small values of trigonal field ∆, resulting in spin liquid ground state. Considering Na3Co2SbO6 as an example, we find that this compound is proximate to a Kitaev spin liquid phase, and can be driven into it by slightly reducing ∆ by ∼ 20 meV, e.g., via strain or pressure control. We argue that due to the more localized nature of the magnetic electrons in 3d compounds, cobaltates offer the most promising search area for Kitaev model physics. The Kitaev honeycomb model [1], demonstrating the key concepts of quantum spin liquids [2] via an elegant exact solution, has attracted much attention (see the recent reviews [3][4][5][6][7]). In this model, the nearest-neighbor (NN) spins S = 1/2 interact via a simple Ising-type coupling S γ i S γ j . However, the Ising axis γ is not global but bond-dependent, taking the mutually orthogonal directions (x, y, z) on the three adjacent NN-bonds on the honeycomb lattice. Having no unique easy-axis and being frustrated, the Ising spins fail to order and form instead a highly entangled quantum many-body state, supporting fractional excitations described by Majorana fermions [1].
Much effort has been made to realize the Kitaev spin liquid (SL) experimentally. From a materials perspective, the Ising-type anisotropy is a hallmark of unquenched orbital magnetism. As the orbitals are spatially anisotropic and bond-directional, they naturally lead to the desired bond-dependent exchange anisotropy via spin-orbit coupling [8]. Along these lines, 5d iridates have been suggested [9] to host Kitaev model; later, 4d RuCl 3 was added [10] to the list of candidates. To date, however, the Kitaev SL remains elusive, as this state is fragile and destroyed by various perturbations, such as small admixture of a conventional Heisenberg coupling [11] caused by direct overlap of the d orbitals. Even more detrimental to Kitaev SL are the longer range couplings [10], unavoidable in weakly localized 5d-and 4d-electron systems with the spatially extended d wavefunctions. We thus turn to 3d systems with more compact d orbitals [13].
While the idea of extending the search area to 3d materials is appealing, and plausible theoretically [4,16], it raises an immediate question crucial for experiment: Is spin-orbit coupling (SOC) in 3d ions strong enough to support the orbital magnetism prerequisite for the Kitaev model design? This is a serious concern, since non-cubic crystal fields present in real materials tend to quench orbital moments and suppress the bond-dependence of the exchange couplings [8]. In this Letter, we give a positive answer to this question. Our quantitative analysis of the crystal field effects on the magnetism of 3d cobaltates shows that the orbital moments remain active and generate a Kitaev model as the leading term in the Hamiltonian. Most importantly, we identify the trigonal crystal field as the key and experimentally tunable parameter, which decides the strength of the non-Kitaev terms in 3d compounds.
Our main results are summarized in Fig. 1, displaying various magnetic phases of spin-orbit entangled pseudospin-1/2 Co 2+ ions on a honeycomb lattice. The phase diagram is shown as a function of trigonal field ∆, in a window relevant for honeycomb cobaltates, and a ra- The calculated magnetic phase diagram of honeycomb cobaltates. The Kitaev SL phase is surrounded by ferromagnetic (FM) states with moments in the honeycomb ab-plane and along the c-axis, zigzag-type states with moments in the ab-plane (zz1), along Co-O bonds (zz2), and in the ac-plane (zz3). Vortex-and stripy-type phases take over at smaller U/∆ pd . The color map shows the second-NN spin correlation strength (divided by S 2 = 1/4), which drops sharply in the SL phase. The star indicates the rough position of Na3Co2SbO6. tio of Coulomb repulsion U and the charge-transfer gap ∆ pd [17]. From the analysis of experimental data, we find that Na 3 Co 2 SbO 6 [15,18,19] is located at just ∼ 20 meV "distance" from the Kitaev SL phase (see Fig. 1), and could be driven there by a c-axis compression that reduces ∆. This seems feasible, given that ∆ variations within a window of ∼ 70 meV were achieved by strain control in a cobalt oxide [21].
We now describe our calculations resulting in Fig. 1. In short, we first derive the pseudospin exchange interactions from a microscopic theory, as a function of various parameters, and then obtain the corresponding ground states numerically by exact diagonalization.
Exchange interactions.-In an octahedral environment, Co 2+ ion with t 5 2g e 2 g configuration possesses spin S = 3/2 and effective orbital moment L = 1, which form, via spin-orbit coupling, a pseudospin S = 1/2 [1]. Over decades, cobaltates served as a paradigm for quantum magnetism, providing a variety of pseudospin-1/2 models ranging from the Heisenberg model in perovskites with corner-sharing octahedra [22,23] to the Ising model when the CoO 6 octahedra share their edges [24].
A microscopic theory of Co 2+ interactions in the edgesharing geometry has been developed just recently [4,16], assuming an ideal cubic symmetry. Here we consider a realistic case of trigonally distorted lattices, where t 2g orbitals split as shown in Fig. 2(a). Our goal is to see if such distortions leave enough room for the Kitaev model physics in real compounds. This is decided by the spinorbital structure of the pseudospin S = 1/2 wavefunctions. In terms of |S Z , L Z states, where the trigonal axis Z c is perpendicular to the honeycomb plane, they read as: The coefficients C 1,2,3 depend on a relative strength ∆/λ of the trigonal field ∆(L 2 Z − 2 3 ) and SOC λL · S [2,26]. At ∆ = 0, one has (C 1 , , and all the three components of L are equally active. A positive (negative) ∆ field tends to quench L Z (L X/Y ). The next step is to project various spin-orbital exchange interactions in cobaltates [4] onto the above pseudospin-1/2 subspace. The calculations are standard but very lengthy, and the readers interested in details are referred to the Supplemental Material [26]. At the end, we obtain the S = 1/2 Kitaev model K S γ i S γ j , supplemented by Heisenberg J and off-diagonal anisotropy Γ, Γ terms; for γ = z type NN-bonds, they read as: Interactions H (γ) ij for γ = x, y type bonds follow from a cyclic permutation among S x j , S y j , and S z j , as dictated by lattice symmetry. While the Hamiltonian (S6) is of the same form as in d 5 Ir/Ru systems [3,5], the microscopic origin of its parameters K, J, Γ, Γ is completely different in d 7 Co compounds. This is due to the spin-active e g electrons of Co(t 5 2g e 2 g ) ions, which generate new spin-orbital exchange channels t 2g -e g and e g -e g , shown in Fig. 2(b), in addition to the t 2g -t 2g ones operating in d 5 systems with t 2g -only electrons. In fact, the new terms make a major contribution to the exchange parameters, as illustrated in Figs. 2(c)-2(f). In particular, Kitaev coupling K comes almost entirely from the t 2g -e g process. It is also noticed that t 2g -e g and e g -e g contributions to J, Γ, and Γ are of opposite signs and largely cancel each other, resulting in only small overall values of these couplings. Figure 2 shows that the trigonal field ∆, which acts via modification of the pseudospin wavefunction (S2), has an especially strong impact on the non-Kitaev couplings J, Γ, and Γ . As a result, the relative strength (J/K, etc) of these "undesired" terms is very sensitive to ∆ variations. This suggests the orbital splitting ∆ as an efficient (and experimentally accessible) parameter that controls the proximity of cobaltates to the Kitaev-model regime.
Another important parameter in the theory is the U/∆ pd ratio. In contrast to Ir/Ru based Mott insulators with small U/∆ pd ∼ 0.5, cobaltates are charge-transfer insulators [17], with typical values of U/∆ pd ∼ 2 − 3 depending on the material chemistry. Including both Mott-Hubbard U and charge-transfer ∆ pd excitations, we have calculated [26] the exchange couplings as a function of U/∆ pd and ∆/λ. Figure 3(a) shows that Kitaev coupling K is not much sensitive to U/∆ pd variations. On the other hand, the non-Kitaev terms, especially Heisenberg coupling J, are quite sensitive to U/∆ pd , see Figs. 3(b)-3(d). However, their values relative to K remain small over a broad range of parameters.
Phase diagram.-Having quantified the exchange parameters in Hamiltonian (S6), we are now ready to address the corresponding ground states. As Kitaev coupling is the leading term, the model is highly frustrated. We therefore employ exact diagonalization which has been widely used to study phase behavior of the extended Kitaev-Heisenberg models (see, e.g., Refs. [11][12][13][28][29][30]). In particular, by utilizing the method of coherent spin states [12,13], we can detect and identify the magnetically ordered phases (including easy-axis directions for the ordered moments). When non-Kitaev couplings are small (roughly below 10% of the FM K value), a quantum spin-liquid state is expected. Reflecting the unique feature of the Kitaev model [1], this state is characterized by short-range spin correlations that are vanishingly small beyond nearest-neighbors [11].
The resulting phase diagram, along with the data quantifying spin correlations beyond NN distances, is presented in Fig. 3(e,f). The main trends in the phase map are easy to understand considering the variations of non-Kitaev couplings with ∆/λ and U/∆ pd . As we see in Fig. 3(c,d), Γ exactly vanishes at the ∆ = 0 line, and Γ is very small too. Thus, in the cubic limit, the model (S6) essentially becomes the well studied K − J model, with large FM Kitaev K term, and J correction changing from AF J > 0 to FM J < 0 as a function of U/∆ pd . Consequently, the ground state changes from stripy AF (at small U/∆ pd ) to FM order at large U/∆ pd , through the Kitaev SL phase in between [28]. In the SL phase, spin correlations are indeed short-ranged and bond-selective: for z-type NN bonds, we find S z S z / S 2 0.52 (as in the Kitaev model), while they nearly vanish at farther distances, see Fig. 3(e,f).
As we switch on the trigonal field ∆, the Γ term comes into play confining the SL phase to the window of |∆|/λ < 1 (where |Γ /K| < 0.1). In the FM phases, the sign of Γ decides the direction of the FM moments, as can be expected from classical considerations. On the left-top (left-bottom) part of the phase map, where Heisenberg coupling J is AF, we observe that the stripy state gives way to a vortex-type [3] (zigzag-type) ordering, stabilized by the combined effect of Γ and Γ terms.
To summarize up to now, the nearest-neighbor pseudospin Hamiltonian is dominated by the FM Kitaev model, which appears to be robust against trigonal split- ting of orbitals. Subleading terms, represented mostly by J and Γ couplings, shape the phase diagram, which includes a sizeable SL area. While these observations are encouraging, it is crucial to inspect how the picture is modified by longer range interactions, especially by the third-NN Heisenberg coupling J 3 S i · S j , which appears to be the major obstacle on the way to a Kitaev SL in 5d and 4d compounds [5,10]. We have no reliable estimate for J 3 , since long-range interactions involve multiple exchange channels and are thus sensitive to material chemistry details. As such, they have to be determined experimentally. We note that |J 3 /K| 0.1 was estimated [33,34] in the 4d compound RuCl 3 ; in cobaltates with more localized 3d orbitals [13], this ratio is expected to be smaller.
Adding a J 3 term to the model (S6), we have re-examined the ground states and found that the Kitaev SL phase is stable up to |J 3 /K| ∼ 0.06 [26]. The modified phase diagram, obtained for a representative value of J 3 = 0.15t 2 /U 0.04|K|, is shown in Fig. 1 [35]. Its comparison with Fig. 3 tells that the main effect of J 3 is to support the zigzag-type states (with different orientation of moments) at the expense of other phases. Note also that the SL area is shifted to the right, where FM J and AF J 3 tend to frustrate each other. The phase diagram in Fig. 1 should be generic to Co 2+ honeycomb systems, and will be used in the following discussion.
Traditionally, experimental data in Co 2+ compounds is analysed in terms of an effective S = 1/2 models of XXZ-type [41,43,47,[50][51][52]. As S = 1/2 magnons (∼ 10 meV) are well separated from higher lying spinorbit excitations (∼ 30 meV), the pseudospin picture itself is well justified; however, a conventional XXZ model neglects the bond-directional nature of pseudospin S = 1/2 interactions. A general message of our work is that a proper description of magnetism in cobaltates should be based on the model of Eq. (S6), supplemented by longer-range interactions. We note in passing that the XXZ model also follows from Eq. (S6) when the Kitaev-type anisotropy is suppressed [3]; however, such an extreme limit is unlikely for realistic trigonal fields, given the robustness of the K coupling, see Fig. 3.
As an example, we consider Na 3 Co 2 SbO 6 which has low Néel temperature and a reduced ordered moment [15]. Analysing the magnetic susceptibility data [15] including all spin-orbit levels [26], we obtain a positive trigonal field ∆ 38 meV and λ 28 meV; these values are typical for Co 2+ ions in an octahedral environment (see, e.g., Ref. [47]). With ∆/λ 1.36, we evaluate S = 1/2 doublet g-factors g ab 4.6 and g c 3, from which a saturated moment of 2.3µ B , consistent with the magnetization data [15], follows.
Zigzag-ordered moments in Na 3 Co 2 SbO 6 are confined to the ab-plane [15]; this corresponds to the zz1 phase in Fig. 1. The easy-plane anisotropy is due to the Γ term, which is positive for ∆ > 0, see Fig. 3(d). Regarding the location of Na 3 Co 2 SbO 6 on the U/∆ pd axis of Fig. 1, we believe it is close to the FM//ab phase, based on the following observations. First, a sister compound Li 3 Co 2 SbO 6 has ab-plane FM order [37] (most likely due to smaller Co-O-Co bond angle, 91 • versus 93 • , slightly enhancing the FM J value). Second, zigzag order gives way to fully polarized state at small magnetic fields [15,18]. These facts imply that zz1 and FM//ab states are closely competing in Na 3 Co 2 SbO 6 . Based on the above considerations, we roughly locate Na 3 Co 2 SbO 6 in the phase diagram as shown in Fig. 1. In this parameter area, the exchange couplings are K −3.6 t 2 /U , J/|K| ∼ −0.14, Γ/|K| ∼ −0.03, and Γ /|K| ∼ 0.16, see Figs. 3(a)-3(d). The large K implies the proximity to the Kitaev-model regime, explaining thereby a strong reduction of the ordered moments from their saturated values [15]. As a crucial test for our theory, we show in Fig. 4 the expected magnon dispersions and spectral weights; neutron scattering experiments on Na 3 Co 2 SbO 6 are desired to verify these predictions. In addition to magnons, a continuum of (spinon) excitations as in RuCl 3 [54,55] is also expected.
If the above picture is confirmed by experiments, the next step should be to drive Na 3 Co 2 SbO 6 into the Kitaev SL regime. As suggested by Fig. 1, this requires a reduction of the trigonal field by ∼ 20 meV, e.g. by means of strain or pressure control. At this point, the relative smallness of SOC for 3d Co ions comes as a great advantage: while strong enough to form the pseudospin moments, it makes the lattice manipulation of the S = 1/2 wavefunctions (and hence magnetism) far easier than in iridates [56]. Monitoring the magnetic behavior of Na 3 Co 2 SbO 6 and other honeycomb cobaltates under uniaxial pressure would be thus very interesting.
To conclude, we have presented a comprehensive theory of exchange interactions in honeycomb cobaltates, and studied their magnetic phase behavior. The analysis of Na 3 Co 2 SbO 6 data suggests that this compound is proximate to a Kitaev SL phase and could be driven there by a c-axis compression. A broader message is that as one goes from 5d Ir to 4d Ru and further to 3d Co, magnetic d orbitals become more localized, and this should improve the conditions for realization of the nearest-neighbor-only interaction model designed by Kitaev. We thank A. Yaresko, T. Takayama, and A. Smerald for useful discussions. G.Kh. acknowledges sup- The d 7 Co 2+ ions in an octahedral crystal field have predominantly t 5 2g e 2 g configuration with a high spin S = 3/2 [1]. A trigonal distortion along Z-axis splits the t 2g manifold into an orbital singlet a 1g and a doublet e g by energy ∆, see Fig. S1(a,b). In the electron representation, it is captured by the Hamiltonian H ∆ = 1 3 ∆(2n a1g − n e g ). In terms of the effective angular momentum L = 1 of the Co 2+ ions, the a 1g -hole configuration corresponds to L Z = 0, while the e g doublet hosts the L Z = ±1 states. Consequently, the trigonal field Hamiltonian translates into H ∆ = ∆(L 2 Z − 2 3 ). The following relations between the L-states and orbitals hold: where shorthand notations a = d yz , b = d zx , and c = d xy are used.
) and H λ = λL · S results in a level structure shown in Fig. S1(c). The states are labeled according to the total angular momentum J eff = 1 2 , 3 2 , and 5 2 . The ground state Kramers doublet hosts a pseudospin S = 1/2; its wavefunctions, written in the basis of |S Z , L Z , read as: The coefficients obey a relation C 1 : . The ground state energy is The exchange Hamiltonian between the pseudospins S = 1/2 is obtained by projecting the Kugel-Khomskii type spin-orbital Hamiltonians onto the ground state doublet (S2). We also specify the excited states, which will be needed in Sec. IV to calculate the magnetic susceptibility. The wavefunctions and energies for 3 2 , ± 1 2 and 5 2 , ± 1 2 states share the same form as of Eq. S2 and Eq. S3, but with different r 1 . Namely, the above equation ∆ λ = r1+3 2 − 3 r1 − 4 r1+2 has three roots with r 1 > 0, 0 > r 1 > −2, and r 1 < −2. The root r 1 > 0 corresponds to the ground state. The other two roots correspond to 3 2 , ± 1 2 and 5 2 , ± 1 2 , respectively. The wavefunctions and energies of the remaining states are: Here, c ϕ = cos ϕ, s ϕ = sin ϕ, and tan 2ϕ = 2 √ 6λ/(2∆ + λ).

II. Pseudospin S = 1/2 Hamiltonian and calculation of its parameters
In the cubic reference frame, pseudospin-1/2 interactions on z-type bonds have a general form The interactions on x and y type bonds are obtained by cyclic permutations among S x j , S y j , and S z j . The Hamiltonian in Eq. S5 can also be written in global XY Z reference frame [3]: with c γ ≡ cos φ γ and s γ ≡ sin φ γ . The angles φ γ = 0, 2π 3 , 4π 3 refer to the z, x, and y type bonds, respectively. The transformations between the two sets of parameters entering Eq. S5 and Eq. S6 are: Since the pseudospin wavefunctions (S2) are defined in the trigonal XY Z basis, it is technically simpler to derive S = 1/2 Hamiltonian in a form of (S6), and then convert the results onto a cubic xyz reference frame via Eqs. S7.
As discussed in the main text, there are three basic exchange channels in d 7 systems, which we consider now in detail. General form of the Kugel-Khomskii type spin-orbital Hamiltonians were obtained earlier [4]; for completeness, they will be reproduced below. Here, the major task is to derive the corresponding pseudospin-1/2 Hamiltonians in a realistic case of finite trigonal splitting of t 2g orbitals. As the S = 1/2 wavefunctions (S2) are somewhat complicated, the calculations are tedious but can still be done analytically.

Intersite U processes
The spin-orbital Hamiltonian for these exchange processes is given by equations (A2) and (3) of Ref. [4]: Here n a = a † a, etc. denote the orbital occupations, t is the hopping between a and b orbitals via ligand ions, t is the direct overlap of c orbitals. The Mott-Hubbard excitation energies are Now, we need to express various combinations of the spin and orbital operators above in terms of the pseudospins S = 1/2 defined by Eq. S2. To this end, we have derived a general projection table, presented in subsection 4 below. Using this table, we obtain the pseudospin Hamiltonian in the form of Eq. S6, with the following exchange constants: Coefficients u i (i = 1, 2, . . . , 7) are given by Eqs. S33 below; they depend on the spatial shape of the pseudospin wavefunctions (S2), and thus decide how the relative values of the pseudospin interactions vary as a function of trigonal field ∆.

Charge-transfer processes
The spin-orbital Hamiltonian is (Eq. 9 in Ref. [4]): where ∆ pd is charge-transfer gap. U p and U p = U p − 2J p H are the intra-and inter-orbital Coulomb repulsion of the ligand p orbitals, respectively, and J p H is the Hund's coupling. Using the projection table of subsection 4, we find the exchange constants in the form of Eq. S6: (S11)

Cyclic exchange processes
The spin-orbital Hamiltonian is (Eq. 11 in Ref. [4]): After projection, we obtain the exchange constants as: (S13) The total contribution from t 2g -t 2g hopping channel to Eq. S6 is given by The corresponding K, J, Γ, and Γ values can be obtained using Eqs. S7.

Intersite U processes
The corresponding spin-orbital exchange Hamiltonian is (Eq. A5 in Ref. [4]): Here, t e = t 2 pdσ /∆ e , with t pdσ representing hopping between p and e g orbitals via the charge-transfer gap ∆ e = ∆ pd +D. Parameter D is the splitting between t 2g and e g levels. The constants α 1 and 1/ U are: After projection onto pseudospin-1/2 doublet (S2), we get the exchange constants in the form of Eq. S6: (S17)

Charge-transfer processes
The spin-orbital Hamiltonian describing these processes is (Eq. 19 in Ref. [4]): where The corresponding pseudospin exchange constants are:

Cyclic exchange processes
The corresponding spin-orbital Hamiltonian is (Eq. 22 in Ref. [4]): with α 4 = 1 − 1 2 D ∆ pd +D . After projection onto pseudospin-1/2 doublet, we obtain: The total contribution from t 2g -e g exchange channel to Eq. S6 is given by The corresponding Hamiltonian is very simple (see Eq. 27 in Ref. [4]): After projecting this onto pseudospin subspace, we find while A 3 = B 3 = 0. The latter implies that e g -e g interaction channel supports the XXZ-type model. Total values of the exchange constants are obtained by summing up t 2g -t 2g , t 2g -e g , and e g -e g contributions [Eqs. S14, S23, and S25, respectively], and converted into K, J, Γ, and Γ using Eqs. S7.

Microscopic parameters used in the calculations
Apart from an overall energy scale t 2 /U , a number of microscopic parameters appeared in the above expressions for exchange constants. Hund's coupling J H ∼ 0.8 eV follows from optical data in CoO [5]; cubic splitting D for 3d ions is of the order of 1.0 − 1.5 eV. With the ab initio estimates of U ∼ 5.0 − 7.8 eV [6][7][8] [9], while U p is about ∼ 4 eV, so we use the representative values of J p H /U p = 0.3 and U p /U = 0.7. We set a direct hopping t = 0.2t (i.e. smaller than in 5d/4d compounds [10]), but this value is nearly irrelevant here since t 2g -t 2g exchange is of minor importance anyway, see Fig. 2 of the main text. A ratio t pdσ /t pdπ = 2 [11] is used. Regarding ∆/λ and U/∆ pd values, we vary them rather broadly, as they most sensitively control the exchange interactions. With the above input parameters, we arrive at K, J, Γ, and Γ values presented in the main text. We have verified that while variations of the input parameters result in some changes of the exchange constants, they do not affect the overall picture and conclusions.

III. Exact diagonalization: Phase diagrams based on static correlations and coherent-state analysis
We consider the nearest-neighbor (NN) interaction model (Eq. 2 of the main text or S5 in the previous section), supplemented by the third-NN Heisenberg exchange J 3 that appears as the major one among the long-range interactions in ab-initio studies [10]. In this section we show the full evolution of the phase diagram with the parameter J 3 and also demonstrate the robustness of our picture with respect to variations of the Hund's exchange J H . The data presented here complements Fig. 1 and Fig. 3(e,f) of the main manuscript.
To determine the magnetic state, we have performed exact diagonalization using the values of exchange parameters derived in Sec. II. Utilizing the Lanczos method, we have obtained exact ground states of the exchange Hamiltonian for a symmetric, hexagon-shaped cluster containing 24 sites. Periodic boundary conditions were applied, corresponding to a periodic tiling of an infinite lattice. Since the small cluster does not allow for spontaneous symmetry breaking, we inspect its magnetic state by analyzing the static spin correlations and by employing the method of coherent spin states introduced in Ref. [12].
We focus on real-space correlations that enable us to judge the extent of the Kitaev spin liquid phase which should be characterized by vanishing correlations beyond nearest neighbors. By evaluating the static spin correlations in momentum space, we would be able to detect the magnetically ordered states that show peaks at the characteristic momenta of the particular ordering pattern. Here, however, it is favorable to utilize the method of coherent spin states that provides a better access to the magnetic order encoded in the complex cluster wavefunction. In essence, it constructs "classical" states (coherent spin states) with spins pointing in prescribed directions and identifies a "classical" state having maximum overlap with the exact cluster ground state. Thanks to its full flexibility in the individual spin directions, the method can precisely determine both collinear patterns as well as non-collinear ones. The "classical" trial state is a product state of spins pointing in prescribed directions (in the sense of finding spin up with 100% probability when measuring in that particular direction) and as such it excludes quantum fluctuations. The maximum overlap is therefore a useful indicator of the amount of quantum fluctuations. For a fluctuation-free state and nondegenerate cluster ground state, the corresponding probability reaches the value 1/(number of degenerate patterns). In contrast, Kitaev spin liquid is highly fluctuating and does not contain a pronounced "classical" state which leads to a tiny maximum overlap (see [12,13] for details). Figures S3 and S4 show phase diagram data as functions of U/∆ pd and ∆/λ for several values of J 3 . The static correlations up to fourth NN presented in upper three rows of panels clearly localize the Kitaev spin liquid phase spreading in the area with dominant K. It is surrounded by several phases with long-range correlations that are identified by the method of coherent spin states. For J 3 = 0, these include two types of FM orders with the magnetic moments lying in the honeycomb plane and perpendicular to it, respectively, stripy phase, zigzag phase zz3, and finally a vortex phase of the type depicted in Fig. S2.
The effect of nonzero antiferromagnetic J 3 may be estimated by considering the correlations of third NN in the individual phases. Strongly supported by J 3 is the zigzag phase that is characterized by AF oriented spins on all third-neighbor bonds. Similarly, a large suppression may be expected for FM and stripy phases that have FM aligned third NN spins. The effect on the vortex phase is weak as each spin has one FM aligned third neighbor and two third neighbors at an angle of 120 • , leading to a cancellation of J 3 in energy on classical level. Finally, in the Kitaev spin liquid phase the third neighbors are not correlated at all, so that small J 3 has a moderate negative impact when trying to align them in AF fashion. The consequences of the above energetics are well visible in Figs. S3 and S4. Once including nonzero J 3 , the Kitaev spin liquid phase slightly grows first, at the expense of FM and stripy phases. At the same time, the Kitaev spin liquid phase is also being expelled from the bottom left corner by the expanding zz3 phase. With increasing J 3 between J 3 = 0.05 and 0.15 in t 2 /U units, two new zigzag phases zz1 and zz2 around Kitaev SL are successively formed. Once J 3 reaches 0.25t 2 /U , the zigzag order quickly takes over, suppressing the Kitaev SL phase completely.
In the large area covered by the zigzag order, various ratios and combinations of signs of the nearest-neighbor interactions are realized. This is the origin of three distinct zigzag phases zz1, zz2, and zz3, differing in their moment directions as seen in bottom panels of Figs. S3 and S4. Negative Γ and positive Γ found in zz1 phase space [see Fig. 3(c,d) of the main text] lead to the ab-plane moment direction. The zz3 phase is characterized by opposite signs of Γ and Γ interactions which stabilizes the zigzag order as in Na 2 IrO 3 [12,14]. Finally, in the zz2 phase, Γ and Γ terms maintain only small values and moment directions pointing along cubic axes x, y, z are selected by order-from-disorder mechanism [12].
To check the robustness of our picture, we have also performed the exact diagonalization for a different J H value. The trends discussed above remain quite similar as demonstrated in Figs. S5 and S6 calculated for J H /U = 0.2. Roughly speaking, when we increase the J H /U value, the whole scenario merely shifts to smaller U/∆ pd region.

IV. Trigonal crystal field ∆ in Na3Co2SbO6
The parameter ∆ determines the effective magnetic moment values µ α eff (α = ab or c), and thus can be obtained from paramagnetic susceptibility χ α (T ). One has to keep in mind that extracting the moments from a standard Curie-Weiss fit χ(T ) = C/(T − Θ) + χ 0 assumes that the excited levels are high in energy (as compared to k B T ) and hence thermally unpopulated. The Curie constant C is then indeed temperature independent, providing the ground state g-factors and moments. For Co 2+ ions, where the excited level at ∼ 30 meV is thermally activated already at the room temperature, we have to use instead a general expression for a single-ion susceptibility: Here, n and m run over all the 12 states (6 doublets in Fig. S1), with the wavefunctions and energies calculated in Sec. I. The partition function Z(T ) = n e −βEn , and β = 1/k B T . M α nm = n|M α |m is matrix element of the magnetic moment operator M = (2S − 3 2 κL) (in units of Bohr magneton µ B ). We use the covalency reduction factor κ = 0.8 typical for Co 2+ ion [1]. χ α ion includes both the Curie and Van-Vleck contributions and depends on two parameters, ∆ and λ.
We have fitted the data of Ref. [15] with χ α (T ) = χ α ion + χ α 0 , and obtained a fair agreement with experiment for both χ ab and χ c , using ∆ = 38 meV and λ = 28 meV, see Fig. S7(a,b). In particular, the characteristic changes in the slopes of both 1/χ ab and 1/χ c data are well reproduced by the calculations. In fact, this behavior is common for layered cobaltates and deserves some discussion.