Derivation of the spin Hamiltonians for Fe in MgO

A method to calculate the effective spin Hamiltonian for a transition metal impurity in a non- magnetic insulating host is presented and applied to the paradigmatic case of Fe in MgO. In a first step we calculate the electronic structure employing standard density functional theory (DFT), based on generalized-gradient approximation (GGA), using plane waves as a basis set. The corresponding basis of atomic-like maximally localized Wannier functions is derived and used to represent the DFT Hamiltonian, resulting in a tight-binding model for the atomic orbitals of the magnetic impurity. The third step is to solve, by exact numerical diagonalization, the N electron problem in the open shell of the magnetic atom, including both effect of spin-orbit and Coulomb repulsion. Finally, the low energy sector of this multi-electron Hamiltonian is mapped into effective spin models that, in addition to the spin matrices S, can also include the orbital angular momentum L when appropriate. We successfully apply the method to Fe in MgO, considering both, the undistorted and Jahn-Teller (JT) distorted cases. Implications for the influence of Fe impurities on the performance of magnetic tunnel junctions based on MgO are discussed.


I. INTRODUCTION
Understanding the electronic properties of magnetic transition metals embedded in diamagnetic hosts plays a central role in several branches of condensed matter physics and materials science. The presence of transition metal impurities is known to modify the electronic properties of insulators [1], semiconductors [2] and molecular crystals [3]. Thus, diluted semiconductors become paramagnetic and their optoelectronic properties, such as the photoluminescence spectrum become extremely sensitive to the application of magnetic fields, resulting in the so called giant Zeeman splitting [2]. In turn, the electronic and spin properties of the magnetic atoms are very sensitive to their environment [1]. This permits inferring local information about the host by means of spin probing techniques such as electron paramagnetic resonance [1].
Very often, the spin properties of a magnetic system are described in terms of effective single spin Hamiltonians [1,3] built in terms of atomic spin operators only. Whereas the symmetry of a given system determines which terms are possible in an effective spin Hamiltonian, prediction of the values of the various parameters can be a difficult problem. Extraordinary progress in instrumentation techniques makes it now possible to probe individual magnetic atoms in a solid state environment [4,5] using a variety of techniques, such as scanning tunneling microscope (STM) inelastic electron spectroscopy (IETS) [6,7], and single quantum dot photoluminescence * Permanent address: Departamento de Física Aplicada, Universidad de Alicante. [8,9]. These techniques permit assessing the delicate interplay between spin properties of the transition metal and electronic and structural properties of the local environment at the atomic scale [7,10] and motivate the quest of quantitative methods to address this interplay.
Conventional density functional theory [11,12] provides an accurate description of the electronic properties of the ground state of solids but it does not provide a direct route to describe the fine details of the low energy spin excitations inherent in magnetic atoms in insulating hosts. For instance, the ground state of the effective Hamiltonian in conventional functionals in DFT is a unique Slater determinant with broken spin symmetry, which is fundamentally different from the multiplet nature of the real system. In this context, we find it convenient to have a constructive theoretical approach to derive the effective spin Hamiltonian, starting from an atomistic DFT description of the electronic properties of the system, but describing the electronic properties of the system with a multi-electron approach that captures the multiplet nature of the relevant electronic states.
Here we propose a method to obtain an effective spin Hamiltonian for a magnetic atom in an insulating host, starting from density functional calculations, in four well defined steps. First, a density functional calculation of the electronic properties of the magnetic atom inside the non-magnetic host. The second step is to represent the effective DFT Hamiltonian with a basis of localized atomic orbitals, which allows obtaining the crystal and ligand fields terms of the atomic orbitals of the relevant open shell of the magnetic atom, defining thereby a multiorbital Hubbard Hamiltonian. Since our DFT approach makes use of a plane-wave basis, we implement this step by means of the wannierization [13] technique. Up to this point, the methodology is very similar to previous work [14][15][16][17][18][19][20][21]. In a third step we add to the Hubbard model the intra-atomic Coulomb repulsion and the spin-orbit coupling for the electrons in the open-shell. The final step is a symmetry analysis of the spectrum and wave functions, obtained by numerical diagonalization of the effective Hubbard model. The resulting multi-electron states analysis permits constructing an effective spin Hamiltonian for the system.
Below we describe in more detail the method and apply it to the paradigmatic case of Fe 2+ as a substitutional impurity of Mg in MgO [1], a band insulator. The spin properties of this system have been studied in detail by means of several techniques, including far infrared spectroscopy [22], acoustic paramagnetic resonance [23], infrared spectroscopy [24], and XPS [25]. The interplay between atomic structure and spin properties is beautifully illustrated in this system: we consider both the case of undistorted Fe/MgO, where the octahedral symmetry does not quench completely the orbital angular momentum L of Fe 2+ as well as the system with a Jahn Teller distortion in which case L is quenched, resulting in a very different type of effective spin Hamiltonian. Our findings might shed some light on recent results [7] probing a single Cobalt atom on MgO, that indicate that the orbital moment is partially quenched, in contrast with magnetic adatoms deposited on other surfaces. In addition, we are interested in Fe as a possible impurity in MgO tunnel barriers in magnetic tunnel junctions with Fe based electrodes [26,27] and we discuss how it could reduce the spin-filter properties, when compared to the ideal system.
The rest of this manuscript is organized as follows. In Sec. II we study the Electronic Structure of Fe 2+ as a substitutional impurity of Mg in MgO using DFT calculations. In Sec. III we discuss the derivation of the single-particle part of the magnetic atom Hamiltonian from the DFT calculation using the wannierization approach. In Sec. IV, we build and solve by numerical diagonalization the generalized Hubbard model and derive the effective spin Hamiltonians for two different geometries. In Sec. V we summarize the advantages and shortcomings of this work and discuss the effect of Fe impurities in MgO barriers on the magnetoresistance of magnetic tunnel junctions.

II. ELECTRONIC STRUCTURE: DFT CALCULATIONS
In this section we describe our DFT calculations, for pristine MgO as well calculations for super cells of Mg 31 O 32 Fe. For the super cells, we consider two geometries, with and without Jahn-Teller distortion of the Fe atom. In addition, and for reasons discussed below, we did both spin-polarized and spin-unpolarized calculations.
Our calculations were done using the generalizedgradient approximation (GGA) for exchange-correlation energy [28], using plane-wave basis sets and ultrasoft pseudopotential method for Mg and O, and Projector Augmented-Wave (PAW) [29] for Fe as implemented in QUANTUM ESPRESSO (QE) code [30]. Since we are interested in the spin-unpolarized calculation, there is no need to include the DFT+U correction. Although a proper DFT calculation of magnetism will require DFT+U corrections [34], here we do not include the Uterm at this level since, as discussed in Secs. III and IV, our approach to derive an effective spin Hamiltonian requires to start with a spin unpolarized DFT calculation to which we should add the many body Coulomb repulsion between the d-orbital electrons in the Fe.
For the case of the super cells, the number of k points was taken to be 8 × 8 × 8 and we used a Fermi-Dirac smearing with a broadening parameter of 0.0035 Ry. Finally we fixed the cutoff energies for the wave function and charge density at 65 Ry and 700 Ry respectively. The calculation of the magnetic atom in MgO is done using a 64 atoms supercell of Mg 31 O 32 Fe, with lattice parameter 2a. The supercell, including the Fe atom is shown in Fig.  1

(a).
MgO is an insulator with a NaCl type crystal lattice (see Fig. 1). Using the experimental lattice constant (a = 4.22Å), our DFT calculations give a band gap of 5.4 eV , below the actual value 7.8 eV [31][32][33]. This discrepancy is very common and quite close to other DFT calculations for MgO (see for example Ref. [34] where the calculated gap is 5.85 eV ). Our calculations show that the valence band of the MgO is mainly formed by O 2p orbitals and the conduction band by Mg 3p and 2s orbitals.
Our calculations show that the main effect of the Fe impurity on the MgO band structure is the appearance of 10 in-gap very narrow bands that, as we show below, are associated to the d orbitals of the Fe atom, see Fig. 2(a). Six of these levels are below the Fermi energy and, for spin-polarized calculations, correspond to 5 levels spin ↑ and 1 ↓ resulting in a spin S = 2.
In the ideal MgO-like bulk crystal, where the Fe substitutes a Mg atom, the Fe is in an octahedral environment surrounded by 6 oxygen neighbors. In this undistorted geometry, the Fe-d levels are expected to split in a lower energy triplet, t 2g , and a higher energy doublet e g , due both to the interaction with the charged neighbor oxygens (crystal field contribution) and the hybridization with the oxygen atomic orbitals (ligand field contributions).
In the undistorted octahedral environment, the Fermi energy lies exactly at the t 2g orbital triplet of the minority spin, so that the system has a orbital degeneracy that leads to Jahn-Teller instability [1,24,25], which we model by letting the system relax from an initial configuration in which Fe is slightly off the center of the octahedron. The distorted solution so found has lower energy than the undistorted one. In both cases, distorted and undistorted, the relaxation was performed until the forces acting on atoms were smaller than 10 −3 a.u.. In the undistorted octahedral environment, the O surround- ing atoms are all at 2.135Å from the Fe atom. In order to characterize the deformed configuration it is convenient to set Fe as the origin of coordinates and label oxygens as in the right panel of Fig. 1 [35]. Distortions happen to be symmetric, i.e., with δ R 1 = −δ R 4 , and we express them in terms of the normal modes of the octahedron. It turns out [35] that the computed distortion can be expressed as a linear combination of the breathing mode , which clearly preserves the octahedral symmetry, and the q 3 = 1 which singles out the z axis symmetrywise and preserves the planar square symmetry of the xy plane. The obtained distortion can be written as 0.01q 1 + 0.03q 3 , where q 1 are expressed inÅ and is said to be tetragonal [35]. It should be emphasized that we have not made a systematic attempt to study all possible Jahn-Teller distortions in this system. Instead, we are testing our method for a particular distortion, which is in line with previous work [25].
The effect of the tetragonal distortion is apparent in both, the spin-polarized, Fig. 2, and the spin-unpolarized cases, Fig. 3. The finite width ∼ 50 meV of the DOS peaks, much smaller than the crystal field splitting, is a finite size effect due hybridization of d orbitals between Fe atoms located at different unit cells. In both cases the t 2g triplet degeneracy is split into a doublet and a singlet, and the e g doublet is also split. Importantly, the tetragonal distortion does preserve the 4µ B magnetic moment (S = 2). However, the very different orbital arrangement will result in important differences in the spin Hamiltonian, discussed below.
The spin-polarized calculations discussed so far provide a mean-field-like description of the magnetism of Fe in MgO. However, in order to determine the parameters for the Hamiltonian in the multiplet configuration interaction calculation, presented in section IV, we start from a spin-unpolarized calculation, a strategy used as well in previous work [21]. It must be noted that, for spin-polarized calculations, the crystal field splitting ∆ is spin dependent, which is clearly a feature of a mean field solution that breaks spin-rotational symmetry. In the distorted case the sign of the splitting of the t 2g levels, as well as the magnitude of the splitting of the e g levels are spin-dependent. Since it is convenient to have a spin-independent crystal field Hamiltonian, we have performed a spin-unpolarized calculation of Fe in MgO. For the undistorted case we obtain a ground state configuration (0e g , 6t 2g ) where all d−electrons of Fe occupy the degenerate (orbital and spin) states d xy , d yz and d xz [see Fig. 3a),b)]. The computed crystal field splitting is ∆ = 1.45 eV. For the tetragonal distortion, the spin unpolarized calculation still shows that 6 electrons occupy the t 2g levels, but the d xy level is now split from d xy and d yz , as shown in Fig. 3c).

III. CALCULATION OF THE CRYSTAL FIELD HAMILTONIAN USING WANNIER FUNCTIONS
The discussion of the previous section shows that it is possible to describe Fe in terms of 6 electrons occupying the in-gap levels, which are predominantly formed by Fe d orbitals. To do so, we would like to extract from the DFT calculation a one-body Hamiltonian projected over these d orbitals that includes their interaction with the host crystal. However, the DFT Hamiltonian is expressed in terms of Bloch waves that in our calculations are expressed as linear combinations of plane waves. In order to go to an atomic like description, we compute the so called maximally localized Wannier functions (MLWF) [13,[36][37][38][39] associated to the Bloch states of the DFT calculation, using the package Wannier90. The Wannier functions form an orthogonal and complete basis set that we use to express the Hamiltonian. Interestingly, we find 5 atomic-like MLWF with the same symmetry than the real ℓ = 2 spherical harmonics. Therefore, we take the representation of the DFT Hamiltonian in this subspace, as the crystal field Hamiltonian H CF (although it also contains ligand field contributions).
This wannierization [13,[36][37][38][39] procedure is implemented as follows. First, we select the group of Bloch bands from the spin unpolarized calculation for which the MLWF are calculated. For Fe/MgO, we take the valence bands as well as the 10 (counting spin) in-gap states. These groups of bands do not overlap in energy with others, so that it is not necessary to perform the disentanglement procedure [13,38]. Second, the Bloch states |ψ kn are projected over a set of localized functions. Based on the population analysis of the DFT calculation, we project over the atomic like d orbitals centered around Fe and p orbitals centered around oxygens. In total, there are 96 p orbitals (32 Oxygen atoms) and 5 d orbitals. After an iterative procedure, the MLWF are determined. Expectedly, the calculation yields five MLWF localized around the Fe atom that, in the neigh-borhood of the atom, have the same symmetry that the real spherical harmonics with ℓ = 2, as shown in the left panels of Fig. 4 (a,c and e). It is important to point out that the MLWF are not strictly identical to the atomic orbitals, because the tail of the wave functions have a different symmetry, as shown in the right panel Fig. 4 (b,d and f).
The representation of the DFT Hamiltonian in the basis of the MLWF yields a tight-binding Hamiltonian whose energy bands are identical to the valence and ingap bands obtained from DFT. For the purpose of this work, we are interested in the intra-cell Hamiltonian: where H dd has dimension 5, corresponding to the d orbitals of Fe, and H pp has dimension 96, corresponding to the 3 p orbitals of the 32 oxygen atoms in the unit cell.
The H dd part describes the crystal field splitting of the d levels. For the undistorted case, it describes the t 2g triplet and e g doublet, separated by a crystal field splitting ∆ CF . Interestingly, diagonalization of H dd yields, for the undistorted case, ∆ CF = 0.83 eV, much smaller than the DFT splitting 1.45 eV, that is only recovered if the whole H DF T matrix is diagonalized. Thus, we see that this approach permits us to quantify the ligand and crystal field contributions to the splitting. We find that almost half of the t 2g − e g splitting comes from the so called [1] ligand field contribution, described by V dp , the hybridization between the d like orbitals and the p states that form the valence band of MgO.
In order to preserve a small dimension of the Hilbert space, so that the number of multi-electron configurations can be handled numerically, it is convenient to work with a truncated Hamiltonian for the d electrons only, but that includes their hybridization to the p levels. Such a Hamiltonian could be produced using degenerate second order perturbation theory for the different degenerate manifolds within the 5 d levels, discussed below: This second order Hamiltonian yields eigenvalues within 10 percent of the exact ones. It is possible to do better by realizing that the projection of the exact eigenstates of H DF T over the the d like MLWF is always higher than 80 percent, and in most cases higher than 90 percent. More important, the spectrum and wave functions projected over the MLWF of the 5 in-gap states can be described with: where l a are the ℓ = 2 angular momentum matrices, and d 2 , d 4 and a are obtained by fitting. Here we approximate the MLWF by the real spherical harmonics centered in the Fe ion. The same approximation is used in the calculation of spin-orbit and on-site Coulomb integrals later on. This methodology has been used before [40] with good qualitative results. In order to fit a, d 2 and d 4 we employ the analytical expressions of the eigenvalues of H CF : 18a+d 2 +d 4 , 18a+ d 2 + d 4 , 18a + 4d 2 + 16d 4 , 24a, 24a + 4d 2 + 16d 4 . For the undistorted case, the in-gap d levels obtained from diagonalization of Eq. (3) feature a triplet (t 2g ) and a doublet (e g ), and are fitted with d 2 = d 4 = 0, as expected from the octahedral symmetry. The t 2g − e g splitting is thus given by 6a, which yields a = 0.241 eV. For the JT distorted case, the t 2g triplet is split into a singlet and a doublet, see Fig. 3a), while the e g doublet is also split. The fitting yields a = 0.250 eV, d 2 = 0.461 eV and d 4 = −0.1 eV. The difference between the fitted and computed energy levels are always smaller than 2 meV.

IV. EFFECTIVE FEW ELECTRON HAMILTONIAN
In the previous section we have demonstrated that, starting from a DFT calculation for Fe in a supercell of MgO we are able to derive a crystal field Hamiltonian for the in-gap d levels, including both crystal and ligand fields contributions, expressed in a basis of localized atomic-like orbitals provided by the maximally localized Wannier Functions.
In this section we derive an effective spin Hamiltonian that accounts for the low energy spectrum of a magnetic impurity within the MgO. This is done in two stages. We first build and solve a Hamiltonian for the 6 electrons in the d levels or Fe, including the effect of crystal and ligand field as described at the DFT level, and adding the Coulomb and spin-orbit coupling interactions. This few electron problem can be diagonalized numerically. In the second stage we analyze the symmetry and properties of the low energy levels and propose an effective spin Hamiltonian that accounts for them. This is done for the undistorted and distorted configurations studied in Sec. II. By so doing, we arrive at effective spin models in agreement with the literature [1,24,25].

A. Multiplet calculation
We consider a Hamiltonian for the N = 6 electrons in the d orbitals of Fe in MgO that includes four terms, electron-electron, crystal-field and ligand field, spin-orbit and Zeeman interactions: The Coulomb term reads: where d † mσ (d mσ ) denotes the creation (annihilation) operator of an electron with spin σ in the ℓ = 2, ℓ z = m state of the magnetic atom, denoted by φ m ( r), assumed to be equal to the product of a radial hydrogenic function ( with effective charge Z and a effective Bohr radius a µ ) and a spherical harmonic. Thus, we are only considering d 6 configurations and leaving out pd 7 configurations.
It turns out that the Coulomb integrals V mnm ′ n ′ scale linearly with the value of V 0000 ≡ U . Explicit expressions for the on-site Coulomb integrals are given in the appendix. Specifically, U could be computed using Eq. (A7). Another option, followed here, is to consider U as an adjustable parameter. In this work we take U = 19.6 eV, which yields the correct splitting between the 3 P 2 excited state and the 5 D ground state of the free ion, measured [42] to be 2.41 eV.
The second term in Eq. (4) corresponds to the crystal and ligand fields Hamiltonian discussed in the previous section: with m|H CF |m ′ derived from DFT using the procedure described above and, which is a very good approximation, is given by Eq. (3).
The last term in the Hamiltonian describes spin-orbit coupling: where ζ is the single particle spin-orbit coupling of the d-electrons. It is also very frequently expressed as λ L · S with L the total angular momentum. For the case of Fe 2+ , both parameters ζ and λ are related by λ = −ζ/2S [1], with ζ = 50.1 meV and S = 2. The last term in Eq. (4) corresponds to the Zeeman Hamiltonian: where g = 2. So, if we assume that the CF term is given by Eq. (3), the multiplet Hamiltonian (4) depends on five energy scales: U , a, d 2 , d 4 and ζ as well as the magnetic field. For N = 6 electrons, the total number of d 6 configurations is 210. Therefore, numerical diagonalization of the Hamiltonian is straightforward. In agreement with Hund's rules, we obtain a ground state multiplet that maximizes S and L. Thus, the ground state, denoted by 5 D, has a degeneracy of (2L + 1) * (2S + 1) = 25, with L = S = 2. This low energy many body spectrum is fully independent of U provided that the crystal field is not high enough to mix the 5 D with the 3 P2 multiplet. This could change if d 7 p configurations are included.
In order to analyze the results, it is convenient to add the different terms in the Hamiltonian one by one, in order of importance: Coulomb U , the crystal field (a, d 2 ,d 4 ), and spin-orbit coupling (ζ). Thus, in a first step the problem is solved considering only H Coul . In this case the Hamiltonian commutes with S 2 and L 2 , the square of total spin and orbital angular momentum.

B. Undistorted Fe/MgO
We discuss first the case of Fe 2+ in MgO without Jahn-Teller distortion. The effect of the octahedral component (a) of the crystal field on the L = S = 2 multiplet is shown in Fig. 5. As a result of the breaking of the orbital rotational symmetry, L is no longer a good quantum number and the 2L + 1 degeneracy is partially lifted. As we turn on a, see Eq. (3), the 5 D levels of iron splits into two, an orbital Γ 5 triplet ground state, with total degeneracy 15, and an orbital doublet excited state, 10 times degenerate see Fig. 5a).
The 15-fold degeneracy of the ground state multiplet of the Fe in the octahedral environment of MgO can be interpreted as if the ground state multiplet had a L = 1 orbital momentum. Actually, the representation of the ℓ operator on the subspace of the t 2g orbitals is isomorphic to the ℓ = 1 operators multiplied by −1 [1]. When SOC is added to the Hamiltonian, the 15-fold degenerate ground state splits into a triplet, a quintuplet and a septuplet, in ascending energy order [see Fig. 5b)]. This pattern can be rationalized in terms of the following effective Hamiltonian where the total spin is coupled to the pseudo angular momentum L = −1 [1]: where λ = −ζ/2S. The first term naturally leads to a spectrum with multiplets associated toJ 2 = ( L + S) 2 : J = 1 (ground),J = 2 (first excited) andJ = 3 (third excited), with degeneracies 2J + 1 = 3, 5 and 7 respectively. The result of the calculation of the expectation values ψ|S z |ψ for ψ in the ground state triplet with J = 1, backs up the idea that the S = 2 spin is coupled to an pseudo-angular momentum with L = −1. In Figs. 5d) and e) we plot the expectation values of S z and L z for the 3 states of the ground state triplet as a function of the magnetic field. Notice that theJ z = ±1 and J z = 0 values are recovered by subtracting ψ|S z |ψ and ψ|L z |ψ , in contrast with the common case of a total angular momentum. The CI calculation for the 15 lowest energy states for Fe 2+ in the undistorted environment has some fine structure not captured by the first term in Eq. (9). In particular, the multiplets withJ = 1 and 2 have some fine structure [see Fig.5b)], that can be accounted by the a second term in the effective Hamiltonian This operator does not break the triple degeneracy of the ground state, breaks theJ = 2 into a triplet and a doublet (being isomorphic to the problem of the octahedral crystal field splitting of the ℓ = 2 orbitals), and theJ = 3 into a singlet and two triplets. In summary, in the undistorted case, our calculation portraiys Fe 2+ as a system with S = 2 and pseudo-orbital momentum L = 1 [1]. Spin-orbit coupling leads to a ground state triplet withJ = 1. The energy splitting to the first excited state, withJ = 2, is approximately linear in the atomic spin-orbit coupling, reflecting the fact that the octahedral symmetry quenches only in part the orbital angular momentum. Thereby, the effective model has to take into account L, and not only S. The Jahn-Teller distortion, that we discuss next, changes this situation.

C. Jahn-Teller distorted Fe/MgO
We now discus the effect of the tetragonal distortion on the multiplet structure of Fe 2+ in MgO. As discussed in Sec. II, this distortion introduces the uniaxial terms d 2 l 2 z + d 4 l 4 z in Eq. (3). The effect of the uniaxial terms on the many-body 15-fold degenerate ground state of the Hamiltonian with ζ = 0 and a = 0.250 meV is shown in Fig. 6a). It is apparent that the JT distortion splits these 15 states into a ground state quintuplet, corresponding to a S = 2 spin with quenched orbital momentum, and a excited manifold with 10 states. Thus, it takes a JT distortion on top of the octahedral crystal field to eliminate the extra 2L + 1 = 15 degeneracy of the Γ 5 orbital triplet. When spin-orbit coupling is added [ Fig. 6b)] the 2S + 1 = 5 degeneracy is broken into a singlet, a doublet, and a split doublet [see Fig. 7a)]. Finally, the Zeeman splitting breaks the remaining degeneracies, as observed in Fig. 6c).
Interestingly, the five low energy states, corresponding tol z = 0, can be described by an effective S = 2 Hamiltonian of the form The comparison of the spectra as a function of a magnetic field B z , calculated both with the full CI Hamiltonian and the effective spin model, is shown in Fig. 7a). The parameters of the effective Hamiltonian are obtained by fitting to the multiplet calculation. We obtain D = 0.734  Fig. 7b),c) as a function of B z . It is apparent that the ground state (black line) has S z = 0, as a result of the dominant uniaxial term DS 2 z favoring the minimum spin projection as ground state, S z = 0. The first excited doublet, split by B z , has S z = ±1. The S 4 x +S 4 y term couples the otherwise degenerate doublet S z = ±2, resulting in a quantum spin tunneling splitting. The mixing of the wave functions is apparent in the non-linear evolution of the expectation value of S z as a function of B z . At small field the magnetic moment is quenched. At higher field the Zeeman contribution overcomes the quantum spin tunneling. We note in passing that, in contrast with the S = 2 spin with C 2 in plane symmetry [3], in our case there is no quantum spin tunneling splitting within the S z = ±1 doublet, that remains degenerate.

V. DISCUSSION AND CONCLUSIONS
The results of the previous sections illustrate how, for the cases of Fe 2+ in MgO with and without Jahn Teller distortion, we have been able to derive effective spin Hamiltonians [Eq. (9) and Eq. (11)] that reproduce the spectra obtained from the few-electron Hamiltonian. The parameters are derived directly from a DFT calculation of the electronic structure of this system. We now list possible improvements for the method. In addition, we briefly discuss the implications for a technologically relevant system, MgO tunnel barriers with Fe electrodes [ 26,27,43,44], and present our conclusions.
A. Improvements for the method There are several ways in which the method presented in this manuscript could be improved. First, the approximation that the Wannier functions are atomic orbitals in the evaluation of the matrix elements of both the spinorbit coupling and Coulomb interaction could be avoided at the price of performing the numeric integration using the actual Wannier functions. This would also allow extending the method to situations in which the localized atomic orbital lives in interstitial sites, such as the technologically relevant [45,46] cases of NV centers in diamond [47], or Mg vacancies in the MgO [48,49]. Second, a more accurate quantitative description would require us to correct the double counting of some of the Coulomb interactions [19,[50][51][52]. Third, the Hilbert space in the multiplet calculation could be expanded in two ways, either including more intra-atomic configurations [7] , such as pd 5 , as well as configurations where the charge is transferred to the neighbor oxygen atoms [17,53]. Fourth, the GGA calculation underestimates the gap of insulators, which most likely has some influence on the d levels as well. The use of a hybrid functional, or of an approximation adequate to compute energy gaps, such as the GW method [54][55][56][57][58] would be an improvement, but the computational overhead for unit cells with tens of atoms is far from small. Finally, the method presented here could be improved obtaining U from first principles calculations [20,[57][58][59].

B. Influence of Fe impurities in the barrier MgO of a magnetic tunnel junction
We now briefly discuss some relevant consequences drawn from our calculation in the context of spin dependent transport in MgO magnetic tunnel junctions with Fe based electrodes such as Fe or CoFeB [26,27]. A key figure of merit of magnetic tunnel junctions is the magnetoresistance, defined as M R = 100 × (R AP − R P )/R P , where R P and R AP are the resistance for parallel and antiparallel orientation of the electrode magnetizations. A very large MR, exceeding 1000, was predicted for Fe/MgO/Fe MTJ [43,44]. Actual experiments in this system have found room temperature MR above 600 [60]. that have permitted a tremendous boost of this technology, but quite below the theoretical limit.
The very likely presence of substitutional impurities of Fe in the MgO barrier would affect transport in two ways, opening two additional tunneling channels in the magnetic tunnel junction. On one side, electrons could tunnel through the in-gap d levels (see Fig. 2). Elastic tunneling through these states is possible at large bias (≃ 1eV ), when the Fermi energy of one of the electrodes is set on resonance with the in-gap d levels. This would yield characteristic resonance line shapes at finite bias, not much different from those observed experimentally [61]. At small bias, electrons could still tunnel through these d levels through second order cotunneling processes, in which the transport electron would excite a spin transition between the low energy states of the Fe, within a range of a few meV, see Fig. 7 (a) . Whereas this process will give a much smaller contribution to transport, they are known to be an efficient [62] source of spin-flip. These problems will be addressed qualitatively elsewhere.

C. Summary
In summary we have presented a method to derive effective spin Hamiltonians for magnetic atoms inside insulators, starting from a DFT calculation based on plane waves. This is achieved by post-processing the DFT calculation to obtain the maximally localized Wannier functions, which, in the system considered here, happen to be atomic-like orbitals in the magnetic atom. Expressed in the basis of the Wannier functions, we can build a manybody Hamiltonian [Eq. (4)] that includes the effect of crystal and ligand fields, as given by DFT, and the effect of spin-orbit interaction and on-site Coulomb repulsion at the magnetic atom. This model is solved by numerical diagonalization. An analysis of the symmetry of the spectrum and the multi-electron wave functions allows us to postulate a much simpler effective spin Hamiltonian [Eq. (9) and Eq. (11)] that accurately describe the low energy sector of the spectrum. We apply this method to the case of Fe 2+ in MgO, considering both the undistorted and distorted geometries. In the former the orbital momentum is not quenched which results in very different type of effective Hamiltonian, featuring both S and L operators. In the Jahn Teller distorted case, orbital momentum is quenched, and a spin S = 2 Hamiltonian is enough to describe the lowest energy states of Fe 2+ . The method can be implemented to study a variety of systems, including diluted magnetic semiconductors, magnetic adsorbates on insulating surfaces, and magnetic atoms migrated from the electrodes into the barrier in magnetic tunnel junctions.

VI. ACKNOWLEDGEMENTS
We acknowledge J. W. González for technical assistance with the use of Quantum Espresso and Wan-nier90. We also acknowledge J. L. Lado for fruitful discussions and assistance with technical aspects of DFT. AF acknowledges funding from the European Union's Seventh Framework Programme for research, technological development and demonstration, under the PEO-PLE programme, Marie Curie COFUND Actions, grant agreement number 600375 and CONICET. JFR acknowledges financial support by Generalitat Valenciana (ACOMP/2010/070), Prometeo.
This integral is solved numerically for l = 0, 2, and 4. From Eqs. (A7) and (A3) it is clear that all matrix elements V ijkl scales proportional to z. For convenience, instead of using z as a free parameter, we use U = V 0000 as the free parameter. In particular, z = 1.95/a 0 , with a 0 the Bohr radius, for U = 19.6 eV.