# Magnetic Couplings in Edge-Sharing $d^{7}$ Compounds

Stephen M. Winter  Department of Physics and Center for Functional Materials,
Wake Forest University, NC 27109, USA

###### Abstract

High-spin $d^{7}$ Co(II) compounds have recently been identified as possible
platforms for realising highly anisotropic and bond-dependent couplings
featured in quantum-compass models such as the celebrated Kitaev model. In
order to evaluate this potential, we consider all symmetry-allowed
contributions to the magnetic exchange for ideal edge-sharing bonds. Though a
combination of ab-initio and cluster many-body calculations we conclude that
bond-dependent couplings are generally suppressed in favor of Heisenberg
exchange for realistic materials. Consequences for several prominent materials
including Na2Co2TeO6 and BaCo2(AsO4)2 are discussed.

##  I Introduction

Pursuit of strongly anisotropic $d$-block magnets has been motivated by the
possibility of material realization of so-called quantum compass
modelsNussinov and Van Den Brink (2015), such Kitaev’s celebrated honeycomb
modelKitaev (2006). In these materials, competition between different bond-
dependent magnetic interactions produces an extensive classical degeneracy
conducive to quantum spin liquid ground statesHermanns  _et al._ (2018);
Broholm  _et al._ (2020); Zhou  _et al._ (2017). Realising these conditions in
real materials requires precise tuning and suppression of the usual isotropic
magnetic exchange. This can be accomplished, in principle, in edge-sharing
compounds with $d^{5}$ filling and strong spin-orbital coupling. Remarkably,
for ideal considerations, the specific spin-orbital composition of the local
moments suppresses all couplings except those bond-dependent Ising couplings
precisely prescribed by Kitaev’s modelJackeli and Khaliullin (2009). This
revelation led to a flurry of studiesSingh and Gegenwart (2010); Plumb  _et
al._ (2014); Banerjee  _et al._ (2016); Winter  _et al._ (2017a); Trebst
(2017); Banerjee  _et al._ (2017) in $5d^{5}$ Ir(IV) compounds such as A2IrO3
(A = Na, Li) and the $4d^{5}$ Ru(III) compound $\alpha$-RuCl3. These studies
have revealed clear evidence of dominant bond-dependent anisotropic couplings
in these compoundsHwan Chun  _et al._ (2015); Suzuki  _et al._ (2021), leading
to a variety of anomalous behaviors from the breakdown of conventional magnon
excitationsWinter  _et al._ (2017b) to the possibility of a field-induced
spin-liquidBanerjee  _et al._ (2018); Kasahara  _et al._ (2018); Yokoi  _et
al._ (2021). However, while Kitaev couplings are thought to be the largest
interaction, other couplings of similar magnitude always lift the classical
degeneracy leading to magnetic order at zero-field.

In this context, the seminal work of Liu et al.Liu and Khaliullin (2018); Liu
(2021); Liu  _et al._ (2020) and Sano et al. Sano  _et al._ (2018) renewed
hope for realizing Kitaev’s spin liquid, by showing that the magnetic exchange
in high-spin $3d^{7}$ Co(II) ions may also produce dominant Kitaev
interactions for ideal considerations. In particular, these studies assumed
the dominant hopping between metals occurs via metal-ligand hybridization,
which suppresses other couplings. While this condition is satisfied for
$5d^{5}$ Ir(IV) compounds such as A2IrO3, the presence of significant direct
metal-metal hopping in $4d^{5}$ $\alpha$-RuCl3 is the primary source of non-
Kitaev interactionsRau  _et al._ (2014); Winter  _et al._ (2016). It is not
clear that these assumptions are satisfied in $3d$ systems.

The possibility of strong bond-dependent Kitaev interactions also challenges
the conventional view that Co(II) compounds typically have bond-independent
XXZ anisotropy largely driven by the effects of local crystal field
distortions on the $j_{1/2}$ doublets Lines (1963); Oguchi (1965). For
example, CoNb2O6 (CNO) is considered to be a prototypal 1D Ising
ferromagnetScharf  _et al._ (1979); Maartense  _et al._ (1977); Kobayashi  _et
al._ (1999), and has been studied in the context of transverse-field Ising
criticality Lee  _et al._ (2010); Coldea  _et al._ (2010); Morris  _et al._
(2014). The structure consists of zigzag chains of edge-sharing CoO6
octahedra, which can be considered as alternating X- and Y-bonds per Fig.
1(a). While this bonding geometry might be expected to produce large bond-
dependent couplings, the dominant nearest neighbor interaction is known to
have an Ising form $S_{i}^{\alpha}S_{j}^{\alpha}$ with a common $\alpha$-axis
for all bonds. Recent studies have highlighted the importance of small
deviationsFava  _et al._ (2020); Morris  _et al._ (2021), but it is
nonetheless evident that the Kitaev coupling is not dominant.

More recently, the pursuit of 2D honeycomb materials with large bond-dependent
couplings has drawn attention to Na3Co2SbO6 (NCSO), and Na2Co2TeO6 (NCTO).
Both materials show zigzag antiferromagnetic orderLefrançois  _et al._ (2016);
Bera  _et al._ (2017); Wong  _et al._ (2016); Bera  _et al._ (2017); Chen  _et
al._ (2021). This ground state is natural for strong bond-dependent
couplingsChaloupka  _et al._ (2013); Rau  _et al._ (2014), although longer-
range Heisenberg $J_{2}$ and $J_{3}$ may also be invokedFouet  _et al._
(2001); Kimchi and You (2011) Indeed, analysis of inelastic neutron scattering
has led to a wide variety of proposed models for the couplingsChen  _et al._
(2021); Songvilay  _et al._ (2020); Kim  _et al._ (2021); Lin  _et al._
(2021), which span the entire range from dominant Heisenberg to dominant
Kitaev. Overall, the relative role of nearest neighbor bond-dependent coupling
vs. longer range Heisenberg exchange remains unclear.

Figure 1: (a) Three types of edge-sharing bonds, with definition of global
$(xyz)$ coordinates. (b) Energy level diagram showing the splitting of the
local electronic levels in the absence of spin-orbit coupling.

Two more honeycomb materials of recent interest are BaCo2(AsO4)2 (BCAO) and
BaCo2(PO4)2 (BCPO). Of these, BCPO displays only short-range incommensurate
correlations, suggesting strong frustrationNair  _et al._ (2018). BCAO orders
in a state intermediary between zigzag antiferromagnetic and ferromagnetic
states with unconventional magnon dispersionRegnault  _et al._ (2006, 2018),
which has been discussed as an incommensurate helimagnetRegnault  _et al._
(1977) or double stripe zigzagRegnault  _et al._ (2018). Under applied field
in-plane, BCAO undergoes a series of phase transitionsRegnault  _et al._
(1977) between magnetization plateaus, and was proposed to host a field-
induced spin-liquidZhong  _et al._ (2020); Zhang  _et al._ (2021). However,
this was recently called into question due to the appearance of sharp magnon
modes in each of the phasesShi  _et al._ (2021). As with NCTO, the relative
role of different couplings is a subject of much discussion; the first ab-
initio studiesDas  _et al._ (2021); Maksimov  _et al._ (2022) favored a nearly
XXZ model, in contrast with the assumption of large Kiteav interactions.

All of these findings call for a reinvestigation of the magnetic couplings in
edge-sharing Co(II) materials. In this work, we find, in contrast to the
assumptions of the initial theoretical analysis, that ligand-mediated hopping
is not large in these compounds. For this reason the character of the magnetic
couplings is significantly altered from the expected Kitaev form. In
particular, ferromagnetic Heisenberg $J$ typically dominates, while the myriad
of smaller anisotropic couplings may appear depending on the specific details
of the hopping and crystal field distortions.

The paper is organized as follows: We first review the single-ion ground
state, and the effect of crystal field distortions on the spin-orbital
composition of the $j_{1/2}$ moments. We then analyze the full set of relevant
symmetry-allowed hoppings in edge-sharing bonds. On the basis of these
hoppings, we then compute the resulting magnetic couplings. Finally, we
discuss the results in the context of materials of recent interest.

##  II Single-Ion Considerations

###  II.1 Local Electronic State

At each metal atom, we consider a Hamiltonian that is a sum, respectively, of
Coulomb interactions, crystal-field splitting, and spin-orbit coupling:

$\displaystyle\mathcal{H}_{i}=\mathcal{H}_{U}+\mathcal{H}_{\rm
CFS}+\mathcal{H}_{\rm SOC}$ (1)

The Coulomb interactions are most generally written:

$\displaystyle\mathcal{H}_{U}=\sum_{\alpha,\beta,\delta,\gamma}\sum_{\sigma,\sigma^{\prime}}U_{\alpha\beta\gamma\delta}\
c_{i,\alpha,\sigma}^{\dagger}c_{i,\beta,\sigma^{\prime}}^{\dagger}c_{i,\gamma,\sigma^{\prime}}c_{i,\delta,\sigma}$
(2)

where $\alpha,\beta,\gamma,\delta$ label different $d$-orbitals. 111The
coefficients $U_{\alpha\beta\gamma\delta}$ may be grouped according to the
number of unique orbital indices, from one to four. For example, the intra-
orbital Hubbard terms $n_{i,\alpha,\uparrow}n_{i,\alpha,\downarrow}$ have one
unique index $\alpha$, while the inter-orbital Hubbard terms
$n_{i,\alpha,\sigma}n_{i,\beta,\sigma^{\prime}}$ have two unique indices
$\alpha,\beta$. In the spherically symmetric approximation Sugano (2012), the
Coulomb coefficients with three and four indices vanish unless at least one of
the orbitals is an $e_{g}$ orbital. For this reason, $t_{2g}$-only (and
$e_{g}$-only) models reduce to the familiar Kanamori formgeorges2013strong;
Pavarini (2014), which includes only Hubbard density-density repulsion, Hund’s
exchange, and pair-hopping contributions. However, when both $e_{g}$ and
$t_{2g}$ orbitals are considered together, it is important to include the full
rotationally symmetric Coulomb terms. This is particularly true when computing
anisotropic magnetic exchange, because any approximations to the Coulomb
Hamiltonian are likely to explicitly break rotational symmetry, leading to
erroneous sources of anisotropy. In the spherically symmetric approximation
Sugano (2012), the coefficients $U_{\alpha\beta\gamma\delta}$ are all related
to the three Slater parameters $F_{0},F_{2},F_{4}$. In terms of these, the
familiar $t_{2g}$ Kanamori parameters are, for example:

$\displaystyle U_{t2g}=F_{0}+\frac{4}{49}\left(F_{2}+F_{4}\right)$ (3)
$\displaystyle J_{t2g}=\frac{3}{49}F_{2}+\frac{20}{441}F_{4}$ (4)

We take the approximate ratio $F_{4}/F_{2}=5/8$, following Ref. Pavarini,
2014. The full parameterization is described in Ref. Sugano, 2012. Unless
otherwise stated, we use $U_{t2g}=3.25$ eV, and $J_{t2g}=0.7$ eV to model
Co(II) compounds, following Ref. Das  _et al._ , 2021.

For the crystal-field Hamiltonian, we consider an ideal trigonal distortion
within $D_{3d}$ site symmetry. The Hamiltonian can be written:

$\displaystyle\mathcal{H}_{\rm
CFS}=\sum_{\sigma}\mathbf{c}_{i,\sigma}^{\dagger}\ \mathbb{D}\
\mathbf{c}_{i,\sigma}$ (5)

where:

$\displaystyle\mathbf{c}_{i,\sigma}^{\dagger}=\left(c_{i,yz,\sigma}^{\dagger}\
c_{i,xz,\sigma}^{\dagger}\ c_{i,xy,\sigma}^{\dagger}\
c_{i,z^{2},\sigma}^{\dagger}\ c_{i,x^{2}-y^{2},\sigma}^{\dagger}\right)$ (6)

In terms of the global $(xyz)$ coordinates defined in Fig. 1(a), the CFS
matrix can be written:

$\displaystyle\mathbb{D}=\left(\begin{array}[]{ccccc}0&\Delta_{2}&\Delta_{2}&0&0\\\
\Delta_{2}&0&\Delta_{2}&0&0\\\ \Delta_{2}&\Delta_{2}&0&0&0\\\
0&0&0&\Delta_{1}&0\\\ 0&0&0&0&\Delta_{1}\end{array}\right)$ (12)

where $\Delta_{1}$ is the $t_{2g}$-$e_{g}$ splitting, and $\Delta_{2}$ is the
trigonal term. Generally, $\Delta_{2}>0$ corresponds to trigonal elongation,
as shown in Fig. 1(b), although the actual sign is further influenced by the
details of the ligand environments and longer ranged Coulomb potentials.
Without SOC, the $t_{2g}$ levels are split into a doubly degenerate $e$ pair
and a singly degenerate $a$ level, with $E_{a}-E_{e}=3\Delta_{2}$. As
discussed below, the trigonal splitting has a strong impact on the nature of
the local moments.

For the “high-spin” $d^{7}$ case, the ground state has nominal configuration
$(t_{2g})^{5}(e_{g})^{2}$, with three unpaired electrons ($S=3/2$), as shown
in Fig. 1(b). In the absence of trigonal splitting ($\Delta_{2}=0$), there is
a three-fold orbital degeneracy associated with the $t_{2g}$ levels, leading
to an effective orbital momentum $L_{\rm eff}=1$. Spin-orbit coupling
$\mathcal{H}_{\rm SOC}=\lambda\mathbf{L}\cdot\mathbf{S}$ splits the resulting
multiplets into $J_{\rm eff}=1/2$, 3/2, and $5/2$ states. The $j_{1/2}$
doublet is always the ground state, and furnishes the effective spin-orbital
moment relevant at low energiesLines (1963). The resulting $J_{\rm eff}$
multiplets are composed of many configurations belonging to different orbital
occupancies and spin values. However, third row metals typically satisfy
$J_{H},\Delta_{1}\gg\lambda,\Delta_{2}$, such that configurations belonging
precisely to the $(t_{2g})^{5}(e_{g})^{2}$, $S=3/2$ manifold carry the
dominant weight. When projected into this manifold, the ground state doublet
can be written in terms of $|m_{L},m_{S}\rangle$ states asLines (1963):

$\displaystyle\left|j_{1/2},+\frac{1}{2}\right\rangle=$ $\displaystyle\
c_{1}\left|-1,\frac{3}{2}\right\rangle+c_{2}\left|0,\frac{1}{2}\right\rangle+c_{3}\left|1,-\frac{1}{2}\right\rangle$
(13) $\displaystyle\left|j_{1/2},-\frac{1}{2}\right\rangle=$ $\displaystyle\
c_{1}\left|1,-\frac{3}{2}\right\rangle+c_{2}\left|0,-\frac{1}{2}\right\rangle+c_{3}\left|-1,\frac{1}{2}\right\rangle$
(14)

where the coefficients $c_{n}$ vary with $\Delta_{2},\lambda$. The pure $L,S$
multiplets $|m_{L},m_{S}\rangle$ can be conveniently expressed in terms of the
single-particle levels with precise orbital momentum:

$\displaystyle|e_{a,\sigma}\rangle=$ $\displaystyle\ |d_{z^{2},\sigma}\rangle$
(15) $\displaystyle|e_{b,\sigma}\rangle=$ $\displaystyle\
|d_{x^{2}-y^{2},\sigma}\rangle$ (16) $\displaystyle|t_{+,\sigma}\rangle=$
$\displaystyle\
-\frac{1}{\sqrt{2}}\left(|d_{yz,\sigma}\rangle+i|d_{xz,\sigma}\rangle\right)$
(17) $\displaystyle|t_{0,\sigma}\rangle=$ $\displaystyle\
|d_{xy,\sigma}\rangle$ (18) $\displaystyle|t_{-,\sigma}\rangle=$
$\displaystyle\
\frac{1}{\sqrt{2}}\left(|d_{yz,\sigma}\rangle-i|d_{xz,\sigma}\rangle\right)$
(19)

This leads to:

$\displaystyle\left|-1,\frac{3}{2}\right\rangle=$ $\displaystyle\
\left|e_{a,\uparrow}e_{b,\uparrow}t_{+,\uparrow}t_{0,\uparrow}t_{0,\downarrow}t_{-,\uparrow}t_{-,\downarrow}\right\rangle$
(20) $\displaystyle\left|0,\frac{1}{2}\right\rangle=$ $\displaystyle\
\frac{1}{\sqrt{3}}\left|e_{a,\uparrow}e_{b,\uparrow}t_{+,\uparrow}t_{+,\downarrow}t_{0,\downarrow}t_{-,\uparrow}t_{-,\downarrow}\right\rangle$
(21) $\displaystyle\
+\frac{1}{\sqrt{3}}\left|e_{a,\uparrow}e_{b,\downarrow}t_{+,\uparrow}t_{+,\downarrow}t_{0,\uparrow}t_{-,\uparrow}t_{-,\downarrow}\right\rangle$
$\displaystyle\
+\frac{1}{\sqrt{3}}\left|e_{a,\downarrow}e_{b,\uparrow}t_{+,\uparrow}t_{+,\downarrow}t_{0,\uparrow}t_{-,\uparrow}t_{-,\downarrow}\right\rangle$
$\displaystyle\left|1,-\frac{1}{2}\right\rangle=$ $\displaystyle\
\frac{1}{\sqrt{3}}\left|e_{a,\uparrow}e_{b,\downarrow}t_{+,\uparrow}t_{+,\downarrow}t_{0,\uparrow}t_{0,\downarrow}t_{-,\downarrow}\right\rangle$
(22) $\displaystyle\
+\frac{1}{\sqrt{3}}\left|e_{a,\downarrow}e_{b,\uparrow}t_{+,\uparrow}t_{+,\downarrow}t_{0,\uparrow}t_{0,\downarrow}t_{-,\downarrow}\right\rangle$
$\displaystyle\
+\frac{1}{\sqrt{3}}\left|e_{a,\downarrow}e_{b,\downarrow}t_{+,\uparrow}t_{+,\downarrow}t_{0,\uparrow}t_{0,\downarrow}t_{-,\uparrow}\right\rangle$

The time-reversed partners can be similarly obtained. For $\Delta_{2}=0$, the
coefficients are $c_{1}=1/\sqrt{2}$, $c_{2}=1/\sqrt{3}$, $c_{3}=1/\sqrt{6}$.
In this same limit, the multiplet energies satisfy:

$\displaystyle E_{3/2}-E_{1/2}=\frac{1}{2}\lambda$ (23) $\displaystyle
E_{5/2}-E_{1/2}=\frac{4}{3}\lambda$ (24)

With $\lambda_{\rm Co}\approx 60$ meV, the $j_{1/2}\to j_{3/2}$ excitation is
expected to appear in the range of $\sim 30$ meV, as has been seen
experimentally in numerous compoundsSarte  _et al._ (2018); Ross  _et al._
(2017); Songvilay  _et al._ (2020); Kim  _et al._ (2021).

###  II.2 Local Effects of Trigonal Distortion

Figure 2: Evolution of the single-ion properties with trigonal CFS
$\Delta_{2}$. (a) Energy spectrum. (b) Wavefunction coefficients $c_{n}$. (c)
$g$-tensor components. $g_{||}$ refers to the direction parallel to the
trigonal distortion axis, $\hat{x}+\hat{y}+\hat{z}$.

For finite $\Delta_{2}$, the composition of the doublet is significantly
altered. Here, we review similar discussions in Ref. Lines, 1963; Liu, 2021.
In Fig. 2, we show the evolution of the local spectrum as a function of
$\Delta_{2}/\lambda$, as well as the coefficients $c_{n}$ and $g$-tensor for
the lowest doublet.

In the limit of large trigonal elongation $\Delta_{2}>0$, the unpaired hole in
the $t_{2g}$ levels occupies the singly degenerate $a$ level, thus quenching
the orbital moment completely. This corresponds to $c_{2}\to 1$ and
$c_{1},c_{3}\to 0$. The energetic splitting between the lowest two doublets
becomes small, thus restoring the fourfold degeneracy of the nearly pure
$S=3/2$ moment. The $m_{s}=\pm 1/2$ states lie slightly below the $m_{s}=\pm
3/2$ states, due to residual easy-plane single-ion anisotropy. As such, the
$g$-tensor for the lowest doublet satisfies $g_{\perp}>g_{||}$, where $g_{||}$
refers to the component along the trigonal axis. However, a model
incorporating only the lowest doublet remains valid only as long as the
single-ion anisotropy remains large compared to the intersite magnetic
exchange (roughly, if $\Delta_{2}<\lambda/2$).

For the opposite case of trigonal compression $\Delta_{2}<0$, the unpaired
hole in the $t_{2g}$ levels occupies the doubly degenerate $e$ levels, thus
retaining some orbital degeneracy consistent with $L_{\rm eff}=1/2$. This
corresponds to $c_{1}\to 1$, $c_{2},c_{3}\to 0$. The effect of SOC is then to
split the $S=3/2$, $L_{\rm eff}=1/2$ manifold into four doublets. Since the
lowest doublet corresponds to pure $m_{s}=\pm 3/2$, this may be considered as
strong easy-axis single-ion anisotropy. Consistently, the $g$-tensor satisfies
$g_{||}\gg g_{\perp}$ in this limit. The gap between the lowest doublets
converges to $\lambda/3\sim 20$ meV, which should typically exceed the
intersite magnetic coupling. For this reason, a model incorporating only the
lowest doublet may remain valid for large $\Delta_{2}<0$.

##  III Edge-Sharing Bond Hoppings

###  III.1 General Form

The effective $d$-$d$ hopping between metal sites is described by:

$\displaystyle\mathcal{H}_{\rm
hop}=\sum_{ij,\sigma}\mathbf{c}_{i,\sigma}^{\dagger}\ \mathbb{T}_{ij}\
\mathbf{c}_{j,\sigma}$ (25)

For an ideal edge-sharing bond, $C_{2v}$ symmetry restricts the form of the
hopping matrices. In terms of the global $(x,y,z)$ coordinates defined in Fig.
3, the matrices are constrained to take the following form, for the Z-bond:

$\displaystyle\mathbb{T}_{Z}=\left(\begin{array}[]{ccccc}t_{1}&t_{2}&0&0&0\\\
t_{2}&t_{1}&0&0&0\\\ 0&0&t_{3}&t_{6}&0\\\ 0&0&t_{6}&t_{4}&0\\\
0&0&0&0&t_{5}\end{array}\right)$ (31)

Of these, $t_{1},t_{3},t_{4}$, and $t_{5}$ are primarily direct hopping
between metal atoms, as shown in Fig. 3. Only $t_{2}$ and $t_{6}$ have
significant contributions from both direct hopping and hybridization with the
ligands.

Figure 3: Summary of symmetry allowed hoppings for ideal Z-bonds with $C_{2v}$
symmetry. $t_{1},t_{3},t_{4}$, and $t_{5}$ arise from direct metal-metal
hopping, while $t_{2}$ and $t_{6}$ have contributions from both direct and
ligand-assisted processes. The global $(xyz)$ and local
$(\hat{e}_{1}\hat{e}_{2}\hat{e}_{3})$ coordinates are shown.

###  III.2 Survey of Materials

There are two main factors affecting the balance of direct vs. ligand-assisted
hopping: (i) the degree of hybridization with the ligands, and (ii) the Co-Co
bond lengths. In general, metal-ligand hybridization is typically lower in
third row metal compounds than their $4d$ and $5d$ counterparts, particularly
for $t_{2g}$ orbitals. It is precisely this effect that reduces $t_{2g}-e_{g}$
splitting $\Delta_{1}$ for $3d$ metals, which is required for stability of the
high-spin state in Co $3d^{7}$ compounds. For this reason, ligand-assisted
hopping is expected to be suppressed overall. Real materials span a wide range
Co-Co distances in edge-sharing Co(II) compounds, e.g. from $\sim 2.9$ Å in
BaCo2(AsO4)2Dordević (2008) to $\sim 3.9$ Å in CoI2Wyckoff and Wyckoff (1963).
While we leave complete discussion of individual materials for later work, it
is useful to establish realistic ranges of hoppings. In order to do so, we
employed fully relativistic density functional theory calculations performed
with FPLO at the GGA (PBE) level. Hopping integrals were extracted by
formulating Wannier orbitals via projection onto atomic $d$-orbitals and/or
$p$-orbitals.

Figure 4: Evolution of the relevant hoppings as a function of Co-Co distance.
Solid lines correspond to hypothetical stretched cubic CoO (see text). Points
correspond to real materials; BCAO = BaCo2(AsO4)2, NCTO = Na2Co2TeO6, NCSO =
Na3Co2SbO6, CNO = CoNb2O6. (a) Hoppings in the $d$-only scheme. (b) Hoppings
in the $p+d$ scheme.

In order to get a general idea of the of the bond-length dependence, we first
considered hypothetical cubic CoO (NaCl type; $Fm\bar{3}m$) structures with a
symmetrically stretched unit cells. Hoppings for the 5-band $3d$-only fitting
are shown in Fig. 4(a). This construction maintains 90∘ Co-O-Co bond angles,
which deviates slightly from real materials, but nonetheless provides insight.
In particular, we find in the entire range of Co-Co bond lengths, that direct
hopping is the largest, leading to $|t_{3}|>|t_{2}|,|t_{6}|$. This trend is
also true for estimates of real materials. In particular, we show in Fig. 4(a)
results for several prominent materials based on literature structures:
CoNb2O6 (Ref. Sarvezuk  _et al._ , 2011), BaCo2(AsO4)2 (Ref. Dordević, 2008),
Na3Co2SbO6 (Ref. Songvilay  _et al._ , 2020), and Na2Co2TeO6 (Ref. Xiao  _et
al._ , 2019)222The Na2Co2TeO6 structure contains disorder in the Na position,
in which each Na position has occupancy 2/3. To perform calculations, we
artificially increased the occupancy to 1, which corresponds to Na3Co2TeO6. It
is expected this change in the filling should have minimal impact on the
computed hoppings.. For each case, the cubic projection coordinates were
defined to be orthogonal but minimize the difference with the corresponding
Co-O bond vectors in the (distorted) octahedra. From these results, it is
evident that the physical region corresponds to large $t_{3}$ and subdominant
$t_{6}$. By contrast, $t_{2}$ is suppressed, such that
$|t_{2}|\sim|t_{1}|,|t_{4}|,|t_{5}|\lesssim 0.05$ eV. A similar situationWellm
_et al._ (2021) was recently proposed for Na2BaCo(PO4)2. For materials with
Co-Co bond lengths $\sim$ 3.0 Å, direct hopping almost certainly dominates.
This differs from the previous theoretical works predicting large Kitaev
couplingsLiu and Khaliullin (2018); Liu (2021); Liu  _et al._ (2020); Sano
_et al._ (2018), which considered ligand-mediated hopping $t_{2}$ and $t_{6}$
to be the largest. This discrepancy calls for a reexamination of the magnetic
couplings.

Finally, in Fig. 4(b), we show Slater-Koster hoppings extracted from the CoO
calculations by fitting with an 8-band $(3d+2p)$ model including explicitly
the O orbitals. These are relevant for considering some exchange processes
(see below). In terms of these, the $d$-only hoppings are given approximately
by:

$\displaystyle t_{2}\approx$ $\displaystyle\
-\frac{1}{2}t_{dd}^{\pi}+\frac{(t_{pd}^{\pi})^{2}}{\Delta_{pd}}$ (32)
$\displaystyle t_{3}\approx$ $\displaystyle\ t_{dd}^{\sigma}$ (33)
$\displaystyle t_{6}\approx$ $\displaystyle\
\frac{\sqrt{3}}{4}t_{dd}^{\sigma}-\frac{t_{pd}^{\pi}t_{pd}^{\sigma}}{\Delta_{pd}}$
(34)

where $\Delta_{pd}=4.5$ eV is the charge-transfer energy from Co $d$ to O $p$
orbitals. For $t_{2}$ and $t_{6}$, the contributions from ligand-assisted
hopping is positive, while the direct hopping contribution is negative.

Figure 5: (a-c) Evolution of the magnetic couplings in the $d$-only model for
ideal edge-sharing bond with no trigonal distortion ($\Delta_{2}=0$) and
$\Delta_{1}=1.1$ eV, $U=3.25$ eV, $J_{H}=0.7$ eV,
$t_{1}=|t_{3}|/4,t_{4}=t_{5}=-|t_{3}|/4,t_{6}=+0.1$ eV. The ferromagnetic
correction $\delta J$ due to ligand exchange processes is not included (see
text). (d) Computed couplings along the path indicated in (a-c), interpolating
between the direct and ligand-assisted hopping regimes. A correction $\delta
J=-2$ meV has been added.

##  IV Magnetic Couplings

###  IV.1 General Form

For ideal edge-sharing bonds with $C_{2v}$ symmetry, the magnetic couplings
may be written in the familiar formRau  _et al._ (2014):

$\displaystyle\mathcal{H}_{ij}=$ $\displaystyle\ J\
\mathbf{S}_{i}\cdot\mathbf{S}_{j}+K\
S_{i}^{\gamma}S_{j}^{\gamma}+\Gamma\left(S_{i}^{\alpha}S_{j}^{\beta}+S_{i}^{\beta}S_{j}^{\alpha}\right)$
$\displaystyle+\Gamma^{\prime}\left(S_{i}^{\alpha}S_{j}^{\gamma}+S_{i}^{\gamma}S_{j}^{\alpha}+S_{i}^{\beta}S_{j}^{\gamma}+S_{i}^{\gamma}S_{j}^{\beta}\right)$
(35)

where $(\alpha,\beta,\gamma)=(x,y,z)$ for the Z-bonds, $(y,z,x)$ for the
X-bonds, and $(z,x,y)$ for the Y-bonds, in terms of the global $xyz$
coordinates. In order to estimate the couplings in the following sections, we
exactly diagonalize the full $d$-only Hamiltonian
$\mathcal{H}_{U}+\mathcal{H}_{\rm CFS}+\mathcal{H}_{\rm SOC}+\mathcal{H}_{\rm
hop}$ for two neighboring sites. The couplings are extracted by projecting
onto the ideal $j_{1/2}$ doublets defined in eq’n (13, 14). This procedure is
analogous to Ref. Winter  _et al._ , 2016, and is guaranteed to yield
couplings that converge to the results of perturbation theory with respect to
$\mathcal{H}_{\rm hop}$.

As discussed in Ref. Liu and Khaliullin, 2018, there are several different
electronic processes that contribute to the magnetic couplings at low orders
in the full $p+d$ model. These can be grouped into two categories: (i) those
involving excited states with up to one hole occupying the ligand orbitals,
and (ii) those involving multiple excited ligand holes simultaneously. The
majority of these processes are captured, in principle, in the downfolded
$d$-only hopping model, assuming suitably renormalized hopping and on-site
Coulomb terms. We assume that the hopping integrals extracted from DFT
incorporate this renormalization already. However, our approach does not
capture the subset of processes in category (ii) in which two holes meet on a
ligand in different $p$-orbitals, and interact via Hund’s coupling. Such
processes effectively renormalize the nearest neighbor Coulomb terms when
downfolded, which we have not considered explicitly. We therefore estimate the
effects of the additional contributions. From Ref. Liu and Khaliullin, 2018,
there is a correction to both $J$ and $K$ given approximately by:

$\displaystyle\delta J\approx$ $\displaystyle\
-\frac{\gamma}{\Delta_{pd}^{2}}\left(\frac{5}{2}(t_{pd}^{\sigma})^{4}+\frac{3}{2}(t_{pd}^{\sigma})^{2}(t_{pd}^{\pi})^{2}+(t_{pd}^{\pi})^{4}\right)$
(36) $\displaystyle\delta K\approx$ $\displaystyle\
-\frac{\gamma}{\Delta_{pd}^{2}}\left(\frac{1}{2}(t_{pd}^{\sigma})^{2}(t_{pd}^{\pi})^{2}-\frac{1}{2}(t_{pd}^{\pi})^{4}\right)$
(37) $\displaystyle\gamma=$ $\displaystyle\
\frac{40J_{H}^{p}}{81(\Delta_{pd}+U_{p}/2)^{2}}$ (38)

where $J_{H}^{p}$ is the Hund’s coupling at the ligand, $U_{p}$ is the excess
Coulomb repulsion at the ligand, $\Delta_{pd}\approx 4.5$ eV is the charge-
transfer gap. We take the same approximations as Ref. Liu and Khaliullin, 2018
($U_{p}=0.7\ U_{t2g},\ J_{H}^{p}=0.3\ U_{p}$), and consider
$t_{pd}^{\sigma}\approx 1$ eV, $t_{pd}^{\pi}\approx-0.5$ eV, according to Fig.
4(b). From this, we estimate the correction to the Kitaev coupling to be
negligible $\delta K\sim 0.1$ meV, while the corrections to the Heisenberg
coupling may be typically in the range $\delta J\sim-2$ to $-6$ meV. The
remaining contributions to the exchange are investigated in the next sections.

###  IV.2 General Hopping Dependence

In the following, we focus on the contributions to the magnetic exchange from
(downfolded) $d$-$d$ hopping. We first consider the case $\Delta_{2}=0$. With
the choice, $\Gamma^{\prime}=0$ strictly. Up to second order in hopping, the
couplings may be written:

$\displaystyle J=$ $\displaystyle\
\mathbf{t}\cdot\mathbb{M}_{J}\cdot\mathbf{t}^{T}+\delta J$ (39) $\displaystyle
K=$ $\displaystyle\ \mathbf{t}\cdot\mathbb{M}_{K}\cdot\mathbf{t}^{T}+\delta K$
(40) $\displaystyle\Gamma=$ $\displaystyle\
\mathbf{t}\cdot\mathbb{M}_{\Gamma}\cdot\mathbf{t}^{T}$ (41)

where:

$\displaystyle\mathbf{t}=\left(t_{1}\ t_{2}\ t_{3}\ t_{4}\ t_{5}\
t_{6}\right)$ (42)

and $\mathbb{M}$ is a function of $F_{n},\lambda,\Delta_{n}$. We use
$\Delta_{1}=1.1$ eV and $\lambda=0.06$ eV, which is consistent with estimates
from DFT in the previous sections and $U_{t2g}=3.25$ eV, and $J_{t2g}=0.7$ eV,
following Ref. Das  _et al._ , 2021. To estimate $\mathbb{M}$ for these
parameters, we computed the magnetic couplings for a grid of hoppings
$-0.05<t_{n}<+0.05$ and fit the resulting couplings. This provides an estimate
of the couplings in the perturbative regime:

$\displaystyle\mathbb{M}_{J}=$ $\displaystyle\
\left(\begin{array}[]{c|cccccc}&t_{1}&t_{2}&t_{3}&t_{4}&t_{5}&t_{6}\\\
\hline\cr t_{1}&-55&0&143&10&2&0\\\ t_{2}&0&-76&0&0&0&-77\\\
t_{3}&0&0&-33&2&2&0\\\ t_{4}&0&0&0&260&0&0\\\ t_{5}&0&0&0&0&259&0\\\
t_{6}&0&0&0&0&0&165\end{array}\right)$ (50) $\displaystyle\mathbb{M}_{K}=$
$\displaystyle\
\left(\begin{array}[]{c|cccccc}&t_{1}&t_{2}&t_{3}&t_{4}&t_{5}&t_{6}\\\
\hline\cr t_{1}&128&0&-119&-9&35&0\\\ t_{2}&0&-108&0&0&0&86\\\
t_{3}&0&0&-8&-2&5&0\\\ t_{4}&0&0&0&-4&0&0\\\ t_{5}&0&0&0&0&1&0\\\
t_{6}&0&0&0&0&0&-147\end{array}\right)$ (58)
$\displaystyle\mathbb{M}_{\Gamma}=$ $\displaystyle\
\left(\begin{array}[]{c|cccccc}&t_{1}&t_{2}&t_{3}&t_{4}&t_{5}&t_{6}\\\
\hline\cr t_{1}&0&-34&0&0&0&49\\\ t_{2}&0&0&-116&-2&1&0\\\
t_{3}&0&0&0&0&0&-31\\\ t_{4}&0&0&0&0&0&-67\\\ t_{5}&0&0&0&0&0&0\\\
t_{6}&0&0&0&0&0&0\end{array}\right)$ (66)

in units of 1/eV. Recall, for real materials we generally anticipate
$|t_{3}|>|t_{6}|>|t_{1}|\sim|t_{2}|\sim|t_{4}|\sim|t_{5}|$. Furthermore,
$t_{1}>0$, $t_{3}<0$, $t_{4}<0$, $t_{6}>0$. These results highlight several
key aspects:

Heisenberg $J$: For $J$, there are various contributions of different sign.
Those arising from hopping between $t_{2g}$ orbitals ($t_{1},t_{2},t_{3}$) are
exclusively ferromagnetic. Processes involving hopping between $e_{g}$
orbitals ($t_{4},t_{5}$) are exclusively antiferromagnetic. The terms related
to $e_{g}$-$t_{2g}$ hopping tend to be antiferromagnetic $\propto t_{6}^{2}$
and $t_{2}t_{6}$ given that ab-initio tends to yield $t_{2}<0$ and $t_{6}>0$.

Kitaev $K$: There are also different contributions to $K$ of varying sign.
Hopping between $e_{g}$ orbitals ($t_{4},t_{5}$) makes little contribution to
the anisotropic couplings overall. The sign of the contribution from
$t_{2g}$-$t_{2g}$ hopping depends on the balance of transfer integrals: terms
$\propto t_{2}^{2}$ are ferromagnetic, while terms $\propto t_{1}^{2}$ and
$t_{1}t_{3}$ are antiferromagnetic. Contributions related to $t_{2g}$-$e_{g}$
hopping may take both signs: terms $\propto t_{6}^{2}$ are ferromagnetic,
while terms $\propto t_{2}t_{6}$ depend on the sign of $t_{2}$.

Off-diagonal $\Gamma$: For the off-diagonal couplings, the primary
contribution arises at order $t_{2}t_{3}$, and as a result
$\text{sign}(\Gamma)\approx-\text{sign}(t_{2}t_{3})$. There are no
contributions that are diagonal with respect to the hopping pairs. A similar
result appears in Ref. Liu and Khaliullin, 2018.

The appearance of hopping combinations such as $t_{2}t_{6}$, which do not
conserve $t_{2g}$ and $e_{g}$ occupancy, may appear surprising at first. If
the ground state doublets have approximate configuration
$(t_{2g})^{5}(e_{g})^{2}$, one might expect terms mixing the occupancy to be
forbidden at low orders, because they do not connect ground states. However,
in reality, neither occupancy is preserved by either spin-orbit coupling or
the full Coulomb terms, which are treated exactly (not perturbatively) in this
approach.

For general parameters, we expect all three couplings to be finite. The
computed hopping-dependence of $K,J,\Gamma$ are shown in Fig. 5 for the choice
$t_{1}=|t_{3}|/4,t_{4}=t_{5}=-|t_{3}|/4,t_{6}=+0.1$ eV, which is compatible
with the ab-initio estimates. With this choice, we interpolate between the
limits of dominant ligand vs. direct hopping.

In the hypothetical regime of pure ligand-mediated hopping ($t_{2}$ and
$t_{6}$), we find $\Gamma=0$, while contributions from $d$-$d$ hopping satisfy
$J>0$ is antiferromagnetic and $K<0$ is ferromagnetic. These findings verify
expectations from perturbation theory for this limitLiu and Khaliullin (2018);
Sano  _et al._ (2018). The Kitaev coupling is the largest, with values
$|K/J|\sim 1-10$ depending on the precise balance of hoppings. If we consider
also the ferromagnetic correction $\delta J\sim-2$ to $-6$ meV discussed in
the previous section, the overall sign of $J$ should reverse, and the
magnitude may be suppressed, such that dominant Kitaev coupling is possible
with some tuning.

By contrast, for the physically relevant region of large $t_{3}$ and finite
values of all hoppings, we anticipate that ferromagnetic $J<0$ is the dominant
coupling, particularly due to contributions $\propto t_{1}t_{3}$ and the
correction $\delta J$. In fact, $\delta J$ (which is just the regular
ferromagnetic exchange for 90∘ bondsGoodenough (1963)) is the largest
contribution. All possible combinations of signs of $K$ and $\Gamma$ are
possible depending on the balance of hoppings, but their magnitude is
suppressed relative to $J$. For very short Co-Co bond lengths, where the
direct hopping contribution to $t_{2}$ is the largest ($t_{2}<0$), the
tendency is for $K,\Gamma<0$. For longer bond lengths, where ligand-mediated
contributions are the largest ($t_{2}>0$), then $K,\Gamma>0$.

###  IV.3 Effect of Trigonal Distortion

We next consider the effects of trigonal distortion. Given the relatively
small value of the atomic SOC constant $\lambda_{\rm Co}\approx 60$ meV, small
distortions may be relevant for Co(II) compounds. This make clarifying the
size and sign of $\Delta_{2}$ important for modelling such materials.
Following Ref. Lines, 1963, the alterations to the nature of the local moments
are expected to induce significant uniaxial anisotropy along the trigonal
axis. While the $K,J,\Gamma,\Gamma^{\prime}$ notation is convenient for
discussing the couplings in the limit $\Delta_{2}\to 0$, the axial anisotropy
is more apparent in alternative local XXZ coordinates shown in Fig. 3. In
particular, for each bond we define local coordinates: $\hat{e}_{1}$ is
parallel to the bond and $\hat{e}_{3}=(\hat{x}+\hat{y}+\hat{z})/\sqrt{3}$ is
along the global trigonal axis. Thus, the couplings may be written Ross  _et
al._ (2011); Maksimov  _et al._ (2019):

$\displaystyle\mathcal{H}_{ij}=$ $\displaystyle\
J_{xy}\left(S_{i}^{1}S_{j}^{1}+S_{i}^{2}S_{j}^{2}\right)+J_{z}S_{i}^{3}S_{j}^{3}$
(67) $\displaystyle\
+2J_{\pm\pm}\left(S_{i}^{1}S_{j}^{1}-S_{i}^{2}S_{j}^{2}\right)+J_{z\pm}\left(S_{i}^{3}S_{j}^{2}+S_{i}^{2}S_{j}^{3}\right)$

where the superscript numbers refer to the local directions. The two
parameterizations may be relatedMaksimov  _et al._ (2019) via:

$\displaystyle J_{xy}=$ $\displaystyle\
J+\frac{1}{3}\left(K-\Gamma-2\Gamma^{\prime}\right)$ (68) $\displaystyle
J_{z}=$ $\displaystyle\ J+\frac{1}{3}\left(K+2\Gamma+4\Gamma^{\prime}\right)$
(69) $\displaystyle J_{\pm\pm}=$ $\displaystyle\
-\frac{1}{6}\left(K+2\Gamma-2\Gamma^{\prime}\right)$ (70) $\displaystyle
J_{z\pm}=$ $\displaystyle\
-\frac{\sqrt{2}}{3}\left(K-\Gamma+\Gamma^{\prime}\right)$ (71)

A similar parameterization was also suggested in Ref. Liu, 2021; Liu  _et al._
, 2020. In general, for $\Delta_{2}<0$, as the moments become more axial,
components of the exchange along the $\hat{e}_{3}$ direction are expected to
be enhanced compared to the $\hat{e}_{1}$ and $\hat{e}_{2}$ directionsLines
(1963). As a result $J_{xy}$ and $J_{\pm\pm}$ should be suppressed relative to
$J_{z}$ and $J_{z\pm}$. Trigonal elongation $\Delta_{2}>0$ should have the
opposite effect. As the moments become more planar, $J_{xy}$ and $J_{\pm\pm}$
should be relatively enhanced.

Figure 6: Anisotropic ferromagnetic corrections as a function of trigonal
distortion, for $\delta J_{0}=-2$ meV. (a) Local XXZ scheme. (b) Global
$J,K,\Gamma,\Gamma^{\prime}$ scheme.  Figure 7: Magnetic couplings for ideal
edge-sharing bond with finite trigonal distortion. The parameters are
otherwise the same as Fig. 5. (a-f): $\Delta_{2}/\lambda=-0.5$, (g-l):
$\Delta_{2}/\lambda=+0.5$. The ferromagnetic corrections $\delta
J,\delta\Gamma,\delta\Gamma^{\prime}$ due to ligand exchange processes are not
included (see text). (e,k): Couplings along the path interpolating between
direct and ligand-assisted hoppings depicted in (a,g) in the
$J,K,\Gamma,\Gamma^{\prime}$ scheme. (f,l): Couplings along the path in the
$J_{z},J_{xy},J_{z\pm},J_{\pm\pm}$ scheme. For (e,f,k,l), solid lines indicate
the results of $d$-$d$ exchange only, and dashed lines indicate corrected
values $J+\delta J$, according to $\delta J_{0}=-2$ meV.

Following Ref. Lines, 1963; Liu  _et al._ , 2020, the ferromagnetic
corrections to $J$ resulting from ligand exchange processes are rendered
anisotropic, with:

$\displaystyle\delta J_{xy}\approx$ $\displaystyle\ u_{xy}^{2}\ \delta J_{0}$
(72) $\displaystyle\delta J_{z}\approx$ $\displaystyle\ u_{z}^{2}\ \delta
J_{0}$ (73) $\displaystyle u_{xy}=$ $\displaystyle\
\frac{3}{5}\left(2\sqrt{3}c_{1}c_{3}+2c_{2}^{2}\right)$ (74) $\displaystyle
u_{z}=$ $\displaystyle\ \frac{3}{5}\left(1+2(c_{1}^{2}-c_{3}^{2})\right)$ (75)

with $c_{n}$ given in Fig. 2(b) and $\delta J_{0}$ defined according to eq’n
(36). In the $J,K,\Gamma,\Gamma^{\prime}$ parameterisation, these corrections
correspond to:

$\displaystyle\delta J=$ $\displaystyle\
\frac{1}{3}\left(u_{z}^{2}+2u_{xy}^{2}\right)\delta J_{0}$ (76)
$\displaystyle\delta\Gamma=$ $\displaystyle\
\delta\Gamma^{\prime}=\frac{1}{3}\left(u_{z}^{2}-u_{xy}^{2}\right)\delta
J_{0}$ (77)

As with the undistorted case, the corrections to the anisotropic couplings
$J_{\pm\pm}$ and $J_{z\pm}$ are predicted to be small. The evolution of the
ferromagnetic corrections with distortion are shown in Fig. 6. For
$\Delta_{2}<0$, the ratio $\delta J_{z}/\delta J_{xy}>1$ is, in principle,
unbounded (and should increase continuously with trigonal distortion). For
$\Delta_{2}>0$ the degree of anisotropy is restricted, because the distortion-
induced effects are bounded $1/4\lesssim J_{z}/J_{xy}<1$. The lower bound is
reached for large $\Delta_{2}$, where the orbital moment is quenched, thus
restoring the full degeneracy of the $S=3/2$ states. However, a low-energy
model including only the lowest doublet would no longer be sufficient, so this
limit does not represent a physically sensible model. In terms of global
coordinates, the trigonal distortion primarily introduces off-diagonal
couplings, where $\Delta_{2}<0$ tends to be associated with
$\delta\Gamma,\delta\Gamma^{\prime}<0$, and vice versa.

To explore the exchange contributions from (downfolded) $d$-$d$ hopping, we
recomputed the couplings using exact diagonalization with significant
distortion $\Delta_{2}/\lambda=\pm 0.5$ to emphasize the effects. Results are
shown in Fig. 7 for the choice
$t_{1}=|t_{3}|/4,t_{4}=t_{5}=-|t_{3}|/4,t_{6}=+0.1$ eV, which is compatible
with the ab-initio estimates. In Fig. 7(e,f) and (k,l), we also show the
effect of corrections $\delta J$. The results are as follow:

Trigonal compression: For $\Delta_{2}<0$, as shown in Fig. 7(a-e), we find all
four of the couplings $J,K,\Gamma,\Gamma^{\prime}$ may be of similar
magnitude. This is particularly true in the region of large ligand-assisted
hopping. For the physically relevant region of large direct hopping ($t_{3}\gg
t_{2}$), we find that $K$ is still relatively suppressed (same as for
$\Delta_{2}=0$), but large $\Gamma,\Gamma^{\prime}$, with
$\text{sign}(\Gamma,\Gamma^{\prime})\sim\text{sign}(J)$ are induced. These
results are more easily interpreted in the alternative XXZ parameterization
shown in Fig. 7(e). In particular, as the local moments become more axial with
larger trigonal distortion, the coupling becomes dominated by a ferromagnetic
Ising exchange $J_{z}$. Overall, the estimated ferromagnetic correction
$\delta J_{z}$ is quite large compared to the regular $d$-$d$ contributions.
For the physically relevant region, we anticipate $J_{xy}=-2$ to $0$ meV,
$J_{z}=-3$ to $-10$ meV, $J_{\pm\pm}=-0.5$ to $+0.5$ meV, and $J_{z\pm}=-0.5$
to $+1.5$ meV for significant trigonal distortion of $\Delta=-\lambda/2$. . As
a result, we expect such materials to be described mostly by Ising couplings
with a common axis for every bond.

Trigonal elongation: For $\Delta_{2}>0$, we find that $K$ is less suppressed.
The distortions induce off-diagonal couplings following roughly
$\text{sign}(\Gamma,\Gamma^{\prime})\sim-\text{sign}(J)$. In the XXZ
parameterization, this corresponds to an enhancement of $J_{xy}$. In the
hypothetical ligand-assisted hopping region, we find that $J_{z}$ may be
almost completely suppressed due to different values of the ferromagnetic
shifts $\delta J_{z}$ and $\delta J_{xy}$. While $J_{xy}$ appears to be the
largest coupling in this limit, the bond-dependent couplings $J_{z\pm}$ and
$J_{\pm\pm}$ may also remain significant. For the physically relevant region,
we find that $J_{xy}$ is typically the dominant coupling, with
$J_{xy}/J_{z}\sim 4$, which is the hypothetical limit. Overall, we anticipate
$J_{xy}=-2$ to $-10$ meV, $J_{z}=-0.5$ to $-4$ meV, $J_{\pm\pm}=-2$ to $+1$
meV, and $J_{z\pm}=0$ to $+1$ meV for significant trigonal distortion of
$\Delta=+\lambda/2$.

###  IV.4 Longer Range Couplings

While we have discussed above that $t_{2g}$-ligand hybridization should
generally be small in $3d$ metal oxides (as reflected by small
$t_{pd}^{\pi}$), the $e_{g}$-ligand hybridization may still play a significant
role through the large $t_{pd}^{\sigma}$. This is particularly relevant for
third neighbor bonds in honeycomb materials, because it gives rise to a large
hopping between $d_{x^{2}-y^{2}}$ orbitals shown in Fig. 8 at order
$(t_{pd}^{\sigma})^{2}t_{pp}^{\sigma}/\Delta_{pd}^{2}\sim 0.05$ to 0.1 eV.
This is equivalent to a 3rd neighbor $t_{5}$, which allows the associated
coupling to be readily estimated from the matrices $\mathbb{M}$. In
particular, we estimate (for $\Delta_{2}=0$):

$\displaystyle J_{3}\approx+0.5\text{ to}+2.5\text{ meV}$ (78) $\displaystyle
K_{3}\approx\Gamma_{3}\approx 0$ (79)

This is the only major third neighbor hopping pathway, so there are no
additional terms to compete, and a relatively large antiferromagnetic $J_{3}$
should be expected for all honeycomb materials with partially occupied $e_{g}$
orbitals.

Figure 8: 3rd neighbor hopping relevant to $J_{3}$.

##  V Conclusions

In this work, we have considered the magnetic couplings in edge-sharing
$d^{7}$ compounds. On this basis, we make several observations:

(1) All of the edge-sharing Co(II) oxides considered in this work appear to
fall outside the regime of primary focus in previous theoretical studiesLiu
and Khaliullin (2018); Liu (2021); Liu  _et al._ (2020); Sano  _et al._
(2018). In particular, direct hopping likely dominates over ligand-assisted
hopping ($t_{3}\gg t_{2}$). In the realistic regime, we find that $K$ is
generally suppressed compared to $J$, which calls into question models with
dominant $K$ proposed for these materials.

(2) Compared to heavy $d^{5}$ Kitaev materials such as iridates A2IrO3 and
$\alpha$-RuCl3, the weak spin-orbit coupling of Co increases the relative
importance of local distortions. The presence of the $e_{g}$ spins also opens
additional exchange pathways, whose balance depends sensitively on local
parameters such as $J_{H},U$, and $\Delta_{1}$. This makes anticipating the
magnetic Hamiltonian somewhat challenging. For oxides, fortuitous fine-tuning
may result in a different balance of couplings, but we anticipate that
ferromagnetic $J$ (or equivalently $J_{z},J_{xy}$) is likely always the
largest coupling. The signs and magnitudes of the other couplings
$K,\Gamma,\Gamma^{\prime}$ are influenced by the crystal field splitting and
specific details of the hoppings. We find regions with all possible signs and
relative magnitudes. Real materials with small trigonal distortions are likely
described by $|K/J|\sim 0.2$ to 0.5, and $K\approx\Gamma$; specifically:
$J\sim-8$ to $-2$ meV, $K\sim-2$ to $+2$ meV, and $\Gamma\sim-1$ to $+3$ meV.
$\Gamma^{\prime}$ is likely small unless there are significant departures from
ideal symmetry of the bonds. These findings are compatible with the overall
scale of those reported in the literatureRegnault  _et al._ (1977); Nair  _et
al._ (2018); Regnault  _et al._ (2018); Fava  _et al._ (2020). It is not clear
that a uniquely dominant $K$ is possible.

(3) For systems with significant crystal field distortions, our findings are
compatible with the historical description of Co(II) magnetic couplings in
terms of XXZ models by M. E. Lines (Ref. Lines, 1963). This is true
particularly because of the importance of ligand exchange processes, which are
responsible for ferromagnetic couplings in materials with 90∘ bond angles in
the Goodenough-Kanamori description Goodenough (1963). We estimate that these
are at least as important as processes involving $d$-$d$ hopping. In this
case, the considerations discussed in Ref. Lines, 1963; Liu  _et al._ , 2020
become equivalent to the classic results of Lines. For trigonal crystal fields
with $\Delta_{2}<0$ (corresponding to $g_{||}>g_{\perp}$), the Ising
anisotropy induced by the crystal field may be very large, such that the
couplings are dominated by a ferromagnetic $J_{z}$ with a common Ising axis
for every bond. For positive crystal field $\Delta_{2}>0$ (corresponding to
$g_{\perp}>g_{||}$), XXZ anisotropy is more limited, but may still be large
for significant distortions $\Delta_{2}\sim\lambda/2$. Such materials are
generally more desirable for realising strongly bond-dependent couplings.

(4) Regarding NCSO and NCTO: Some constraints can be placed on the
interactions on the basis of the ordered moment directions in the zigzag
state. Related discussions appear in Ref. Sanders  _et al._ , 2021. For NCTO,
it is generally agreed that the ordered moments lie nearly in the honeycomb
plane, oriented along the direction of the ferromagnetic chainsLefrançois  _et
al._ (2016). For the zigzag domain with magnetic wavevector parallel to the
Z-bond, this is $\hat{e}_{2}^{Z}=(1,1,-2)/\sqrt{6}$ in cubic coordinates. This
orientation is generally expectedChaloupka and Khaliullin (2016) for
$K,\Gamma>0$, which is compatible with large $t_{3}$ and small $t_{2}>0$, as
we find in ab-initio for NCTO (see Fig. 4). The antiferromagnetic sign of $K$
is driven by a combination of hopping processes $\propto t_{2}t_{6},\
t_{1}^{2}$ and $t_{1}t_{3}$. Most of these were not previously considered in
the literature. For moments precisely along $\hat{e}_{2}^{Z}$, the magnetic
state is left invariant under a 180∘ rotation around (111) followed by time
reversal; it is thus reasonable to assume the couplings bear the same
symmetry. This places the constraint on the couplings
$\Gamma+\delta\Gamma=K+\Gamma^{\prime}+\delta\Gamma^{\prime}$. Experimental
estimatesKim  _et al._ (2021); Liu  _et al._ (2020); Liu (2021) for NCSO and
NCTO suggest $\Delta_{2}\sim+4$ to $+13$ meV, which corresponds to a
$\delta\Gamma,\delta\Gamma^{\prime}\sim 0.2-0.5$ meV. With these suggestions,
one may then consider the small magnitude of the magnon gap ($\sim 1$ meV)
observed in experiment at both the $\Gamma$-pointLin  _et al._ (2021); Chen
_et al._ (2021) and the ordering wavevectorSongvilay  _et al._ (2020); Kim
_et al._ (2021); Chen  _et al._ (2021); Lin  _et al._ (2021). This would be
anomalous for large departures from XXZ-symmetry. It may further be remarked
that the field-evolution of the ESR modesLin  _et al._ (2021) follow
expectations for moderate easy-plane XXZ anisotropy. Taken together, we
suggest $J_{xy}=-3.25,J_{z}=-2.25,J_{\pm\pm}=-0.125,J_{z\pm}=0$ meV as an
appropriate starting point for analysis. These correspond to
$J_{1}=-3,J_{3}=+2.5,K=\Gamma^{\prime}=+0.25,\Gamma=+0.5$ meV, which are
essentially consistent with the model of Ref. Lin  _et al._ , 2021.

(5) Regarding BCAO: The breadth of experimental data on BCAO, in terms of the
progression of field-induced phases and inelastic neutron data provide a
number of clues towards the magnetic model. While we leave full elaboration
for future studyMaksimov  _et al._ (2022), some comments can be made. A recent
reinvestigationRegnault  _et al._ (2018) of the zero-field structure suggested
it might better be described by a double stripe
$\uparrow\uparrow\downarrow\downarrow$ analogue of the zigzag antiferromagnet,
with moments oriented nearly along the in-plane $\hat{e}_{2}^{Z}$ direction,
as with NCTO. This orientation points to $K,\Gamma>0$. The large discrepancy
between in-plane critical fields (0.2, 0.5 T) and the out-of-plane critical
field (4T)Zhang  _et al._ (2021) suggests significant anisotropy. Indeed, the
$g$-tensor appears to satisfyRegnault  _et al._ (2018) $g_{z}\sim 0.5\
g_{xy}$, and within an XXZ model, $J_{z}\sim 0.4\ J_{xy}$. These findings
point to significant crystal field effects, with $\Delta_{2}\sim 0.2$ to
$0.25\ \lambda$, i.e. $\Delta_{2}\sim 15$ meV, implying significant $\Gamma$
and $\Gamma^{\prime}$ in the global coordinate scheme. In the XXZ scheme, the
nearly in-plane moments suggest small $J_{z\pm}$, while an apparently small
anisotropy between in-plane field directionsRegnault  _et al._ (2018) may
place restrictions on $J_{\pm\pm}$. In contrast, the authors of Ref. Zhang
_et al._ , 2021 have advocated for small $J$, large $K<0$, and small average
off-diagonal coupling $\bar{\Gamma}=(\Gamma+2\Gamma^{\prime})/3$ on the basis
of THz spectroscopy experiments. It should be emphasized that these conditions
are not mutually compatible: small $J_{z\pm}$ and $J_{\pm\pm}$ implies small
$K$, and large anisotropy between $J_{xy}$ and $J_{z}$ implies large
$\Gamma+2\Gamma^{\prime}>0$. If we consider BCAO to be in the physical regime
of hoppings, our findings tend to contradict the Kitaev-dominant model. We
propose a model similar to NCTO: $K$ is small and likely antiferromagnetic,
$J<0$ is the dominant coupling, and $\Gamma,\Gamma^{\prime}>0$ reflect a
planar XY-anisotropy $|J_{xy}|>|J_{z}|$. The anomalous aspects of the ground
state are then understood as a competition between $J_{1},J_{2}$, and $J_{3}$,
as previously suggestedRegnault  _et al._ (1977); Nair  _et al._ (2018);
Regnault  _et al._ (2018). These suggestions are compatible with the recent
ab-initio estimatesDas  _et al._ (2021); Maksimov  _et al._ (2022).

(6) From the perspective of chemistry, it is unclear how to access the
desirable ligand-assisted hopping regime where Kitaev coupling is largest. It
is necessary to increase the metal-ligand hybridization relative to direct
hopping between metal atoms. This may typically be achieved by matching the
electronegativity of the metals and ligands, such that $\Delta_{pd}$ is small.
As a general trend, electronegativity of transition metals increases for
heavier atoms, while the opposite is true for $p$-block ligands. The
combination of heavy metals with heavy ligands results in the most covalent
metal-ligand bonds. However, with increased covalency comes increased
$t_{2g}$-$e_{g}$ splitting, and heavier metals tend to have reduced Coulomb
terms $J_{H}$. As a result, the high-spin $(t_{2g})^{5}(e_{g})^{2}$ state is
only typically achievable in Co(II) compounds. By contrast, Rh(II) and Ni(III)
tend to adopt low-spin $(t_{2g})^{6}(e_{g})^{1}$ ground states. The most
promising avenue then appears to be the combination of Co(II) with heavier
ligands. With this in mind, we computed the hopping integrals for the
triangular lattice compounds CoCl2, CoBr2 and CoI2 according to crystal
structures from Ref. Wilkinson  _et al._ , 1959; Wyckoff and Wyckoff, 1963.
For these compounds, we still estimate $|t_{3}/t_{2}|\sim 4$. Thus, it is not
clear that Kitaev-dominant exchange can be achieved without extreme fine
tuning.

##  VI Acknowledgements

We acknowledge useful discussions with D. Smirnov, Y. Jiang, S. Streltsov, P.
Maksimov, and R. Valenti. We also thank P. Dai for bringing the problem to our
attention. This work was supported by a pilot grant from the Center for
Functional Materials. Computations were performed on the Wake Forest
University DEAC Cluster, a centrally managed resource with support provided in
part by Wake Forest University.

## References

  * Nussinov and Van Den Brink (2015) Z. Nussinov and J. Van Den Brink, Rev. Mod. Phys. 87, 1 (2015).
  * Kitaev (2006) A. Kitaev, Ann. Phys. 321, 2 (2006).
  * Hermanns  _et al._ (2018) M. Hermanns, I. Kimchi, and J. Knolle, Annu. Rev. Condens. Matter Phys. 9, 17 (2018).
  * Broholm  _et al._ (2020) C. Broholm, R. Cava, S. Kivelson, D. Nocera, M. Norman, and T. Senthil, Science 367, eaay0668 (2020).
  * Zhou  _et al._ (2017) Y. Zhou, K. Kanoda, and T.-K. Ng, Rev. Mod. Phys. 89, 025003 (2017).
  * Jackeli and Khaliullin (2009) G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 102, 017205 (2009).
  * Singh and Gegenwart (2010) Y. Singh and P. Gegenwart, Phys. Rev. B 82, 064412 (2010).
  * Plumb  _et al._ (2014) K. Plumb, J. Clancy, L. Sandilands, V. V. Shankar, Y. Hu, K. Burch, H.-Y. Kee, and Y.-J. Kim, Phys. Rev. B 90, 041112 (2014).
  * Banerjee  _et al._ (2016) A. Banerjee, C. Bridges, J.-Q. Yan, A. Aczel, L. Li, M. Stone, G. Granroth, M. Lumsden, Y. Yiu, J. Knolle, S. Bhattacharjee, D. L. Kovrizhin, R. Moessner, D. A. Tennant, D. G. Mandrus, and S. E. Nagler, Nat. Mater. 15, 733 (2016).
  * Winter  _et al._ (2017a) S. M. Winter, A. A. Tsirlin, M. Daghofer, J. van den Brink, Y. Singh, P. Gegenwart, and R. Valenti, J. Phys. Condens. Matter 29, 493002 (2017a).
  * Trebst (2017) S. Trebst, arXiv preprint arXiv:1701.07056 (2017).
  * Banerjee  _et al._ (2017) A. Banerjee, J. Yan, J. Knolle, C. A. Bridges, M. B. Stone, M. D. Lumsden, D. G. Mandrus, D. A. Tennant, R. Moessner, and S. E. Nagler, Science 356, 1055 (2017).
  * Hwan Chun  _et al._ (2015) S. Hwan Chun, J.-W. Kim, J. Kim, H. Zheng, C. C. Stoumpos, C. Malliakas, J. Mitchell, K. Mehlawat, Y. Singh, Y. Choi, T. Gog, A. Al-Zein, M. M. Sala, M. Krisch, J. Chaloupka, G. Jackeli, G. Khaliullin, and B. J. Kim, Nat. Phys. 11, 462 (2015).
  * 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. Yavas, M. Daghofer, J. Chaloupka, G. Khaliullin, H. Gretarsson, and B. Keimer, Nat. Commun. 12, 1 (2021).
  * Winter  _et al._ (2017b) S. M. Winter, K. Riedl, P. A. Maksimov, A. L. Chernyshev, A. Honecker, and R. Valentí, Nat. Commun. 8, 1 (2017b).
  * Banerjee  _et al._ (2018) A. Banerjee, P. Lampen-Kelley, J. Knolle, C. Balz, A. A. Aczel, B. Winn, Y. Liu, D. Pajerowski, J. Yan, C. A. Bridges, A. T. Savici, B. C. Chakoumakos, M. D. Lumsden, D. A. Tennant, R. Moessner, D. G. Mandrus, and S. E. Nagler, npj Quantum Mater. 3, 1 (2018).
  * Kasahara  _et al._ (2018) Y. Kasahara, T. Ohnishi, Y. Mizukami, O. Tanaka, S. Ma, K. Sugii, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, T. Shibauchi, and Y. Matsuda, Nature 559, 227 (2018).
  * Yokoi  _et al._ (2021) T. Yokoi, S. Ma, Y. Kasahara, S. Kasahara, T. Shibauchi, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, C. Hickey, S. Trebst, and Y. Matsuda, Science 373, 568 (2021).
  * Liu and Khaliullin (2018) H. Liu and G. Khaliullin, Phys. Rev. B 97, 014407 (2018).
  * Liu (2021) H. Liu, Int. J. Mod. Phys. B 35, 2130006 (2021).
  * Liu  _et al._ (2020) H. Liu, J. Chaloupka, and G. Khaliullin, Phys. Rev. Lett. 125, 047201 (2020).
  * Sano  _et al._ (2018) R. Sano, Y. Kato, and Y. Motome, Phys. Rev. B 97, 014408 (2018).
  * Rau  _et al._ (2014) J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Phys. Rev. Lett. 112, 077204 (2014).
  * Winter  _et al._ (2016) S. M. Winter, Y. Li, H. O. Jeschke, and R. Valentí, Phys. Rev. B 93, 214431 (2016).
  * Lines (1963) M. Lines, Phys. Rev. 131, 546 (1963).
  * Oguchi (1965) T. Oguchi, J. Phys. Soc. Japan 20, 2236 (1965).
  * Scharf  _et al._ (1979) W. Scharf, H. Weitzel, I. Yaeger, I. Maartense, and B. Wanklyn, J. Magn. Magn. Mater. 13, 121 (1979).
  * Maartense  _et al._ (1977) I. Maartense, I. Yaeger, and B. Wanklyn, Solid State Commun. 21, 93 (1977).
  * Kobayashi  _et al._ (1999) S. Kobayashi, S. Mitsuda, M. Ishikawa, K. Miyatani, and K. Kohn, Phys. Rev. B 60, 3331 (1999).
  * Lee  _et al._ (2010) S. Lee, R. K. Kaul, and L. Balents, Nat. Phys. 6, 702 (2010).
  * Coldea  _et al._ (2010) R. Coldea, D. Tennant, E. Wheeler, E. Wawrzynska, D. Prabhakaran, M. Telling, K. Habicht, P. Smeibidl, and K. Kiefer, Science 327, 177 (2010).
  * Morris  _et al._ (2014) C. Morris, R. V. Aguilar, A. Ghosh, S. Koohpayeh, J. Krizan, R. Cava, O. Tchernyshyov, T. McQueen, and N. Armitage, Phys. Rev. Lett. 112, 137403 (2014).
  * Fava  _et al._ (2020) M. Fava, R. Coldea, and S. Parameswaran, Proc. Natl. Acad. Sci. 117, 25219 (2020).
  * Morris  _et al._ (2021) C. Morris, N. Desai, J. Viirok, D. Hüvonen, U. Nagel, T. Room, J. Krizan, R. Cava, T. McQueen, S. Koohpayeh, R. K. Kaul, and N. P. Armitage, Nat. Phys. 17, 832 (2021).
  * 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, Phys. Rev. B 94, 214416 (2016).
  * Bera  _et al._ (2017) A. K. Bera, S. M. Yusuf, A. Kumar, and C. Ritter, Phys. Rev. B 95, 094424 (2017).
  * Wong  _et al._ (2016) C. Wong, M. Avdeev, and C. D. Ling, J. Solid State Chem. 243, 18 (2016).
  * Chen  _et al._ (2021) W. Chen, X. Li, Z. Hu, Z. Hu, L. Yue, R. Sutarto, F. He, K. Iida, K. Kamazawa, W. Yu, X. Lin, and Y. Li, Phys. Rev. B 103, L180404 (2021).
  * Chaloupka  _et al._ (2013) J. Chaloupka, G. Jackeli, and G. Khaliullin, Phys. Rev. Lett. 110, 097204 (2013).
  * Fouet  _et al._ (2001) J. Fouet, P. Sindzingre, and C. Lhuillier, Eur. Phys. J. B 20, 241 (2001).
  * Kimchi and You (2011) I. Kimchi and Y.-Z. You, Phys. Rev. B 84, 180407 (2011).
  * 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, Phys. Rev. B 102, 224429 (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, J. Phys. Condens. 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, J. Mei, Z. Qu, H. Zhou, J.-G. Park, Y. Wan, and J. Ma, Nat. Commun. 12, 1 (2021).
  * Nair  _et al._ (2018) H. S. Nair, J. Brown, E. Coldren, G. Hester, M. Gelfand, A. Podlesnyak, Q. Huang, and K. Ross, Phys. Rev. B 97, 134409 (2018).
  * Regnault  _et al._ (2006) L. Regnault, C. Boullier, and J. Henry, Physica B Condens. Matter 385, 425 (2006).
  * Regnault  _et al._ (2018) L.-P. Regnault, C. Boullier, and J. Lorenzo, Heliyon 4, e00507 (2018).
  * Regnault  _et al._ (1977) L. Regnault, P. Burlet, and J. Rossat-Mignod, Physica B Condens. Matter 86, 660 (1977).
  * Zhong  _et al._ (2020) R. Zhong, T. Gao, N. P. Ong, and R. J. Cava, Sci. Adv. 6, eaay6953 (2020).
  * Zhang  _et al._ (2021) X. Zhang, Y. Xu, R. Zhong, R. Cava, N. Drichko, and N. Armitage, arXiv preprint arXiv:2106.13418 (2021).
  * Shi  _et al._ (2021) L. Shi, X. Wang, R. Zhong, Z. Wang, T. Hu, S. Zhang, Q. Liu, T. Dong, F. Wang, and N. Wang, Phys. Rev. B 104, 144408 (2021).
  * Das  _et al._ (2021) S. Das, S. Voleti, T. Saha-Dasgupta, and A. Paramekanti, Phys. Rev. B 104, 134425 (2021).
  * Maksimov  _et al._ (2022) P. A. Maksimov, A. V. Ushakov, Z. V. Pchelkina, Y. Li, S. M. Winter, and S. V. Streltsov, arXiv preprint (2022).
  * Note (1) The coefficients $U_{\alpha\beta\gamma\delta}$ may be grouped according to the number of unique orbital indices, from one to four. For example, the intra-orbital Hubbard terms $n_{i,\alpha,\delimiter 52568952}n_{i,\alpha,\delimiter 52573049}$ have one unique index $\alpha$, while the inter-orbital Hubbard terms $n_{i,\alpha,\sigma}n_{i,\beta,\sigma^{\prime}}$ have two unique indices $\alpha,\beta$. In the spherically symmetric approximation Sugano (2012), the Coulomb coefficients with three and four indices vanish unless at least one of the orbitals is an $e_{g}$ orbital. For this reason, $t_{2g}$-only (and $e_{g}$-only) models reduce to the familiar Kanamori formgeorges2013strong; Pavarini (2014), which includes only Hubbard density-density repulsion, Hund’s exchange, and pair-hopping contributions. However, when both $e_{g}$ and $t_{2g}$ orbitals are considered together, it is important to include the full rotationally symmetric Coulomb terms. This is particularly true when computing anisotropic magnetic exchange, because any approximations to the Coulomb Hamiltonian are likely to explicitly break rotational symmetry, leading to erroneous sources of anisotropy.
  * Sugano (2012) S. Sugano,  _Multiplets of transition-metal ions in crystals_ (Elsevier, 2012).
  * Pavarini (2014) E. Pavarini, in  _Many-Electron Approaches in Physics, Chemistry and Mathematics_ (Springer, 2014) pp. 321–341.
  * Sarte  _et al._ (2018) P. M. Sarte, R. A. Cowley, E. E. Rodriguez, E. Pachoud, D. Le, V. García-Sakai, J. W. Taylor, C. D. Frost, D. Prabhakaran, C. MacEwen, A. Kitada, A. J. Browne, M. Songvilay, Z. Yamani, W. J. L. Buyers, J. P. Attfield, and C. Stock, Phys. Rev. B 98, 024415 (2018).
  * Ross  _et al._ (2017) K. A. Ross, J. Brown, R. Cava, J. Krizan, S. E. Nagler, J. Rodriguez-Rivera, and M. B. Stone, Phys. Rev. B 95, 144414 (2017).
  * Dordević (2008) T. Dordević, Acta Crystallogr. E 64, i58 (2008).
  * Wyckoff and Wyckoff (1963) R. W. G. Wyckoff and R. W. Wyckoff,  _Crystal structures_ , Vol. 1 (Interscience publishers New York, 1963).
  * Sarvezuk  _et al._ (2011) P. W. C. Sarvezuk, E. J. Kinast, C. Colin, M. Gusmão, J. Da Cunha, and O. Isnard, Int. J. Appl. Phys. 109, 07E160 (2011).
  * Xiao  _et al._ (2019) G. Xiao, Z. Xia, W. Zhang, X. Yue, S. Huang, X. Zhang, F. Yang, Y. Song, M. Wei, H. Deng, and D. Jiang, Crystal Growth & Design 19, 2658 (2019).
  * Note (2) The Na2Co2TeO6 structure contains disorder in the Na position, in which each Na position has occupancy 2/3. To perform calculations, we artificially increased the occupancy to 1, which corresponds to Na3Co2TeO6. It is expected this change in the filling should have minimal impact on the computed hoppings.
  * Wellm  _et al._ (2021) C. Wellm, W. Roscher, J. Zeisner, A. Alfonsov, R. Zhong, R. J. Cava, A. Savoyant, R. Hayn, J. van den Brink, B. Büchner, O. Janson, and V. Kataev, Phys. Rev. B 104, L100420 (2021).
  * Goodenough (1963) J. B. Goodenough,  _Magnetism and the chemical bond_ , Vol. 1 (Interscience publishers, 1963).
  * Ross  _et al._ (2011) K. A. Ross, L. Savary, B. D. Gaulin, and L. Balents, Phys. Rev. X 1, 021002 (2011).
  * Maksimov  _et al._ (2019) P. Maksimov, Z. Zhu, S. R. White, and A. Chernyshev, Phys. Rev. X 9, 021017 (2019).
  * Sanders  _et al._ (2021) A. L. Sanders, R. A. Mole, J. Liu, A. J. Brown, D. Yu, C. D. Ling, and S. Rachel, arXiv preprint arXiv:2112.12254 (2021).
  * Chaloupka and Khaliullin (2016) J. Chaloupka and G. Khaliullin, Phys. Rev. B 94, 064435 (2016).
  * Wilkinson  _et al._ (1959) M. Wilkinson, J. Cable, E. Wollan, and W. Koehler, Phys. Rev. 113, 497 (1959).