Derivation of the spin Hamiltonians for Fe in MgO

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 [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 is performed. The second step is to represent the effective DFT Hamiltonian with a basis of localized atomic orbitals, which allows us to obtain the crystal and ligand fields terms of the atomic orbitals of the relevant open shell of the magnetic atom, defining thereby a multi-orbital 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–21]. In the 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 state analysis permits the construction of an effective spin Hamiltonian for the system.

Below we describe in more detail the method and apply it to the paradigmatic case of Fe²⁺ 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²⁺ 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] showing partial quench of the orbital moment of a Co adatom on MgO, in contrast to the full quenching taking place on other surfaces like Cu₂N [26].

In addition, we are interested in Fe as a possible impurity in MgO tunnel barriers in magnetic tunnel junctions with Fe based electrodes [27, 28] 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 section 2 we study the Electronic Structure of Fe²⁺ as a substitutional impurity of Mg in MgO using DFT calculations. In section 3 we discuss the derivation of the single-particle part of the magnetic atom Hamiltonian from the DFT calculation using the wannierization approach. In section 4, we build and solve by numerical diagonalization the generalized Hubbard model and derive the effective spin Hamiltonians for two different geometries. In section 5 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.

In this section we describe our DFT calculations, for pristine MgO as well as calculations for super cells of Mg₃₁O₃₂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 generalized gradient approximation (GGA) for exchange-correlation energy [29], using plane–wave basis sets and ultrasoft pseudopotential method for Mg and O, and projector augmented-wave (PAW) [30] for Fe as implemented in Quantum ESPRESSO (QE) code [31]. 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 [35], here we do not include the U-term at this level since, as discussed in sections 3 and 4, our approach to deriving an effective spin Hamiltonian requires us 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\times 8\times 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₃₁O₃₂Fe, with lattice parameter $2a$. The supercell, including the Fe atom is shown in figure 1(a).

Zoom In Zoom Out Reset image size

MgO is an insulator with a NaCl-type crystal lattice (see figure 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 [32–34]. This discrepancy is very common and quite close to other DFT calculations for MgO (see for example [35] 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 with the d orbitals of the Fe atom; see figure 2(a). Six of these levels are below the Fermi energy and, for spin-polarized calculations, correspond to 5 levels spin $\uparrow$ and 1 $\downarrow$ resulting in a spin S = 2.

Zoom In Zoom Out Reset image size

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 neighbour 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}_{2\,g}}$ orbital triplet of the minority spin, so that the system has an 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 surrounding 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 figure 1, with coordinates $\vec{R}=({{X}_{i}},{{Y}_{i}},{{Z}_{i}})$ , i = 1,.., 6 [36]. Distortions happen to be symmetric, i.e., with $\delta {{\vec{R}}_{1}}=-\delta {{\vec{R}}_{4}}$, and we express them in terms of the normal modes of the octahedron. It turns out [36] that the computed distortion can be expressed as a linear combination of the breathing mode ${{q}_{1}}=\frac{1}{\sqrt{3}}({{X}_{1}}+{{Y}_{2}}+{{Z}_{3}})$, which clearly preserves the octahedral symmetry, and the ${{q}_{3}}=\frac{1}{\sqrt{6}}(2{{Z}_{3}}-{{X}_{1}}-{{Y}_{2}})$, 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.01{{q}_{1}}+0.03{{q}_{3}}$, where q₁ are expressed in Å and is said to be tetragonal [36]. 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, figure 2, and the spin-unpolarized figure 3, cases. The finite width ∼50 meV of the DOS peaks, much smaller than the crystal field splitting, is a finite size effect due to 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{{\mu }_{B}}$ magnetic moment (S = 2). However, the very different orbital arrangement will result in important differences in the spin Hamiltonian, as discussed below.

Zoom In Zoom Out Reset image size

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 4, 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}_{2\,g}}$ 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 $(0{{e}_{g}},6{{t}_{2\,g}})$ where all d-electrons of Fe occupy the degenerate (orbital and spin) states d_xy, d_yz and d_xz (see figures 3(a) and (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}_{2\,g}}$ levels, but the d_xy level is now split from d_xy and d_yz, as shown in figure 3(c).

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, 37–40] associated with 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 as the real $\ell =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, 37–40] 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 the others, so that it is not necessary to perform the disentanglement procedure [13, 39]. Second, the Bloch states $|{{\psi }_{{\bf k}n}}\rangle$ 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 MLWFs localized around the Fe atom that, in the neighbourhood of the atom, have the same symmetry as the real spherical harmonics with $\ell =2$, as shown in figures 4(a), (c) and (e). It is important to point out that the MLWF are not strictly identical to the atomic orbitals, because the tails of the wave functions have a different symmetry, as shown in figures 4(b), (d) and (f).

Zoom In Zoom Out Reset image size

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 in-gap bands obtained from DFT. For the purpose of this work, we are interested in the intra-cell Hamiltonian:

$$\begin{eqnarray}{{H}_{DFT}}=\left( \begin{array}{ccccccccccccccc} {{H}_{dd}} & {{V}_{dp}} \\ {{V}_{pd}} & {{H}_{pp}} \\ \end{array} \right),\end{eqnarray} \tag{ 1 }$$

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 ${{\Delta }_{CF}}$. Interestingly, diagonalization of H_dd yields, for the undistorted case, ${{\Delta }_{CF}}=0.83$ eV, much smaller than the DFT splitting 1.45 eV, that is only recovered if the whole H_DFT 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:

$$\begin{eqnarray}&&{{\mathcal{H}}_{d{{d}^{\prime }}}}={{H}_{d{{d}^{\prime }}}}+\mathop{\sum }\limits_{p}\frac{{{V}_{dp}}{{V}_{p{{d}^{\prime }}}}}{{{E}_{d}}-{{E}_{p}}}.\end{eqnarray} \tag{ 2 }$$

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_DFT 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:

$$\begin{eqnarray}&&{{H}_{CF}}=a(\ell _{X}^{4}+\ell _{Y}^{4}+\ell _{z}^{4})+{{d}_{2}}l_{z}^{2}+{{d}_{4}}l_{z}^{4},\end{eqnarray} \tag{ 3 }$$

where l_a are the $\ell =2$ angular momentum matrices, and d₂, d₄ 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 [41] with good qualitative results.

In order to fit a, d₂ and d₄ we employ the analytical expressions of the eigenvalues of H_CF: $18a+{{d}_{2}}+{{d}_{4}},18a+{{d}_{2}}+{{d}_{4}},18a+4{{d}_{2}}+16{{d}_{4}},24\;a,24a+4{{d}_{2}}+16{{d}_{4}}$. For the undistorted case, the in-gap d-levels obtained from diagonalization of equation (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 figure 3(a), while the e_g doublet is also split. The fitting yields a = 0.250 eV, d₂ = 0.461 eV and ${{d}_{4}}=-0.1$ eV. The difference between the fitted and computed energy levels are always smaller than 2 meV.

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 section 2. By so doing, we arrive at effective spin models in agreement with the literature [1, 24, 25].

4.1. 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:

$$\begin{eqnarray}&&H={{H}_{{\rm Coul}}}+{{H}_{{\rm CF}}}+{{H}_{{\rm SO}}}+{{H}_{{\rm Zeem}}},\end{eqnarray} \tag{ 4 }$$

The Coulomb term reads:

$$\begin{eqnarray}&&{{H}_{{\rm Coul}}}=\frac{1}{2}\mathop{\sum }\limits_{\frac{m,m^{\prime} }{n,n^{\prime} }}{{V}_{mnm^{\prime} n^{\prime} }}\mathop{\sum }\limits_{\sigma \sigma ^{\prime} }d_{m\sigma }^{\dagger }d_{n\sigma ^{\prime} }^{\dagger }{{d}_{n^{\prime} \sigma ^{\prime} }}{{d}_{m^{\prime} \sigma }},\end{eqnarray} \tag{ 5 }$$

where $d_{m\sigma }^{\dagger }$ (${{d}_{m\sigma }}$) denotes the creation (annihilation) operator of an electron with spin σ in the $\ell =2,{{\ell }_{z}}=m$ state of the magnetic atom, denoted by ${{\phi }_{m}}(\vec{r})$, assumed to be equal to the product of a radial hydrogenic function (with effective charge Z and a effective Bohr radius ${{a}_{\mu }}$) and a spherical harmonic. Thus, we are only considering d⁶ configurations and leaving out $p{{d}^{7}}$ configurations.

It turns out that the Coulomb integrals ${{V}_{mnm^{\prime} n^{\prime} }}$ scale linearly with the value of ${{V}_{0000}}\equiv U$. Explicit expressions for the on-site Coulomb integrals are given in the appendix. Specifically, U could be computed using equation (A.7). 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 ^3P2 excited state and the ^5D ground state of the free ion, measured [43] to be 2.41 eV. Although this un-screened value of the Coulomb interaction is certainly reduced in the solid, the only role of U in the low energy spectra of our single electronic configuration, Fe-d⁶, is to avoid the mixing of the ⁵D and ³P₂ states.

The second term in equation (4) corresponds to the crystal and ligand fields Hamiltonian discussed in the previous section:

$$\begin{eqnarray}&&{{H}_{CF}}=\frac{1}{2}\mathop{\sum }\limits_{m,m^{\prime} ,\sigma }\langle m|{{H}_{{\rm CF}}}|m^{\prime} \rangle d_{m\sigma }^{\dagger }{{d}_{m^{\prime} \sigma }},\end{eqnarray} \tag{ 6 }$$

with $\langle m|{{H}_{{\rm CF}}}|m^{\prime} \rangle$ derived from DFT using the procedure described above and a very good approximation, is given by equation (3).

The last term in the Hamiltonian describes spin–orbit coupling:

$$\begin{eqnarray}&&{{H}_{{\rm SO}}}=\zeta \mathop{\sum }\limits_{mm^{\prime} ,\sigma \sigma ^{\prime} }\langle m\sigma |\vec{\ell }\cdot \vec{S}|m^{\prime} \sigma ^{\prime} \rangle d_{m\sigma }^{\dagger }{{d}_{m^{\prime} \sigma ^{\prime} }},\end{eqnarray} \tag{ 7 }$$

where ζ is the single particle spin–orbit coupling of the d-electrons. It is also very frequently expressed as $\lambda \vec{L}\cdot \vec{S}$ with $\vec{L}$ the total angular momentum. For the case of Fe²⁺, both parameters ζ and λ are related by $\lambda =-\zeta /2S$ [1], with $\zeta =50.1$ meV and S = 2.

The last term in equation (4) corresponds to the Zeeman Hamiltonian:

$$\begin{eqnarray}&&{{H}_{{\rm Zeem}}}={{\mu }_{B}}\vec{B}\cdot \mathop{\sum }\limits_{mm^{\prime} ,\sigma \sigma ^{\prime} }\langle m,\sigma |\left( \vec{l}+g\vec{S} \right)|m^{\prime} \sigma ^{\prime} \rangle d_{m\sigma }^{\dagger }{{d}_{m^{\prime} \sigma ^{\prime} }},\end{eqnarray} \tag{ 8 }$$

where g = 2. So, if we assume that the CF term is given by equation (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⁶ 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 ⁵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 ⁵D with the ³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₂, d₄), and spin–orbit coupling (ζ). Thus, in a first step the problem is solved considering only ${{H}_{{\rm Coul}}}$. In this case the Hamiltonian commutes with S² and L², the square of total spin and orbital angular momentum.

4.2. Undistorted Fe/MgO

We discuss first the case of Fe²⁺ 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 figure 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 equation (3), the ⁵D levels of iron splits into two, an orbital ${{\Gamma }_{5}}$ triplet ground state, with total degeneracy 15, and an orbital doublet excited state, 10 times degenerate (see figure 5(a)).

Zoom In Zoom Out Reset image size

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 $\vec{\ell }$ operator on the subspace of the ${{t}_{2g}}$ orbitals is isomorphic to the $\ell =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 figure 5(b)). This pattern can be rationalized in terms of the following effective Hamiltonian where the total spin is coupled to the pseudo-angular momentum $\mathcal{L}=-1$ [1]:

$$\begin{eqnarray}&&{{\mathcal{H}}_{{\rm eff}1}}=\lambda \vec{\mathcal{L}}\cdot \vec{S}+\delta {{\mathcal{H}}_{{\rm eff}1}},\end{eqnarray} \tag{ 9 }$$

where $\lambda =-\zeta /2S$. The first term naturally leads to a spectrum with multiplets associated with ${{\tilde{J}}^{2}}={{(\vec{\mathcal{L}}+\vec{S})}^{2}}$: $\tilde{J}=1$ (ground), $\tilde{J}=2$ (first excited) and $\tilde{J}=3$ (third excited), with degeneracies $2\tilde{J}+1=3,5$ and 7 respectively. The result of the calculation of the expectation values $\langle \psi |{{S}_{z}}|\psi \rangle$ for ψ in the ground state triplet with $\tilde{J}=1$, backs up the idea that the S = 2 spin is coupled to an pseudo-angular momentum with $\mathcal{L}=-1$. In figures 5(d) 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 the ${{\tilde{J}}_{z}}=\pm 1$ and ${{\tilde{J}}_{z}}=0$ values are recovered by subtracting $\langle \psi |{{S}_{z}}|\psi \rangle$ and $\langle \psi |{{L}_{z}}|\psi \rangle$, in contrast with the common case of a total angular momentum.

The CI calculation for the 15 lowest energy states for Fe²⁺ in the undistorted environment has some fine structure not captured by the first term in equation (9). In particular, the multiplets with $\tilde{J}=1$ and 2 have some fine structure (see figure 5(b)), which can be accounted for with the second term in the effective Hamiltonian

$$\begin{eqnarray}&&\delta {{\mathcal{H}}_{{\rm eff}1}}=\tilde{a}(\tilde{J}_{x}^{4}+\tilde{J}_{y}^{4}+\tilde{J}_{z}^{4}).\end{eqnarray} \tag{ 10 }$$

This operator does not break the triple degeneracy of the ground state, but breaks the $\tilde{J}=2$ into a triplet and a doublet (being isomorphic to the problem of the octahedral crystal field splitting of the $\ell =2$ orbitals), and the $\tilde{J}=3$ into a singlet and two triplets.

In summary, in the undistorted case, our calculation portrays Fe²⁺ as a system with S = 2 and pseudo-orbital momentum $\mathcal{L}=-1$ [1]. Spin–orbit coupling leads to a ground state triplet with $\tilde{J}=1$. The energy splitting to the first excited state, with $\tilde{J}=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, which we discuss next, changes this situation.

4.3. Jahn-Teller distorted Fe/MgO

We now discuss the effect of the tetragonal distortion on the multiplet structure of Fe²⁺ in MgO. As discussed in section 2, this distortion introduces the uniaxial terms ${{d}_{2}}l_{z}^{2}+{{d}_{4}}l_{z}^{4}$ in equation (3). The effect of the uniaxial terms on the many-body 15-fold degenerate ground state of the Hamiltonian with $\zeta =0$ and a = 0.250 meV is shown in figure 6(a). 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 ${{\Gamma }_{5}}$ orbital triplet. When spin–orbit coupling is added (figure 6(b)) the $2S+1=5$ degeneracy is broken into a singlet, a doublet, and a split doublet (see figure 7(a)). Finally, the Zeeman splitting breaks the remaining degeneracies, as observed in figure 6(c).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Interestingly, the five low energy states, corresponding to ${{\tilde{l}}_{z}}=0$, can be described by an effective $\mathcal{S}=2$ Hamiltonian of the form

$$\begin{eqnarray}&&{{\mathcal{H}}_{{\rm eff}2}}=D\mathcal{S}_{z}^{2}+\frac{a}{6}\left[ \mathcal{S}_{x}^{4}+\mathcal{S}_{y}^{4}+\mathcal{S}_{z}^{4} \right]+{{g}^{*}}{{\mu }_{B}}\vec{B}\cdot \vec{\mathcal{S}}.\end{eqnarray} \tag{ 11 }$$

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 figure 7(a). The parameters of the effective Hamiltonian are obtained by fitting to the multiplet calculation. We obtain D = 0.734 meV, a = 0.130 meV and ${{g}^{*}}=2.03$. The expectation values of S_z and L_z, computed with the eigenstates of the full CI Hamiltonian, are shown in figures 7(b), (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 $D\mathcal{S}_{z}^{2}$ favouring the minimum spin projection as ground state, ${{\mathcal{S}}_{z}}=0$. The first excited doublet, split by B_z, has ${{S}_{z}}=\pm 1$. The $\mathcal{S}_{x}^{4}+\mathcal{S}_{y}^{4}$ term couples the otherwise degenerate doublet ${{S}_{z}}=\pm 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 $\langle {{S}_{z}}\rangle$ 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₂ in plane symmetry [3], in our case there is no quantum spin tunneling splitting within the ${{S}_{z}}=\pm 1$ doublet, which remains degenerate.

The results of the previous sections illustrate how, for the cases of Fe²⁺ in MgO with and without Jahn–Teller distortion, we have been able to derive effective spin Hamiltonians (equations (9) and (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 [27, 28, 44, 45], and present our conclusions.

5.1. Improvements for the method

5.2. Influence of Fe impurities in the barrier MgO of a magnetic tunnel junction

5.3. Summary

We acknowledge J W González for technical assistance with the use of Quantum ESPRESSO and Wannier90. 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 PEOPLE programme, Marie Curie COFUND Actions, grant agreement number 600375 and CONICET. JFR acknowledges financial support by Generalitat Valenciana (ACOMP/2010/070), Prometeo, and MEC-Spain (FIS2013-47328-C2-2-P).

The Coulomb parameters V_ijkl are calculated assuming nd hydrogen-like wavefunctions

$$\begin{eqnarray}&&{{V}_{ijkl}}=\frac{{{{\rm e}}^{2}}}{4\pi {{\epsilon }_{0}}}\int {\rm d}{{\vec{r}}_{1}}{\rm d}{{\vec{r}}_{2}}\frac{\phi _{i}^{*}({{{\vec{r}}}_{1}})\phi _{j}^{*}({{{\vec{r}}}_{2}}){{\phi }_{k}}({{{\vec{r}}}_{2}}){{\phi }_{l}}({{{\vec{r}}}_{1}})}{|{{{\vec{r}}}_{1}}-{{{\vec{r}}}_{2}}|},\end{eqnarray} \tag{ A.1 }$$

where

$$\begin{eqnarray}&&{{\phi }_{i}}(\vec{r})={{R}_{n,2}}(r)Y_{2}^{i}(\Omega ),\end{eqnarray} \tag{ A.2 }$$

with $Y_{l}^{m}(\Omega )$ the spherical harmonic and ${{R}_{n,2}}(r)$ the hydrogen wavefunction corresponding to quantum numbers n and l = 2 for an effective nuclear charge Z and effective radius ${{a}_{\mu }}$

$$\begin{eqnarray}&&{{R}_{n2}}(r)=\frac{4}{27\sqrt{10}}{{\left( z/3 \right)}^{3/2}}{{\left( zr \right)}^{2}}{{{\rm e}}^{-zr/3}}\end{eqnarray} \tag{ A.3 }$$

where $z=Z/{{a}_{\mu }}$. Using the spherical harmonic expansion

$$\begin{eqnarray}&&\frac{1}{|{{{\vec{r}}}_{1}}-{{{\vec{r}}}_{2}}|}=4\pi \mathop{\sum }\limits_{l=0}^{\infty }\mathop{\sum }\limits_{m=-l}^{+l}\frac{1}{2\;l+1}\frac{r_{\lt }^{l}}{r_{\gt }^{l+1}}Y_{l}^{m*}({{\Omega }_{1}})Y_{l}^{m}({{\Omega }_{2}}),\end{eqnarray} \tag{ A.4 }$$

where ${{r}_{\lt }}={\rm min} ({{r}_{1}},{{r}_{2}})$ and ${{r}_{\gt }}={\rm max} ({{r}_{1}},{{r}_{2}})$ , one can divide the integral in equation (A.2) into an angular part and a radial part, writing then

$$\begin{eqnarray}&&{{V}_{ijkl}}=\mathop{\sum }\limits_{\ell }{{U}_{\ell }}\chi _{ijkl}^{\ell ,m},\end{eqnarray} \tag{ A.5 }$$

where ${{U}_{\ell }}$ and $\chi _{ijkl}^{\ell ,m}$ contains the radial and the angular information respectively. The angular integrations over the solid angles ${{\Omega }_{1}}$ and ${{\Omega }_{2}}$ factorizes, $\chi _{ijkl}^{\ell ,m}={{(-1)}^{i+j+m}}\Phi _{2,\ell ,2}^{-i,m,l}\Phi _{2,\ell ,2}^{-j,-m,k}$, where each part can be written in terms of the Wigner 3-j symbols [64]

$$\begin{eqnarray}&&\begin{array}{l} \Phi _{{{l}_{1}},{{l}_{2}},{{l}_{3}}}^{{{m}_{1}},{{m}_{2}},{{m}_{3}}} \\ \quad =\int {\rm d}\Omega Y_{{{l}_{1}}}^{{{m}_{1}}}(\Omega )Y_{{{l}_{2}}}^{{{m}_{2}}}(\Omega )Y_{{{l}_{3}}}^{{{m}_{3}}}(\Omega ) \\ \quad =\sqrt{\frac{(2{{l}_{1}}+1)(2{{l}_{2}}+1)(2{{l}_{3}}+1)}{4\pi }}\left( \begin{array}{ccccccccccccccc} {{l}_{1}} & {{l}_{2}} & {{l}_{3}} \\ 0 & 0 & 0 \\ \end{array} \right) \\ \quad \times \left( \begin{array}{ccccccccccccccc} {{l}_{1}} & {{l}_{2}} & {{l}_{3}} \\ {{m}_{1}} & {{m}_{2}} & {{m}_{3}} \\ \end{array} \right). \\ \end{array}\end{eqnarray} \tag{ A.6 }$$

The radial part is given by

$$\begin{eqnarray}&&{{U}_{\ell }}\equiv \frac{{{{\rm e}}^{2}}}{{{\epsilon }_{0}}(2\ell +1)}\int {\rm d}r{\rm d}r^{\prime} {{r}^{2}}r{{^{\prime} }^{2}}\frac{r_{\lt }^{\ell }}{r_{\gt }^{\ell +1}}R_{n2}^{2}(r)R_{n2}^{2}(r^{\prime} ).\end{eqnarray} \tag{ A.7 }$$

This integral is solved numerically for l = 0, 2, and 4. From equations (A.7) and (A.3) it is clear that all matrix elements V_ijkl scale 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₀ the Bohr radius, for U = 19.6 eV.