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

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:

$\displaystyle\mathcal{H}=\mathcal{H}_{\rm hop}+\sum_{i}\mathcal{H}_{i}$

(1)

where $\mathcal{H}_{\rm hop}$ denotes the hopping between Co sites. $\mathcal{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:

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

(2)

The SOC term is written:

$\displaystyle\mathcal{H}_{\rm SOC}=\lambda\mathbf{L}_{i}\cdot\mathbf{S}_{i}$

(3)

where the atomic SOC constant for Co is roughly $\lambda_{\rm Co}\approx 60$ meV. 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}$

(4)

where $\alpha,\beta,\gamma,\delta$ label different $d$-orbitals. ¹¹1The 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 formGeorges et al. (2013); 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)$		(5)
	$\displaystyle J_{t2g}=\frac{3}{49}F_{2}+\frac{20}{441}F_{4}$		(6)

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. We also take the approximate ratio $F_{4}/F_{2}=5/8$, which is appropriate for $3d$ elementsPavarini (2014).

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}$

(7)

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)$

(8)

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:

$\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)$

(14)

where $\Delta_{1}$ is the $t_{2g}$-$e_{g}$ splitting, and $\Delta_{2}$ is the trigonal term. Trigonal distortions are relevant for BCAO, BCPO, NCTO, and NCSO. Generally, $\Delta_{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\Delta_{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:

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

(15)

For an ideal edge-sharing bond with $90^{\circ}$ metal-ligand-metal 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:

$\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)$

(21)

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.3.

Previous worksLiu and Khaliullin (2018); Liu (2021); Liu et al. (2020); Sano et al. (2018) 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 ($\Delta_{2}\ll\lambda$) and dominant ligand-assisted hopping ($t_{2}\gg 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.

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 $\mathcal{H}_{i}$ typically satisfy $J_{H},\Delta_{1}\gg\lambda,\Delta_{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 $\mathcal{H}_{i}$ into this manifold, and diagonalizing. This procedure was employed in Ref. Lines, 1963, and is repeated in this section for the purpose of discussion. More detailed discussions can be found in Ref. Lines, 1963; Liu, 2021.

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 low energy multiplets into $J_{\rm eff}=1/2$, 3/2, and $5/2$ states. For $\Delta_{2}=0$, the multiplet energies satisfy:

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

Given $\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).

For finite $\Delta_{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_{\rm eff}=1/2$ doublet as $\Delta_{2}$ vanishes. 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$		(24)
	$\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$		(25)

where the coefficients $c_{n}$ vary with $\Delta_{2},\lambda$ as shown in Fig. 4(b). More detailed expressions for the $|m_{L},m_{S}\rangle$ 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. Lines, 1963; for $\Delta_{2}=0$, the coefficients are $c_{1}=1/\sqrt{2}$, $c_{2}=1/\sqrt{3}$, $c_{3}=1/\sqrt{6}$. From Fig. 4(b), it is evident that the composition of the lowest doublet changes dramatically with finite $\Delta_{2}$.

In the limit of large trigonal elongation $\Delta_{2}>0$, the wavefunction coefficients converge to $c_{2}\to 1$ and $c_{1},c_{3}\to 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}=\pm 1/2$ states lie slightly below the $m_{s}=\pm 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_{\perp}>g_{||}$, where $g_{||}$ refers to the component along the trigonal $(111)$ axis.

For the opposite case of trigonal compression $\Delta_{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_{\rm eff}=1/2$. As depicted in Fig. 4(a), the effect of SOC is then to split the $S=3/2$, $L_{\rm 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}\to 1$, $c_{2},c_{3}\to 0$, implying that it is composed of pure $m_{s}=\pm 3/2$ states, according to Eq. (10, 11). Consistently, the $g$-tensor satisfies $g_{||}\gg g_{\perp}$ in this limit (Fig. 4(c)).

Finally, it should be emphasized that an effective $J=\frac{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 $\Delta_{2}<0$, 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 even for large $\Delta_{2}<0$. For $\Delta_{2}>0$, which is expected to be relevant for BCAO, BCPO, NCTO, and NCSO, the range of validity is roughly constrained to $\Delta_{2}<\lambda/2$. For larger values of $\Delta_{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 $\Delta_{2}$ for various Co(II) compounds may be found in Ref. Liu, 2021. All of these materials are expected to be well described by a $J=\frac{1}{2}$ model, which can be verified experimentally by the observation of the well-defined local $j_{1/2}\to j_{3/2}$ excitations in an energy range larger than the dispersion of magnetic excitations Sarte et al. (2018); Ross et al. (2017); Songvilay et al. (2020); Kim et al. (2021). 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.

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 metal-metal vs. ligand-assisted hopping processes. In terms of the usual Slater-Koster hoppings, the downfolded $d$-only hoppings are given approximately by:

	$\displaystyle t_{1}\approx$	$\displaystyle\ \frac{1}{2}t_{dd}^{\pi}$		(26)
	$\displaystyle t_{2}\approx$	$\displaystyle\ -\frac{1}{2}t_{dd}^{\pi}+\frac{(t_{pd}^{\pi})^{2}}{\Delta_{pd}}$		(27)
	$\displaystyle t_{3}\approx$	$\displaystyle\ \frac{3}{4}t_{dd}^{\sigma}$		(28)
	$\displaystyle t_{4}\approx$	$\displaystyle\ -\frac{1}{2}t_{dd}^{\sigma}$		(29)
	$\displaystyle t_{5}\approx$	$\displaystyle\ t_{dd}^{\pi}$		(30)
	$\displaystyle t_{6}\approx$	$\displaystyle\ \frac{\sqrt{3}}{4}t_{dd}^{\sigma}-\frac{t_{pd}^{\pi}t_{pd}^{\sigma}}{\Delta_{pd}}$		(31)

where $\Delta_{pd}$ is the charge-transfer energy from Co $d$ to the ligand $p$ orbitals. Here, we have ignored contributions from $t_{dd}^{\delta}$. 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 $\propto t_{pd}/\Delta_{pd}$).

In order to investigate the hopping trends with bond-length dependence, we first considered hypothetical cubic CoO (NaCl type; $Fm\bar{3}m$) 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 FPLOKoepernik and Eschrig (1999); Opahle et al. (1999) at the GGA (PBEPerdew et al. (1996)) 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}^{\sigma}$, which quantifies the hopping between ligand $p$ and metal $e_{g}$ orbitals. This is followed in magnitude by $t_{pd}^{\pi}$, which quantifies the hopping between ligand $p$ and metal $t_{2g}$ orbitals. The difference in magnitude between $t_{pd}^{\sigma}$ and $t_{pd}^{\pi}$ 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 $\Delta_{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 $\Delta_{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 $\Delta_{1}$ for $3d$ transition metal compounds Gray (1965), 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 compatible with (but not exactly given by) Eq’ns (26) - (31) with $\Delta_{pd}\sim 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: CoNb₂O₆ (Ref. Sarvezuk et al., 2011), BaCo₂(AsO₄)₂ (Ref. Dordević, 2008), Na₃Co₂SbO₆ (Ref. Songvilay et al., 2020), and Na₂Co₂TeO₆ (Ref. Xiao et al., 2019)²²2The Na₂Co₂TeO₆ 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 Na₃Co₂TeO₆. It is expected this change in the filling should have minimal impact on the computed hoppings., as well as triangular lattice compounds CoCl₂, CoBr₂ and CoI₂ according to crystal structures from Ref. Wilkinson et al., 1959; Wyckoff and Wyckoff, 1963. These edge-sharing Co(II) compounds span a wide range Co-Co distances, e.g. from $\sim 2.9$ Å in BaCo₂(AsO₄)₂Dordević (2008) to $\sim 3.9$ Å in CoI₂Wyckoff and Wyckoff (1963). For each case, we employed fully relativistic density functional theory calculations performed with FPLOKoepernik and Eschrig (1999); Opahle et al. (1999) at the GGA (PBEPerdew et al. (1996)) 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}|\sim|t_{1}|,|t_{4}|,|t_{5}|\lesssim 0.05$ eV. We find that this applies equally to Co oxides and halides. A similar situationWellm et al. (2021) was recently proposed for Na₂BaCo(PO₄)₂. On this basis, we conclude: for edge-sharing 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.

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)$		(32)

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. This definition turns out to be convenient for the case of trigonal splitting $\Delta_{2}=0$. For finite $\Delta_{2}$, following Ref. Lines, 1963, 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: $\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}$		(33)
		$\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)$		(34)
	$\displaystyle J_{z}=$	$\displaystyle\ J+\frac{1}{3}\left(K+2\Gamma+4\Gamma^{\prime}\right)$		(35)
	$\displaystyle J_{\pm\pm}=$	$\displaystyle\ -\frac{1}{6}\left(K+2\Gamma-2\Gamma^{\prime}\right)$		(36)
	$\displaystyle J_{z\pm}=$	$\displaystyle\ -\frac{\sqrt{2}}{3}\left(K-\Gamma+\Gamma^{\prime}\right)$		(37)

A similar parameterization was also employed in Ref. Liu, 2021; Liu et al., 2020. Both parameterizations are employed below, where appropriate.

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 $\mathcal{O}[\ (t_{pd}^{2}/\Delta_{pd})^{2}/U_{dd}\ ]$ to the magnetic exchange, where $\Delta_{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_{\rm eff}=t_{pd}^{2}/\Delta_{pd}$ as an effective ligand-mediated $d-d$ hopping, as discussed above. In this case, the magnetic couplings arise at order $\mathcal{O}[\ t_{\rm eff}^{2}/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 $\mathcal{O}[\ (t_{pd}^{2}/\Delta_{pd})^{2}/(2\Delta_{pd})\ ]$. In the downfolded $d$-only basis, these processes are not distinguished from those in Fig. 6(a). Rather, they reflect the fact that the Coulomb repulsion within the $d$-orbital Wannier functions is renormalized to $\tilde{U}^{-1}\sim U_{dd}^{-1}+(2\Delta_{pd})^{-1}$ via hybridization of the $d$- and $p$-orbitals. Thus both the “antiferromagnetic” superexchange and cyclic exchange processes are captured, in principle, by exchange contributions at order $\mathcal{O}[\ t_{\rm eff}^{2}/\tilde{U}\ ]$, assuming the hoppings and Coulomb repulsion terms are suitably renormalized.

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. 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)$		(38)
	$\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)$		(39)
	$\displaystyle\gamma=$	$\displaystyle\ \frac{40J_{H}^{p}}{81(\Delta_{pd}+U_{p}/2)^{2}}$		(40)

where $J_{H}^{p}$ is the Hund’s coupling at the ligand, $U_{p}$ is the excess Coulomb repulsion at the ligand ion, $\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. 5(a). 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 reason for the discrepancy in magnitudes between $\delta J$ and $\delta K$ is that the dominant contribution to $\delta J$ comes from hopping between the $e_{g}$ and $p$ orbitals (depicted in Fig. 6(c)), as parameterized by the large hopping $t_{pd}^{\sigma}$. This process does not contribute to $K$. The small correction $\delta 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 (24, 25). This procedure is described in more detail in 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}$ for small values of $t_{n}/U$ (with all other contributions to the Hamiltonian, which appear in $\mathcal{H}_{i}$, being treated exactly). The reason for adopting this method is that analytical perturbation theory is generally challenging for finite $\Delta_{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 $\delta J$ must be added to the computed couplings to estimate the omitted ligand exchange processes depicted in Fig. 6(c).

In this section, we first consider the magnetic couplings for the case $\Delta_{2}=0$, for which $\Gamma^{\prime}=0$ strictly. Up to second order in $d-d$ hopping (and fourth order in $d-p$ hopping), the couplings may be written:

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

where:

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

(44)

and $\mathbb{M}$ parameterizes the contributions that are captured by (downfolded) $d-d$ hopping. The matrices are a function of $F_{n},\lambda,\Delta_{n}$. To estimate $\mathbb{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 $\mathcal{H}_{U}+\mathcal{H}_{\rm CFS}+\mathcal{H}_{\rm SOC}+\mathcal{H}_{\rm hop}$ for two metal sites, and fit the resulting couplings. For this purpose, 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. 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)$

(52)

$\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)$

(60)

$\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)$

(68)

in units of $10^{-3}$/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$ (see section II.3). 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)$ is typically the same as $-\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. 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.3. 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 $\Gamma=0$. If we ignore the ferromagnetic correction $\delta J$, then $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 super-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$.

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.

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.

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}$		(69)
	$\displaystyle\delta J_{z}\approx$	$\displaystyle\ u_{z}^{2}\ \delta J_{0}$		(70)
	$\displaystyle u_{xy}=$	$\displaystyle\ \frac{3}{5}\left(2\sqrt{3}c_{1}c_{3}+2c_{2}^{2}\right)$		(71)
	$\displaystyle u_{z}=$	$\displaystyle\ \frac{3}{5}\left(1+2(c_{1}^{2}-c_{3}^{2})\right)$		(72)

with $c_{n}$ given in Fig. 4(b) and $\delta J_{0}$ defined according to eq’n (38). Thus, 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}$		(73)
	$\displaystyle\delta\Gamma=$	$\displaystyle\ \delta\Gamma^{\prime}=\frac{1}{3}\left(u_{z}^{2}-u_{xy}^{2}\right)\delta J_{0}$		(74)

Here, we have continued to neglect the small corrections $\delta K_{0}$. As with the undistorted case, the corrections to the anisotropic couplings $J_{\pm\pm}$ and $J_{z\pm}$ 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 $\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. 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 $\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. 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 $\delta J$. The results are as follow:

Trigonal compression: For $\Delta_{2}<0$, as shown in Fig. 9(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}$ are induced, with $\text{sign}(\Gamma,\Gamma^{\prime})$ typically being the same as $\text{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 $\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})$ typically the same as $-\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$.

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}^{\pi}$), the $e_{g}$-ligand hybridization may still play a significant role in the magnetic exchange due to the large $t_{pd}^{\sigma}$. This fact is what leads to significant ferromagnetic corrections $\delta 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}^{\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}$		(75)
	$\displaystyle K_{3}\approx\Gamma_{3}\approx 0$		(76)

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 $\Delta_{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₃CoSb₂O₉, for which the closest Co²⁺ neighbors have the same geometry as depicted in Fig. 10. In this case, it is well establishedMa et al. (2016); Ito et al. (2017); Ghioldi et al. (2018); Kamiya et al. (2018) 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.