First principles Heisenberg models of 2D magnetic materials: the importance of quantum corrections to the exchange coupling

We will consider the Heisenberg model with nearest neighbor interactions given by

$$\begin{equation}H=-\frac{J}{2}\sum _{\langle ij\rangle }{\mathbf{S}}_{i}\cdot {\mathbf{S}}_{j},\end{equation} \tag{ 1 }$$

where the sum is over nearest neighbor (magnetic) atoms and J is the exchange constant. For a bipartite lattice there is a unique [up to an SO(3) rotation in spin space] anti-ferromagnetic state where all sites are anti-aligned with neighboring sites. If J < 0 this comprises the classical ground state (the Neel state) and for simplicity we will assume the spins to be aligned along the z-direction in this state. However, this is not an eigenstate of the quantum model and standard spinwave analysis shows that there is a state with lower energy referred to as the non-interacting magnon (NIM) state [19]. We will regard this as an approximate eigenstate of the Heisenberg Hamiltonian in the following. The NIM state has the same spin symmetry as the Neel state, but it is not an eigenstate of ${S}_{i}^{z}$. Instead the expectation value is given by

$$\begin{equation}{\langle {S}_{z}\rangle }_{\mathrm{N}\mathrm{I}\mathrm{M}}=S\left(1-\alpha /S\right),\end{equation} \tag{ 2 }$$

where the constant S is the largest eigenvalue of S_z and α is given by

$$\begin{equation}\alpha =\frac{1}{2}\left[{\left\langle \frac{1}{\sqrt{1-\vert {\gamma }_{\mathbf{q}}{\vert }^{2}}}\right\rangle }_{\mathrm{B}\mathrm{Z}}-1\right]\end{equation} \tag{ 3 }$$

where

$$\begin{equation}{\gamma }_{\mathbf{q}}=\frac{1}{{N}_{\mathrm{n}\mathrm{n}}}\sum _{{\Delta}}{\mathrm{e}}^{-\mathrm{i}\mathbf{q}\cdot {\mathbf{R}}_{{\Delta}}}.\end{equation} \tag{ 4 }$$

Here ⟨...⟩_BZ denotes a Brillouin zone average over q and R_Δ are the N_nn lattice vectors connecting nearest neighbor sites. The constant α is larger for low-dimensional models and thus becomes more important for 2D materials than for 3D materials. In the case of the honeycomb and square lattices the value of α can be evaluated numerically yielding 0.258 and 0.197 respectively.

Although the NIM state is obviously useful to describe properties of anti-ferromagnets (J < 0), the derivation does not require that J is negative. In particular, for a ferromagnetic Heisenberg model (J > 0) the NIM state can be regarded as the approximate eigenstate of highest energy. With an accurate exchange–correlation functional, it is expected that such a state should be represented by a configuration where the spins share the symmetry of the NIM state. This coincides with the spin symmetry of the Neel state and the spin configuration thus corresponds qualitatively to the classical anti-ferromagnetic configuration of interest. Since the associated spin densities are in principle accurately described by DFT, the ratio m_AFM/m_FM should yield (1 − α/S) provided that the magnetic moments are localized. It is, however, far from obvious that standard approximations for the exchange–correlation functional will capture the intricate correlations in the anti-ferromagnetic state. In section 3 we will provide evidence that a proper renormalization of the spin is captured in simple generalized gradient approximations exemplified by the Perdew–Burke–Ernzerhof (PBE) [26] functional.

For an N-site periodic bipartite lattice, the expectation value of the Hamiltonian (1) using the Neel state is NJS²N_nn/2 where N_nn is the number of nearest neighbors. For an anti-ferromagnetic lattice (J < 0) this comprises the classical ground state and for a ferromagnetic lattice (J > 0) it is the classical state of highest energy. However, the NIM state has a lower (higher) energy for anti-ferromagnetic (ferromagnetic) models of the form (1). It is given by

$$\begin{equation}{E}^{\mathrm{N}\mathrm{I}\mathrm{M}}=\frac{N}{2}\left({N}_{\mathrm{n}\mathrm{n}}{S}^{2}J\right)\left[1+\beta /S\right],\end{equation} \tag{ 5 }$$

where

$$\begin{equation}\beta =1-{\left\langle \sqrt{1-\vert {\gamma }_{\mathbf{q}}{\vert }^{2}}\right\rangle }_{\mathrm{B}\mathrm{Z}}.\end{equation} \tag{ 6 }$$

For the honeycomb and square lattices the value of β are given by 0.202 and 0.158 respectively.

Exchange coupling constants are routinely evaluated from DFT using ferromagnetic and anti-ferromagnetic spin configurations in the simulations [20–25]. But it is usually assumed that such configurations can be mapped to the ferromagnetic state (E^FM = −NJS²N_nn/2) as well as the Neel state (E^Neel = NJS²N_nn/2), which always leads to an overestimation of J. In exact DFT, the total energy of a given spin configuration should be mapped to the eigenstate of the Heisenberg model of the same spin symmetry. In particular, anti-ferromagnetic configurations have to be mapped to the NIM state, which provides a much better description of the anti-ferromagnetic state than the Neel state. For bipartite lattices this yields the expression

$$\begin{equation}J=\frac{{\Delta}E}{{N}_{\mathrm{n}\mathrm{n}}{S}^{2}\left(1+\beta /2S\right)},\end{equation} \tag{ 7 }$$

where ΔE = E^NIM − E^FM is the energy difference per magnetic atom obtained from DFT. Again, we remark that the value of β is in general smaller in 3D systems due to the 3D BZ average and it is often a better approximation to neglect the correction when evaluating exchange constants in 3D. However, as we will see below the inclusion of correlation effects in the energy mapping analysis can lead to significant corrections to the predictions of exchange constants and critical temperatures in 2D. Most spectacular are the effect on spin-1/2 systems on a honeycomb lattice where the exchange constants will be reduced by 17% compared to the classical prediction.

The model (1) does not allow for magnetic order in 2D due to the Mermin–Wagner theorem [9] and one needs to consider models with terms that explicitly break the spin-rotational symmetry. Such terms originate from spin–orbit coupling and here we will assume that the most important effect on the magnetic order comes from single-ion anisotropy and nearest neighbor anisotropic exchange. The Hamiltonian then takes the form

$$\begin{equation}H=-\frac{J}{2}\sum _{\langle ij\rangle }{\mathbf{S}}_{i}\cdot {\mathbf{S}}_{j}-\frac{\lambda }{2}\sum _{\langle ij\rangle }{S}_{i}^{z}{S}_{j}^{z}-A\sum _{i}{\left({S}_{i}^{z}\right)}^{2},\end{equation} \tag{ 8 }$$

where we have assumed isotropy in the xy-plane, which we take to comprise the atomic plane of the material. From hereon we restrict ourselves to the case of J > 0. If one assumes that the easy-axis is along the z-direction, a simple spin-wave analysis then shows that the magnetic excitation spectrum has a gap given by

$$\begin{equation}{\Delta}=A\left(2S-1\right)+\lambda S{N}_{\mathrm{n}\mathrm{n}}.\end{equation} \tag{ 9 }$$

A finite gap in the spectrum implies a broken spin-rotational symmetry and the model is expected to exhibit magnetic order at finite temperatures. If, however, the z-axis is not the easy axis one will obtain a negative spinwave gap (Δ < 0) and the spin-wave analysis is faulty because it has not been carried out on the magnetic ground state. If one assumes in-plane magnetic isotropy a negative spin-wave gaps thus implies the presence of a residual continuous rotational symmetry and the material cannot exhibit magnetic order by the Mermin–Wagner theorem. The sign of the spin-wave gap [evaluated from equation (9)] thus represents a descriptor for magnetic order at finite temperatures—even if it only corresponds to a physical quantity in the case of Δ > 0.

The anisotropy constants λ and A can be evaluated from DFT including spin–orbit coupling by considering energy differences between in-plane and out-of-plane spin configurations [16]. However, the anisotropy terms in equation (8) will alter the expression for the exchange coupling and if spin–orbit coupling is strong the spin–orbit corrections to J will typically be larger than the quantum corrections. For S ≠ 1/2. We thus propose to evaluate the Heisenberg parameters according to

$$\begin{equation}A=\frac{{\Delta}{E}_{\text{FM}}\left(1-\frac{{N}_{\text{FM}}}{{N}_{\text{AFM}}}\right)+{\Delta}{E}_{\text{AFM}}\left(1+\frac{{N}_{\text{FM}}}{{N}_{\text{AFM}}}\right)}{\left(2S-1\right)S},\end{equation} \tag{ 10 }$$

$$\begin{equation}\lambda =\frac{{\Delta}{E}_{\text{FM}}-{\Delta}{E}_{\text{AFM}}}{{N}_{\text{AFM}}{S}^{2}},\end{equation} \tag{ 11 }$$

$$\begin{equation}J=\frac{{E}_{A\mathrm{F}\mathrm{M}}^{\parallel }-{E}_{\text{FM}}^{\parallel }}{{N}_{\text{AFM}}{S}^{2}\left(1+\beta /2S\right)},\end{equation} \tag{ 12 }$$

where ${\Delta}{E}_{\text{FM}\left(\mathrm{A}\mathrm{F}\mathrm{M}\right)}={E}_{\text{FM}\left(\mathrm{A}\mathrm{F}\mathrm{M}\right)}^{\parallel }-{E}_{\text{FM}\left(\mathrm{A}\mathrm{F}\mathrm{M}\right)}^{\perp }$ are the DFT energy differences between in-plane and out-of-plane magnetization for ferromagnetic (anti-ferromagnetic) spin configurations and N_FM(AFM) is the number of nearest neighbors with aligned (anti-aligned) spins in the anti-ferromagnetic configuration. Except for the quantum correction in equation (12) and the factor of (2S − 1) in equation (10) this comes out of an exact energy mapping to the classical Heisenberg model and the quantum corrections are simply added ad hoc to the expression for the exchange coupling. In addition we take 2S → (2S − 1) in the denominator of equation (10) to ensure that the sign of the spinwave gap is strcitly determined by the sign of ΔE_FM. In the case of S = 1/2 we take A = 0. It should be stressed that in contrast to equation (7) this is no longer rigorously derived since only two of the four spin configurations needed for the evaluation of equations (10)–(12) are eigenstates of the quantum mechanical Heisenberg model. For example, with J > 0 and out-of-plane easy-axis, the ground state of the Heisenberg model will be mapped ${E}_{\text{FM}}^{\perp }$ and the eigenstate of highest energy will be mapped to ${E}_{\text{AFM}}^{\parallel }$, but the remaining two spin configurations (which have energies obtainable from DFT) are not expected to be represented by eigenstates of the Heisenberg model.

The Curie temperature of the model (8) may be obtained by either classical Monte Carlo simulations [18] or a renormalized spin-wave analysis [12, 17, 18]. We have previously shown that renormalized spin-wave theory breaks down in the case of large single-ion anisotropy [18] and in the present work we will evaluate Curie temperatures from classical Monte Carlo simulations. It may seem odd to rely on a classical analysis since we have argued that it is crucial to include quantum corrections when mapping DFT calculations to Heisenberg models. However, at temperatures in the vicinity of the critical temperature quantum fluctuations tend to be quenched by thermal fluctuations and a classical analysis becomes reliable even if they cannot be trusted at low temperatures. We note that spin-1/2 systems may comprise an important exception to this, since the single-ion anisotropy term becomes proportional to the identity in that case. It follows that magnetic order cannot exist at finite temperatures in spin-1/2 systems unless λ ≠ 0, which is in stark contrast to the predictions of classical Monte Carlo simulations where the value of S simply introduces a rescaling of the Heisenberg parameters. For S ≠ 1/2, Monte Carlo simulations of the model (8) can be accurately fitted to a function of the form [18]

$$\begin{equation}{T}_{\text{C}}={T}_{\text{C}}^{\text{Ising}}f\left(\frac{{\Delta}}{J\left(2S-1\right)}\right)\end{equation} \tag{ 13 }$$

where

$$\begin{equation}f\left(x\right)={\mathrm{tanh}}^{1/4}\left[\frac{6}{{N}_{\mathrm{n}\mathrm{n}}}\;\mathrm{log}\left(1+\gamma x\right)\right]\end{equation} \tag{ 14 }$$

and γ = 0.033. ${T}_{\text{C}}^{\text{Ising}}$ is the critical temperature of the corresponding Ising model, which (in units of JS²/k_B) is given by 1.52, 2.27, and 3.64 for honeycomb, square and hexagonal lattices respectively.

In order to assess the performance of semi-local functionals for the correlated anti-ferromagnetic configuration of real materials, we have calculated the integral of the norm of the magnetization density

$$\begin{equation}{m}_{\mathrm{A}\mathrm{b}\mathrm{s}}=\int \vert {m}_{{\uparrow}}\left(\mathbf{r}\right)-{m}_{{\downarrow}}\left(\mathbf{r}\right)\vert \mathrm{d}\mathbf{r}\end{equation} \tag{ 15 }$$

for a wide variety of 2D materials in both the ferromagnetic and anti-ferromagnetic state. We have included all insulating magnetic materials present in the Computational 2D Materials Database [27] (C2DB) where the magnetic atoms form either a honeycomb, square or triangular lattice. In this context, we define an insulator by a threshold of 0.2 eV for the band gap in both FM and AFM configurations to ensure that the basic electronic structure is not altered too much between the different spin configurations. The computational details can be found in reference [27], but here we just mention that all calculations were carried out with the electronic structure code GPAW [28–30] using the PBE approximation [26] for the exchange–correlation energy. Since the triangular lattice is not bipartite the classical state of lowest/highest energy is not collinear, but comprises a 120° non-collinear structure [31]. However, as the far majority of DFT calculations in the literature (including energy mapping studies) are based on collinear structures we choose to base the energy mapping analysis of triangular lattices on a 'striped' anti-ferromagnetic configuration where each site has two aligned and four anti-aligned nearest neighbors [16]. Since this state does not comprise an extremum in the classical energy landscape it is not possible to perform a spinwave analysis to obtain the quantum corrections to the spin, but due to the four anti-aligned nearest neighbors the correction is expected to be similar to the square lattice.

In figure 1 we display the ratio m_AFM/m_FM [calculated from equation (15)] per magnetic atom along with the theoretical predictions given by equation (2). Although the calculations exhibit clear deviations from the Heisenberg prediction, there is a significant trend toward increasing quantum corrections with decreasing spin. It may be argued that a reduction of magnetic moments is expected in the anti-ferromagnetic state due to incomplete localization of the magnetic moments, which introduces cancellation of magnetization densities in the interstitial regions. However, such cancellation effects are expected to yield a ratio m_AFM/m_FM, which is independent of spin when averaged over a large class of materials. In contrast, figure 1 shows a clear tendency toward decreasing staggered magnetization with decreasing spin, which is in accordance with the correlated state predicted by the Heisenberg model.

Zoom In Zoom Out Reset image size

The reduction of the ratio m_AFM/m_FM with decreasing spin provides some confidence that DFT is able to capture (at least partly) the intricate correlations of the anti-ferromagnetic state at the level of generalized gradient calculations. This implies that exchange parameters calculated from energy mapping to the Heisenberg model should be corrected according to equation (12). This in turn may have a strong influence on the calculation of critical temperatures from the expression (13). In figure 2 we show the quantum corrections to the exchange coupling parameters as well as the corrected critical temperatures relative to the classical estimates for all ferromagnetic insulators in the C2DB with Δ > 0. In the case of exchange parameters the reduction only depends on the value of S, whereas the corrected critical temperatures also depend on the spinwave gap Δ. However, if Δ/J ≪ 1 the change in critical temperature becomes δT_C/T_C = 3δJ/4J and the change in critical temperature thus largely follows the exchange parameter, which is also evident from figure 2. In general the reduction in critical temperatures is on the order of 5–7%. We note here that the non-corrected values of J and T_C calculated here differ slightly from previously published results calculated with the same method [16] since the effect of anisotropy is included in the evaluation of exchange constants, which were not the case in reference [16]. All the results are compiled in table 1.

Zoom In Zoom Out Reset image size

Although we have argued that PBE (partly) captures quantum nature of anti-ferromagnetic configurations, a correct prediction of ground state energies for magnetic insulators pose a major challenge for DFT, since these materials typically involve strongly correlated electrons. This implies a large uncertainty in the prediction of exchange constants and critical temperatures. To quantify this we have repeated the calculations of the materials in the BiI₃ and CdI₂ prototypes using PBE + U where we have adopted the values of the Hubbard corrections (U) used in reference [32]. Table 2 shows a detailed list of materials obtained with the Hubbard correction and the calculated parameters. The overall trend in the corrections is very similar to figure 2 and are not shown. However the list of ferromagnetic candidates in table 2 is rather different from table 1. Several of the materials in table 2 are lacking from table 1 and several materials present in table 2 are absent from table 1. Part of the reason is that we only include ferromagnetic insulators and the inclusion of the Hubbard correction may open a gap in a material that is predicted to be metallic without the Hubbard correction. In addition, the Hubbard correction may change the sign of both the spinwave gaps and exchange constants when U is included in the calculations. For example, with PBE + U the out-of-plane direction becomes a hard axis in NiI₂ and the material is predicted to have T_C = 0, but NiCl₂ acquires an out-of-plane easy with PBE + U and is thus predicted to have a finite critical temperature. In addition the Mn halides becomes ferromagnetic with PBE + U whereas they are anti-ferromagnetic with bare PBE. In table 2 we have highlighted all the materials present in both cases with bold face. We note that the value of T_C = 47 K for CrI₃ obtained with PBE + U is in better agreement with the experimental value of 45 K compared to the bare PBE result of 29 K. This could, however, be coincidental since several other approximations enter the evaluation of T_C in the present framework. In addition, it has previously been shown that the value depends rather strongly on the choice of functional [33]. In particular, the results of an LDA + U calculation gives a reduction of J by 20% and there is no a priori reason to believe that PBE should be better than LDA for the calculation of exchange constants. Nevertheless, the inclusion of Hubbard corrections typically provides a better description of the electronic structure in correlated materials, but it is not clear which value of U one should use.