Onset of Magnetic Order in Strongly-Correlated Systems from ab initio Electronic Structure Calculations: Application to Transition Metal Oxides

We describe an ab initio theory of finite temperature magnetism in strongly-correlated electron systems. The formalism is based on spin density functional theory, with a self-interaction corrected local spin density approximation (SIC-LSDA). The self-interaction correction is implemented locally, within the KKR multiple-scattering method. Thermally induced magnetic fluctuations are treated using a mean-field `disordered local moment' (DLM) approach and at no stage is there a fitting to an effective Heisenberg model. We apply the theory to the 3d transition metal oxides, where our calculations reproduce the experimental ordering tendencies, as well as the qualitative trend in ordering temperatures. We find a large insulating gap in the paramagnetic state which hardly changes with the onset of magnetic order.


Introduction
Many materials characterised by strong electron-electron correlations are of technological interest. In the case of dilute magnetic semiconductors and the colossal magnetoresitive manganites, amongst others, this interest stems from their potential application in spintronic devices. As such, it is important to understand the interplay between charge and spin degrees of freedom of the many interacting electrons in this class of materials. This requires a full quantum description of the electron spin. For low temperatures, a magnetic material has an electronic structure which has a fixed spin-polarisation, e.g. a uniform spin-polarisation for a ferromagnet or fixed sublattice spin polarisations for an antiferromagnet. With increasing temperature, spin fluctuations are induced which eventually destroy the long-range magnetic order and hence the overall spin polarization. These collective electron modes interact as the temperature increases, depending upon and affecting the underlying electronic structure. The implications can be explored by invoking a time scale separation between the fast electronic motions and the, much slower, spin fluctuations. For intermediate times, τ , the spin orientations of electrons leaving an atomic site are sufficiently correlated with those arriving that the magnetisation, averaged over τ , is nonzero. 'Local-moments' are thus established. Although set up by the collective motion of the interacting electrons these can be described as classical spin like variables,{ê i }, provided the temperature is not so low that the dynamics of the spin fluctuations become important. So, for most finite temperatures, ensemble averages over the static orientational configurations of the local moments determine the magnetic properties of a system. In the so-called disordered local moment (DLM) approach [1], these averages are carried out using a mean-field technique.
A first-principles implementation of the DLM method [2], based on density functional theory (DFT) in the local spin density approximation (LSDA) and a multiple scattering effective medium method to handle the local moment disorder, has proved to be extremely successful at describing the magnetic properties of many transition metal systems [3,4,5]. In these systems, however, charge correlations are weak and the electrons are fully itinerant. This is certainly not the case in strongly-correlated systems, such as the transition metal oxides (TMO), which often contain highly localised electron states for which the LSDA fails.
The self-interaction corrected local spin density approximation (SIC-LSDA) [6], can account for correlation effects such that the ground state properties of the TMOs are well reproduced [7,8]. Recent papers by Lueders et al. [10] and Daene et al. [11] set out how to implement this approach using a multiple scattering Kohn-Korringa-Rostoker (KKR) method. This implementation of the SIC, the so-called local SIC (LSIC) [10] offers an immediate generalisation to disordered systems and opens up the possibility of implementing the DLM picture of magnetism within such a SIC based electronic structure scheme.
Indeed recently [12], we showed how such a DLM-SIC approach provides an accurate description of finite temperature magnetism in the heavy rare earths. In these systems, the magnetic moments arise from highly localised 4f states and are coupled via an indirect RKKY exchange mechanism [13], mediated by the sd conduction electrons. In this paper we apply the DLM-SIC to the TMOs, where a superexchange mechanism, mediated by the oxygen ions, operates between local moments set up by the 3d states of the transition metal ions. We investigate the spin fluctuations that characterise the paramagnetic state of these systems above their magnetic transition temperatures, gaining information about the type of ordering that is likely to occur as the temperature is lowered through a phase transition. This finite temperature study is complementary to calculations of the ground state energies of the same systems in different magnetic states at T = 0K carried out by several groups [7,8,9]. The work by Daene et al. [11] is particularly relevant.
The paper is organised as follows. In section 2 we outline the DLM formalism and its first-principles SIC-LSDA implementation within the KKR method. We examine in detail the electronic structure of MnO in the paramagnetic state in section 3, and, by evaluating the paramagnetic spin susceptibility, investigate the presence of any underlying ordering tendencies. We extend our study to the whole TMO series in section 4, where we find good agreement with experiment, with the notable exception of the Nèel temperature of NiO, which we find to be too small by a third. The paper is summarised in section 5, where we include some discussions as to what additional aspects, not accounted for in our present formalism, may be needed in order to capture fully the physics of NiO.

Formalism
We start with a key assumption of a time scale separation between that associated with fast electronic motions, i.e. the electron hopping time scale, and that characteristic of, typically much slower, spin fluctuations. At an intermediate time scale well defined moments exist at all lattice sites, the orientations of which we describe using a set of unit vectors, {ê i }. The local moment phase space specified by {ê i } is assumed to be ergodic and hence long time averages can be replaced by ensemble averages. These averages use the Gibbsian measure where the partition function, Z, is given by and β = (k B T ) −1 . Ω({ê i }) is a 'generalised' grand potential, the term 'generalised' meaning that Ω({ê i }) is not associated with a thermal equilibrium state.
In the DLM approach [2,4] a mean-field approximation for Ω({ê i }) is constructed by expanding it about a single-site reference spin Hamiltonian, Ω 0 ({ê i }) = − i h i ·ê i . Here, the parameters h i play the role of a Weiss field and are determined using a Feynman variational approach [14], whereby the free energy of the system, F = −β −1 logZ, is minimised. This free energy takes into account both the entropy associated with transverse spin fluctuations and also the production of electron-hole pairs associated with Stoner excitations.
The probability function, P 0 , associated with Ω 0 can be written as a product of single-site measures: where and Specifying these single-site probabilities, P i , means that a class of mean-field theories, developed originally for the study of substitutionally disordered alloys, becomes available to us to treat the local moment disorder. In particular, we deploy the coherent potential approximation (CPA) [15], which is known to be the best single-site approximation.
Here, an effective medium is specified. The motion of an electron through this medium approximates, on the average, the motion of an electron in the disordered lattice. This effective medium is determined self-consistently and in the context of multiple-scattering theory, i.e. the so-called KKR-CPA [16,17], the self-consistency condition states that there should be no further scattering of an electron, on the average, when a single site in the effective medium is replaced by one of the constituent 'alloy' potentials. For the local moment disorder, this condition reads mathematically as where quantities with a tilde superscript (˜) are 2 × 2 matrices in spin space andD i (ê i ) are the so-called CPA projectors, defined by (in L(≡ l, m) representation) Here, the single-site matrixt(ê i ) describes the scattering from a site with local moment orientated in the directionê i such that whereσ x ,σ y andσ z are the three Pauli spin matrices defined according to the global z-axis. In the local reference frame, where the z-axis is aligned withê i , we evaluate the matrices t + /t − , representing the scattering of an electron with spin parallel/antiparallel to the local moment directionê i . These matrices are calculated according to where the phaseshifts δ L (ǫ) are computed using effective DFT potentials. These effective potentials, v + and v − , differ on account of the 'local exchange splitting' [2], which is the cause of the local moment formation. Unlike the conventional LSDA implementation, the potentials v + /v − are orbital dependent in our new SIC-LSDA approach. This dependency comes about by our SI-correcting certain L channels, the details of which can be found in reference [11]. Importantly, the SI-corrected channels of v + and v − may differ. Indeed the valence channels to which we apply the self-interaction correction are those with a resonant phase shift [10,11]. Such resonant behaviour is characteristic of well localised electron states, which will establish quasi-atomic like moments. Through the influence they exert on the electron motions, these moments will be reinforced by the spins of more itinerant-like electrons. It thus follows that resonant states will tend to define the local moment orientation and, as such, we expect to SI-correct a greater number of channels of v + than we do for v − . t c i describes a site of the CPA effective medium. The scattering-path matrix, τ (where underlined matrices are in site representation), is related to the single-site scattering matricest viã where ǫ is a complex energy and g is the structural Green's function, which describes the free propagation of an electron between scattering centres. In the paramagnetic regime,t c = t c1 andτ c,00 = τ c,001 , the local moments have no preferred orientation and the P i 's become site independent. Moreover using Eq. 8 we findD where The superscript 0 signifies that the CPA projector is evaluated in the paramagnetic state. Substituting P i (ê i ) = P 0 = 1 4π into Eq. 6 we obtain 1 4π which becomes, on carrying out the integration, Equation 14 is evidently just the CPA equation for a system with 50% of moments pointing 'up' and 50% pointing 'down', i.e. an Ising-like system. The electronic structure problem is thus reduced to that of an equiatomic binary alloy, where the two 'alloy' components have anti-parallel local moments. Treating this 'alloy' problem with the KKR-CPA, in conjunction with the LSIC charge self-consistency procedure outlined in references [10,11], we arrive at a fully self-consistent LSIC-CPA description of the DLM paramagnetic state.
It should be noted that the equivalence of the DLM electronic structure problem to that of an Ising-like system is purely a consequence of the symmetry of the paramagnetic state, and is not the result of our imposing any restriction on the moment directions.
Indeed, in the formalism for the paramagnetic spin susceptibility, which we outline now, we maintain and consider the full 3D orientational freedom of the moments.
Within the DLM method the magnetisation at a site i, M i , is given by M i = µm i , where m i = dê i P i (ê i )ê i and µ is the local moment magnitude, determined selfconsistently. In the paramagnetic regime, where P i is independent ofê i , M i = 0. Using a perturbative approach, we investigate the onset of magnetic order, where M i = 0 becomes finite. In particular, we consider the response of the paramagnetic state to the application of an external, site-dependent magnetic field. Focusing on the dominant response of the system to line up the moments with the applied field, we obtain the following expression for the static spin susceptibility: where S (2) is the direct correlation function for the local moments, defined by In the paramagnetic state, S ik depends only on the vector difference between the positions of sites i and k. A lattice Fourier transform can hence be taken of Eq. 15, giving An expression for S (2) (q), involving scattering quantities obtained from the electronic structure of the paramagnetic state, can be found in reference [18]. By investigating the wavevector dependence of the susceptibility we gain information about the spin fluctuations that characterise the high temperature paramagnetic state. Conversely most first-principles calculations for the TMOs have, up to now, concentrated on the groundstate. These T = 0 calculations have demonstrated the importance of including strong electron correlation effects. In particular for conventional LDA calculations, in which such effects are neglected, the size of the groundstate magnetic moments are found to be substantially underestimated. When the effects of strong Coulomb repulsion between electrons occupying the partially filled d states are taken into account, such as through the self-interaction correction [7,8,9] or LDA+U method [42], the resulting moments are found to be in much better agreement with experiment.
In some respects, the incorporation of correlation effects is even more important when dealing with finite temperatures. In particular, it is essential to describe the localised nature of the d states in order that local moments survive into the paramagnetic phase. Indeed, when strong correlation effects were neglected in our DLM calculations we were not able to stabilise a local moment for NiO. For the other TMOs a local moment could be stabilised but, contrary to experimental observations, a large magnetovolume effect was exhibited.
Recently, using the LDA+DMFT method [20], it was shown how, by taking suitable account of strong electron-electron correlations, a reasonable description of the electronic spectrum of NiO in its paramagnetic state can be obtained [22]. In another LDA+DMFT study [43], in addition to NiO, calculations for other TMOs have also been performed. Our DLM-SIC investigation of the TMO series, which we go on to describe now, focuses on the onset of magnetic order and can be considered equivalent to these studies, but without quantum fluctuations. In particular, we concentrate on the energy scales associated with magnetic (spin) fluctuations, enabling us to obtain estimates of the Néel temperatures, and we investigate the importance of the dynamical quantum fluctuations.

MnO
In this section we discuss in detail our calculations for MnO. Using the LSIC-CPA approach outlined in section 2, we perform an electronic structure calculation for the paramagnetic state. We use the atomic sphere approximation (ASA) with a unit cell that has, in addition to a manganese and an oxygen ion, two empty spheres so as to obtain a better representation of the charge density and space filling. Daene et al. [11] compare the energies of the transition metal oxides for specific magnetically ordered states at T = 0K between different SIC configurations corresponding to different numbers of localised states. The configuration of lowest energy is determined by the balance between localisation (self-interaction) and band formation (hybridisation) energy. We follow the same procedure here but for the high temperature paramagnetic state. For MnO, we find the energy to be minimised when all Mn d states with spins parallel to the local moment direction are SI-corrected, with none corrected in the opposite spin channel. This picture is in keeping with Hund's first rule.
In Figure 1 we show the local DOS, where an exchange-splitting is evident. Of course, when an average is taken over all spin orientations, the electronic structure does not have an overall spin-polarisation. Nevertheless, it is possible to identify such 'local exchange splitting' experimentally, using photoemission [27] and inverse photoemission [28] techniques. The local moment obtained for the Mn sites in the DLM paramagnetic state differs from that in the groundstate by ≈ 0.03µ B (see Table 1). Such a small change between the ordered (zero temperature) and disordered (high temperature) states justifies our DLM picture, where the magnitude of the 'local exchange splitting' is independent of the orientational configurations of the moments. Furthermore this feature causes the spin-summed electronic structure to show little difference between the paramagnetic state above T N and the magnetically ordered ground state. Notably a sizeable band gap at the Fermi energy persists into the paramagnetic phase. Figure 2 shows the electronic structure of the paramagnetic (DLM) state in more detail and demonstrates the wave-vector, k, dependence of the local exchange splitting. This figure is constructed from calculations of the Bloch spectral function, A B (k, E). For non-interacting electrons in an ordered system at T = 0K, A B (k, E) comprises a set of Dirac delta functions which trace out the electronic band structure. When electron interaction, finite temperature or disorder effects are included the spectral function is a set of broadened peaks describing quasi-particle excitations. The broadening of the peaks shown in Figure 2 is a consequence of the local moment disorder in the paramagnetic state [29].
In figure 3 we present the results of our paramagnetic spin susceptibility calculations for MnO. These show the paramagnetic state to be dominated by spin fluctuations with wavevector q max = (0.5, 0.5, 0.5) (in units of 2π/(lattice constant)). This indicates that the system will order into an antiferromagnetic type II (AF2) structure, where moments within a <111> layer are aligned but are antiparallel in successive layers. This concurs with the experimentally observed groundstate of this system and also T = 0K calculations [7,8,11], where the most stable structure was determined by comparing the total energies of different magnetic configurations.
We examine the temperature dependence of the static spin susceptibility χ(q max ), in particular looking for a divergence which would signify that the paramagnetic state becomes unstable with respect to the formation of a spin density wave, characterised by the wavevector q max . For the theoretical (experimental) lattice parameters, such a divergence occurs at 102K (103K). This mean field estimate of the Néel temperature is in good agreement with the experimental value of 118K (see Table 1).

TMO series
In this section we extend our study to the other TMOs which also occur in the cubic rocksalt structure above the Néel temperature, namely FeO, CoO and NiO. We find the energies to be minimised when all TM d states with spins parallel to the local moment direction are SI-corrected together with one, two and three t 2g d-states corrected in the opposite spin channel for FeO, CoO and NiO respectively. This is again in accord with Hund's rule. The magnetic moments obtained are listed in Table 1 and are found to agree closely with their groundstate values. This concurs with experimental data for NiO [30], where the local moment size remains essentially unchanged as the Néel temperature is passed through. Indeed, more generally, the transition between the magnetically ordered and disordered phases is known to have little effect on the valence band photoemission spectra [31], and this is reflected in our spin-summed density of states which hardly changes between the paramagnetic state above T N and the magnetically ordered ground state at T = 0K. The insulating gaps of the paramagnetic DLM states are also very close in magnitude to those found in our calculations of the magnetically ordered states [11]. We find the gaps' sizes, in eV, of the paramagnetic states to be 3.3 (3.7) for MnO (see Figs.1,2), 3.5 (2.5) for FeO, 2.8 (2.4) for CoO and 3.8 (4.1) for NiO, where the experimental values shown in brackets are taken from the summary in ref. [7] and also [21] Our paramagnetic susceptibility calculations indicate that, like MnO, the other members of the TMO series have a tendency to order into the AF2 structure. The temperatures at which we predict this ordering to occur are listed in Table 1, where we find good agreement with experiment, with the exception of NiO where we underestimate the temperature by about a third. Interestingly, the Heisenberg mean-field results from a recent investigation of the TMOs [32], where a Hartree-Fock treatment of the electrons was used, show a similar discrepancy for the Néel temperature of NiO. Also, recent quasiparticle self-consistent GW [33] results for MnO and NiO show similar behaviour. This suggests that some additional physics, not at work in the other transition metal oxides, may be of relevance to the determination of the ordering temperatures of NiO. We return to discussing this point in section 5. In Table 1 we also give the Néel temperatures obtained from a muffin-tin (MT), as opposed to an ASA, implementation of our electronic structure scheme. Qualitatively, the trend in temperatures is the same, although the MT implementation gives generally better agreement with experiment.

Conclusion
We have used a locally self-interaction corrected implementation of the disordered local moment model to study the onset of magnetic order in MnO, FeO, CoO and NiO. Specifically, by taking into account the strong intrasite Coulomb correlations between the 3d electrons through SIC, we obtained an accurate first-principles, finite temperature description of the paramagnetic state. We found, in agreement with recent DMFT Table 1. Local magnetic moments and Néel temperatures for the series of transition metal oxides, obtained from our DLM-SIC calculations. The magnetic moments are compared to those obtained in reference [11] for the groundstate (AF2) configuration. Values not enclosed (enclosed) in round brackets were calculated using theoretical (experimental) lattice parameters. The variation of the Néel temperatures with respect to changes in the lattice spacing is also given. a Taken from ref. [7], for detailed references see references therein. b Reference [34]. c Reference [35].

Compound
studies of Ren et al [22] and Wan et al [43], that the paramagnetic state from our DLM-SIC study is characterized by a large insulating gap (see Figs. 1 and 2 for MnO). We conclude that this gap is the consequence of the strong electronic correlations in the paramagnetic state which the DLM-SIC describes particularly well without the need of including dynamical fluctuations. Also the local moment obtained for the transition metal sites in the DLM-SIC paramagnetic state differs little from that in the ground state (e.g. by ≈ 0.03µ B for MnO). We see also from the present work that the magnetic order in the AF2 ground state has no significant influence on the spin-summed electronic structure of NiO. This is consistent with our previous results that the band gap of AF1, AF2 and F magnetically ordered NiO differ only slightly [9] mostly on account of changes in hybridization. With the successful application of the DLM-SIC to the transition metal oxides we have been able to demonstrate that it is a generalized framework, capable of computing finite temperature magnetic interactions for correlated systems as well as the itinerant systems [2,4], where the SIC reduces to LSDA. The underlying magnetic ordering tendencies that our DLM-SIC study revealed were found to be in agreement with the experimentally observed magnetic groundstates. The corresponding Néel temperatures were, with the notable exception of NiO, also in good agreement with experiment. In order to begin to understand our theoretical underestimate of the Néel temperature of NiO, it is informative to look at the Néel temperatures of each of the TMOs as a function of lattice spacing. We found the temperatures to vary approximately linearly with lattice spacing and, as shown in Table 1, this dependency becomes more pronounced as the TMO series is crossed. This in turn implies that, as the series is crossed, the Néel temperatures become more sensitive to the underlying electronic structure. In particular, the hybridisation between the strongly-correlated d states and the exchange-mediating oxygen (s and p) sites becomes more delicate. This suggests that some of the error in the Néel temperature might be due to the incorrect positioning of these states. Indeed, although the LSIC does a much better job than the LDA at reproducing the insulating gaps of TMOs, the LSIC band gaps reported in [11] and here, although qualitatively correct, are in disagreement with experiment by up to 40% (FeO). In the recent LDA+DMFT study of Ren et al [22], the band gap of NiO was reproduced with very high accuracy, of course subject to an appropriate choice of the Hubbard U parameter. In that investigation a dynamic self-energy was used to incorporate correlation effects in contrast to the static self-energy used in the SIC approach. Our DLM treatment of essentially transverse magnetic fluctuations can also be considered as the static limit of some, as yet undeveloped, dynamical mean field theory of spin fluctuations. Since, however, we deal with temperatures where the dynamics of transverse spin fluctuations should not be important, our static treatment of these electronic degrees of freedom is a reasonable approximation. The effects of the faster dynamical charge correlations and longitudinal spin fluctuations, on the other hand, may be significant even at high temperatures and a more sophisticated account of these might therefore improve our estimate of the Néel temperature in due course.
In this paper we have focussed on the paramagnetic states of the TM-oxides. Here the crystal structure is that of cubic rock salt. It has been well-documented [23,24] how the establishment of AF2 magnetic order prompts a distortion of the crystalline lattices in these materials below T N . By including relativistic effects, such as spinorbit coupling, into our description of the TMOs in principal it should be possible to deduce how these magnetostrictive effects arise at T N as the symmetry is broken. The inclusion of spin-orbit coupling would also determine how the magnetic moments orient with respect to the crystal lattice as magnetic order develops. [25,26] This fundamental approximation made in our investigation of the neglect of spin-orbit coupling needs some further comment. For transition metals this is a good approximation since crystal field effects quench the 3d orbital moments [36]. In TMOs, however, strong Coulomb correlations can lead to a reduction of these crystal field effects and hence a preservation of the orbital moment. Indeed, the presence of an orbital moment in CoO has long been established [37]. Recent experimental [38] and theoretical [39] work has suggested that a significant orbital moment is also exhibited by NiO, with an orbital to spin angular momentum ratio as high as L/2S=0.17 [38]. There are important implications associated with the presence of an orbital moment, in particular with regards to which orbitals are occupied. More specifically, for a doublet groundstate the orbital angular momentum is completely quenched [40] and hence the presence of an orbital moment means that the d states cannot be in a pure (t 2g ) 6 (e g ) 2 configuration. A recent experimental study of NiO [41] has suggested that, in the groundstate, there is a fractional occupation of orbitals by unpaired spin electrons. Due to a very slight rhombohedral distortion, brought about by antiferromagnetic ordering, the symmetry changes from cubic to rhombohedral below the Néel temperature. As a result, the t 2g orbital splits into a single a g and doubly degenerate e g level, with the original e g level retaining its symmetry characteristics but renamed e ′ g . In reference [41] it was reported that the a g level has a partial filling of 1.69, giving rise to an orbital moment µ L = 0.31µ B . In our calculations, although we allow the energetics to determine the orbital configurations, we are restricted to integer occupancies of the SI-corrected, localised, orbitals. In order to consider fractional occupancies it is feasible to use the CPA to describe a random distribution of orbital configurations, analogous to our treatment of random moment orientations described in section 2. Through careful choice of the probabilities associated with the different orbital configurations, various fractional occupancies can be mimicked to be followed by an investigation of orbital ordering in the TMOs. This type of study is being planned and holds promise for a further contribution to the understanding of these materials.