Higher-order exchange interactions in two-dimensional magnets

Magnetism in recently discovered van der Waals materials has opened new avenues in the study of fundamental spin interactions in truly two-dimensions. A paramount question is what effect higher-order interactions beyond bilinear Heisenberg exchange have on the magnetic properties of few-atom thick compounds. Here we demonstrate that biquadratic exchange interactions, which is the simplest and most natural form of non-Heisenberg coupling, assume a key role in the magnetic properties of layered magnets. Using a combination of nonperturbative analytical techniques, non-collinear first-principles methods and classical Monte Carlo calculations that incorporate higher-order exchange, we show that several quantities including magnetic anisotropies, spin-wave gaps and topological spin-excitations are intrinsically renormalized leading to further thermal stability of the layers. We develop a spin Hamiltonian that also contains antisymmetric exchanges (e.g. Dzyaloshinskii-Moriya interactions) to successfully rationalize numerous observations currently under debate, such as the non-Ising character of several compounds despite a strong magnetic anisotropy, peculiarities of the magnon spectrum of 2D magnets, and the discrepancy between measured and calculated Curie temperatures. Our results lay the foundation of a universal higher-order exchange theory for novel 2D magnetic design strategies.

calculations that incorporate higher-order exchange, we show that several quantities including magnetic anisotropies, spin-wave gaps and topological spin-excitations are intrinsically renormalized leading to further thermal stability of the layers. We develop a spin Hamiltonian that also contains antisymmetric exchanges (e.g. Dzyaloshinskii-Moriya interactions) to successfully rationalize numerous observations currently under debate, such as the non-Ising character of several compounds despite a strong magnetic anisotropy, peculiarities of the magnon spectrum of 2D magnets, and the discrepancy between measured and calculated Curie temperatures. Our results lay the foundation of a universal higher-order exchange theory for novel 2D magnetic design strategies.

Main
The finding of magnetism in atomically thin van der Waals (vdW) materials 1,2 has attracted an increasing amount of interest in the investigation of novel magnetic phenomena at the nanoscale 3 . Important in the understanding and applications of 2D vdW magnets in real technologies is the elucidation of fundamental interactions that determine their magnetic properties, such as the exchange interactions between the spins. These interactions govern the magnetic ordering and can be either symmetric, which determine collinear ferromagnets and anti-ferromagnets, or antisymmetric, which promote topological non-trivial spin textures, e.g. skyrmions. Higher-order exchange terms involving the hopping of two or more electrons play a pivotal role in the spin-ordering of low-dimensional nanostructures.
For instance, biquadratic (BQ) exchange interactions are critical in the elucidation of the magnetic features of several systems, such as multilayer materials 4 , perovskites 5 , iron-based superconductors 6,7 , iron tellurides 8 and oxides 9 . Indeed, in materials where the exchange is for some reason weak 10 (e.g. low-temperature magnets) BQ exchange has a particularly strong influence being the case of several 2D magnets discovered so far 11,12 .
Here we identify that several families of 2D magnets including metals, insulators and small band gap semiconductors, develop substantially large BQ exchange interactions. The delicate interplay between the superexchange process through the non-magnetic atom and the Coulomb repulsion at neighboring spin-sites involving more than one electron induces sizable corrections to the total energy of the systems beyond bilinear spin models, e.g. Kitaev, Heisenberg, Ising. We developed a generalized non-Heisenberg spin Hamiltonian that includes both BQ and Dzyaloshinskii-Moriya interactions (DMI), providing a universal picture of the spin properties of 2D vdW magnets. We show that while BQ exchange interactions give substantial contributions to the thermal properties, such as in critical temperatures (T c ), non-collinear spin effects via DMI are negligible. The model also offers insights via analytical equations into the spin excitations of 2D materials as a general formalism is developed for materials that hold honeycomb crystal structure, ferromagnetism and out-of-plane easy axis in a universal basis. Our results provide the conceptual framework to understand a variety of 2D magnetic materials in a consistent way explaining numerous experimental observations, including controversies on the magnetic properties of CrI 3 or CrBr 3 .

Biquadratic exchange interactions in 2D magnets
In order to calculate the different exchange contributions to the total energy at the level of first-principles methods (see Supplementary Sections S1-S4 for details), we map the angular dependence of the spins S j = µ s s j , where µ s is the magnetic moment and |s j |= 1, in the unit cell of the layered material (Fig. 1a). We rotate the spins by θ between two known spin configurations: from ferromagnetic (FM) at θ = 0, to anti-ferromagnetic (AFM) at θ = 180 o .
Small steps in θ generate a path of quasi-continuously configurations where both energy and magnetization are allowed to relax self-consistently without any fixed constraint on the direction. The resulting curve in energy includes contributions from bilinear (BL) exchanges up to higher-order terms, e.g. BQ exchange interactions. We quantify the contribution of each kind of interactions using the following non-Heisenberg spin Hamiltonian: where J ij and λ ij are the isotropic and anisotropic BL exchanges between spins S i and S j on atomic sites i and j; D i is the on-site anisotropy with easy axis e i ; and K ij is the BQ exchange interactions which is due to electron hopping between two adjacent sites 13 . We restrict the discussions to first nearest-neighbor BQ interactions, that is, K ij = K bq , and K ij = 0 otherwise. This assumption has been shown to be sufficient to study a variety of nanomagnetic systems with higher-order exchanges 5,6,8,9 . We can write Eq.1 in a similar form separating the terms with angular and non-angular dependence as E tot  1T ). Surprisingly, as the spins are spatially rotated away a sizable deviation from a cos θ-like behavior characteristic of BL exchange interactions (shaded area) is observed (Fig. 1b-e). The overall trend seems independent of the stoichiometric formula or the atomic elements composing the structure. The universal character of the BQ exchange interactions can be appreciated clearer in Fig. 1f where a quadratic regression (Q-Regression) can be undertaken over the computed dataset. Materials with alike chemical environment (e.g. bond lengths, electron affinity, binding energy) such as CrI 3 , CrBr 3 and CrCl 3 present close variation of the energy and consequently similar magnitudes of the BQ exchange (Table 1). Other compounds tend to gain energy with the spin rotation and stabilize in a different magnetic coupling, for instance, CrF 3 (Fig. 1b) becomes AFM ordered. This is also the case of Mn-based compounds, 2H-MnS 2 , 2H-MnSe 2 , MnPS 3 , MnPSe 3 (Fig. 1d-e), which agreed with magnetometry measurements [14][15][16]  These results suggest that BL models are insufficient to describe the magnetic features of 2D vdW magnets.

A Hubbard-based microscopic model
To understand the microscopic mechanism of the BQ exchange in 2D materials, we will use the following 2D half-filled Hubbard Hamiltonian applied to the honeycomb lattice ( Fig. 2ab): where indices i and j denote the lattice sites, the d-spin states on the metal atoms (M A,B ) are labeled as σ =↑, ↓, the sum < ij > is over the nearest neighbors and t ef f is the effective nearest neighbor hoping between M A,B ions. t ef f is due to the hybridization between 3d n and np orbitals (n = 2, 3, 4, 5 depending on the atomic element involved) at M A,B and X atoms, respectively (Fig. 2b). The direct hopping between M A,B scales with t dd ∼ r −5 and therefore for relatively big distances can be neglected 13 . U > 0 is the non-negative on-site Coulomb is the creation (annihilation) operator for a fermion with spin σ at site where J (2) FM and J (2) AFM are the exchange energies for FM and AFM coupling at second-order, and the BL exchange is defined as AFM . U ex is an energy correction for the internal spin exchange whether a spin flip is required during the hopping between M A,B sites ( Fig. 2b) 22 . Such term can be used to stabilize or destabilize spin transfer through the M A,B −X covalent bonds as the exchange occurs 13,22 . If the electron-hopping is to an occupied orbital of the neighboring site, U ex will favor AFM alignment of the two spins via a superexchange interaction (see Supplementary Section S9 for details). However, if the electron-hopping is to an unoccupied or a virtual state, a FM alignment will be favored by U ex 23 . The competition between the amount of energy U ex to stabilize a specific coupling and the Coulomb repulsion U between the electrons at an energy state may compensate each other leading to a small value of J bl . Indeed, several 2D magnets have shown low-temperature magnetism 1,2,12 which is directly related with the small magnitude of the exchange interactions. In this case, it is necessary to extend the perturbation in t ef f /U to fourth-order in Eq.2 involving at least one electron from both M A,B sites which resulted in 13 : AFM being the BQ exchange energy. Both Eqs.3−4 show that a competition between FM and AFM couplings takes place once the electrons are hopping between different spin sites. The stabilization of one or another magnetic order is determined by several factors such as the ligand-field splitting ∆ 0 between t 2g and e g states in the metal atom in the honeycomb lattice. As the filling of both type of states determines the magnitude of the 3d − sp hybridization between metals and ligands, we can approach ∆ 0 ≈ U − U ex being proportional to the bandgap of the material 10 . It has been shown that the role of high-order exchange interactions increases on reduction of the bandgap mediated by the non-magnetic atoms which is related to the ratio of BQ and BL exchanges 21,24 . If we divide Eq.4 and Eq.3 and use the definition of ∆ 0 , we can write a direct relationship between the exchange interactions and the bandgap for 2D magnets as: This equation can be understood as a direct interplay between the Coulomb repulsion and the hopping of electrons between different sites subjected to the crystal field or bandgap of the material. We can consider two situations in Eq.5: ∆ 0 → 0 and ∆ → 2U which correspond to small and large bandgap materials, respectively. This resulted in: (7) Figure 2c shows the variation of K bq /J bl as a function of ∆ 0 for the core of materials displaying BQ exchange interactions. Strikingly, both Eqs.6−7 correctly describes the overall behavior observed in our simulations. Materials with similar bonding environment, for instance, either in terms of Cr (Mn) atoms follow an increase (decrease) of K bq /J bl with the bandgap, respectively. It is worth mentioning that as the compounds tend to a AFM spin alignment 10,23 , i.e. CrF 3 , they increase the value of K bq /J bl being the case within the Crbased trihalide family (CrX 3 , X=F, Cl, Br, I). There is also an abrupt change in behavior for narrow bandgap materials and metals with no dependence on ∆ 0 (inset in Fig. 2c). These results indicate that one can design the amount of BQ exchange in a 2D magnet tuning its bandgap, for instance using an electric bias as recently used in CrI 3 25 .

Thermal effects in 2D magnets
To verify whether BQ exchange interactions will have any effect on the thermal properties of For atomistic spin dynamics simulations we calculate the effective magnetic field B i by taking the first-derivative of Eq.8 on the different spin components: The effective field is then included within the total effective field describing the time evolution of each atomic spin using the stochastic Landau-Lifshitz-Gilbert equation 27 . We calculate up to third nearest-neighbors λ ij and J ij in Eq.1 on a representative set of 2D magnets, e.g. Cr-based trihalide family (Table 2). Supplementary Section S10 gives a thorough discussion on the calculation of the exchange parameters. Figure  closely those measured for both compounds with almost no difference (Fig. 3a-d).

Enhancement of magnetic stability
To understand the intrinsic effect of higher-order exchange interactions in the thermal features of 2D magnetic materials and provide a consistent description of the underlying spin interactions, we have developed an analytical model based on spin-wave theory 31,32 . We generalized Eq.1 to second order contributions on the magnetic anisotropies which resulted in: where D bq i and λ bq ij are the BQ on-site anisotropy and BQ anisotropic exchange, respectively.
We replace the spin operators S i,j in Eq.10 by bosonic creation (a † i , b † i ) and annihilation (a i , b i ) operators over the honeycomb sub-lattices A and B using Holstein-Primakoff transfor-mations 32 : For i ∈ A sublattice: For i ∈ B sublattice: where we assume that a † i a i << S and b † i b i << S for small deviations of the spins from their ground state orientations. This gives (see Supplementary Section S13 for details): where the sum over i and j runs over the sub-lattices and first nearest neighbors, respectively, and Z is the number of first nearest-neighbors. This procedure outlines a strong implication of the inclusion of higher-order exchange interactions in the description of the magnetic properties of 2D magnets. That is, the enhancement of several magnetic quantities: where λ bq and D bq are the BQ anisotropic exchanges, and the BQ on-site magnetic anisotropy, respectively, for first nearest neighbors. Incidentally Eqs.12−14 yield several implications on the magnetic properties of the sheets, in particular, on the stabilization of magnetism in truly 2D. It is well known that in order to overcome thermal fluctuations that could destroy any magnetic order 33 , sizable magnetic anisotropies need to be developed to gap the lowenergy modes in the magnon spectra. That is, a spin wave gap needs to appear at the energy dispersion to act as a barrier to excitations of long-wavelength spin waves. We can show that the relation between the spin wave gap taken into account BQ exchange interactions (∆ bq ) and that at the level BL exchange (∆ bl ) for first-nearest neighbors is given by (see Supplementary Section S13): where ∆ bl = 2S D + Zλ 2 . By using some parameters for monolayer CrI 3 from Table 2, and approaching the second term in Eq.15 as D bq +3λ bq ≈ 10.72 µeV (see Supplementary Section S7) we can estimate ∆ bl = 0.81 meV and ∆ bq = 1.0 meV. The magnitude of ∆ bq is consistent with recent measurements of the magnon dispersion for bulk CrI 3 , which a spin wave gap of approximately 1.3 meV was measured 34 . Furthermore, the increment of the on-site and anisotropic magnetic anisotropies (Eqs. 13−14) indicates that not only spin-orbit mechanisms are behind the substantial anisotropy in CrI 3 but rather higher-order exchange processes.
Such BQ exchange-driven large magnetic anisotropy mechanism has been proposed for ironbased superconductors 7,17,35 which successfully described their magnetic properties. A direct consequence Eqs.12−14 is the increment of Curie temperatures by factor r given by (see Supplementary Section S13): whereT C and T C are the Curie temperatures with and without BQ interactions, respectively. Including few values in Eq.16, we can roughly estimate an enhancement of r ≈ 39% for monolayer CrI 3 which follows the Monte Carlo calculations (Fig. 3a-b). It is worth noticing that the model in Eq. 11, i) takes into account only first-nearest neighbors in the exchange interactions, and ii) we assume a mean-field approach in the solution of the nonlinear Holstein-Primakoff transformation (e.g. magnon-magnon interactions) to simplify the complex mathematical terms, i.e. four-operator product. Supplementary Section S14 provides a full discussion on the details involved.

Moriya interactions in topological spin excitations
An intriguing question that raised by the presence of BQ exchange interactions is whether they play an important role in the description of magnetic quasiparticles such as magnons and non-trivial spin textures in 2D vdw magnets. It has recently been shown using neutron scattering 34  the delicate balance between them, that need to be considered. In order to account for all these quantities, we extended the model in Eq.1 with the addition of DMI: where A ij is the DMI between spins S i and S j . For a honeycomb ferromagnetic layer, with an easy axis perpendicular to the surface (z-direction), there is no breaking of the inversion symmetry of the lattice at first nearest-neighbors, e.g. A 1st ij = 0. However, contributions from the second nearest-neighbors become non-negligible as space inversion is not present.
Therefore, we consider the DMI vector as A = ν ij A z z, where A z is the magnitude of the DMI along of the easy-axis, and ν ij = ±1 represents the hopping of spins at second nearest-neighbors from sites i to j and vice-versa, respectively (Fig. 2a). Similarly as shown above, we can use Holstein-Primakoff transformations for J i > 0 to write Eq.17 in terms of bosonic creation and annihilation operators (see details in Supplementary Section S15): where we separate the terms due to BL exchange (H BL ), from those due to BQ (H BQ ) and DMI (H DM I ) which can be written as: where H D is the on-site anisotropy term, and H Isotropic i and H Anisotropic i are respectively the isotropic and anisotropic parts of the BL exchange Hamiltonian taken into account up to third nearest-neighbors (i = 1, 2, 3). The sum in ij runs over the first nearest-neighbors at both sublattices A and B while that on ij runs over the second nearest-neighbors specifically on either A or B lattice. We notice that the second nearest neighbors not only break the inversion symmetry of the honeycomb lattice but also generate a magnetic flux φ involving A z and J 2 given by: The magnetic flux (circular arrow) can be appreciated in Fig. 2a  fermions. The topological features of the magnon bands can be described in the k−space by using the Fourier transform of the creation (a † k , b † k ) and annihilation (a k , b k ) operators in Eq.18 as: where the different terms can be written as (see Supplementary Section S15 for details): where C(k) = cos(φ) 3 j=1 cos(k · u j ) and S(k) = sin(φ) 3 j=1 sin(k · u j ). The vectors τ j and u j are respectively between 1 st and 2 nd nearest neighbors (Fig. 2a).
The eigenvalues of Eq.23 can be written in terms of the lower and upper energy bands as: Note that Eqs.24 and Eq.25 are general for any honeycomb material with ferromagnetic order, easy-axis perpendicular to the surface and develop DMI and BQ exchange interactions. For instance, we can use them to predict the energy dispersion of the magnon bands over the first Brillouin zone (BZ) of any 2D magnet, e.g. CrI 3 . Figure 4a-e shows different levels of theory either using a simple XXZ model or including more sophisticated terms through the BQ exchange, DMI or simultaneously all of them. It is clear that XXZ models without any contribution from DMI (Fig. 4a-c) does not describe the gap opening at the Dirac points due to the breaking of the inversion symmetry. Moreover, a model at the level of XXZ+DMI as initially used to understand the magnon dispersion of CrI 3 34 does not capture entirely the full profile of the bands (Fig. 4d). The upper branch E + becomes nearly flat with the increment of J 2 at the path K − M − K while the lower magnon branch E − turns more curved. This picture modifies substantially when BQ exchange interactions are included in the XXZ+BQ+DMI model (Fig. 4e) as the magnon bands follow a similar curvature throughout the variation of J 2 , although the upper branch at Γ increased to E + Γ =27.8 meV.
A sound comparison between the XXZ+BQ+DMI model and the experimental results for CrI 3 is obtained (Fig. 4f) when J 1 is varied similarly (Fig. 4a) indicating that the first nearestneighbors are important on the stabilization of E + Γ . It is worth mentioning that the value of the exchange interactions taken into account in the fitting of the neutron scattering spectra 34 do not separate BL contributions from BQ as shown in Eq.12-14. Hence, it is not known from the fitting procedure 34 what is the contribution of K bq to the magnon dispersion. However, such separation can be clearly stated in our model as indicated in Fig. 4f. Furthermore, even though DMI is important for the gap opening at the Dirac point, it does not contribute to the magnitudes of the magnetization or critical temperatures for any 2D magnet with an out-of-plane easy-axis and ferromagnetic aligned spins (see details in Supplementary Section S15).

Implications and prospects
The effects of the BQ exchange interactions proposed here should manifest in experimentally accessible temperature range. One indication is already the accurate reproduction of experimental Curie temperatures 1,29 including higher-order exchange which could not be obtained at the level of Ising, Heisenberg or Kitaev models. The magnitudes of BQ exchange can be in principle extracted from accurate hysteresis loops using phenomenological models 41,42 .
Importantly in such analyses are potential temperature variations of BL and BQ exchanges with layer thickness which may indicate tunable interlayer exchanges still to be explored in 2D magnets.

Supplementary sections
Materials and Methods.

Data Availability
The data that support the findings of this study are available within the paper and its Supplementary Information.

Competing interests
The Authors declare no conflict of interests.   Table 2: Computed values of several magnetic quantities for CrX 3 (X=I, Br, Cl, F) at different number of nearest neighbors: isotropic (J 1 , J 2 , J 3 ) and anisotropic (λ 1 , λ 2 , λ 3 ) BL exchanges. The on-site magnetic anisotropy D is also included. See Supplementary Section S10 for details.  , chromium based ternary tellurides (Cr 2 X 2 Te 6 , X=Ge, P, Si), manganese based ternary chalcogenides (Mn 2 P 2 X 6 , X=S, Se, Te), transition metal dichalcogenides (MnX 2 , X=S, Se, Te) of different phases (2H, 1T ) and an iron-based dichalcogenide (2H-FeS 2 ). The reference energy is taken at 0 o as the spins oriented at the same direction. Rotations can occur inplane or out-of-plane with similar behavior. Symbols are calculated energies. Dashed lines correspond to a quadratic fitting using E tot bq (θ) = A bq 0 + A bq 1 · S 2 cos(θ) + A bq 2 · S 4 cos 2 (θ), while color filling areas indicate the deviation between the quadratic E tot bq and a linear fitting using E tot bl (θ) = A bl 0 + A bl 1 · S 2 cos(θ). Materials that show large deviation, such as CuBr 3 or 2H-FeS 2 , develop large BQ exchange interactions. f, Logarithm of the total energies for the dataset in b-e, as a function of θ (dots). A quadratic regression (Q-Regression) is evaluated over the calculated DFT energies (solid line) indicating a universal behavior of BQ exchange interactions in 2D magnets. , which breaks the inversion symmetry of the lattice. The dashed lines show the magnon hopping between second NN as the magnons gain a phase given by φ (see text). The lattice vectors u i (i = 1, 2) and τ j (i = 1, 2, 3) show the first and second NN on the lattice, respectively. b, Zoom-in on the BQ exchange process involving two electrons between sites M A and M B with 3d m electrons in the valence. The BQ exchange K bq is mediated by non-magnetic atoms X with a valence given by np electrons, where n will depend on the atomic elements involved. An on-site Hubbard U term is at the M A,B sites with U ex representing the potential spin-splitting when the spins of the electrons involved in the BQ exchange align ferro-or anti-ferromagnetically (Supplementary Section S9). The difference between up and down spin coupling is given by 2U ex . c, K bq /J bl versus ∆ 0 (eV) for all materials displaying BQ exchange interactions. Magnetic atoms with similar chemical environment in terms of Coulomb repulsion, exchange interactions and valence follow alike behavior for K bq /J bl . For instance, Cr in blue and Mn in green (inset). Orange dots show materials with dissimilar electronic configurations but with ∆ 0 = 0.