Magnetic Couplings in Edge-Sharing High-Spin $d^7$ Compounds

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 real materials. Consequences for several prominent materials including Na$_2$Co$_2$TeO$_6$ and BaCo$_2$(AsO$_4$)$_2$ are discussed.


I. INTRODUCTION
Pursuit of strongly anisotropic d-block magnets has been motivated by the possibility of material realization of quantum compass models 1 , such as Kitaev's celebrated honeycomb model 2 . In these materials, competition between different bond-dependent magnetic interactions produces an extensive classical degeneracy conducive to quantum spin-liquid ground states [3][4][5] . Realising these conditions in real materials requires precise tuning and suppression of the usual isotropic magnetic exchange. This can be accomplished, in principle, in edgesharing 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 the Kitaev model 6 . This revelation initiated a number of studies 7,8 in 5d 5 Ir(IV) compounds such as A 2 IrO 3 (A = Na, Li) 9 and the 4d 5 Ru(III) compound 10-12 α-RuCl 3 . These studies have revealed clear evidence of dominant bond-dependent anisotropic couplings in these compounds 13,14 , leading to a variety of anomalous behaviors from the breakdown of conventional magnon excitations 15 to the possibility of a field-induced spin liquid [16][17][18] . 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. [19][20][21] and Sano et al. 22 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 transition metal ions occurs via hybridization with ligand ions, which suppresses other couplings. While this condition is satisfied for 5d 5 Ir(IV) compounds such as A 2 IrO 3 , the presence of significant direct metal-metal hopping in 4d 5 α-RuCl 3 is the primary source of non-Kitaev interactions 23,24 . 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 25,26 . For example, CoNb 2 O 6 (CNO) is considered to be a prototypical 1D Ising ferromagnet [27][28][29] , and has been studied in the context of transverse-field Ising criticality [30][31][32] . The structure consists of zigzag chains of edge-sharing CoO 6 octahedra. 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 S α j with a common α-axis for all bonds. Recent studies have highlighted the importance of small deviations 33,34 , 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 Na 3 Co 2 SbO 6 (NCSO), and Na 2 Co 2 TeO 6 (NCTO). Both materials show zigzag antiferromagnetic order [35][36][37][38] . This ground state is natural for strong bonddependent couplings 23,39 , although longer-range Heisenberg J 2 and J 3 may also be invoked 40,41 Indeed, analysis of inelastic neutron scattering has led to a wide variety of proposed models for the couplings 38,[42][43][44] , which span the entire range from dominant Heisenberg to dominant Kitaev. Overall, the relative role of nearest neighbor bonddependent coupling vs. longer range Heisenberg exchange remains unclear.
Two more honeycomb materials of recent interest are BaCo 2 (AsO 4 ) 2 (BCAO) and BaCo 2 (PO 4 ) 2 (BCPO). Of these, BCPO displays only short-range incommensurate correlations, suggesting strong frustration 45 . BCAO orders in a state intermediary between zigzag antiferromagnetic and ferromagnetic states with unconventional magnon dispersion 46,47 , which has been discussed as an incommensurate helimagnet 48 or double stripe zigzag 47 . Under applied field in-plane, BCAO undergoes a series of phase transitions 48 between magnetization plateaus, and was proposed to host a field-induced spin liquid 49,50 . However, this was recently called into question due to the appearance of sharp magnon modes in each of the phases 51 . As with NCTO, the relative role of different couplings is a subject of much discussion; the first abinitio studies 52,53 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 the nearest neighbor couplings, while a variety of smaller anisotropic couplings may appear depending on the specific details of the hopping and crystal field distortions, which are discussed in this work.
The paper is organized as follows: We first define the electronic Hamiltonian, and 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, and present calculations of such hoppings for real materials to establish their physically relevant ranges. On the basis of these hoppings, we then compute the resulting magnetic couplings. Finally, we summarize the results, and discuss them in the context of materials of recent interest.

A. General Form
We first define the electronic Hamiltonian relevant for the edge-sharing bond depicted in Fig. 1(a), which is conventionally called a Z-bond to indicate that the global cubic z-axis is perpendicular to the shared edge. There are generally two different orbital bases that can be employed in the description of magnetic couplings. One may consider explicitly either (i) the full basis of metal d-orbitals and ligand p-orbitals, or (ii) the d-orbitals only, for which the effects of the p-orbitals are downfolded into the d-orbital Wannier function basis. The latter Wannier functions are anti-bonding combinations of d-and p-orbitals depicted schematically in Fig. 1(b). In the following, we mostly consider the downfolded basis, but discuss the deficiencies of this approach below.
In the downfolded d-only basis, the total Hamiltonian is given by: where H hop denotes the hopping between Co sites. H i is the local Hamiltonian at each Co site i, which is a sum, respectively, of Coulomb interactions, crystal-field splitting, and spin-orbit coupling: The SOC term is written: where the atomic SOC constant for Co is roughly λ Co ≈ 60 meV. The Coulomb interactions are most generally written: where α, β, γ, δ label different d-orbitals. 54 In the spherically symmetric approximation 55 , the coefficients U αβγδ 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: Unless otherwise stated, we use U t2g = 3.25 eV, and J t2g = 0.7 eV to model Co(II) compounds, following Ref. 52. We also take the approximate ratio F 4 /F 2 = 5/8, which is appropriate for 3d elements 56 . For the crystal-field Hamiltonian, we consider an ideal trigonal distortion within D 3d site symmetry. The Hamiltonian can be written: where: creates an electron in one of the d-orbital (Wannier) functions at site i. In terms of the global cubic (xyz) coordinates defined in Fig. 2, the CFS matrix can be written: where ∆ 1 is the t 2g -e g splitting, and ∆ 2 is the trigonal term. Trigonal distortions are relevant for BCAO, BCPO, NCTO, and NCSO. Generally, ∆ 2 > 0 corresponds to trigonal elongation, as shown in Fig. 2, 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∆ 2 . As discussed below, the trigonal splitting has a strong impact on the nature of the local moments.
The effective d-d hopping between transition metal sites is described by: For an ideal edge-sharing bond with 90 • metal-ligandmetal bond angles, 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: Of these, t 1 , t 3 , t 4 , and t 5 are primarily direct hopping between the d-orbitals on the metal atoms, as shown in Fig. 3. By contrast, t 2 and t 6 have significant contributions from both direct hopping and ligand-assisted hopping, which arises from the hybridization between the metal and ligand orbitals. Realistic ranges for these hoppings are discussed in section II C. Previous works [19][20][21][22] on the possibility of dominant Kitaev magnetic exchange in high-spin 3d 7 Co(II) ions have highlighted the ideal regime is found for small trigonal splitting (∆ 2 λ) and dominant ligand-assisted hopping (t 2 t 1 , t 3 , t 4 , t 5 , t 6 ). In principle, if either of these conditions is violated, the interactions may deviate significantly from the Kitaev form. We consider these possible deviations in the following sections.

B. Composition of 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. 2. For 3d transition metals, the constants in the Hamiltonian H i typically satisfy J H , ∆ 1 λ, ∆ 2 . This ensures that configurations belonging precisely to the (t 2g ) 5 (e g ) 2 , S = 3/2 manifold carry the dominant weight in the low energy multiplets, even when SOC and trigonal splitting are considered. For this reason, the low-energy single-ion states may be conveniently described by projecting H i into this manifold, and diagonalizing. This procedure was employed in Ref. 25, and is repeated in this section for the purpose of discussion. More detailed discussions can be found in Ref. 20 and 25. In the absence of trigonal splitting (∆ 2 = 0), there is a three-fold orbital degeneracy associated with the t 2g levels, leading to an effective orbital momentum L eff = 1. Spin-orbit coupling H SOC = λL · S splits the low energy multiplets into J eff = 1/2, 3/2, and 5/2 states. For ∆ 2 = 0, the multiplet energies satisfy: Given λ Co ≈ 60 meV, the j 1/2 → j 3/2 excitation is expected to appear in the range of ∼ 30 meV, as has been seen experimentally in numerous compounds 42,43,57,58 . For finite ∆ 2 , the evolution of the multiplet energies is shown in Fig. 4(a). For all values of trigonal splitting, the single-ion ground state is a doublet; it is smoothly connected to the pure J eff = 1/2 doublet as ∆ 2 vanishes. The ground state doublet can be written in terms of |m L , m S states as 25 : where the coefficients c n vary with ∆ 2 , λ as shown in Fig. 4(b). More detailed expressions for the |m L , m S states in terms of the occupation of the single-particle d-orbitals is given in Appendix A. Analytical expression for c n can be found in Ref. 25; for ∆ 2 = 0, the coefficients are c 1 = 1/ √ 2, c 2 = 1/ √ 3, c 3 = 1/ √ 6. From Fig. 4(b), it is evident that the composition of the lowest doublet changes dramatically with finite ∆ 2 .
In the limit of large trigonal elongation ∆ 2 > 0, the wavefunction coefficients converge to c 2 → 1 and c 1 , c 3 → 0. This reflects the fact that the unpaired hole in the t 2g levels would occupy the singly degenerate a level, thus quenching the orbital moment completely (see Fig. 2). As a result, spin-orbit coupling effects are suppressed, and the fourfold degeneracy of the nearly pure S = 3/2 moment is restored, such that the energetic splitting between the lowest two doublets becomes small ( Fig. 4(a)). In this limit, the m s = ±1/2 states lie slightly below the m s = ±3/2 states due to the residual effects of SOC, and thus the splitting of the lowest two doublets may be considered as a residual easy-plane single-ion anisotropy. As such, the g-tensor (Fig. 4(c)) for the lowest doublet satisfies g ⊥ > g || , where g || refers to the component along the trigonal (111) axis.
For the opposite case of trigonal compression ∆ 2 < 0, the unpaired hole in the t 2g levels occupies the doubly degenerate e levels (see Fig. 2), thus retaining some orbital degeneracy consistent with L eff = 1/2. As depicted in Fig. 4(a), the effect of SOC is then to split the S = 3/2, L eff = 1/2 manifold into four doublets. The lowest doublet remains energetically separated from the higher states. For increasing trigonal distortion, the wavefunction coefficients for the lowest doublet approach c 1 → 1, c 2 , c 3 → 0, implying that it is composed of pure m s = ±3/2 states, according to Eq. (10, 11). Consistently, the g-tensor satisfies g || g ⊥ in this limit ( Fig. 4(c)).
Finally, it should be emphasized that an effective J = 1 2 model that incorporates only the lowest doublet remains valid only as long as the energetic splitting between the lowest doublet and first excited state remains large compared to the intersite magnetic exchange. For large ∆ 2 < 0, the gap between the lowest doublets converges to λ/3 ∼ 20 meV, which should typically exceed the intersite magnetic coupling. For this reason, a model incorporating only the lowest doublet may remain valid even for large ∆ 2 < 0. For ∆ 2 > 0, which is expected to be relevant for BCAO, BCPO, NCTO, and NCSO, the range of validity is roughly constrained to ∆ 2 < λ/2. For larger values of ∆ 2 , it would be more appropriate to consider an S = 3/2 model with easy-plane single-ion anisotropy. A detailed analysis of the experimental values of ∆ 2 for various Co(II) compounds may be found in Ref. 20. All of these materials are expected to be well described by a J = 1 2 model, which can be verified experimentally by the observation of the well-defined local j 1/2 → j 3/2 excitations in an energy range larger than the dispersion of magnetic excitations 42,43,57,58 . For the materials of interests, the trigonal splitting is therefore not sufficiently large to invalidate the j 1/2 picture, but may nonetheless have a considerable effect on the magnetic couplings, as discussed in the following sections.

C. Survey of Hopping in Real Materials
As noted in previous sections, the key question for the realization of dominant Kitaev coupling in high-spin d 7 compounds is the relative importance of direct metalmetal vs. ligand-assisted hopping processes. In terms of the usual Slater-Koster hoppings, the downfolded d-only hoppings are given approximately by: where ∆ pd is the charge-transfer energy from Co d to the ligand p orbitals. Here, we have ignored contributions from t δ dd . For t 2 and t 6 , the second terms represent ligand-assisted hopping. There are two main factors that can affect these: (i) the Co-Co bond lengths (which primarily affect the direct overlap of metal d-orbitals on adjacent metal sites) and (ii) the degree of hybridization with the ligands (as quantified by the relative weight of the p-orbitals in the d-orbital Wannier functions, which scales as ∝ t pd /∆ pd ).
In order to investigate the hopping trends with bondlength dependence, we first considered hypothetical cubic CoO (NaCl type; F m3m) structures with a symmetrically stretched unit cells, allowing us to investigate a continuous range of geometries. This construction maintains 90 • Co-O-Co bond angles, which deviates slightly from real materials, but allows us to consider ideal trends. In order to estimate the hoppings, we employed fully relativistic density functional theory calculations performed with FPLO 59,60 at the GGA (PBE 61 ) level. Hopping integrals were extracted by formulating Wannier orbitals via projection onto atomic p and/or d-orbitals. For example, the Slater-Koster hoppings were extracted from the CoO calculations by fitting with an 8-band (3d + 2p) model including explicitly the O orbitals. Results are shown in Fig. 5(a) as a function of Co-Co distance.
As expected, the largest magnitude corresponds to t σ pd , which quantifies the hopping between ligand p and metal e g orbitals. This is followed in magnitude by t π pd , which quantifies the hopping between ligand p and metal t 2g orbitals. The difference in magnitude between t σ pd and t π pd results in a larger degree of metal-ligand hybridization for the Wannier e g orbitals than the t 2g orbitals. This is, of course, the origin of octahedral crystal field splitting ∆ 1 , as described by standard ligand field theory. From the signs of the Slater-Koster hoppings, one can see that the contribution to t 2 and t 6 from ligand-assisted hopping is always positive, while the direct hopping contribution is always negative. The balance of this competition is therefore governed by the charge-transfer energy ∆ pd . For 3d transition metal compounds, this quantity is typically large, which reduces ligand-metal hybridization relative to their 4d and 5d counterparts. It is precisely this effect that leads to smaller t 2g − e g splitting ∆ 1 for 3d transition metal compounds 62 , which is a key requirement for stability of the high-spin state in Co 3d 7 compounds. For this reason, ligand-assisted hopping is expected to be suppressed overall. This leads to strong competition between hopping processes and overall suppression of t 2 and t 6 . In order to demonstrate this effect, we show in Fig. 5(b) (solid lines) the downfolded hoppings extracted for the hypothetical CoO structures by Wannier fitting with a 5-band (3d-only) basis. The results are compati-ble with (but not exactly given by) Eq'ns (16) -(21) with ∆ pd ∼ 4.5 eV. For all lattice constants, we find that t 3 is the dominant hopping rather than the ligand-assisted t 2 and t 6 . While we leave complete discussion of individual materials for later work, it is useful to confirm realistic ranges of hoppings. In order to do so, we performed additional calculations on several prominent oxide materials based on literature structures:  68 . For each case, we employed fully relativistic density functional theory calculations performed with FPLO 59,60 at the GGA (PBE 61 ) level (as above). Hoppings in the downfolded d-only basis were estimated by projection onto d-orbitals, for which the cubic projection coordinates were defined to be orthogonal but minimize the difference with the corresponding Co-X (X = O, Cl, Br, I) bond vectors in the (distorted) octahedra.
Computed hoppings for the real materials are shown in Fig. 5(b). In general, t 1 , t 4 and t 5 are small, and are thus omitted from the figure. From these results, it is evident that the physical region corresponds to large direct hopping t 3 and subdominant ligand-mediated hopping t 6 . Similarly, the ligand-mediated t 2 is significantly suppressed, such that |t 2 | ∼ |t 1 |, |t 4 |, |t 5 | 0.05 eV. We find that this applies equally to Co oxides and halides. A similar situation 69 was recently proposed for Na 2 BaCo(PO 4 ) 2 . On this basis, we conclude: for edgesharing materials with Co-Co bond lengths ∼ 3.0Å, direct hopping almost certainly dominates. This differs from the previous theoretical works predicting large Kitaev couplings [19][20][21][22] , which considered ligand-mediated hopping t 2 and t 6 to be the largest. This discrepancy calls for a reexamination of the magnetic couplings.

A. General Form
For ideal edge-sharing bonds with C 2v symmetry, the magnetic couplings may be written in the familiar form 23 : where (α, β, γ) = (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. This definition turns out to be convenient for the case of trigonal splitting ∆ 2 = 0. For finite ∆ 2 , following Ref. 25, the alterations to the nature of the local moments are expected to induce significant uniaxial anisotropy along the trigonal axis. This axial anisotropy is more apparent in alternative local XXZ coordinates shown in Fig. 3. In particular, for each bond we define local coordinates:ê 1 is parallel to the bond andê 3 = (x +ŷ +ẑ)/ √ 3 is along the global trigonal axis. Thus, the couplings may be written 70,71 : where the superscript numbers refer to the local directions. The two parameterizations may be related 71 via: A similar parameterization was also employed in Ref. 20 and 21. Both parameterizations are employed below, where appropriate.

B. Method and Downfolding Schemes
In previous works, the magnetic couplings have been discussed in terms of perturbation theory in both the full d + p basis, and the downfolded d-only basis. It is useful to discuss how terms in each scheme are related.
As an example of this downfolding, we may consider the ligand-mediated exchange process shown in Fig. 6(a), in which a hole hops from a metal d-orbital to a ligand p-orbital, and subsequently to a d-orbital on an adjacent metal site, thus yielding an excited (d 8 , d 6 ) configuration. From this excited configuration, the reverse process returns the system to the ground (d 7 , d 7 ) configuration. Such processes contribute terms of order O[ (t 2 pd /∆ pd ) 2 /U dd ] to the magnetic exchange, where ∆ pd is the charge transfer energy, and U dd is the Coulomb repulsion between d-electrons. The same process can also be captured within a d-only picture, by identifying t eff = t 2 pd /∆ pd as an effective ligand-mediated d − d hopping, as discussed above. In this case, the magnetic couplings arise at order O[ t 2 eff /U dd ]. Such terms are naturally captured in the downfolded d-only scheme.
We next consider the "cyclic" exchange process depicted in Fig. 6(b). In this case, a hole first hops from one metal to a ligand, and then a second hole hops from the second metal to another ligand. The system is returned to the ground (d 7 , d 7 ) configuration by two additional hops of the holes from the ligands to opposite metal sites from where they started. Such processes contribute to the magnetic exchange at order O[ (t 2 pd /∆ pd ) 2 /(2∆ pd ) ]. In the downfolded d-only basis, these processes are not distinguished from those in  Finally, we consider the process depicted in Fig. 6(c). In this case, two holes hop from their parent metal sites to the same ligand, but occupy different orbitals. They interact with eachother via Hund's coupling at the ligand site, and then hop back to their parent metal sites. Such processes constitute the usual Goodenough-Kanamori ferromagnetic superexchange. When downfolded into the d-only basis, these processes reflect the fact that the tails of the Wannier orbitals extend onto the ligands, which results in an enhancement of the effective intersite Hund's coupling. Such terms are not explicitly included in our d-only model Eq'n (1). We therefore estimate the effects of the additional contributions using previously derived perturbative expressions. From Ref. 19, there is a cor-rection to both J and K given approximately by: where J p H is the Hund's coupling at the ligand, U p is the excess Coulomb repulsion at the ligand ion, ∆ pd ≈ 4.5 eV is the charge-transfer gap. We take the same approximations as Ref. 19 , and consider t σ pd ≈ 1 eV, t π pd ≈ −0.5 eV, according to Fig. 5(a). From this, we estimate the correction to the Kitaev coupling to be negligible δK ∼ 0.1 meV, while the corrections to the Heisenberg coupling may be typically in the range δJ ∼ −2 to −6 meV. The reason for the discrepancy in magnitudes between δJ and δK is that the dominant contribution to δJ comes from hopping between the e g and p orbitals (depicted in Fig. 6(c)), as parameterized by the large hopping t σ pd . This process does not contribute to K. The small correction δK is henceforth neglected.
In the following, we compute all exchange contributions that are captured within the downfolded d-only basis using exact diagonalization of the model Eq'n (1) for two Co sites. The couplings are extracted by projecting onto the ideal j 1/2 doublets defined in eq'n (14,15). This procedure is described in more detail in Ref. 24, and is guaranteed to yield couplings that converge to the results of perturbation theory with respect to H hop for small values of t n /U (with all other contributions to the Hamiltonian, which appear in H i , being treated exactly). The reason for adopting this method is that analytical perturbation theory is generally challenging for finite ∆ 2 and the full set of allowed hoppings t 1−6 . The expressions below are intended to be used with hopping integrals extracted from DFT, which are assumed to be suitably renormalized. The correction δJ must be added to the computed couplings to estimate the omitted ligand exchange processes depicted in Fig. 6(c).

C. General Hopping Dependence
In this section, we first consider the magnetic couplings for the case ∆ 2 = 0, for which Γ = 0 strictly. Up to second order in d − d hopping (and fourth order in d − p hopping), the couplings may be written: where: and M parameterizes the contributions that are captured by (downfolded) d − d hopping. The matrices are a function of F n , λ, ∆ n . To estimate M, we computed the magnetic couplings for a grid of hoppings −0.05 < t n < +0.05 eV using exact diagonalization of the full d-orbital model H U +H CFS +H SOC +H hop for two metal sites, and fit the resulting couplings. For this purpose, we use ∆ 1 = 1.1 eV and λ = 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. 52. This provides an estimate of the couplings in the perturbative regime: in units of 10 −3 /eV. Recall, for real materials we generally anticipate |t 3 | > |t 6 | > |t 1 | ∼ |t 2 | ∼ |t 4 | ∼ |t 5 |. Furthermore, t 1 > 0, t 3 < 0, t 4 < 0, t 6 > 0 (see section II C). 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 gt 2g hopping tend to be antiferromagnetic ∝ t 2 6 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 ∝ t 2 2 are ferromagnetic, while terms ∝ t 2 1 and t 1 t 3 are antiferromagnetic. Contributions related to t 2g -e g hopping may take both signs: terms ∝ t 2 6 are ferromagnetic, while terms ∝ t 2 t 6 depend on the sign of t 2 .
Off-diagonal Γ: For the off-diagonal couplings, the primary contribution arises at order t 2 t 3 , and as a result sign(Γ) is typically the same as −sign(t 2 t 3 ). There are no contributions that are diagonal with respect to the hopping pairs. A similar result appears in Ref. 19.
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, Γ 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. These ratios of hoppings are compatible with the ab-initio estimates in section II C.
With this choice, we interpolate between the limits of dominant ligand-mediated vs. direct hopping.
In the hypothetical regime of pure ligand-mediated hopping (finite t 2 , t 6 > 0 ), we find Γ = 0. If we ignore the ferromagnetic correction δJ, then J > 0 is antiferromagnetic and K < 0 is ferromagnetic. These findings verify expectations from perturbation theory for this limit 19,22 . The Kitaev coupling is the largest, with values |K/J| ∼ 1 − 10 depending on the precise balance of hoppings. If we consider also the ferromagnetic correction δJ ∼ −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 ∝ t 1 t 3 and the correction δJ. In fact, δJ (which is just the regular ferromagnetic superexchange for 90 • bonds 72 ) is the largest contribution. All possible combinations of signs of K and Γ 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, Γ < 0. For longer bond lengths, where ligand-mediated contributions are the largest (t 2 > 0), then K, Γ > 0.

D. Effect of Trigonal Distortion
We next consider the effects of trigonal distortion. Given the relatively small value of the atomic SOC constant λ Co ≈ 60 meV, small distortions may be relevant for Co(II) compounds. This make clarifying the size and sign of ∆ 2 important for modelling such materials.
In general, for ∆ 2 < 0, as the moments become more axial, components of the exchange along theê 3 direction are expected to be enhanced compared to theê 1 andê 2 directions 25 . As a result J xy and J ±± should be suppressed relative to J z and J z± . Trigonal elongation ∆ 2 > 0 should have the opposite effect. As the moments become more planar, J xy and J ±± should be relatively enhanced.
Following Ref. 21 and 25, the ferromagnetic corrections to J resulting from ligand exchange processes are rendered anisotropic, with: with c n given in Fig. 4(b) and δJ 0 defined according to eq'n (28). Thus, in the J, K, Γ, Γ parameterisation, these corrections correspond to: Here, we have continued to neglect the small corrections δK 0 . As with the undistorted case, the corrections to the anisotropic couplings J ±± and J z± are predicted to be small. They have also been neglected in the following. The evolution of the ferromagnetic corrections with distortion are shown in Fig. 8. For ∆ 2 < 0, the ratio δJ z /δJ xy > 1 is, in principle, unbounded (and should increase continuously with trigonal distortion). For ∆ 2 > 0 the degree of anisotropy is restricted, because the distortion-induced effects are bounded 1/4 J z /J xy < 1. The lower bound is reached for large ∆ 2 , where the orbital moment is quenched, thus restoring the full degeneracy of the S = 3/2 states. As discussed above, 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 cubic coordinates [ Fig. 8(b)], the trigonal distortion primarily introduces off-diagonal couplings, where ∆ 2 < 0 tends to be associated with δΓ, δΓ < 0, and vice versa.
To explore the exchange contributions from (downfolded) d-d hopping, we recomputed the couplings using exact diagonalization with significant distortion ∆ 2 /λ = ±0.5 to emphasize the effects. Results are shown in Fig. 9  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. 9(e,f) and (k,l), we also show the effect of corrections δJ. The results are as follow: Trigonal compression: For ∆ 2 < 0, as shown in Fig. 9(a-e), we find all four of the couplings J, K, Γ, Γ 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 t 2 ), we find that K is still relatively suppressed (same as for ∆ 2 = 0), but large Γ, Γ are induced, with sign(Γ, Γ ) typically being the same as sign(J). These results are more easily interpreted in the alternative XXZ parameterization shown in Fig. 9(f). 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 δ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 ±± = −0.5 to +0.5 meV, and J z± = −0.5 to +1.5 meV for significant trigonal distortion of ∆ = −λ/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 ∆ 2 > 0, we find that K is less suppressed. The distortions induce off-diagonal couplings following roughly sign(Γ, Γ ) typically the same as −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 δJ z and δJ xy . While J xy appears to be the largest coupling in this limit, the bond-dependent couplings J z± and J ±± may also remain significant. For the physically relevant region, we find that J xy is typically the dominant coupling, with J xy /J z ∼ 4, which is the hypothetical limit. Overall, we anticipate J xy = −2 to −10 meV, J z = −0.5 to −4 meV, J ±± = −2 to +1 meV, and J z± = 0 to +1 meV for significant trigonal distortion of ∆ = +λ/2.

E. Longer Range Couplings
While we have discussed above that t 2g -ligand hybridization should generally be small in 3d transition metal oxides (as reflected by small t π pd ), the e g -ligand hybridization may still play a significant role in the magnetic exchange due to the large t σ pd . This fact is what leads to significant ferromagnetic corrections δJ for nearest neighbor bonds. For longer range bonds, such as third neighbors in edge-sharing honeycomb lattices, it gives rise to a large hopping between d x 2 −y 2 orbitals shown in Fig. 10 at order (t σ pd ) 2 t σ pp /∆ 2 pd ∼ 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 M. In particular, we estimate (for ∆ 2 = 0): 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. For finite ∆ 2 , the interactions are rendered anisotropic, with XXZ anisotropy reflecting the composition of the local moments. These estimates may be confirmed by considering the triangular lattice compound Ba 3 CoSb 2 O 9 , for which the closest Co 2+ neighbors have the same geometry as depicted in Fig. 10. In this case, it is well established 73-76 that the interactions are described by J xy = +1.7 meV, J z = +1.5 meV. Third neighbor interactions of similar magnitude should be expected for all honeycomb materials BCAO, BCPO, NCSO, and NCTO.

IV. DISCUSSION
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 studies [19][20][21][22] . In particular, direct hopping likely dominates over ligand-assisted hopping (t 3 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) The presence of the e g spins also opens additional exchange pathways. The e g orbitals mix relatively strongly with the ligand p-orbitals in contrast to the t 2g orbitals. This enhances the importance of certain ligand exchange processes, which are responsible for ferromagnetic couplings in materials with 90 • bond angles in the Goodenough-Kanamori description 72 . Analogous contributions are typically an order of magnitude smaller in low-spin d 5 compounds, such as A 2 IrO 3 and α-RuCl 3 , where the e g orbitals are empty. The e g spins also participate in longer range antiferromagnetic couplings, corresponding to third neighbors for honeycomb lattice materials. It may be commonly expected that first and third neighbor interactions have similar magnitude in materials such as BCAO, BCPO, NCTO, and NCSO. Thus, while it was initially hoped that the smaller spatial extent of the 3d orbitals would suppress longer range coupling relative to those found in 4d 5 and 5d 5 compounds, the presence of e g spins likely enhances long-range interactions.
(3) Compared to heavy d 5 Kitaev materials such as iridates A 2 IrO 3 and α-RuCl 3 , the weak spin-orbit coupling of Co increases the relative importance of local distortions. There are also many separate contributions to the magnetic couplings whose balance depends sensitively on local parameters such as J H , U , and ∆ 1 . This makes anticipating the magnetic Hamiltonian somewhat challenging. For Co(II) 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, Γ, Γ are influenced by the crystal field splitting and specific details of the hoppings. We find regions with all possible signs and relative magnitudes.
(4) Real materials with small trigonal distortions (discussed in Section III C) are likely described by |K/J| ∼ 0.2 to 0.5, and K ≈ Γ; specifically: J ∼ −8 to −2 meV, K ∼ −2 to +2 meV, and Γ ∼ −1 to +3 meV. These magnitudes are compatible with the overall scale of interactions reported in the experimental literature 33,45,47,48 . Γ is predicted to be small unless there are significant departures from ideal symmetry of the bonds (i.e. ∆ 2 = 0, C 2h bond symmetry, 90 • M-L-M bond angles). It is not clear that a uniquely dominant K is possible even for weak trigonal distortions.
(5) For systems with significant crystal field distortions (section III D), our findings are compatible with the historical description of Co(II) magnetic couplings in terms of XXZ models by M. E. Lines (Ref. 25). This is true particularly because these considerations apply to the ferromagnetic ligand exchange processes, which we estimate are at least as important as processes involving d-d hopping. In this case, the considerations discussed in Ref. 21 become equivalent to the classic results of Lines.
For trigonal crystal fields with ∆ 2 < 0 (corresponding to g || > g ⊥ ), the Ising anisotropy induced by the crystal field may be very large (J z /J xy 1), such that the couplings are dominated by a ferromagnetic J z with a common Ising axis for every bond. For positive crystal field ∆ 2 > 0 (corresponding to g ⊥ > g || ), XXZ anisotropy is more limited (J xy /J z ≤ 4), but may still be quite significant. Such materials are generally more desirable for realising strongly bond-dependent couplings, because the XXZ anisotropy is not as dominant.
(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 transition metal ions. This may typically be achieved by matching the electronegativity of the metals and ligands, such that ∆ pd is small. As a general trend, electronegativity of transition metals increases for heavier atoms, while the opposite is true for p-block ligand ions. The combination of heavy metals with heavy ligands results in the most covalent metal-ligand bonds. However, with increased covalency comes increased t 2ge 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 CoCl 2 , CoBr 2 and CoI 2 . For these compounds, we still estimate |t 3 /t 2 | ∼ 4. Thus, it is not clear that Kitaev-dominant exchange can be achieved without extreme fine tuning.
(7) Regarding NCSO and NCTO: Experimental estimates 20,21,43 for NCSO and NCTO suggest ∆ 2 ∼ +4 to +13 meV, which corresponds to a significant easyplane XXZ anisotropy with δΓ, δΓ ∼ 0.2 − 0.5 meV. Large departures from XXZ symmetry are also unlikely, given the small magnitude of the magnon gap (∼ 1 meV) observed in experiment at both the Γ-point 38,44 and the ordering wavevector 38,[42][43][44] . It may further be remarked that the field-evolution of the ESR modes 44 follow expectations for moderate easy-plane XXZ anisotropy. Thus, we may anticipate J 1 < 0, and Γ, Γ > 0, with the sign of K being dependent on the microscopic details.
It may be possible to place some constraints on the interactions on the basis of the ordered moment directions in the zigzag state. Related discussions appear in Ref. 77. For NCTO, it was initially reported that the ordered moments lie nearly in the honeycomb plane, oriented along the direction of the ferromagnetic chains 35 . This orientation was later disputed, with a significant tilting of the ordered moments out of the plane being reported in Ref. 78. We may consider each case separately.
If we consider the zigzag domain with magnetic wavevector parallel to the Z-bond, then moments oriented strictly in-plane would point alongê Z 2 = (1, 1, −2)/ √ 6 in cubic coordinates. The magnetic state would be left invariant under a 180 • rotation around (111) followed by time reversal; it is thus reasonable to assume the couplings possess the same symmetry. This orientation generally requires 79 K, Γ > 0, and places the constraint on the couplings Γ + δΓ = K + Γ + δΓ (i.e. J z± = 0). Taken together, we suggest J 1 = −3, J 3 = +2.5, K = Γ = +0.25, Γ = +0.5 meV as an appropriate starting point for analysis. These considerations favor the model proposed in Ref. 44, and correspond to J xy = −3.25, J z = −2.25, J ±± = −0.125, J z± = 0 meV. The antiferromagnetic sign of K is possible with large t 3 and small t 2 > 0.
In contrast, if the moments are tilted out of plane, then fewer constraints may be placed on the couplings a-priori. However, such tilting likely requires K < 0. In this case, we suggest J 1 = −3, J 3 = +2.5, K = −0.25, Γ = +0.25, Γ = +0.5 meV as an appropriate starting point for analysis. The ferromagnetic sign of K is possible with large t 3 and small t 2 < 0.
(8) Regarding BCAO: The breadth of experimental data on BCAO, in terms of the progression of fieldinduced phases 80 and inelastic neutron data, provide a number of clues towards the magnetic model. While we leave full elaboration for future study 53 , some comments can be made. A recent reinvestigation 47 of the zero-field structure suggested it might better be described by a double stripe ↑↑↓↓ analogue of the zigzag antiferromagnet, with moments oriented nearly along the in-planê e Z 2 direction, as with NCTO. This orientation points to K, Γ > 0. The large discrepancy between in-plane critical fields (0.2, 0.5 T) and the out-of-plane critical field (4T) 50 suggests significant anisotropy. Indeed, the g-tensor appears to satisfy 47 g z ∼ 0.5 g xy , and within an XXZ model, J z ∼ 0.4 J xy . These findings point to significant crystal field effects, with ∆ 2 ∼ 0.2 to 0.25 λ, i.e. ∆ 2 ∼ 15 meV, implying significant Γ and Γ in the global coordinate scheme. In the XXZ scheme, the nearly in-plane moments suggest small J z± , while an apparently small anisotropy between in-plane field directions 47 may place restrictions on J ±± . In contrast, the authors of Ref. 50 have advocated for small J, large K < 0, and small average off-diagonal couplingΓ = (Γ + 2Γ )/3 on the basis of THz spectroscopy experiments. It should be emphasized that these conditions are not mutually compatible: small J z± and J ±± implies small K, and large anisotropy between J xy and J z implies large Γ+2Γ > 0. If we consider BCAO to be in the physical regime of hoppings, our findings contradict the Kitaev-dominant model. We propose a model similar to NCTO with in-plane moments: K is small and likely antiferromagnetic, J < 0 is the dominant coupling, and Γ, Γ > 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 suggested 45,47,48 . These suggestions are compatible with the recent ab-initio estimates 52,53 .

V. CONCLUSIONS
In high-spin Co(II) 3d 7 , the spin-orbital nature of the j eff = 1/2 ground state results in complex bonddependent magnetic couplings in edge-sharing compounds. The theoretical description of these interactions is complicated by the existence of a wide number of different contributions, the balance of which is difficult to predict without detailed knowledge of the microscopic hopping, Coulomb, and crystal field parameters. Based on DFT studies of a range of Co oxides and halides, we find particular discrepancies with previous assumptions regarding these microscopic parameters. Specifically, direct metal-metal hopping is more prominent than ligand-assisted hopping, which results in significant suppression, and possible sign-reversal of the Kitaev inter-actions. Ferromagnetic exchange resulting from e g spins likely strongly dominates the nearest neighbor exchange. Off-diagonal couplings are likely enhanced by crystal field effects, resulting in significant XXZ anisotropy that is consistent with previous description of materials like the Ising ferromagnet CoNb 2 O 6 . Further neighbor exchange is also promoted by the presence of the e g spins. Taken together, Co(II) materials will continue to represent interesting platforms for highly complex anisotropic magnetism. However, the realisation of Kitaev magnetism specifically appears to be unrealistic from the microscopic perspective.