The anisotropic Heisenberg model

The Heisenberg model can be derived as an approximation to the many-body electronic Hamiltonian in a localized basis where only the spin degrees of freedom are retained. In that framework, the interactions between localized spins are referred to as direct exchange and arise as a consequence of Coulomb interactions and Pauli exclusion.^{[Reference Yosida8]} A rather different type of magnetic interactions may be derived by assuming strongly localized spins and include hybridization by second order perturbation theory. This will typically be mediated by non-magnetic ligand atoms and is referred to as superexchange.^{[Reference Anderson21,Reference Anderson22]} Remarkably, the superexchange interaction between localized spins has precisely the same form as the direct exchange interaction and can be included in the Heisenberg model by a simple redefinition of the interaction parameters.

A rather general form of the Heisenberg model, which includes magnetic anisotropy, can be written as

(1)$$H = -\displaystyle{1 \over 2}\sum\limits_{ij\alpha \beta} {S_i^\alpha {\rm {\frak J}}_{ij}^{\alpha \beta} S_i^\beta} -\sum\limits_{i\alpha} {A_i^\alpha {\lpar {S_i^\alpha} \rpar }^2}, $$

where $S_i^\alpha$ is the α component of the spin operator for the magnetic atom at site i and ${\rm {\frak J}}_{ij}^{\alpha \beta}$ is a 3 × 3 interaction matrix that couple spin states at sites i and j. The second sum contains single-ion anisotropy constants that determine the energy cost of a global spin rotation. In general, the off-diagonal elements of ${\rm {\frak J}}_{ij}^{\alpha \beta}$ lead to Kitaev interactions^{[Reference Xu, Feng, Xiang and Bellaiche23,Reference Banerjee, Bridges, Yan, Aczel, Li, Stone, Granroth, Lumsden, Yiu, Knolle, Bhattacharjee, Kovrizhin, Moessner, Tennant, Mandrus and Nagler24]} and Dzyaloshinskii–Moriya interactions^{[Reference Heide, Bihlmayer and Blügel25–Reference Liu, Shi, Lu and Anantram27]} that may give rise to highly intriguing physical properties. However, the off-diagonal components are typically an order of magnitude smaller than the isotropic part^{[Reference Xu, Feng, Xiang and Bellaiche23]} and have marginal influence on the ordering temperatures. In the following, we will thus stick to a simplified version where the off-diagonal components of ${\rm {\frak J}}_{ij}^{\alpha \beta}$ is neglected and we assume in-plane isotropy such that ${\rm {\frak J}}_{ij}^{xx} = {\rm {\frak J}}_{ij}^{yy}$ and $A_i^x = A_i^y$. The Heisenberg model can then be written as

(2)$$H = -\displaystyle{{\widetilde{{J_{ij}}}} \over 2}\sum\limits_{ij} {{\bi S}_i\cdot {\bi S}_j-\displaystyle{{\widetilde{{J_{ij}^z}}} \over 2}} \sum\limits_{ij} {\lpar {2S_i^z S_j^z -S_i^x S_j^x -S_i^y S_j^y} \rpar } -\sum\limits_i {A_i{\lpar {S_i^z} \rpar }^2}, $$

where $A_i = A_i^z -A_i^x$, $\widetilde{{J_{ij}^z}} = \lpar {{\rm {\frak J}}_{ij}^{zz} -{\rm {\frak J}}_{ij}^{xx}} \rpar /3$, and $\widetilde{{J_{ij}}} = {\rm Tr}\lsqb {{\rm {\frak J}}_{ij}^{\alpha \beta}} \rsqb /3$ is the isotropic part of the exchange tensor. Equation (2) is then equivalent to equation (1) except for a constant term given by $\sum\nolimits_i {A_i^x S_i\; \lpar {S_i + 1} \rpar }$, where S _i is the maximum eigenvalue of $S_i^z$. In equation (2), we have explicitly split the exchange tensor into an isotropic and a traceless part. However, it is often more convenient to rewrite it in the form

(3)$$H = -\displaystyle{1 \over 2}\sum\limits_{ij} {J_{ij}{\bi S}_i\cdot {\bi S}_j} -\displaystyle{1 \over 2}\sum\limits_{ij} {\lambda _{ij}S_i^z S_j^z} -\sum\limits_i {A_i{\lpar {S_i^z} \rpar }^2}, $$

with $J_{ij} = {\rm \;} \widetilde{{J_{ij}}}-\widetilde{{J_{ij}^z}}$, and $\lambda = 3\widetilde{{J_{ij}^z}}$. This is the version of the Heisenberg Hamiltonian that will be used in the following. We note that different conventions for the parameters are used in the literature and it is always crucial to specify a convention when referring to a set of parameters.

Heisenberg parameters from first principles

Before delving into an analysis of equation (3), we will briefly show how to obtain the Heisenberg parameters J _ij, λ _ij, and A _i from first principles simulations. It is, in principle, possible to calculate the exchange coupling constants from microscopic expressions involving either an exchange integral in the case of direct exchange or from a combination of hopping matrix elements and on-site Coulomb interactions in the case of superexchange.^{[Reference Besbes, Nikolaev, Meskini and Solovyev28,Reference Wang, Nikolaev, Ren and Solovyev29]} This approach is, however, problematic for several reasons: (1) a quantitative calculation has to be carried out in a localized basis set (for example, Wannier functions^{[Reference Marzari, Mostofi, Yates, Souza and Vanderbilt30]}), and the results will inevitably depend on a specific choice of basis. (2) The exchange mechanism is often a combination of exchange and superexchange, and there is no way of a priori determining how to balance the contributions in a quantitative treatment of the parameters. (3) Methods based on DFT yield accurate total energies (if a good xc-functional is used), but “Kohn–Sham parameters” such as band energies, hopping parameters, or exchange integrals do not comprise genuine physical properties.

A much better approach is based on an energy mapping analysis, where DFT total energies corresponding to different spin symmetries^{[Reference Görling31,Reference Görling32]} are mapped to the model (3) and used to extract the parameters.^{[Reference Xiang, Lee, Koo, Gong and Whangbo33–Reference Bose and Kudrnovský38]} This requires no knowledge of the underlying exchange mechanism, and the accuracy is only limited by the quality of the applied exchange-correlation functional. In particular, if we restrict ourselves to a single type of magnetic atom and nearest neighbor interactions only, a classical treatment of the Heisenberg model leads to^{[Reference Torelli, Thygesen and Olsen15]}

(4)$$J = \displaystyle{{E_{{\rm AFM}}^\parallel -E_{{\rm FM}}^\parallel} \over {N_{{\rm AFM}}S^2}},$$

(5)$$\lambda = \displaystyle{{E_{{\rm FM}}^\parallel -E_{{\rm FM}}^\bot -E_{{\rm AFM}}^\parallel + E_{{\rm AFM}}^\bot} \over {N_{{\rm AFM}}S^2}},$$

(6)$$A = \displaystyle{{\lpar {E_{{\rm FM}}^\parallel -E_{{\rm FM}}^\bot} \rpar \delta _--\lpar {E_{{\rm AFM}}^\parallel -E_{{\rm AFM}}^\bot} \rpar \delta _ +} \over {2S^2}},$$

where E _FM and E _AFM are the energies of the ferromagnetic and an antiferromagnetic configuration, respectively. The superscripts denotes whether the magnetization is in-plane $(\parallel )$ or out-of-plane (⊥) and δ _± = 1 ± N _FM/N _AFM, where N _FM and N _AFM are the numbers of nearest neighbors with aligned and anti-aligned spins, respectively, in the anti-ferromagnetic configuration. In the case of a bipartite lattice of magnetic atoms, there is a fully anti-ferromagnetic configuration with N _FM = 0 and δ _± = 1. However, for non-bipartite lattices (a triangular lattice for example), there is no natural choice of the antiferromagnetic state and one has to choose a frustrated configuration with N _FM ≠ 0. Typically, the Heisenberg parameters are on the order of meV or less, and it is crucial to use the same basis sets for all configurations in order to obtain well converged results. This approach is easily generalized to yield additional Heisenberg parameters such as next-nearest and third-nearest neighbor interactions, but for insulators these are often much smaller than the nearest neighbor interactions, which usually govern the critical temperatures. In general, however, one has to check whether it is sufficient to restrict the analysis to nearest neighbor interactions.

If the ground state is antiferromagnetic, it may not be obvious, which value to use for S in the energy mapping analysis. However, a ferromagnetic configuration without spin–orbit coupling is bound to yield a half integer for any insulator and that defines the value of S to be used in the Heisenberg model. We note that DFT codes often provide an estimate of local magnetic moments based on an integrated spin-density in some vicinity of the atomic nuclei, but this comprises a rather arbitrary measure and should not be used for quantitative analysis. Moreover, spin–orbit effects may yield a non-integer value of the magnetic moments (in unit of Bohr magnetons per cell) in a ferromagnetic configuration, which is at odds with the Heisenberg model that assumes a half-integer spin in the ferromagnetic ground state. This is due to interband spin mixing, which is not captured by the Heisenberg model and such effects cannot be captured by a model like equation (3). While it would be highly interesting to extend the Heisenberg model to include multiple band effects, this is a rather small effect that is unlikely to affect the prediction of critical temperatures. For example, in CrI₃, the inclusion of spin–orbit coupling leads to an LDA magnetic moment per Cr ion of 3.01 Bohr magnetons, which is very close to the nominal value of 3.0. Finally, the inclusion of spin–orbit coupling that is needed to obtain the four DFT energies $E_{{\rm FM/AFM}}^{\bot /\parallel}$ introduces certain subtleties in the calculations. In the case of CrI₃, the ground state is ferromagnetic with a strong out-of-plane easy axis and from this state $E_{{\rm FM}}^\bot$ is obtained. However, the ferromagnetic state with in-plane magnetization comprises a saddle point in spin-configuration space and any unconstrained calculation will converge toward the state with out-of-plane magnetization. To resolve this, one can either perform a constrained DFT calculation or include spin–orbit corrections non-selfconsistently.

It should be emphasized that the energy mapping analysis leading to equations (4)–(6) was based on the classical Heisenberg model. For a ferromagnetic system with easy-axis anisotropy, the ground state energy coincides with the classical energy. But for an anti-ferromagnetic system, the classical ground state energy only provides an upper estimate of the true ground state energy. For example, for a square lattice with nearest neighbor interactions and no anisotropy, the ground state energy per site is approximately 2JS ²(1 + 0.158/S),^{[Reference Yosida8]} whereas the classical energy is simply 2JS ² (J is negative). If one assumes that DFT provides the correct ground state energy, the classical energy mapping analysis will lead to an overestimation of the exchange coupling constant. This problem also pertains to ferromagnetic systems since the antiferromagnetic state that enters the analysis is presumably correctly described in DFT and therefore does not correspond to the classical energy assumed in equations (4)–(6). In principle, it is straightforward to include quantum corrections in the analysis for any isotropic Heisenberg model, but the situation becomes somewhat more complicated when anisotropy is introduced. For example, a ferromagnetic Heisenberg model with out-of-plane easy axis has a lowest and a highest eigenstate, respectively, that must correspond to $E_{{\rm FM}}^\bot$ and $E_{{\rm AFM}}^\parallel$ obtained from DFT. But the energies $E_{{\rm FM}}^\parallel$ and $E_{{\rm AFM}}^\bot$ cannot be mapped onto any eigenstate of the Heisenberg model and it is not obvious how to proceed with the quantum mechanical energy mapping in this case.

The accuracy of commonly applied exchange-correlation functionals for exchange and anisotropy constant has not yet been completely clarified, but there has been a few studies on the performance of hybrids^{[Reference Ködderitzsch, Hergert, Temmerman, Szotek, Ernst and Winter35]} and LDA + U/PBE + U for three-dimensional bulk systems.^{[Reference Olsen37]} It was found that for the anti-ferromagnets NiO and MnO, LDA + U and PBE + U were able to provide good agreement with experimental values for nearest and next nearest neighbor interactions if the “right” value of U was chosen. In Table I, we compare the Heisenberg parameters for CrI₃ calculated with LDA, LDA + U, PBE, PBE + U, and PBEsol, and we have used a value of U = 3.5 eV for LDA + U and PBE + U. All calculations were performed with the electronic structure code GPAW, which is based on plane waves and the projector-augmented wave method.^{[Reference Enkovaara, Rostgaard, Mortensen, Chen, Dułak, Ferrighi, Gavnholt, Glinsvad, Haikola, Hansen, Kristoffersen, Kuisma, Larsen, Lehtovaara, Ljungberg, Lopez-Acevedo, Moses, Ojanen, Olsen, Petzold, Romero, Stausholm-Møller, Strange, Tritsaris, Vanin, Walter, Hammer, Häkkinen, Madsen, Nieminen, Nørskov, Puska, Rantala, Schiøtz, Thygesen and Jacobsen39–Reference Larsen, Mortensen, Blomqvist, Castelli, Christensen, Dułak, Friis, Groves, Hammer, Hargus, Hermes, Jennings, Jensen, Kermode, Kitchin, Kolsbjerg, Kubal, Kaasbjerg, Lysgaard, Maronsson, Maxson, Olsen, Pastewka, Peterson, Rostgaard, Schiøtz, Schütt, Strange, Thygesen, Vegge, Vilhelmsen, Walter, Zeng and Jacobsen41]} All structures were relaxed with the given functional. It is observed that the predicted exchange constants may differ by more than a factor of three depending on the applied functional. In addition, the anisotropy constant is strongly dependent on the whether or not an on-site Coulomb repulsion U is introduced. However, the spin-wave gap Δ = A(2S − 1) + λN _nnS (to be introduced below) largely determines the spin–orbit mediated magnetic properties and is not very sensitive to the functional. This signifies that the spin–orbit effects are rather insensitive to the xc-functional. Including an onsite U term in the calculation merely transfers the single-ion anisotropy energy into the anisotropic exchange energy. However, the large span in exchange constants indicates that first principles predictions of magnetic properties such as critical temperatures may be off by a factor of three, and there is a strong need for a systematic study of the accuracy of exchange-correlation functionals for exchange constants.

Table I.

Exchange and anisotropy constants of a monolayer CrI₃ calculated with a few different functionals.

All values are in meV. Spin–orbit coupling was included non-selfconsistently except for the values in brackets, which were obtained with selfconsistent spin–orbit coupling. LDA + U and PBE + U were performed with U = 3.5 eV.

Easy-axis order

Once a first-principles Heisenberg Hamiltonian has been obtained through the energy mapping analysis, we are left with the pertinent problem of solving it. In three-dimensional magnets, a rough estimate of the critical temperature can be obtained from Weiss mean field theory. In 2D, however, this approach is bound to fail since there is no reference to the dimensionality of the problem and finite critical temperatures will be predicted even in the absence of anisotropy. In 2D, the effect of anisotropy must be included on equal footing with the exchange interactions. In the presence of an easy-axis (here chosen as the z-direction), a rigorous quantum mechanical analysis can be carried out by introducing the Holstein–Primakoff transformation, which replaces the spin operators in equation (3) by bosonic operators that create or annihilate spin-wave excitations (magnons). For a single type of magnetic atom, the transformation is

(7)$$\eqalign{& S_i^- = \sqrt {2S} a_i^\dagger \sqrt {1-\displaystyle{{a_i^\dagger a_i} \over {2S}}}, \cr & \quad S_i^ + = \sqrt {2S} \sqrt {1-\displaystyle{{a_i^\dagger a_i} \over {2S}}} a_i,\quad S_i^z = S-a_i^\dagger a_i,}$$

with $S_i^ \pm = S_i^x \pm iS_i^y$. The bosonic commutation relations for A _i and $a_i^\dagger$ imply that the commutation relations for the spin operators are fulfilled. In order to proceed, the square roots must be Taylor expanded and one obtains

(8)$$H = E_0 + H_2 + H_4 + H_6 + \ldots, $$

where E ₀ is zeroth order in raising and annihilation operators, H ₂ is second order, H ₄ is fourth order and so forth. In equation (8), we assume each term to be normal ordered such that all annihilation operators are to the right. This implies that E ₀ is simply the ground state energy, H ₂ determines the spectrum of single magnon excitations, H ₄ gives rise to two-magnon interactions, and H ₆ gives rise to three-magnon interactions. It should be noted that the anisotropy parameters only enter E ₀, H ₂, and H ₄, and a truncation at fourth order thus includes anisotropy exactly.

If one is interested in the properties at low temperatures, it may be assumed that two-magnon excitations are rare and the spectrum can be obtained from H ₂, which can be solved directly by a Fourier transform. The spectrum then becomes^{[Reference Torelli and Olsen11]}

$$\varepsilon _n\lpar {\bi q} \rpar = \varepsilon _n^0 \lpar {\bi q} \rpar + A\lpar {2S-1} \rpar + \lambda SN_{nn},$$

where n denotes a band index (the range is equal to the number of magnetic atoms in the unit cell) and $\varepsilon _n^0 \lpar {\bi q} \rpar$ is the spectrum without anisotropy, which satisfies $\varepsilon _{n_0}^0 \lpar 0 \rpar = 0$ for the lowest band n ₀. The term Δ₀ ≡ A(2S − 1) + λSN _nn is therefore the spin-wave gap, and magnetic order in 2D relies on Δ₀ ≠ 0. It should be noted that for S = 1/2, the single-ion anisotropy term cannot open a gap and anisotropic exchange is required for magnetic order.

At finite temperatures, the spin-wave approximation breaks down due to the presence of thermally excited magnons. The magnon interactions may be included through a mean field approximation of the fourth order term, which leads to a temperature-dependent (renormalized) spin-wave spectrum.^{[Reference Yosida8, Reference Gong, Li, Li, Ji, Stern, Xia, Cao, Bao, Wang, Wang, Qiu, Cava, Louie, Xia and Zhang9, Reference Tyablikov42]} In particular, the gap is effectively decreased by temperature effects and for a single site in the unit cell, it may be written as^{[Reference Torelli and Olsen11]}

(9)$${\rm \Delta} = {\rm \Delta} _0-\; 4A\langle n \rangle -\lambda N_{nn}\langle n \rangle, $$

where n is the Bose distribution of magnons averaged over the Brillouin zone. Since the Bose distribution depends on temperature as well as the magnon spectrum, one has to calculate the gap and spectrum self-consistently. In Fig. 1, we show the renormalized spectrum of CrI₃ calculated at different temperatures using the parameters of PBE + U from Table I. It is clear that the gap decreases as the temperature is increased. The Bose distribution of magnons will diverge at the point where the gap vanishes and this signals a phase transition in the present approximations. In Fig. 2, we show the spin-wave gap as a function of temperature, and we observe that it vanishes at T _c = 35 K. This appears to be in good agreement with the experimental value of 45 K. However, the agreement might be fortuitous and certainly depends on our choice of the PBE + U functional for obtaining the Heisenberg parameters.

Figure 1.

Renormalized spin-wave spectrum of CrI₃ calculated at different temperatures.

Figure 2.

Temperature dependent gap of CrI₃ calculated from renormalized spin-wave theory equation (9). The gap closes at T = 35 K, which signals a phase transition with loss of magnetic order.

It is far from clear that the mean field approach is a good approximation close to the critical temperature where a high density of magnons is expected. In fact, it is not even obvious that multiple-magnon interactions (terms beyond fourth order in the Hamiltonian) can safely be neglected. On the other hand, quantum effects tend to be quenched by thermal fluctuations at elevated temperatures and for the sole purpose of evaluating critical temperatures, a purely classical approach might be expected to work well. The average energy and magnetization of the classical Hamiltonian (3) can be obtained straightforwardly from Monte Carlo simulations, and in Fig. 3, we show an example of such calculations; again using the Heisenberg parameters obtained with PBE + U. The magnetization decreases monotonously with increasing temperature and drops abruptly to zero at T = 50 K. Similarly, the heat capacity has a sharp peak at T = 50 K, which thus signals a phase transition where magnetic order is lost. The critical temperatures resulting from such calulations can be tabulated for different lattices and values of anisotropy parameters. In Fig. 4, we display such a compilation for honeycomb, square and triangular lattices in case of λ = 0. In the limit of A → ∞, the in-plane spin components are quenched and the critical temperature approaches that of the corresponding Ising model with coupling parameter J. This is an exact condition, which is naturally obtained in a classical treatment, but breaks down in renormalized spin-wave theory.^{[Reference Torelli and Olsen11]}

Figure 3.

Magnetization per Cr atom (M _z in unit of Bohr magnetons) and heat capacity (dE/dT) of a monolayer CrI₃ obtained from Monte Carlo simulations. Both quantities show a clear phase transition at a T = 50 K where the magnetic order is lost.

Figure 4.

Critical temperature as a function of single-ion anisotropy calculated from Monte Carlo simulations of honeycomb, square, and hexagonal (triangular) lattices. The solid lines are fitted functions given by equation (10). Figure from Ref. ^{[Reference Torelli and Olsen11]} © IOP Publishing. Reproduced with permission. All rights reserved.

As it turn out all the MC calculations can be fitted to the function^{[Reference Torelli and Olsen11]}

(10)$$T_c = S^2T_c^{{\rm Ising}} \tanh ^{1/4}\left[ {\displaystyle{6 \over {N_{nn}}}log\left( {1 + \gamma \displaystyle{{{\rm \Delta}_0} \over {J\lpar {2S-1} \rpar }}} \right)} \right],$$

where γ = 0.033, N _nn is the number of nearest neighbors, and $T_c^{{\rm Ising}}$ is the critical temperature of the corresponding Ising model, which is given by 1.52, 2.27, and 3.64 in units of J/k _B for the honeycomb, square, and triangular lattices, respectively. This expression is useful for high-throughput computational screening of ferromagnetic compounds, since it only requires three Heisenberg parameters that are easily obtained from DFT as explained above. Such an approach has recently been applied to the Computational 2D Materials Database (C2DB)^{[Reference Haastrup, Strange, Pandey, Deilmann, Schmidt, Hinsche, Gjerding, Torelli, Larsen, Riis-Jensen, Gath, Jacobsen, Mortensen, Olsen and Thygesen43,Reference Haastrup, Strange, Pandey, Deilmann, Schmidt, Hinsche, Gjerding, Torelli, Larsen, Riis-Jensen, Gath, Jacobsen, Mortensen, Olsen and Thygesen44]} where 3712 2D materials were screened for magnetic properties and yielded a prediction of 17 novel stable 2D ferromagnetic materials.^{[Reference Torelli, Thygesen and Olsen15]} It should also be straightforward to conduct a screening study on experimental databases like the Inorganic Crystal Structure Database (ICSD)^{[Reference Allmann and Hinek45]} and the Crystallography Open Database (COD)^{[Reference Gražulis, Daškevič, Merkys, Chateigner, Lutterotti, Quirós, Serebryanaya, Moeck, Downs and Le Bail46]} by applying a measure that identifies exfoliable materials from three-dimensional parent compounds.^{[Reference Mounet, Gibertini, Schwaller, Campi, Merkys, Marrazzo, Sohier, Castelli, Cepellotti, Pizzi and Marzari47,Reference Larsen, Pandey, Strange and Jacobsen48]}

Easy-plane order

The above analysis assumes an easy-axis that coincides with the out-of-plane direction. If the ground state has a component of magnetization in the plane of the material and in-plane isotropy is assumed, the rotational freedom prohibits long range magnetic order due to the Mermin–Wagner theorem. Mathematically, the above analysis would then result in a negative spin-wave gap that implies the instability of a state with out-of-plane magnetization. For high-throughput calculations, the sign of the spin-wave gap can thus be used as a descriptor that determines if long-range order is possible at finite temperatures.

Although long-range order is forbidden for easy-plane magnets, finite size effects may give rise to a macroscopic magnetization if the anisotropy is large. For example, in the limit of A _i → −∞, the S ^z components of the spins become quenched and the model (3) effectively reduces to the XY model given by

(11)$$H_{{\rm XY}} = -\displaystyle{1 \over 2}\sum\limits_{ij} {J_{ij}\lpar {S_i^x S_j^x + S_i^y S_j^y} \rpar }. $$

Kosterlitz and Thouless have shown that this model exhibits critical behavior for all temperatures below a certain temperature T _KT, meaning that spin correlations decay with a power law dependence on distance. This implies that long range order may be observed in any macroscopic sample although the order strictly speaking vanishes for an infinite system. Assuming small spin deflections such that S_i · S_j = S ²cos(θ) ≈ S ²(1 − θ ²/2) in equation (11), yields the harmonic XY model (HXY) and a classical analysis then shows that the magnetization of a finite sample is given by^{[Reference José, Kadanoff, Kirkpatrick and Nelson49]}

(12)$$M = \left( {\displaystyle{1 \over {2N}}} \right)^{k_{\rm B}T/8\pi J},$$

where N is the number of sites in the system. However, the analysis neglects the possibility of spin vortices that effectively renormalize the magnetization. The energy of a single vortex diverges, but bound vortex/anti-vortex pairs have a finite energy and constitute fundamental excitations of the system at low energies. At the Kosterlitz–Thouless (KT) transition temperature, the vortex/anti-vortex pairs unbind and the system becomes disordered. In Fig. 5, we show the magnetization (12) as a function of temperature for different values of n. We also indicate $T_{{\rm KT}}^{{\rm HXY}} = 1.351\,J/k_{\rm B}$ above which the magnetic order is destroyed by unbound vortex/anti-vortex pairs. As shown by Bramwell and Holdsworth, the KT temperature is largely unaltered by finite size effects, although strictly speaking it only marks a true phase transition in the infinite system. This implies that the magnetization is lost above the KT temperature and we may take the critical temperature for ordering as $T_c\approx T_{{\rm KT}}^{{\rm HXY}}$. Only the magnitude of the magnetization below, $T_{{\rm KT}}^{{\rm HXY}}$ will depend on the size of a sample and will vanish in the limit of an infinite system. The decay of magnetization with system size is, however, very slow. For example, a sample with an area of 1 μm² and a distance between magnetic atoms of 3 Å, the magnetization is reduced by 60% at the KT transition compared to the value at T = 0. For a 1 mm², sample the reduction will be 80% and a macroscopic magnetization will thus be prominent at any temperature below $T_{{\rm KT}}^{{\rm HXY}}$. The phase transition in the XY model (without the harmonic approximation) is a bit lower (T _KT = 0.898 J/k _B), but the conclusion remains the same.

Figure 5.

The magnetization in units of S as a function of temperature in the HXY model. The dashed vertical line indicates the KT transition above which the spins become disordered due to unbound vortex/anti-vortex pairs.

In real materials, the condition of A → −∞ will, of course, not be satisfied, but part of the KT behavior may be expected in easy-plane 2D magnets with finite anisotropy. Macroscopic samples of 2D materials with an easy-plane may thus show magnetic order due to finite size effects even if the Mermin–Wagner theorem prohibits long-range magnetic order. In particular, the easy-plane 2D magnets, Cr₂Ge₂Te₆,^{[Reference Gong, Li, Li, Ji, Stern, Xia, Cao, Bao, Wang, Wang, Qiu, Cava, Louie, Xia and Zhang9]} NiPS₃,^{[Reference Kim, Lim, Lee, Lee, Kim, Park, Jeon, Park, Park and Cheong14]} and VSe₂ could exhibit magnetism of the KT type described above. Only VSe₂^{[Reference Bonilla, Kolekar, Ma, Diaz, Kalappattil, Das, Eggers, Gutierrez, Phan and Batzill5]} have been reported to show macroscopic magnetization at finite temperatures though and those measurements have subsequently been questioned.^{[Reference Wong, Zhang, Bussolotti, Yin, Herng, Zhang, Huang, Vinai, Krishnamurthi, Bukhvalov, Zheng, Chua, N'Diaye, Morton, Yang, Ou Yang, Torelli, Chen, Goh, Ding, Lin, Brocks, de Jong, Castro Neto and Wee19]} Nevertheless, the subtle correlations and intricate magnetic frustration reported for VSe₂ could indicate that this material exhibits rich magnetic phenomena, and the presence of spin vortices could hinder the direct measurements of magnetization and perhaps make it strongly dependent on the precise conditions under which the measurements are carried out. In the case of Cr₂Ge₂Te₆, the KT transitions are expected at ~ J/k _B, which is roughly 75 K,^{[Reference Gong, Li, Li, Ji, Stern, Xia, Cao, Bao, Wang, Wang, Qiu, Cava, Louie, Xia and Zhang9,Reference Xu, Feng, Xiang and Bellaiche23]} but the persistence of in-plane magnetic order is hindered by the exceedingly small anisotropy in this material.

Outlook

In this Prospective Article, we have shown how to obtain accurate predictions for critical temperatures of 2D ferromagnetic insulators using DFT combined with either renormalized spin-wave theory or classical Monte Carlo simulations. Both methods agree well with experiments for the cases of and CrI₃ and CrBr₃, but the renormalized spin-wave method is expected to break down in highly anisotropic systems. We also presented a universal, but simple analytical expression for the critical temperature that were fitted to MC simulations and only depends on the Heisenberg parameters as well the number of nearest neighbors in a given material. Thus, the critical temperature of any 2D ferromagnetic insulator can be estimated from just four DFT calculations.

A major caveat of this method is the fact that insulating magnetic materials often exhibit strongly correlated physics and the calculated Heisenberg parameters can be rather sensitive to the chosen DFT functional. As such DFT must be regarded as inaccurate for the prediction of Heisenberg parameters in these systems. At least there is a need for more systematic assessment of the performance of various functionals for the prediction of magnetic properties. Another—perhaps more fundamental—deficiency of the method is the energy mapping analysis, which is based on the classical Heisenberg model. The classical Monte Carlo simulations of Curie temperatures are justified because the thermal fluctuations quench the quantum fluctuations in the vicinity of the critical temperature, but the parameters are calculated from the ground state at zero K, and should be mapped to the quantum mechanical Heisenberg model. Without anisotropy, it is straightforward to obtain the quantum corrected value of the nearest neighbor exchange coupling constants. It is, however, challenging to generalize that analysis to the case where anisotropy is included, but we hope that the challenge will be taken up in near future such that a universal and rigorous energy mapping scheme can be defined.

In the present Prospective, we have focused on the calculation of transition temperatures for ferromagnetic order in 2D materials. A natural next step is then to extend the analysis to anti-ferromagnetic order in 2D, which has been observed in FePS₃^{[Reference Lee, Lee, Ryoo, Kang, Kim, Kim, Park, Park and Cheong7,Reference Wang, Du, Fredrik Liu, Hu, Zhang, Zhang, Owen, Lu, Gan, Sengupta, Kloc and Xiong50]} and NiPS₃.^{[Reference Kim, Lim, Lee, Lee, Kim, Park, Jeon, Park, Park and Cheong14,Reference Wildes, Simonet, Ressouche, McIntyre, Avdeev, Suard, Kimber, Lançon, Pepe, Moubaraki and Hicks51]} The calculation of Heisenberg parameters and Monte Carlo simulations would proceed exactly as in the ferromagnetic case and for any bipartite lattice (for example, square or honeycomb), the spin-wave analysis can be carried out by means of a Bogoliubov transformation.^{[Reference Yosida8]} However, for non-bipartite lattices (for example, triangular), anti-ferromagnetic exchange coupling constants will give rise to geometric frustration and the magnetic phase diagram may become extremely rich with several competing phases for a given set of Heisenberg parameters.^{[Reference Chernyshev and Zhitomirsky52,Reference Maksimov, Zhu, White and Chernyshev53]} Even in the absence of anisotropy, the classical ground state is non-collinear and the mechanism by which anisotropy opens a gap is much more complicated than in the ferromagnetic case. A thorough analysis of the effect of anisotropic exchange on the magnetic properties of triangular lattices has recently been carried out,^{[Reference Maksimov, Zhu, White and Chernyshev53]} but an analysis of the effect of single-ion anisotropy in non-bipartite lattices with anti-ferromagnetic exchange coupling still seems to be lacking.

The most important remaining question is how to deal with itinerant magnets in 2D. Most of the reported high-temperature 2D magnets like FePS₃,^{[Reference Lee, Lee, Ryoo, Kang, Kim, Kim, Park, Park and Cheong7,Reference Wang, Du, Fredrik Liu, Hu, Zhang, Zhang, Owen, Lu, Gan, Sengupta, Kloc and Xiong50]} VSe₂,^{[Reference Bonilla, Kolekar, Ma, Diaz, Kalappattil, Das, Eggers, Gutierrez, Phan and Batzill5]} and Fe₃GeTe₂^{[Reference Fei, Huang, Malinowski, Wang, Song, Sanchez, Yao, Xiao, Zhu, May, Wu, Cobden, Chu and Xu3]} are indeed metallic, but it is not at all clear if the Heisenberg model (3) can provide even a qualitatively correct description in those cases. Naively, one may expect that the Heisenberg model could provide a decent approach if the exchange coupling constants can be converged with respect to distance. However, there are several issues with such an approach; it is—for example—not clear what value to use for the spin S. Alternatively, one may try to apply Stoner theory of magnetism,^{[Reference Yosida8]} which comprises a more natural starting point for itinerant magnets, but there is no simple way to obtain critical temperatures in that case and the theory requires inclusion of a parametric Hubbard U that measures on-site Coulomb repulsion between electrons. The prediction of thermodynamic properties for itinerant magnets thus poses the biggest and most important theoretical challenge in the field at the moment.

2D materials comprise an extremely versatile class of compounds. There are endless possibilities of optimizing particular properties by bottom up design of van der Waals heterostructures composed of 2D materials. In the case of magnetic 2D materials, van der Waals heterostructures have already been shown to yield new magnetic properties. For example, multilayers of CrI₃ exhibit anti-ferromagnetic interlayer coupling and the coercive field for interlayer alignment range from 0.65 to 1.84 T depending on the number of layers.^{[Reference Huang, Clark, Navarro-Moratalla, Klein, Cheng, Seyler, Zhong, Schmidgall, McGuire, Cobden, Yao, Xiao, Jarillo-Herrero and Xu54–Reference Jiang, Wang, Chen, Zhong, Yuan, Lu and Ji56]} This implies that the magnetic properties of multilayer structures of CrI₃ can be tuned simply by varying the number of layers. The anti-ferromagnetic structure of multilayers also implies that multilayers can act as spin valves with an interlayer conductance that may change by several orders of magnitude depending on the inter-layer magnetic state.^{[Reference Cardoso, Soriano, García-Martínez and Fernández-Rossier57–Reference Song, Tu, Carnahan, Cai, Taniguchi, Watanabe, McGuire, Cobden, Xiao, Yao and Xu60]} An additional convenient feature of 2D materials is the ease at which they can be electrostatically gated, and for 2D magnets, the magnetic properties may change dramatically under the influence of a gate voltage.^{[Reference Huang, Clark, Klein, MacNeill, Navarro-Moratalla, Seyler, Wilson, McGuire, Cobden, Xiao, Yao, Jarillo-Herrero and Xu61–Reference Suárez Morell, León, Miwa and Vargas63]} Finally, bilayers of CrI₃ have symmetry-allowed nearest neighbor Dzyaloshinskii–Moriya interactions that could give rise to Skyrmions^{[Reference Liu, Shi, Lu and Anantram27,Reference Liu, Shi, Mo and Lu64]} or perhaps topological magnons.^{[Reference Mook, Henk and Mertig65,Reference Pershoguba, Banerjee, Lashley, Park, Ågren, Aeppli and Balatsky66]} For other materials, such as Cr₂Ge₂Te₆ and NiPS₃, the bulk intra-layer magnetic order (ferromagnetic and anti-ferromagnetic, respectively) is preserved down to a bilayer, but vanishes in the monolayer limit. This is due to the easy-plane nature of the anisotropy, by which the long-range order is lacking for a monolayer, whereas even the bilayer structure bypasses the Mermin–Wagner theorem by weak interlayer magnetic interactions. It remains to be seen, however, if a remnant of KT physics can be observed in monolayers of such materials.