# Ab initio guided minimal model for the “Kitaev” material BaCo2(AsO4)2:
Importance of direct hopping, third-neighbor exchange and quantum fluctuations

Pavel A. Maksimov  Bogolyubov Laboratory of Theoretical Physics, Joint
Institute for Nuclear Research, Dubna, Moscow region 141980, Russia  M.N.
Miheev Institute of Metal Physics of Ural Branch of Russian Academy of
Sciences, S. Kovalevskaya St. 18, 620990 Ekaterinburg, Russia  Alexey V.
Ushakov  M.N. Miheev Institute of Metal Physics of Ural Branch of Russian
Academy of Sciences, S. Kovalevskaya St. 18, 620990 Ekaterinburg, Russia
Zlata V. Pchelkina  M.N. Miheev Institute of Metal Physics of Ural Branch of
Russian Academy of Sciences, S. Kovalevskaya St. 18, 620990 Ekaterinburg,
Russia  Department of theoretical physics and applied mathematics, Ural
Federal University, Mira St. 19, 620002 Ekaterinburg, Russia  Ying Li
Department of Applied Physics and MOE Key Laboratory for Nonequilibrium
Synthesis and Modulation of Condensed Matter, School of Physics, Xi’an
Jiaotong University, Xi’an 710049, China  Stephen M. Winter  Department of
Physics and Center for Functional Materials, Wake Forest University, NC 27109,
USA  Sergey V. Streltsov<EMAIL_ADDRESS>M.N. Miheev Institute of Metal
Physics of Ural Branch of Russian Academy of Sciences, S. Kovalevskaya St. 18,
620990 Ekaterinburg, Russia  Department of theoretical physics and applied
mathematics, Ural Federal University, Mira St. 19, 620002 Ekaterinburg, Russia

(August 27, 2024)

###### Abstract

By considering two ab initio-based complementary approaches, we analyze the
electronic structure and extract effective spin models of BaCo2(AsO4)2, a
honeycomb material which has been proposed as a candidate for Kitaev physics.
Both methods show that the dominant direct hopping makes the bond-dependent
Kitaev term negligible averting the material away from the sought-after spin-
liquid regime. As a result, we present a simple three-parameter exchange model
to describe the interactions of the lowest doublet of the honeycomb cobaltate
BaCo2(AsO4)2. Remarkably, it is the third-neighbor interactions, both
isotropic and anisotropic, that are responsible for the standout double-zigzag
ground state of BaCo2(AsO4)2, stabilized by quantum fluctuations. A
significantly large third-nearest neighbor hopping, observed in ab initio,
supports the importance of the third-neighbor interactions in the
stabilization of the unique ground state of BaCo2(AsO4)2.

###### pacs:

71.27.+a, 71.20.-b, 71.15.Mb

Introduction.- In the search for material realizations of Kitaev quantum spin
liquids [1], a whole family of honeycomb cobaltates BaCo2(AsO4)2, BaCo2(PO4)2,
CoTiO3 and many others [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20] were proposed to host dominant nearest-neighbor Kitaev
exchange couplings, following the concept that Co2+ ions ($3d^{7}$) in an
octahedral crystal field environment with total spin $S=3/2$ and orbital
angular momentum $l_{\text{eff}}=1$ build spin-orbit coupled $j_{\text{eff}}$=
1/2 doublet states. However, this proposal was recently put into question by a
combination of first-principles-based calculations with single-site exact
diagonalization and two-site perturbation theory [21] where it was found that
the Kitaev term favouring a spin-liquid ground state must be rather small.

In the present work we perform a detailed analysis of the electronic structure
within density functional theory (DFT) employing the projector-augmented wave
[22] basis as implemented in VASP, as well as the full potential local orbital
(FPLO) basis [23]. We further use two ab initio-based complementary approaches
to extract effective spin models for the honeycomb cobaltate BaCo2(AsO4)2. One
is based on total energy calculations of various magnetic configurations
within DFT+U+SOC and mapping into a spin model to extract the exchange tensor
elements, and the second method is based on the recently introduced projED
method [24, 25] by some of the authors which consists of a combination of DFT
calculations, exact diagonalization (ED) of extracted generalized relativistic
Hubbard models on finite clusters, and projection to low-energy spin
Hamiltonians. We find that the Kitaev model is hardly applicable for a
description of the magnetic properties of BaCo2(AsO)2, which is supported by
both ab initio calculations and the phenomenology of BaCo2(AsO)2.

One of the peculiarities of BaCo2(AsO4)2 is the long search for the correct
magnetic ground state. Early works suggested that the ground state is an
unusual long-range ordered spiral state with in-plane ordering vector
$\mathbf{Q}=(0,\pi/3)$ [13]. However, more refined neutron scattering data
later showed that the ground state is instead a collinear double-zigzag state
[15] with the $++--$ pattern of zigzag chains with the magnetic moments
approximately parallel to the chain direction. Since the model including
isotropic exchange between the nearest and 3rd nearest neighbors
($J_{1}$-$J_{3}$ model) does not reproduce this experimentally established
magnetic order, the Hamiltonian requires extra terms beyond the isotropic
third-neighbor interaction in order to describe the ground state of
BaCo2(AsO4)2. Indeed, ab initio studies in this work and in earlier studies
[21] indicate a significant role of anisotropic couplings associated with the
spin-orbital $j_{\rm eff}=1/2$ moments of $3d^{7}$ Co2+. This is further
supported by the finite energy gap observed in inelastic neutron scattering
data and terahertz spectroscopy, although the small magnitude of this gap
($\sim 1.45$ meV [15]) places constraints on the magnitude of anisotropic
terms.

In this paper we extensively study the phase diagram of the general eight-
parameter model to identify which exchange parameters are essential in
stabilization of the double-zigzag state. We find that it is the third-
neighbor isotropic and anisotropic terms that are crucial and make up a
minimal model required to describe the magnetic properties of BaCo2(AsO4)2.
Even though the ab initio parameter sets yield the zigzag ground state, they
are close to the phase boundary. Therefore, we argue that there can be
additional corrections, such as magnetoelastic coupling, which can tune the
calculation towards the experimentally observed double-zigzag ground state.

DFT: crystal-field and hopping parameters.- We start with the analysis of the
on-site Hamiltonian obtained in non-magnetic DFT using both maximally-
localized Wannier functions [26] within VASP [27] and FPLO [23, 28] basis
sets. All calculation details are given in the Supplemental Material (SM) [29]
(see, also, references [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41]
therein).

Figure 1: (a) The band structure and (b) partial densities of states of
BaCo2(AsO4)2 in GGA approximation. Figure 2: Sketch illustrating two
dominating hopping processes according to DFT calculations. a) Direct nearest
neighbor overlap of two $xy$ orbitals. b) Largest contributions for third
nearest neighbors corresponding to hoppings between $x^{2}-y^{2}$ orbitals via
O $2p$ states.

The crystal structure of BaCo2(AsO4)2 is characterized by the $R\bar{3}$ space
group. It consists of stacked honeycomb layers of edge sharing CoO6 octahedra
along the $c$ axis separated by bilayers of opposite-facing AsO4 tetrahedra
with Ba atoms in between. The cubic crystal field of the O6 octahedron around
Co splits the Co $3d$ orbitals onto $e_{g}^{\sigma}$ and $t_{2g}$ states and
the latter are then further split by additional trigonal distortions onto
lower lying $e_{g}^{\pi}$ and $a_{1g}$ states because of the layered
structure, see inset to Fig. 1(b). In Ref. 21, these splittings were estimated
by rescaling values (to fit available experimental data) obtained from N-th
muffin-tin orbital calculations, leading to 735 meV for the
$e_{g}^{\sigma}-a_{1g}$ splitting and to 94 meV for the splitting within
$t_{2g}$ ($a_{1g}-e_{g}^{\pi}$). Our DFT (VASP) calculations give similar
results, 818 meV and 122 meV respectively, without any additional
renormalization. Interestingly, with the FPLO basis, Wannier projection [36]
leads to somewhat larger values of 963 meV and 133 meV. However, it should be
emphasized that in all cases the latter $t_{2g}$ splitting is comparable to
the atomic spin-orbit coupling constant $\lambda_{\rm Co}\approx 60$ meV
defined by $\mathcal{H}_{\rm SOC}=\lambda\ \mathbf{S}\cdot\mathbf{L}$,
suggesting strong deviations from ideal $j_{\rm eff}={1/2}$ moments at the Co
sites. This finding is also consistent with the strong reported anisotropy of
the $g$-tensor ($g_{ab}\sim 2\ g_{c}$) [15].

We utilized several approaches for the nearest-neighbor hopping integrals:
Wannier function projection [26], projected localized orbitals (PLO) [35] as
implemented in VASP, and FPLO [36]. Notably, all approaches we used led to the
same conclusions. For instance, for the bond depicted in Fig. 2(a), the direct
$xy/xy$ hopping $t^{\prime}\sim-300$ meV dominates over hoppings via $xz$,
$yz$ orbitals, $t\sim$ 50 meV, which are associated with electron transfer via
ligand $p$ states, local coordinate system shown in Fig. 2(b) is chosen.
Similar findings were recently reported for Na2BaCo(PO4)2 [42] and other edge-
sharing Co compounds [43]. Already this fact reveals that one of the key
assumptions used in [2], $\kappa=|t^{\prime}/t|<1$, is far from being
fulfilled. This results in a small Kitaev exchange for the first-nearest
neighbors as we will demonstrate in the next section. It is worth mentioning
that a large $xy/xy$ hopping not only affects exchange coupling, but it also
strongly changes the electronic structure. One can clearly observe formation
of two isolated branches of Co $t_{2g}$ bands in Fig. 1, where each of
$t_{2g}$ orbitals has a direct overlap with the corresponding orbital on one
of the neighboring Co sites. Indeed, the gap between these bands closes, if
the direct $xy/xy$ hopping is put to zero as explained in SM [29].

Another very important feature, which becomes evident already on DFT level, is
that hopping between third-nearest neighbors, $t_{3}$, is not small. The most
important contribution comes from hopping between $e_{g}$ orbitals, which
strongly hybridize with ligand $p$ orbitals. This leads, for example, to an
effective hopping $t_{3}^{x^{2}-y^{2}/x^{2}-y^{2}}\approx 124$ meV, see Fig.
2(b). This hopping is associated with an antiferromagnetic $\text{J}_{3}$ that
may be comparable to the nearest neighbor exchange $\text{J}_{1}$. This
finding is in contradiction with the assumption [2, 44, 3] that longer range
couplings should be suppressed in $3d^{7}$ compounds as a result of stronger
Coulomb interactions and smaller $d$ orbitals localizing the moments in
comparison to traditional $4d^{5}$ and $5d^{5}$ Kitaev candidate materials. In
fact, the partial filling of the $e_{g}$ orbitals in $3d^{7}$ Co2+ provides
additional long-range exchange pathways. Large hopping between the third-
nearest neighbors is also reflected in electronic structure of BaCo2(AsO4)2.
It results in formation of the bonding and antibonding $e_{g}$ states, as
discussed in SM [29], which are clearly seen in Fig. 1.

Table 1: The dependence of exchange interaction parameters on on-site Coulomb $U$ computed from DFT+SOC+U total energy in the extended Kitaev model. Intra-atomic Hund’s exchange was chosen to be $J_{H}$=0.9 eV. $U$ | 5 eV | 6 eV | 7 eV
---|---|---|---
$\text{J}_{1}$ (K) | -61.0 | -40.9 | -37.6
$K_{1}$ (K) | 0.3 | 2.2 | 5.3
$\Gamma_{1}$ (K) | -2.2 | -1.7 | -1.8
$\Gamma_{1}^{\prime}$ (K) | 5.1 | 4.0 | 3.2
$\text{J}_{3}$ (K) | 31.4 | 24.6 | 18.7
$K_{3}$ (K) | -0.2 | 0.2 | -0.2
$\Gamma_{3}$ (K) | -4.5 | -6.0 | -4.5
$\Gamma^{\prime}_{3}$ (K) | -3.6 | -2.3 | -1.8
Table 2: The dependence of exchange interaction parameters on on-site Coulomb $U$ in the crystallographic parameterization. Intra-atomic Hund’s exchange was chosen to be $J_{H}$=0.9 eV. $U$ | 5 eV | 6 eV | 7 eV
---|---|---|---
$J_{1}$ (K) | -63.6 | -42.3 | -37.4
$\Delta_{1}$ | 0.87 | 0.85 | 0.88
$J_{\pm\pm}^{(1)}$ (K) | 2.4 | 1.5 | 0.8
$J_{z\pm}^{(1)}$ (K) | 3.5 | 3.7 | 4.9
$J_{3}$ (K) | 35.2 | 28.2 | 21.3
$\Delta_{3}$ | 0.67 | 0.62 | 0.62
$J_{\pm\pm}^{(3)}$ (K) | 0.3 | 1.2 | 0.9
$J_{z\pm}^{(3)}$ (K) | 0.3 | 1.8 | 1.2

Ab initio Exchange Parameters.- Due to the interplay of crystal field and SOC,
the magnetic interactions of effective doublets are defined by the symmetry of
the lattice. The symmetry of edge-sharing CoO6 octahedra results in the
extended Kitaev-Heisenberg model [45], which has been discussed for
BaCo2(AsO4)2 [3]. Explicitly the exchange Hamiltonian is given by

$\displaystyle\hat{\cal H}_{\text{cubic}}=\sum_{\langle ij\rangle_{n}}$
$\displaystyle\text{J}_{n}\
\mathbf{S}_{i}\cdot\mathbf{S}_{j}+K_{n}S^{\gamma}_{i}S^{\gamma}_{j}+\Gamma_{n}\left(S^{\alpha}_{i}S^{\beta}_{j}+S^{\beta}_{i}S^{\alpha}_{j}\right)$
$\displaystyle+$
$\displaystyle\Gamma^{\prime}_{n}\left(S^{\gamma}_{i}S^{\alpha}_{j}+S^{\gamma}_{i}S^{\beta}_{j}+S^{\alpha}_{i}S^{\gamma}_{j}+S^{\beta}_{i}S^{\gamma}_{j}\right),$
(1)

where the sum is taken over the three types of bonds of the honeycomb lattice,
$\\{\alpha,\beta,\gamma\\}=\\{\text{x,y,z}\\}$ for the Z-type bond and
interactions on the X and Y bonds are obtained through a cyclic permutation
[46, 45], see Fig. 3(a). Note that this model uses cubic axes {x,y,z}, which
are shown in Figs. 2 and 3(a) and related to the ion-ligand bonds. In
accordance with neutron studies, and ab initio results below, we anticipate
that first-neighbor ($n=1$) and third-neighbor ($n=3$) couplings are dominant.

For our purpose, it is also convenient to refer the interactions alternatively
to the crystallographic axes $\\{x,y,z\\}$ which are defined by the honeycomb
plane of magnetic ions, shown in Fig. 3(a). The Hamiltonian in that reference
frame is given by

$\displaystyle\hat{\cal H}_{\text{cryst}}=$ $\displaystyle\sum_{\langle
ij\rangle_{n}}J_{n}\Big{(}S^{x}_{i}S^{x}_{j}+S^{y}_{i}S^{y}_{j}+\Delta_{n}S^{z}_{i}S^{z}_{j}\Big{)}$
$\displaystyle-$ $\displaystyle
2J_{\pm\pm}^{(n)}\Big{(}\Big{(}S^{x}_{i}S^{x}_{j}-S^{y}_{i}S^{y}_{j}\Big{)}c_{\alpha}-\Big{(}S^{x}_{i}S^{y}_{j}+S^{y}_{i}S^{x}_{j}\Big{)}s_{\alpha}\Big{)}$
$\displaystyle-$ $\displaystyle
J_{z\pm}^{(n)}\Big{(}\Big{(}S^{x}_{i}S^{z}_{j}+S^{z}_{i}S^{x}_{j}\Big{)}c_{\alpha}+\Big{(}S^{y}_{i}S^{z}_{j}+S^{z}_{i}S^{y}_{j}\Big{)}s_{\alpha}\Big{)},$
(2)

where $c_{\alpha}\equiv\cos\varphi_{\alpha}$ and
$s_{\alpha}\equiv\sin\varphi_{\alpha}$ with the bond-dependent phases
$\varphi_{\alpha}\\!=\\!\\{0,2\pi/3,-2\pi/3\\}$ for the three types of first
and third-neighbor bonds [45, 47]. The exchange parameters of the extended
Kitaev-Heisenberg model (1) and anisotropic-exchange model in the
crystallographic axes (2) are related through a simple linear transformation:

$\displaystyle\text{J}_{n}$
$\displaystyle=\frac{1}{3}\left(2J_{n}+\Delta_{n}J_{n}+2J_{\pm\pm}^{(n)}-\sqrt{2}J_{z\pm}^{(n)}\right),$
$\displaystyle K_{n}$
$\displaystyle=-2J_{\pm\pm}^{(n)}+\sqrt{2}J_{z\pm}^{(n)},$ (3)
$\displaystyle\Gamma_{n}$
$\displaystyle=\frac{1}{3}\left(-J_{n}+\Delta_{n}J_{n}-4J_{\pm\pm}^{(n)}-\sqrt{2}J_{z\pm}^{(n)}\right),$
$\displaystyle\Gamma^{\prime}_{n}$
$\displaystyle=\frac{1}{6}\left(-2J_{n}+2\Delta_{n}J_{n}+4J_{\pm\pm}^{(n)}+\sqrt{2}J_{z\pm}^{(n)}\right).$

We use a variety of approaches to estimate the magnetic exchange couplings.
Conventional DFT calculations underestimate the effect of strong Coulomb
correlations, which must be taken into account for extraction of the exchange
interaction in the transition metal oxides. We therefore first computed the
magnetic interactions based on DFT(GGA)+U+SOC calculation of total energies of
four non-colinear magnetic configurations [41]. These results are complemented
by extraction of the couplings from exact diagonalization of two-site
clusters, described below.

For the total energy DFT(GGA)+U+SOC calculations, the most important exchanges
for the first and third-nearest neighbors are summarized in Tables 1 and 2 in
terms of both extended Kitaev and crystallographic parameterizations (other
constants were found to be small, e.g. $J_{2}\sim 0.2$K). The calculations
were performed for several values of the on-site Hubbard repulsion parameter
$U$, but all of them demonstrate that (i) the Kitaev exchange is small both
for the first and third-nearest neighbors and that (ii) there is a strong
exchange coupling with third-nearest neighbors. Both factors strongly suppress
formation of a spin-liquid state, and are compatible with previous neutron
scattering analysis, which suggested $J_{1}\sim-38$ K, $J_{3}\sim+10$ K [14],
as well as more recent data [48], which estimated $J_{1}\sim-88$ K, and
$J_{3}\sim+29$ K, and $J_{\pm\pm}^{(1)}\sim-0.6$ K. We note, however, that
large XXZ anisotropy estimated from experiment ($\Delta\sim 0.37$ [14] and
$\Delta_{1}\sim 0.16$ [48]), is not reproduced in this approach.

In order to further examine the magnetic couplings, we employed a
complementary approach similar to Ref. 21: exact diagonalization of the five
$d$-orbital model on two sites [24, 25]. For this purpose, we employ hopping
integrals obtained from VASP, and take the fully spherically symmetric form
[49] of the on-site Coulomb interactions with Slater parameters
$F_{4}/F_{2}=0.625$, and $F_{0}\equiv U$ and $F_{2}\equiv 14J/(1+0.625)$ set
according to $U=5$ to 7 eV, and $J_{H,t_{2g}}=0.9$ eV. The results are shown
in Table 3. For the purpose of comparison, results for Coulomb parameters
equivalent to Ref. 21 ($U=3.25$ eV, $J_{H,t_{2g}}=0.7$ eV) are also shown. We
note that this approach neglects an important contribution to the exchange
involving multiple holes on a given ligand, which can be corrected using
expressions from perturbation theory [43, 50, 3]. This leads to shifts of the
nearest neighbor couplings $\text{J}_{1}\to\text{J}_{1}+\delta\text{J}$, and
$\Gamma_{1}\to\Gamma_{1}+\delta\Gamma$ and
$\Gamma_{1}^{\prime}\to\Gamma_{1}^{\prime}+\delta\Gamma$ that can be estimated
from the trigonal crystal field splitting and realistic metal-ligand hopping
parameters (see [43, 3]). In this case, a rough estimate is $\delta J\sim-20$
K, and $\delta\Gamma\sim+7$ K. Both corrected and uncorrected results are
given in Tables 3 and 4.

Table 3: Exchange interaction parameters (in units of K) computed from exact diagonalization of effective $d$-orbital model. Values in brackets include estimated corrections for omitted ligand exchange processes. $J_{H,t_{2g}}$ | 0.7 eV | 0.9 eV
---|---|---
$U$ | 3.25 eV | 5 eV | 6 eV | 7 eV
$\text{J}_{1}$ (K) | -107 (-127) | -37 (-57) | -18 (-38) | -8.8 (-29)
$K_{1}$ (K) | 32 | 13 | 6.5 | 3.4
$\Gamma_{1}$ (K) | 28 (35) | 14 (21) | 8.0 (15) | 4.8 (12)
$\Gamma_{1}^{\prime}$ (K) | 9.4 (16) | 7 (14) | 4.0 (11) | 2.4 (9)
$\text{J}_{3}$ (K) | 43 | 30 | 27 | 24
$K_{3}$ (K) | -0.6 | -0.4 | -0.3 | -0.3
$\Gamma_{3}$ (K) | -20 | -12 | -10 | -8.9
$\Gamma_{3}^{\prime}$ (K) | -21 | -12 | -11 | -9.2
Table 4: Exchange interaction parameters computed from exact diagonalization of effective $d$-orbital model. Values in brackets include estimated corrections for omitted ligand exchange processes. $J_{H,t_{2g}}$ | 0.7 eV | 0.9 eV
---|---|---
$U$ | 3.25 eV | 5 eV | 6 eV | 7 eV
$J_{1}$ (K) | -113 (-140) | -42 (-69) | -21 (-48) | -11 (-38)
$\Delta_{1}$ | 0.58 (0.51) | 0.36 (0.31) | 0.24 (0.23) | 0.11 (0.19)
$J_{\pm\pm}^{(1)}$ (K) | -12 | -4.5 | -2.4 | -1.4
$J_{z\pm}^{(1)}$ (K) | 6.2 | 2.7 | 1.2 | 0.4
$J_{3}$ (K) | 64 | 42 | 37 | 33
$\Delta_{3}$ | $\sim 0$ | 0.14 | 0.16 | 0.18
$J_{\pm\pm}^{(3)}$ (K) | 0.1 | $\sim 0$ | $\sim 0$ | $\sim 0$
$J_{z\pm}^{(3)}$ (K) | -0.60 | -0.36 | -0.30 | -0.26

The ED + perturbation theory results are essentially compatible with the DFT
results above; ferromagnetic $J_{1}$ dominates the couplings, while the
anisotropic couplings $K_{1},\Gamma_{1},\Gamma_{1}^{\prime}>0$ are all of
smaller but similar magnitude. The main difference is that the XXZ anisotropy
is considerably stronger in ED results, with corrected values of $\Delta_{1}$
ranging between $0.2$ and $0.5$. This anisotropy originates from the effects
of local trigonal crystal field on the $j_{\rm eff}=1/2$ multiplet structure,
which may not be completely captured in one-electron methods such as
DFT+U+SOC. On the other hand, cluster approaches, such as our ED method, tend
to have difficulties with long-range interactions. This fact is embodied by
surprisingly large $J_{3}$, which is comparable to $J_{1}$, as can be seen in
Table 4.

Thus, both calculation methods rather guide than provide exact estimates of
exchange constants. It is important that our methods and perturbation theory
of Ref. [21] give consistent results, which is not always the case for Kitaev
materials (for Li2IrO3 see Ref. [51, 52] and for $\alpha$-RuCl3 see Ref.
[53]). However, both approaches place the system in a region of zigzag order
for all of the presented parameter sets, rather than the double-zigzag
observed in BaCo2(AsO4)2. In the next section we show that these parameter
sets are, in fact, close to the boundary of the double-zigzag phase, and
extensively study the full phase diagram in order to establish which of the
parameters of the model would tune the ground state towards the experimentally
established one.

Figure 3: (a) Honeycomb lattice with three types of bonds; cubic axes {x,y,z}
and crystallographic axes $\\{x,y,z\\}$ are also presented. (b) Brillouin zone
of the honeycomb lattice with high-symmetry points indicated. (c)-(f)
Classical phase diagram of the model (2) calculated with Luttinger-Tisza
method for $J_{1}<0$, $J_{3}>0$, $\Delta=\Delta_{1}=\Delta_{3}$. The ordering
vector of the classical spiral state $Q_{y}$ is represented with the color
intensity map. A sketch of the zigzag shown as an inset. The Spiral-II state
has non-zero $Q_{x}$ component of the ordering vector, thus not related to the
ground state of BaCo2(AsO4)2. Figure 4: Luttinger-Tisza
$J_{3}$-$J_{z\pm}^{(3)}$ phase diagrams for $U=5$ eV sets of parameters from
Tables 1 and 3. DFT and ED parameter sets are indicated with yellow dots.
Figure 5: (a) Luttinger-Tisza phase diagram of the $XY$ (i.e. $\Delta=0$)
$J_{1}$-$J_{3}$-$J_{z\pm}^{(3)}$ model (4) calculated with Luttinger-Tisza
method for $J_{1}<0$, $J_{3}>0$. The ordering vector of the classical state
$Q_{y}$ is represented with the color intensity map. (b) Classical phase
diagram of the model (4) with the same notation for the intensity map. (c)
DMRG phase diagram of the model (4) for $S=1/2$. The gray areas indicate
regions of the intermediate phases which are beyond the scope of this work.
(d) Example of DMRG calculation for a representative parameter set
$J_{3}=J_{z\pm}^{(3)}=0.5|J|$ which exhibits a double-zigzag ground state.

Phase Diagram.- We now investigate the minimal set of interactions that is
compatible with the ground state phenomenology of BaCo2(AsO4)2 [16]. We start
with the isotropic $XY$ $J_{1}$-$J_{3}$ model, guided by ab initio sets
presented above, and study phase transitions when the anisotropic-exchange
terms of the full eight-parameter model (2) are added to the Hamiltonian one
at a time, using the classical Luttinger-Tisza (LT) method [54]. We should
note that in frustrated magnets with broken continuous symmetry,
incommensurate spiral states given by Luttinger-Tisza approach cannot satisfy
local spin constraint, thus are not exact classical ground states [55].
However, even though LT method breaks down, it can still point to regions of
exotic ground states on the phase diagram [56, 46, 57]. Therefore, we use LT
method in this section as guide for the search of classical ground states
beyond ferromagnetic and zigzag phases by plotting the ordering vector
selected by LT, and we take the boundary of stability of LT method as an
indicator of the regions of the phase diagram which can host the double-zigzag
state (which is studied in detail in the next section).

Corresponding phase diagrams are shown in Fig. 3 where we plot the ordering
vector $(0,Q_{y})$ as an intensity map, and the regions where LT method breaks
down are shown with white dashed lines. As depicted in Brillouin zone in Fig.
3(b), $Q_{y}=0$ corresponds to a ferromagnetic (FM) state, $Q_{y}=2\pi/3$ to a
zigzag state, and the spiral state originally proposed for BaCo2(AsO4)2
interpolates between these states with $0<Q_{y}<2\pi/3$. The $J_{1}$-$J_{3}$
model is known to host this spiral state for $0.25|J_{1}|<J_{3}<0.4|J_{1}|$
[58], as shown in Fig. 3(c), and the addition of XXZ anisotropy (with $\Delta$
= $\Delta_{1}$=$\Delta_{3}$) does not alter the relative stabilities of the
phases compared to the pure Heisenberg $J_{1}$-$J_{3}$ model.

In contrast, the addition of anisotropic first-neighbor $J_{\pm\pm}^{(1)}$ and
$J_{z\pm}^{(1)}$ exchange terms stabilizes zigzag state, see Fig. 3(d-e).
Whereas, the $J_{\pm\pm}^{(3)}$ interaction promotes non-FM and non-zigzag
states, shown in gray in Fig. 3(f). However, LT suggests states with the
ordering vector near the K point with non-zero $Q_{x}$ component, thus this
region is not relevant to the ground state of BaCo2(AsO4)2 and we do not go
into the specifics of its structure.

Finally, we also considered phase diagrams starting from the DFT and ED
parameters. For this purpose, we use the results for $U=5$ eV in each case, as
the corresponding interactions are closest in magnitude to the experimental
estimates [14, 48]. Results are shown in Fig. 4 with respect to tuning $J_{3}$
and $J_{z\pm}^{(3)}$. The latter 3rd neighbor interaction was selected because
it tends to stabilize the $(0,Q_{y})$ spiral, when added to the both the
isotropic $J_{1}-J_{3}$ and ED models. As can be seen, both ab initio
parameter sets predict zigzag as a ground state, but both lie relatively close
to the boundary of zigzag state. While we find that various modifications of
these parameters may stabilize the spiral phase (such as increasing
$J_{z\pm}^{(3)}$ with respect to the ED couplings), the discrepancy is most
likely attributable to an overestimation of third neighbor couplings. That is,
rescaling of $J_{3}$ shifts both parameter sets towards the spiral state. As
we show in the next section, the spiral is indeed replaced by the unique
double-zigzag state in a more careful study of the classical ground state and
when quantum fluctuations are taken into consideration.

Double-zigzag state.- In order to search for the double-zigzag state and
explore quantum effects on the relevant phases, we focus on a reduced three-
parameter $J_{1}$-$J_{3}$-$J_{z\pm}^{(3)}$ model as a minimal Hamiltonian for
BaCo2(AsO4)2:

$\displaystyle\hat{\cal H}_{\text{min}}=$ $\displaystyle\sum_{\langle
ij\rangle_{1}}J_{1}\Big{(}S^{x}_{i}S^{x}_{j}+S^{y}_{i}S^{y}_{j}\Big{)}+J_{3}\sum_{\langle
ij\rangle_{3}}\Big{(}S^{x}_{i}S^{x}_{j}+S^{y}_{i}S^{y}_{j}\Big{)}$
$\displaystyle-$ $\displaystyle
J_{z\pm}^{(3)}\Big{(}\Big{(}S^{x}_{i}S^{z}_{j}+S^{z}_{i}S^{x}_{j}\Big{)}c_{\alpha}+\Big{(}S^{y}_{i}S^{z}_{j}+S^{z}_{i}S^{y}_{j}\Big{)}s_{\alpha}\Big{)},$
(4)

We take easy-plane anisotropy $\Delta_{1}=\Delta_{3}=0$. This choice of
minimal interactions is sufficient to provide a wide region of stability for
the double-zigzag phase, although it neglects many interactions relevant to
the real material.

First, the LT phase diagram of the model (4) presented in Fig. 5(a), points to
a state with the $(0,Q_{y})$ ordering vector. However, since the
incommensurate spiral state given by LT does not satisfy the strong local spin
length constraint, we perform a more focused study of the classical ground
state. We use a quasi-one-dimensional cluster, periodic in $x$-axis, using the
fact that spin modulation is only along the $y$-axis, as suggested by the
$(0,Q_{y})$-spiral state from LT calculations. The results of classical energy
minimization on a 24-unit cell cluster are shown in Fig. 5(b), where color
intensity illustrates the magnitude of $Q_{y}$ ordering vector. One can see
that, compared to LT method, double-zigzag state is in fact stabilized in the
phase diagram of the $J_{1}$-$J_{3}$-$J_{z\pm}^{(3)}$ model for
$J_{z\pm}^{(3)}\gtrsim 0.2$, while the spiral state is stable for
$J_{z\pm}^{(3)}\lesssim 0.2$.

Moreover, it is known that there are strong renormalizations of phase diagrams
of frustrated quantum $S=1/2$ models relative to classical
$S\rightarrow\infty$ approximation [59, 46, 60, 61, 62, 63, 64, 65, 66, 67,
68]. In order to study the model (4) in the quantum limit in the context of
$j_{\text{eff}}=1/2$ moments of BaCo2(AsO4)2, we employ density matrix
renormalization group (DMRG) [69] using ITensor library [70] on a 192-site
$S=1/2$ cluster with open boundary conditions using 20 sweeps with error
$<10^{-4}$ and a random initial state. (We have also studied clusters of other
sizes and with periodic boundary conditions, they all yielded very similar
results.) The phases were identified by the maximum value of spin-spin
correlator $\mathcal{S}(\mathbf{k})$ calculated at $\Gamma$, $K$, $M$,
$(0,\pi/3)$ where

$\displaystyle\mathcal{S}(\mathbf{k})=\sum_{i,j}\langle\mathbf{S}_{i}\cdot\mathbf{S}_{j}\rangle
e^{i\mathbf{k}\left(\mathbf{r}_{i}-\mathbf{r}_{j}\right)}.$ (5)

The phase diagram, obtained with DMRG in the quantum $S=1/2$ limit, is shown
in Fig. 5(c). The ordering vector is shown with the color intensity, same as
in Figs. 5(a) and 5(b), the phases in gray are intermediate between FM,
double-zigzag and zigzag but their characterization is beyond the scope of
this work. Note that our investigative DMRG calculation is unable to give a
conclusive result in the region of multiple phase competition for
$J_{z\pm}^{(3)}<0.1|J_{1}|$. Nonetheless, one can see that the double-zigzag
state is stable in a wider region of the phase diagram for $S=1/2$, relative
to classical model predictions. This fact implies that quantum fluctuations
play significant role, which is captured in DMRG. This mechanism is generic,
and applies beyond the minimal model. Such fluctuations are known to stabilize
collinear orders in frustrated systems, such as the field-induced up-up-down
state in the triangular lattice antiferromagnet [71, 72], honeycomb
$J_{1}$-$J_{2}$ model [73], and anisotropic-exchange model on a triangular
lattice [74, 59, 57].

An example of the spin orientations obtained by DMRG for the representative
parameter set $J_{3}=J_{z\pm}^{(3)}=0.5|J_{1}|$ is shown in Fig. 5(c). This
observed spin structure is precisely the same as measured in the latest
neutron data [15]; the $++--$ double-zigzag structure. Moreover, an out-of-
plane canting of the spins around $5^{\circ}$ was also reported [15]. We also
observe the out-of-plane canting, induced by the anisotropic $J_{z\pm}^{(3)}$
term, which couples in-plane and out-of-plane spin components. However, the
canting in our DMRG calculation is not only of opposite sign between the
chains of opposite direction but also has different sign between the A and B
sublattices of the honeycomb lattice.

Finally, we remark on an additional mechanism that may stabilize the
mysterious double-zigzag structure. Our DFT+U+SOC calculations ($U=6$ eV) show
that zigzag is the ground state magnetic structure with the double-zigzag
being 1.2 meV/f.u. higher in energy for the experimental crystal structure.
However, relaxation of the atomic positions completely changes the situation:
the double-zigzag order becomes more stable than the spiral by 0.2 meV/f.u.
($\approx 1$K/Co). The task of deciphering the origin of the stabilization of
the double-zigzag structure from ab initio is extremely complicated due to a
tiny total energy difference (corresponding hopping parameters and nearest
neighbor exchange constants can be found in SM [29]), but one can conclude
with certainty that (i) the system is on the border between two magnetic
phases and phase separation or presence of different domains are not excluded;
(ii) the magneto-elastic coupling is important in BaCo2(AsO4)2 [75].

Conclusions.- By means of ab initio band structure, Luttinger-Tisza and DMRG
calculations we studied the electronic and magnetic properties of BaCo2(AsO4)2
\- a candidate material for the realization of the celebrated Kitaev model.
While previous theoretical results [3] and experimental data [20] suggested
dominant Kitaev interaction, promising proximity to the spin-liquid regime, in
this paper we showed that this notion is not supported neither by
phenomenology, nor by various ab initio methods.

Instead, DFT and ED in this paper establish that large direct exchange due to
$t_{2g}$ orbitals strongly suppresses anisotropic contributions to the
exchange interaction between nearest neighbors in agreement with Ref. 21.
Moreover, there is also substantial coupling with third-nearest neighbors,
such that $\text{J}_{3}>0.3|\text{J}_{1}|$. These two findings together make
the formation a spin-liquid state unfavorable, driving the system towards the
long-range ordered state.

However, proposed ab initio models do not yield the unique double-zigzag
ground state of BaCo2(AsO4)2. Through the extensive search over the eight-
parameter phase space we establish a minimal model that hosts double-zigzag
state in a wide range of parameters. Remarkably, quantum fluctuations inherent
to $j_{\text{eff}}=1/2$ magnets strongly affect the ground state of the
proposed model and stabilize the double-zigzag magnetic structure previously
observed experimentally. We show that proposed parameter sets from ab initio
are near the boundary to double-zigzag state and we suggest that
magnetoelastic coupling can play a crucial role since double-zigzag state can
be stabilized by optimization of the crystal structure in DFT+U+SOC
calculations.

Acknowledgments. We cannot describe the immeasurable gratitude towards Roser
Valentí, whose expertise immensely assisted in the launch of this project, and
whose numerous insightful discussions aided in guiding the paper in the
desirable direction. We would like to thank Sasha Chernyshev and Arun
Paramekanti for useful discussions and comments. S.S. is grateful to B. Cava
and I. Solovyev for fruitful communications. We acknowledge support of the
Russian Science Foundation via project 20-62-46047.

## References

  * Kitaev [2006] A. Kitaev, Anyons in an exactly solved model and beyond, Annals of Physics 321, 2 (2006).
  * Liu and Khaliullin [2018] H. Liu and G. Khaliullin, Pseudospin exchange interactions in ${d}^{7}$ cobalt compounds: Possible realization of the kitaev model, Phys. Rev. B 97, 014407 (2018).
  * Liu  _et al._ [2020] H. Liu, J. c. v. Chaloupka, and G. Khaliullin, Kitaev spin liquid in $3d$ transition metal compounds, Phys. Rev. Lett. 125, 047201 (2020).
  * Kim  _et al._ [2021] C. Kim, J. Jeong, G. Lin, P. Park, T. Masuda, S. Asai, S. Itoh, H.-S. Kim, H. Zhou, J. Ma, and J.-G. Park, Antiferromagnetic kitaev interaction inJ eff = 1/2 cobalt honeycomb materials na3co2sbo6 and na2co2teo6, Journal of Physics: Condensed Matter 34, 045802 (2021).
  * Lin  _et al._ [2021] G. Lin, J. Jeong, C. Kim, Y. Wang, Q. Huang, T. Masuda, S. Asai, S. Itoh, G. Günther, M. Russina, Z. Lu, J. Sheng, L. Wang, J. Wang, G. Wang, Q. Ren, C. Xi, W. Tong, L. Ling, Z. Liu, L. Wu, J. Mei, Z. Qu, H. Zhou, X. Wang, J.-G. Park, Y. Wan, and J. Ma, Field-induced quantum spin disordered state in spin-1/2 honeycomb magnet na2co2teo6, Nature Communications 12, 5559 (2021).
  * Songvilay  _et al._ [2020] M. Songvilay, J. Robert, S. Petit, J. A. Rodriguez-Rivera, W. D. Ratcliff, F. Damay, V. Balédent, M. Jiménez-Ruiz, P. Lejay, E. Pachoud, A. Hadj-Azzem, V. Simonet, and C. Stock, Kitaev interactions in the co honeycomb antiferromagnets ${\mathrm{na}}_{3}{\mathrm{co}}_{2}{\mathrm{sbo}}_{6}$ and ${\mathrm{na}}_{2}{\mathrm{co}}_{2}{\mathrm{teo}}_{6}$, Phys. Rev. B 102, 224429 (2020).
  * Hong  _et al._ [2021] X. Hong, M. Gillig, R. Hentrich, W. Yao, V. Kocsis, A. R. Witte, T. Schreiner, D. Baumann, N. Pérez, A. U. Wolter, _et al._ , Unusual heat transport of the kitaev material na $\\_2$ co $\\_2$ teo $\\_6$: putative quantum spin liquid and low-energy spin excitations, arXiv preprint arXiv:2101.12199 (2021).
  * Wildes  _et al._ [2017] A. R. Wildes, V. Simonet, E. Ressouche, R. Ballou, and G. J. McIntyre, The magnetic properties and structure of the quasi-two-dimensional antiferromagnet CoPS3, Journal of Physics: Condensed Matter 29, 455801 (2017).
  * Lefrançois  _et al._ [2016] E. Lefrançois, M. Songvilay, J. Robert, G. Nataf, E. Jordan, L. Chaix, C. V. Colin, P. Lejay, A. Hadj-Azzem, R. Ballou, and V. Simonet, Magnetic properties of the honeycomb oxide ${\mathbf{na}}_{2}{\mathbf{co}}_{2}{\mathbf{teo}}_{6}$, Phys. Rev. B 94, 214416 (2016).
  * Wong  _et al._ [2016] C. Wong, M. Avdeev, and C. D. Ling, Zig-zag magnetic ordering in honeycomb-layered na3co2sbo6, Journal of Solid State Chemistry 243, 18 (2016).
  * Zvereva  _et al._ [2016] E. A. Zvereva, M. I. Stratan, A. V. Ushakov, V. B. Nalbandyan, I. L. Shukaev, A. V. Silhanek, S. V. Streltsov, and A. N. Vasiliev, Orbitally induced hierarchy of exchange interactions in zigzag antiferromagnetic state of honeycomb silver delafossite ag3co2sbo6, Dalton Transactions 45, 7373 (2016).
  * Eymond  _et al._ [1969] S. Eymond, M. C. Martin, and A. Durif, Donnees cristallographiques sur quelques composes isomorphes du monoarseniate de baryum-nickel : Bani2(aso4)2, Materials Research Bulletin 4, 595 (1969).
  * Regnault  _et al._ [1977] L. Regnault, P. Burlet, and J. Rossat-Mignod, Magnetic ordering in a planar x - y model: Baco2(aso4)2, Physica B+C 86-88, 660 (1977).
  * Regnault  _et al._ [1986] L. Regnault, J. Boucher, J. Rossat-Mignod, J. Bouillot, R. Pynn, J. Henry, and J. Renard, Nonlinear excitations in ld and 2d magnetic systems, Physica B+C 136, 329 (1986), neutron Scattering.
  * Regnault  _et al._ [2018] L.-P. Regnault, C. Boullier, and J. E. Lorenzo, Polarized-neutron investigation of magnetic ordering and spin dynamics in baco¡sub¿2¡/sub¿(aso¡sub¿4¡/sub¿)¡sub¿2¡/sub¿ frustrated honeycomb-lattice magnet, Heliyon 4, 10.1016/j.heliyon.2018.e00507 (2018).
  * Regnault and Rossat-Mignod [1990] L. Regnault and J. Rossat-Mignod, Phase transitions in quasi two-dimensional planar magnets, in  _Magnetic Properties of Layered Transition Metal Compounds_ (Springer, 1990) pp. 271–321.
  * Regnault  _et al._ [2006] L. Regnault, C. Boullier, and J. Henry, Investigation by spherical neutron polarimetry of magnetic properties in baco2(aso4)2, Physica B: Condensed Matter 385-386, 425 (2006).
  * Zhong  _et al._ [2020] R. Zhong, T. Gao, N. P. Ong, and R. J. Cava, Weak-Field Induced Nonmagnetic State in a Co-Based Honeycomb, Sci. Adv. 6, eaay6953 (2020).
  * Shi  _et al._ [2021] L. Y. Shi, X. M. Wang, R. D. Zhong, Z. X. Wang, T. C. Hu, S. J. Zhang, Q. M. Liu, T. Dong, F. Wang, and N. L. Wang, Magnetic excitations of the field-induced states in ${\mathrm{baco}}_{2}{({\mathrm{AsO}}_{4})}_{2}$ probed by time-domain terahertz spectroscopy, Phys. Rev. B 104, 144408 (2021).
  * Zhang  _et al._ [2021] X. Zhang, Y. Xu, R. Zhong, R. Cava, N. Drichko, and N. Armitage, In-and out-of-plane field induced quantum spin-liquid states in a more ideal kitaev material: Baco $\\_2$(aso $\\_4$) $\\_2$, arXiv preprint arXiv:2106.13418 (2021).
  * Das  _et al._ [2021] S. Das, S. Voleti, T. Saha-Dasgupta, and A. Paramekanti, Xy magnetism, kitaev exchange, and long-range frustration in the ${J}_{\mathrm{eff}}=\frac{1}{2}$ honeycomb cobaltates, Phys. Rev. B 104, 134425 (2021).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Physical Review B 50, 17953 (1994).
  * Koepernik and Eschrig [1999] K. Koepernik and H. Eschrig, Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme, Physical Review B 59, 1743 (1999).
  * Riedl  _et al._ [2019] K. Riedl, Y. Li, R. Valentí, and S. M. Winter, Ab initio approaches for low-energy spin hamiltonians, physica status solidi (b) 256, 1800684 (2019).
  * Winter  _et al._ [2016a] S. M. Winter, Y. Li, H. O. Jeschke, and R. Valentí, Challenges in design of kitaev materials: Magnetic interactions from competing energy scales, Physical Review B 93, 214431 (2016a).
  * Mostofi  _et al._ [2014] A. A. Mostofi, J. R. Yates, G. Pizzi, Y. S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, An updated version of wannier90: A tool for obtaining maximally-localised Wannier functions, Comput. Phys. Commun. 185, 2309 (2014).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Physical review B 54, 11169 (1996a).
  * Eschrig  _et al._ [2004] H. Eschrig, M. Richter, and I. Opahle, Relativistic solid state calculations, in  _Theoretical and Computational Chemistry_ , Vol. 14 (Elsevier, 2004) pp. 723–776.
  * [29] See Supplemental Material at [URL will be inserted by publisher] for DFT calculation details.
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Computational materials science 6, 15 (1996b).
  * Kresse and Hafner [1994] G. Kresse and J. Hafner, Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium, Physical Review B 49, 14251 (1994).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Physical Review Letters 77, 3865 (1996).
  * Liechtenstein  _et al._ [1995] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Density-functional theory and strong interactions: Orbital ordering in mott-hubbard insulators, Physical Review B 52, 5467 (1995).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Physical Review B 13, 5188 (1976).
  * Schüler  _et al._ [2018] M. Schüler, O. E. Peil, G. J. Kraberger, R. Pordzik, M. Marsman, G. Kresse, T. O. Wehling, and M. Aichhorn, Charge self-consistent many-body corrections using optimized projected localized orbitals, Journal of Physics Condensed Matter 30, 10.1088/1361-648X/aae80a (2018).
  * Koepernik  _et al._ [2021] K. Koepernik, O. Janson, Y. Sun, and J. Brink, Symmetry conserving maximally projected wannier functions, arXiv preprint arXiv:2111.09652 (2021).
  * Foyevtsova  _et al._ [2013] K. Foyevtsova, H. O. Jeschke, I. I. Mazin, D. I. Khomskii, and R. Valentí, Ab initio analysis of the tight-binding parameters and magnetic interactions in na2iro3, Phys. Rev. B 88, 035107 (2013).
  * Blaha  _et al._ [2018] P. Blaha, K. Schwarz, G. K. H. Madsen, D. Kvasnicka, J. Luitz, R. Laskowski, F. Tran, and L. D. Marks,  _WIEN2k: An Augmented Plane Wave plus Local Orbitals Program for Calculating Crystal Properties_ (Vienna University of Technology, Austria, 2018).
  * Blaha  _et al._ [2020] P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. H. Madsen, and L. D. Marks, Wien2k: An apw+lo program for calculating the properties of solids, The Journal of Chemical Physics 152, 074101 (2020), https://doi.org/10.1063/1.5143061 .
  * Dordević [2008] T. Dordević, BaCo2(AsO4)2, Acta Crystallogr. Sect. E 64, i58 (2008).
  * Xiang  _et al._ [2011] H. Xiang, E. Kan, S.-H. Wei, M.-H. Whangbo, and X. Gong, Predicting the spin-lattice order of frustrated systems from first principles, Physical Review B 84, 224429 (2011).
  * Wellm  _et al._ [2021] C. Wellm, W. Roscher, J. Zeisner, A. Alfonsov, R. Zhong, R. Cava, A. Savoyant, R. Hayn, J. van den Brink, B. Büchner, _et al._ , Frustration enhanced by kitaev exchange in a j eff= 1 2 triangular antiferromagnet, Physical Review B 104, L100420 (2021).
  * Winter [2022] S. M. Winter, to be published (2022).
  * Motome  _et al._ [2020] Y. Motome, R. Sano, S. Jang, Y. Sugita, and Y. Kato, Materials design of kitaev spin liquids beyond the jackeli-khaliullin mechanism, J. Phys. Cond. Matt. 32, 404001 (2020).
  * Rau and Kee [2014] J. G. Rau and H.-Y. Kee, Trigonal distortion in the honeycomb iridates: Proximity of zigzag and spiral phases in na2iro3, arXiv preprint arXiv:1408.4811 (2014).
  * Rau  _et al._ [2014] J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Generic spin model for the honeycomb iridates beyond the kitaev limit, Phys. Rev. Lett. 112, 077204 (2014).
  * Chaloupka and Khaliullin [2015] J. c. v. Chaloupka and G. Khaliullin, Hidden symmetries of the extended kitaev-heisenberg model: Implications for the honeycomb-lattice iridates ${A}_{2}{\mathrm{iro}}_{3}$, Phys. Rev. B 92, 024413 (2015).
  * Halloran  _et al._ [2022] T. Halloran, F. Desrochers, E. Z. Zhang, T. Chen, L. E. Chern, Z. Xu, B. Winn, M. Graves-Brook, M. Stone, A. I. Kolesnikov, _et al._ , Geometrical frustration versus kitaev interactions in baco $\\_2$(aso $\\_4$) $\\_2$, arXiv preprint arXiv:2205.15262 (2022).
  * Sugano [2012] S. Sugano,  _Multiplets of transition-metal ions in crystals_ (Elsevier, 2012).
  * Lines [1963] M. Lines, Magnetic properties of co cl 2 and ni cl 2, Phys. Rev. 131, 546 (1963).
  * Nishimoto  _et al._ [2016] S. Nishimoto, V. M. Katukuri, V. Yushankhai, H. Stoll, U. K. Rößler, L. Hozoi, I. Rousochatzakis, and J. V. D. Brink, Strongly frustrated triangular spin lattice emerging from triplet dimer formation in honeycomb li2iro3, Nature Communications 7, 10273 (2016).
  * Winter  _et al._ [2016b] S. M. Winter, Y. Li, H. O. Jeschke, and R. Valenti, Challenges in design of kitaev materials: Magnetic interactions from competing energy scales, Phys. Rev. B 93, 214431 (2016b).
  * Maksimov and Chernyshev [2020] P. A. Maksimov and A. L. Chernyshev, Rethinking $\alpha\text{$-$}{\mathrm{rucl}}_{3}$, Phys. Rev. Research 2, 033011 (2020).
  * Luttinger and Tisza [1946] J. M. Luttinger and L. Tisza, Theory of dipole interaction in crystals, Phys. Rev. 70, 954 (1946).
  * Li  _et al._ [2016] Y.-D. Li, X. Wang, and G. Chen, Anisotropic spin model of strong spin-orbit-coupled triangular antiferromagnets, Phys. Rev. B 94, 035107 (2016).
  * Rousochatzakis  _et al._ [2016] I. Rousochatzakis, U. K. Rössler, J. van den Brink, and M. Daghofer, Kitaev anisotropy induces mesoscopic ${\mathbb{z}}_{2}$ vortex crystals in frustrated hexagonal antiferromagnets, Phys. Rev. B 93, 104417 (2016).
  * Maksimov  _et al._ [2019] P. A. Maksimov, Z. Zhu, S. R. White, and A. L. Chernyshev, Anisotropic-exchange magnets on a triangular lattice: Spin waves, accidental degeneracies, and dual spin liquids, Phys. Rev. X 9, 021017 (2019).
  * Rastelli  _et al._ [1979] E. Rastelli, A. Tassi, and L. Reatto, Non-simple magnetic order for simple hamiltonians, Physica B+C 97, 1 (1979).
  * Zhu  _et al._ [2018] Z. Zhu, P. A. Maksimov, S. R. White, and A. L. Chernyshev, Topography of spin liquids on a triangular lattice, Phys. Rev. Lett. 120, 207203 (2018).
  * Zhu and White [2015] Z. Zhu and S. R. White, Spin liquid phase of the $s=\frac{1}{2}\phantom{\rule{4.0pt}{0.0pt}}{J}_{1}-{J}_{2}$ heisenberg model on the triangular lattice, Phys. Rev. B 92, 041105 (2015).
  * Cônsoli  _et al._ [2020] P. M. Cônsoli, L. Janssen, M. Vojta, and E. C. Andrade, Heisenberg-kitaev model in a magnetic field: $1/s$ expansion, Phys. Rev. B 102, 155134 (2020).
  * Iqbal  _et al._ [2016] Y. Iqbal, W.-J. Hu, R. Thomale, D. Poilblanc, and F. Becca, Spin liquid nature in the heisenberg ${J}_{1}-{J}_{2}$ triangular antiferromagnet, Phys. Rev. B 93, 144411 (2016).
  * Fouet  _et al._ [2001] J. B. Fouet, P. Sindzingre, and C. Lhuillier, An investigation of the quantum j1-j2-j3model on the honeycomb lattice, The European Physical Journal B - Condensed Matter and Complex Systems 20, 241 (2001).
  * Gong  _et al._ [2013] S.-S. Gong, D. N. Sheng, O. I. Motrunich, and M. P. A. Fisher, Phase diagram of the spin-$\frac{1}{2}$ ${J}_{1}$-${J}_{2}$ heisenberg model on a honeycomb lattice, Phys. Rev. B 88, 165138 (2013).
  * Gong  _et al._ [2014] S.-S. Gong, W. Zhu, D. N. Sheng, O. I. Motrunich, and M. P. A. Fisher, Plaquette ordered phase and quantum phase diagram in the spin-$\frac{1}{2}$ ${J}_{1}\text{$-$}{J}_{2}$ square heisenberg model, Phys. Rev. Lett. 113, 027201 (2014).
  * Gong  _et al._ [2015] S.-S. Gong, W. Zhu, L. Balents, and D. N. Sheng, Global phase diagram of competing ordered and quantum spin-liquid phases on the kagome lattice, Phys. Rev. B 91, 075112 (2015).
  * Rusnačko  _et al._ [2019] J. Rusnačko, D. Gotfryd, and J. c. v. Chaloupka, Kitaev-like honeycomb magnets: Global phase behavior and emergent effective models, Phys. Rev. B 99, 064425 (2019).
  * Suzuki  _et al._ [2021] H. Suzuki, H. Liu, J. Bertinshaw, K. Ueda, H. Kim, S. Laha, D. Weber, Z. Yang, L. Wang, H. Takahashi, K. Fürsich, M. Minola, B. V. Lotsch, B. J. Kim, H. Yavaş, M. Daghofer, J. Chaloupka, G. Khaliullin, H. Gretarsson, and B. Keimer, Proximate ferromagnetic state in the kitaev model material $\alpha$-rucl3, Nature Communications 12, 4512 (2021).
  * White [1992] S. R. White, Density matrix formulation for quantum renormalization groups, Physical Review Letters 69, 2863 (1992).
  * Fishman  _et al._ [2020] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor software library for tensor network calculations (2020), arXiv:2007.14822 .
  * Chubukov and Golosov [1991] A. V. Chubukov and D. I. Golosov, Quantum theory of an antiferromagnet on a triangular lattice in a magnetic field, Journal of Physics: Condensed Matter 3, 69 (1991).
  * Yamamoto  _et al._ [2014] D. Yamamoto, G. Marmorini, and I. Danshita, Quantum phase diagram of the triangular-lattice $xxz$ model in a magnetic field, Phys. Rev. Lett. 112, 127203 (2014).
  * Li  _et al._ [2014] P. H. Y. Li, R. F. Bishop, and C. E. Campbell, Phase diagram of a frustrated spin-$\frac{1}{2}$ ${J}_{1}$-${J}_{2}$ $\mathit{X}xz$ model on the honeycomb lattice, Phys. Rev. B 89, 220408 (2014).
  * Liu  _et al._ [2016] C. Liu, R. Yu, and X. Wang, Semiclassical ground-state phase diagram and $\text{multi-}q$ phase of a spin-orbit-coupled model on triangular lattice, Phys. Rev. B 94, 174424 (2016).
  * Huyan  _et al._ [2022] S. Huyan, J. Schmidt, E. Gati, R. Zhong, R. J. Cava, P. C. Canfield, and S. L. Bud’ko, Hydrostatic pressure effect on the co-based honeycomb magnet ${\mathrm{baco}}_{2}{({\mathrm{AsO}}_{4})}_{2}$, Phys. Rev. B 105, 184431 (2022).