Exchange tensors from broken-symmetry DFT

The front end of the magnon pipeline: how to turn a spin-polarised electronic-structure calculation into the full bilinear spin-interaction tensors that a spin model needs — not just a scalar Heisenberg $J_{ij}$, but the antisymmetric Dzyaloshinskii–Moriya vector ${\mathbf{D}}_{ij}$ and the symmetric-traceless anisotropic exchange ${\mathbf{\Gamma }}_{ij}$ as well, plus the on-site single-ion anisotropy. The downstream consumer is LSWT, which takes exactly these tensors and returns the magnon dispersion; this note is the inverse problem — extracting them from DFT — and the systematics (functional, $U$, SOC, supercell) that decide whether the answer is trustworthy.

Prerequisites

Green’s functions & LKAG (the magnetic force theorem) · Dzyaloshinskii–Moriya (the antisymmetric part and Moriya’s rules) · the magnetic-anisotropy concept (single-ion $D(S^z{)}^2$ and anisotropic exchange; plain text — no dedicated note on disk) · spin-polarised / non-collinear SCF

The target: decomposing the $3\times 3$ exchange tensor

Fix the spin model once. The most general bilinear interaction between two localised moments is

$$H=\sum _{ij}{\mathbf{S}}_i^{\mathsf{T}}{\mathbf{J}}_{ij}{\mathbf{S}}_j+\sum _i{\mathbf{S}}_i^{\mathsf{T}}{\mathbf{A}}_i{\mathbf{S}}_i,$$

with ${\mathbf{J}}_{ij}$ a $3\times 3$ real tensor (Cartesian spin indices) and ${\mathbf{A}}_i$ the single-ion anisotropy matrix. Any such tensor splits uniquely into three irreducible pieces, and each piece is a distinct physical interaction:

$${\mathbf{J}}_{ij}=\underset{\text{isotropic}}{\underbrace{J_{ij}1}}+\underset{\text{antisymmetric (DMI)}}{\underbrace{{\mathbf{J}}_{ij}^{\mathrm{A}}}}+\underset{\text{symmetric traceless}}{\underbrace{{\mathbf{\Gamma }}_{ij}}}.$$

• Isotropic Heisenberg — the trace: $$J_{ij}=\frac{1}{3}\mathrm{Tr}{\mathbf{J}}_{ij},$$ contributing $J_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j$. The only part that survives without spin–orbit coupling.
• Antisymmetric → DMI — the antisymmetric part, dualised to the axial vector ${\mathbf{D}}_{ij}$ exactly as in Dzyaloshinskii–Moriya (same convention used throughout these notes): $${\mathbf{J}}_{ij}^{\mathrm{A}}=\frac{1}{2}({\mathbf{J}}_{ij}-{\mathbf{J}}_{ij}^{\mathsf{T}}),D_{ij}^{\gamma }=\frac{1}{2}{ϵ}^{\gamma \alpha \beta }{J}_{ij}^{\alpha \beta },$$ giving ${\mathbf{D}}_{ij}\cdot ({\mathbf{S}}_i\times {\mathbf{S}}_j)$.
• Symmetric traceless → anisotropic / pseudo-dipolar exchange ${\mathbf{\Gamma }}_{ij}$: $${\mathbf{\Gamma }}_{ij}=\frac{1}{2}({\mathbf{J}}_{ij}+{\mathbf{J}}_{ij}^{\mathsf{T}})-\frac{1}{3}\mathrm{Tr}({\mathbf{J}}_{ij})1,$$ the two-site anisotropic exchange (the off-diagonal Kitaev-type and pseudo-dipolar terms). It is the part that, together with single-ion ${\mathbf{A}}_i$, supplies the magnon gap that Mermin–Wagner says a 2D magnet needs.

The hierarchy in SOC (from the DMI mechanism): $J\sim \mathcal{O}({\lambda }^0)$, $\mathbf{D}\sim \mathcal{O}(\lambda )$, $\mathbf{\Gamma }\sim \mathcal{O}({\lambda }^2)$ — which is also the order in which they are hard to converge.

Route (a): LKAG / the magnetic force theorem (Green’s functions, TB2J)

The Liechtenstein–Katsnelson–Antropov–Gubanov (LKAG) approach computes the exchange parameters as the energy cost of infinitesimal rotations of already-converged local moments, using the magnetic force theorem: to first order, the energy change under a small perturbation is the change in the sum of one-electron (band) energies at fixed self-consistent potential, so no expensive re-SCF of rotated states is needed. Tilt moment $i$ by $\delta {\theta }_i$ and moment $j$ by $\delta {\theta }_j$; the second variation of the band energy with respect to the two tilt angles is the exchange coupling. In the multiple-scattering / Green’s-function formulation this becomes a closed expression — an energy integral up to $E_F$ of a trace over intersite Green’s functions and the on-site exchange splitting (the LKAG formula),

$$J_{ij}=\frac{1}{2\pi }\mathrm{Im}{\int }^{E_F}d\epsilon {\mathrm{Tr}}_L[{\Delta }_i{G}_{ij}^{\uparrow }(\epsilon ){\Delta }_j{G}_{ji}^{\downarrow }(\epsilon )],$$

with ${\Delta }_i$ the local exchange-splitting (the on-site spin potential) and $G_{ij}^{\sigma }$ the spin-resolved intersite Green’s function — the construction set out in Green’s functions & LKAG. The strengths are exactly what you want for a magnon pipeline:

• it returns all $J_{ij}$ for all neighbour shells from one calculation of the reference (usually ferromagnetic or the DFT ground) state — no need to design a configuration per coupling;
• with SOC retained, the Green’s functions and the on-site vertices become $2\times 2$ in spin, the trace generalises to the full $3\times 3$ tensor, and the antisymmetric (${\mathbf{D}}_{ij}$) and symmetric-anisotropic (${\mathbf{\Gamma }}_{ij}$) parts emerge directly as the corresponding tensor components — there is a closed LKAG-type expression for ${\mathbf{D}}_{ij}$;
• it is the engine of TB2J, which post-processes Wannier (or localised-orbital) Hamiltonians from a DFT code and emits a SpinW/UppASD-ready tensor set.

Its main weakness is that it is a perturbative, ground-state estimate (infinitesimal rotations about one reference state), so it inherits any error in that state and in the localised basis used to define the moments — see the caveats below.

Route (b): energy mapping (four-state / broken-symmetry total energies)

The complementary route is non-perturbative: compute total energies of several finite spin configurations and fit them to the spin model. The cleanest version is the four-state (energy-mapping) method (Xiang et al.): to extract one component of ${\mathbf{J}}_{ij}$ between sites $i$ and $j$, run four broken-symmetry SCF calculations with the pair $(i,j)$ set to the four orientation combinations $(\uparrow \uparrow ),(\uparrow \downarrow ),(\downarrow \uparrow ),(\downarrow \downarrow )$ along chosen Cartesian axes (all other moments held in a fixed common reference), and take the combination

$$J_{ij}^{\alpha \beta }=\frac{E^{\uparrow \uparrow }+E^{\downarrow \downarrow }-E^{\uparrow \downarrow }-E^{\downarrow \uparrow }}{4S^2},$$

which cancels the single-site and reference-dependent contributions and isolates the two-site coupling. Looping the chosen axes $(\alpha ,\beta )$ over $x,y,z$ reconstructs the full $3\times 3$ tensor; the antisymmetric combination then gives ${\mathbf{D}}_{ij}$ and the symmetric-traceless one gives ${\mathbf{\Gamma }}_{ij}$. Single-ion ${\mathbf{A}}_i$ is obtained from a separate set rotating one moment. Compared with the alternative of fitting a basis of collinear/non-collinear configurations to a least-squares model, the four-state combination is local and avoids contamination from distant couplings — at the price of many SCF runs and a supercell large enough that the manipulated pair does not interact with its periodic images.

When to use which. LKAG is fast, gives the whole shell structure at once, and is natural with a Wannier interface (TB2J); energy mapping is non-perturbative, makes no infinitesimal-rotation assumption, and is the natural choice with a localised-orbital total-energy code (e.g. CRYSTAL23) where non-collinear total energies are cheap and clean. The two should agree where both apply, and a disagreement is itself diagnostic (usually of basis/localisation sensitivity or of an unstable reference state).

The role of SOC, and single-ion vs. two-site anisotropy

The decomposition above is degenerate without spin–orbit coupling: with SOC switched off, the spin and lattice frames decouple, the energy depends only on relative spin angles, and the calculation can only ever return the isotropic $J_{ij}$ — ${\mathbf{D}}_{ij}=0$ and ${\mathbf{\Gamma }}_{ij}=0$ identically, and the single-ion ${\mathbf{A}}_i$ collapses to a constant. Turning on SOC is what makes the energy depend on the orientation of the spins relative to the lattice, and only then do ${\mathbf{D}}_{ij}$ ($\mathcal{O}(\lambda )$) and ${\mathbf{\Gamma }}_{ij}$, ${\mathbf{A}}_i$ ($\mathcal{O}({\lambda }^2)$) appear. This is the computational restatement of Moriya’s “DMI is first order in SOC.”

A point that the extraction must keep straight is the single-ion vs. two-site distinction, because the two are physically and methodologically different anisotropies that are easy to conflate:

• Single-ion anisotropy ${\mathbf{A}}_i$ ($D(S_i^z{)}^2$ and friends) is an on-site term — the energy of one moment as a function of its own direction in the crystal field. It vanishes identically for spin-$\frac{1}{2}$ ions (where $(S^z{)}^2$ is a constant), and it is extracted by rotating a single moment. It dominates the easy-axis physics in many real 2D magnets.
• Two-site anisotropic exchange ${\mathbf{\Gamma }}_{ij}$ is a bond term, extracted from the symmetric off-diagonal tensor components of a pair. It is present even for spin-$\frac{1}{2}$.

Both feed the magnon gap, but they scale differently with $S$ and with the number of neighbours, and a fit that lumps them together will mis-extrapolate. Keeping ${\mathbf{A}}_i$ and ${\mathbf{\Gamma }}_{ij}$ as separate fit targets (single-moment vs. pair rotations) is the clean way to disentangle them.

Caveats — where the numbers go wrong

The tensors inherit every systematic error of the underlying electronic structure; the extraction method is the easy part. The recurring offenders:

• Functional / +U dependence. The mapping is only as good as the description of moment localisation and hybridisation. The on-site exchange splitting ${\Delta }_i$ and the intersite Green’s functions both depend strongly on the DFT+U value (or the exact-exchange fraction in a hybrid): too small a $U$ over-delocalises the moment and inflates the superexchange $J$, too large a $U$ suppresses hybridisation and shrinks it; the antisymmetric $\mathbf{D}$ and anisotropic $\mathbf{\Gamma }$ are even more sensitive because they ride on the SOC-induced corrections to the same hybridisation. $U$ is therefore not a free knob here — it has to be fixed by an independent criterion, not tuned to the magnetism you are trying to predict.
• Projector / basis dependence. Both routes require defining “the moment on site $i$,” i.e. a projection onto local orbitals (Wannier functions, LCAO basis, PAW/atomic spheres). The extracted $J_{ij}$ depends on that choice — different projectors partition the intersite charge/spin differently — and the dependence is worse for the small anisotropic components. Convergence of the result with respect to the localisation scheme is a real check, not a formality.
• Supercell / image convergence. Energy mapping needs the manipulated moments far enough from their periodic images that the spurious image coupling is below the (small) $\mathbf{D}$/$\mathbf{\Gamma }$ being measured; LKAG needs enough $\mathbf{k}$-points to converge the $E_F$ energy integral. Under-converged cells fake anisotropy.
• The (meta)stability assumption. Both routes presuppose that the reference magnetic state is a genuine (meta)stable solution of the SCF and that the spin model is an adequate parameterisation. If the assumed ground state is actually unstable, the extracted tensors are meaningless — and the LSWT consumer will flag it downstream as a non-positive-definite grand dynamical matrix (the Colpa stability test in LSWT). The two checks are complementary: get a stable SCF state and a positive-definite magnon Hamiltonian.

As example systems (qualitative use only): the layered chromium magnets CrI${}_3$ and CrSBr are the standard testbeds for the energy-mapping route — both are insulating, have well-localised Cr moments, and are the 2D magnets whose anisotropy-stabilised order (CrSBr) is the very thing the extracted $\mathbf{D}$, $\mathbf{\Gamma }$, $\mathbf{A}$ are supposed to capture for the Mermin–Wagner gap. I deliberately quote no computed $J$, $D$, or anisotropy numbers here — those are functional- and projector-dependent outputs, exactly the quantities this note argues are not universal.

Builds toward: LSWT

Key references

• A. I. Liechtenstein, M. I. Katsnelson, V. P. Antropov & V. A. Gubanov, “Local spin density functional approach to the theory of exchange interactions in ferromagnetic metals and alloys,” J. Magn. Magn. Mater. 67, 65 (1987) — the LKAG magnetic force theorem.
• X. He, N. Helbig, M. J. Verstraete & E. Bousquet, “TB2J: A python package for computing magnetic interaction parameters,” Comput. Phys. Commun. 264, 107938 (2021) — the LKAG/Green’s-function implementation that emits the full tensor set.
• H. J. Xiang, E. J. Kan, S.-H. Wei, M.-H. Whangbo & X. G. Gong, “Predicting the spin-lattice order of frustrated systems from first principles,” Phys. Rev. B 84, 224429 (2011) — the four-state / energy-mapping method, including DMI extraction.