Anatomy of Magnetic Anisotropy and Voltage-Controlled Magnetic Anisotropy in Metal Oxide Heterostructure from First Principles

Voltage control of magnetic anisotropy (VCMA) is one of the promising approaches for magnetoelectric control of magnetic tunnel junction (MTJ). Here, we systematically calculated the magnetic anisotropy (MA) and the VCMA energies in the well-known MTJ structure consisting of Fe/MgO interface with Cr buffer layer. In this calculation, we investigated an alloying between Fe and Cr and a strain effect. We used a spin density functional approach which includes both contributions from magnetocrystalline anisotropy energy (MCAE) originating from spin–orbit coupling and shape magnetic anisotropy energy from spin dipole–dipole interaction. In the present approach, the MCAE part, in addition to a common scheme of total energy, was evaluated using a grand canonical force theorem scheme. In the latter scheme, atom-resolved and k-resolved analyses for MA and VCMA can be performed. At first, we found that, as the alloying is introduced, the perpendicular MCAE increases by a factor of two. Next, as the strain is introduced, we found that the MCAE increases with increasing compressive strain with the maximum value of 2.2 mJ/m2. For the VCMA coefficient, as the compressive strain increases, the sign becomes negative and the absolute value becomes enhanced to the number of 170 fJ/Vm. By using the atom-resolved and k-resolved analyses, we clarified that these enhancements of MCAE and VCMA mainly originates from the Fe interface with MgO (Fe1) and are located at certain lines in the two dimensional Brillouin zone. The findings on MCAE and VCMA are fully explained by the spin-orbit couplings between the certain d-orbital states in the second-order perturbation theory.


Introduction
Magnetic random access memory (MRAM) is one type of random access memory, and has some advantages such as relatively fast read/write accessibility and unlimited endurance, compared to the other memories such as static and dynamic access memories and the storages such as NAND (not and)-type flash memory and hard disk drive. In principle, the MRAM uses magnetoresistance effect for reading process. The commercialized MRAM has employed a tunnel magnetoresistance As reviewed in the previous paragraph, the phenomena of large VCMA observed in the Fe/MgO interface with Cr layer are interesting in both magnetism and spintronics application. The system may keep the atomic-scale structure of the interface. In the application, it is important to uncover the mechanism of the enhancement of the VCMA coefficient preserving the high perpendicular MA. However, there is no systematic analysis from the first-principles calculation to explain the mechanism. It may be useful to provide results of the magnetic MA and VCMA from a systematic calculation on the effects of alloying and strain.
In this work, we performed a theoretical study of magnetic anisotropy and voltage-controlled magnetic anisotropy based on spin density functional theory for the Fe/MgO interface with Cr buffer layer. In the calculation of MA energy, we included both magnetocrystalline anisotropy energy (MCAE) and shape magnetic anisotropy energy (SMAE). This allows us to make a direct comparison with available experimental measurements. To the interface, we introduced two kinds of perturbation, namely alloying and strain. These perturbations can effectively modify the electronic structure and eventually enhance the MA energy and VCMA coefficient. We systematically elucidate the possible origin of such enhancement by connecting the electronic structure near the Fermi level with the atom-resolved or k-resolved contributions.

Computational Models and Methods
We investigated electronic structures in the metal oxide heterostructure of Fe/MgO interface with Cr buffer layer by means of the approach based on spin density functional theory (SDFT) [49]. We considered two kinds of heterostructure without and with alloying, called structure-I and structure-II. As shown in Figure 1 In the SDFT approach, the total energy of system was evaluated by solving the Kohn-Sham (KS) equation self-consistently with the densities and wavefunctions for electrons. The density for electron is presented by ρ(r) = n(r)σ 0 + ∑ x,y,z α m α (r)σ α , where n(r) and m α (r) are electron and spin densities, respectively, and σ α and σ 0 are Pauli's spin matrix and unit matrix, respectively. Our SDFT calculation employs fully relativistic and scalar relativistic ultrasoft pseudopotentials (USPPs) [50] and a plane wave basis [51] in the generalized gradient approximation (GGA) for the exchange-correlation energy [52]. With the fully relativistic USPP, one can take into account the spin-orbit coupling (SOC). The energy cut-offs of plane wave basis for the wave function and density are 30 Ryd and 300 Ryd, respectively [53].
In this work, the MA energy (MAE) is defined by the relative energy of [100] magnetization direction (in-plane) with respect to the energy of [001] magnetization direction (perpendicular to the plane). It was evaluated by considering MCAE and SMAE (MAE = MCAE + SMAE). In the evaluation of MCAE, we used two schemes; direct total energy (TE) difference and grand canonical force theorem (GCFT). In the GCFT scheme, the calculation of the KS equation requires two steps. At the first step, a density ρ(r) was generated in a self-consistent field calculation without including the SOC (using scalar relativistic USPP). In the second step, we fixed the density ρ(r) and performed two calculations including the SOC for the magnetization directions along the in-plane and perpendicular to the plane, respectively. As a result, a set of eigenvalues for nth band and wave vector k (ε nk ) was obtained and used to evaluate the MCAE. This GCFT scheme enables us to obtain the anisotropy energy of atom-resolved or k-resolved contributions. The details of the GCFT scheme can be found in the previous report [54]. The SMAE was calculated using discrete approach (DA) and spin density approach (SDA). In the SDA approach, we investigated the components from atomic multipole spin-density, such as the quadrupole component [55].
For the structure-II, we calculated the effect of strain (η) on the MA energy and VCMA coefficient. We considered the variation of strain in the range of η = −8.0% to 3.8% for the ratio with respect to Fe lattice constant (a Fe = 2.87Å) (η = (a − a Fe )/a Fe ). For each in-plane lattice constant, we optimized the atomic structure using the scalar-relativistic level computation with a 24 × 24 × 1 k-mesh and fixing the atomic coordinates of O(3). We used a 32 × 32 × 1 k-mesh for the MAE calculations. In order to apply an EF, we used the effective screening medium (ESM) method [56], and to estimate the EF inside the MgO layer, the dielectric constant ε r (9.8 for MgO) was used [57].

Alloying Effect of Magnetic Anisotropy
The MAEs for structure-I and structure-II are summarized in Table 1, including the MCAE and SMAE contributions. In both structures, the MCAEs are positive and contribute a perpendicular MA. Interestingly, in structure-II where the alloying of Fe and Cr are introduced, the MCAE increases by a factor of two. The SMAEs are negative and contribute an in-plane MA. In structure-II, the absolute of SMAE becomes smaller due to the alloying effect. The situation that an effective thickness of Fe layer becomes smaller is in accordance with an existence of magnetic dead layer found in the measurement [39]. The absolute of SDA is slightly lower than the absolute of DA. This reduction can be explained from a positive sign of quadrupole component in the atomic spin density. The spin density m α (r) is projected to a set of radial spin densities, as follows: [55] where r I = |r − R I | is the distance from the Ith atomic site R I ,r I specifies the angular variable of r I = r − R I , and Y lm (r I ) is the spherical harmonic function. Figure 2 shows the radial spin densities µ α Ilm (r I ) for the spherical component (l = 0, m = 0) and the quadrupole component (l = 2, m = 0) on each Fe atoms. All the other components are negligible. A negative sign indicates an oblate ellipsoid of spin density distribution while a positive sign indicates a prolate ellipsoid. The oblate (prolate) shape in atomic spin density leads to an increase (decrease) of the in-plane shape magnetic anisotropy [55]. In both Figures 2a,b, the large positive quadrupole component at Fe 1 is observed, giving a slight decrease of SMAE. After combining the MCAE and SMAE, the MAE in structure-I becomes negative, while in structure-II it remains positive. The MAE of structure-II has the same sign as the experimental measurement with similar structure [39]. As a result, the enhancement of MCAE in structure-II plays a crucial role for realizing the perpendicular anisotropy. The reduction in SMAE may also support such anisotropy. Table 1. The magnetic anisotropy energy (MAE) from magnetocrystalline anisotropy energy (MCAE) and shape magnetic anisotropy energy (SMAE) for structure-I and structure-II and the comparison with the available experimental data. The MCAE is calculated based on TE and grand canonical force theorem (GCFT). The SMAE is evaluated from discrete approach (DA) and spin density approach (SDA). -

Structure
In order to elucidate the possible origin of the enhancement in MCAE, we calculated the atom-resolved and k-resolved MCAEs using the GCFT scheme. The total MCAE agrees well with that obtained by the TE scheme, as shown in Table 1. The atom-resolved MCAEs for structure-I and structure-II are shown in Figure 3a,b, respectively. The positive and negative values indicate the contributions to perpendicular and in-plane MA, respectively. In both structures, Fe 1 and Fe 2 provide the main contribution to their perpendicular MA. This agrees with the previous theoretical investigation [36,58,59]. Furthermore, we can clearly see that the main origin of the enhanced MCAE in structure-II comes from Fe 1 (see Figure 3). This drastic enhancement on the interface atom may be a useful element for a design of magnetic anisotropy. In addition, in structure-I, the Fe atoms inside the layer (Fe 3 and Fe 4 ) contribute to in-plane MA. As a result, the MCAE in structure-I becomes suppressed. We also calculated the orbital-resolved MCAE. Its result shows that the main contribution to the MCAE originates from d-orbital states, as shown in Tables 2 and 3.  Next, we show the atomic k-resolved MCAE for Fe 1 in structure-I and structure-II by combining the band dispersion curves at the Fermi level in Figure 4. The atomic k-resolved MCAE for the other Fe atoms are presented in Appendix A. On the band dispersion curves, we show only d-orbital states for the minority-spin states, since the majority-spin states are fully occupied and located at the energy lower than −1 eV on Fe 1 component. In the contour maps, the energy ranges from blue to red indicate the negative to positive contributions to the MCAE. From both of Figure 4b,d, we notice that the origin of the enhancement mentioned above originates from the enhancement of positive contribution and the reduction of negative contribution located at certain regions in the two-dimensional Brillouin zone (2DBZ).
In order to understand the mechanism/origin of positive and negative contributions in the k-resolved MCAE, we use a second-order perturbation theory. The formula can be written as [60] where ϕ 0↓ ok and ϕ 0↓ uk represent the occupied and unoccupied minority-spin state and ε 0↓ uk and ε 0↓ ok indicate the corresponding energy, respectively. L z I , L x I are the atomic angular momentum operators and ξ I is the SOC constant. Here, we consider only the coupling between d-orbital minority-spin states, since the majority-spin 3d-orbital states are fully occupied and are located far from the Fermi level.  From Equation (2), the MCAE is strongly related to the states near the Fermi level. the MCAE becomes large when the difference between the energies of occupied and unoccupied states is small. Another important property from Equation (2) is the angular momentum matrix element between d-orbital states. For examples, the couplings between (d xz , d yz ) and (d x 2 −y 2 , d xy ) through L z I provide the positive energies of e and 4e (e > 0), respectively, to the MCAE, when assuming the same pair of atomic radial wavefunctions for the given pair of ε 0↓ uk and ε 0↓ ok . In these couplings, the latter is larger by 4 than the former. Similary, the couplings between (d 3z 2 −r 2 , d yz ), (d xy , d xz ), and (d x 2 −y 2 , d yz ) through L x I provide the negative enegies of −3e, −e, and −e, respectively [60].
Based on the discussion above, for the structure-I, we conclude that there are three main positive contributions labeled by 1, 2, and 3 in the contour map (see Figure 4b), while one negative contribution labeled by 4. By connecting with the detail of electronic structures in structure-I, as shown in Figure 4a, we notice that the label 1 (the positive contribution) is addressed to the SOCs of o(u), xz|L z I |u(o), yz , the label 2 is addressed to the SOCs of o, yz|L z I |u, xz , and the label 3 is addressed to the SOCs of o, xy|L z I |u, x 2 − y 2 . On the other hand, the negative contribution, namely, the label 4 is addressed to the SOCs of o(u), 3z 2 − r 2 |L x I |u(o), yz . By introducing the alloying of Fe and Cr (for structure-I), there are some significant changes in band dispersion curves, as shown in Figure 4c. In Figure 4c, some of d-orbital states become splitted and shifted to the lower energy. For examples, it is clearly observed at Γ and X points. As a result, in structure-II, the positive contribution becomes enhanced and slightly shifted from Γ point at region 1 and 2 in the contour map (see Figure 4c). We checked carefully the electronic structure and identified that the label 1 belongs to the SOCs of o, xy|L z I |u, x 2 − y 2 and label 2 belongs to the SOCs of o, yz|L z I |u, xz and o, xy|L z I |u, x 2 − y 2 . In addition, for the reduction of negative contribution the label 4 is caused by the shift of the states d 3z 2 −r 2 and d yz . As a result, the denominator in Equation (2) becomes larger, and hence the negative contribution to the MCAE becomes reduced.

Strain Effect
We performed the calculations of VCMA for structure-II, as the strain dependence in the range of −8% to 3.8%. The SMAE was found not to contribute to the VCMA because the change of magnetization by EF does not appear. Thus, we considered the contribution of MCAE only. The results for the MCAE is shown in Figure 5. The MCAEs from the TE and GCFT schemes are in good agreement. The results show that by increasing tensile strain (increase in lattice constant), the MCAE decreases, while increasing compressive strain (decrease in lattice constant), the MCAE increases with the maximum value of 2.2 mJ/m 2 . This value is increased by a factor of two, compared to the zero strain.
To clarify the mechanism of the enhancement of MCAE by compressive strain, we calculated the k-resolved MCAE as the strain dependence, as shown in Figure 6. The atomic k-resolved MCAE is presented in Appendix B for one strain value as a reference. From the sequence of k-resolved MCAE, we can see that there are three main changes of MCAE. These are at Γ point, Y point, and along the X-Y line. Furthermore, the atomic k-resolved MCAE shows that the enhancement at X-Y line originating from Fe 1 , while the enhancement at Γ point and Y point originating from the Fe atoms inside the layer (Fe 2 , Fe 3 , and Fe 4 ) (see Figure A2 in Appendix B).
In order to explain such enhancement, especially at X-Y line, we show the band dispersion curves as the strain dependence with the orbital angular components of Fe 1 in Figure 7. From Figure 7a-c, we found that the most significant orbital reconstruction appears at X-Y line. Initially, at η Fe = 3.8% [ Figure 7a], the states of d xy , d x 2 −y 2 , d yz , and d xz are unoccupied at 0.2 eV. By increasing the compressive strain, the states of d x 2 −y 2 , d yz , and d xz are gradually shifted down [ Figure 7b], and finally become occupied at −0.1 eV [ Figure 7c]. Based on Equation (2), due to this orbital reconstruction, the SOCs of u, xy|L z I |o, x 2 − y 2 and o(u), yz|L z I |u(o), xz increase. As a result, the MCAE is enhanced. Next, we show the result of the VCMA coefficient (γ) for structure-II as the strain dependence. The γ was calculated as follows. At first, the magnitude of EF is estimated from the slope of the electrostatic potential at the front of ESM, as seen in Figure 8a. The effective EF in the MgO layer (E) was determined by dividing the above EF with the dielectric constant (ε r = 9.8). As an example, Figure 8b shows the MCAE as a function of EF for η Fe = −5.9%. The VCMA coefficient was eventually assigned from the slope under the positive EF. The summary of γ is shown as the strain dependence in Figure 9. From this result, we can see that at zero strain the sign of γ is positive. By increasing the compressive strain larger than 2%, the sign becomes negative and the absolute value is enhanced to the number of 170 fJ/Vm. To analyze γ, we calculated the atom-resolved MCAE at zero EF and a finite EF. Figure 10a,b are the results for η Fe = 0.3% and η Fe = −5.9%, respectively. It is clear that the main change of MCAE induced by EF appears at the interface Fe (Fe 1 ). This is due to the screening effect of the metal layer and this result agrees with the previous theoretical reports [21,23,61].    Figure 11, we can see that the EF induces MCAE only in certain regions in 2DBZ. In the case of η Fe = 0.3%, positive EF increases the MCAE, while negative EF decreases the MCAE. The reduction and increase in MCAE are located at (4/8) M X-line. The origin of such change can be explained as follows. In zero EF, there is a large negative contribution as mentioned in the previous discussion (see Figure 6) due to the SOCs of o(u), 3z 2 − r 2 |L x I |u(o), yz . By applying a positive EF, corresponding to the reduction of the number of electrons, the Fermi level shifts down. As a result, the SOCs of o(u), 3z 2 − r 2 |L x I |u(o), yz (negative contribution to MCAE) decrease and hence the MCAE increases. For the negative EF is vice versa.
On the other hand, in the case of η Fe = −5.9% [see Figure 11d-f], positive EF reduces the MCAE, while negative EF increases the MCAE. Those behaviors located at XY-line were observed and a large γ was obtained. The explanations of such behavior and large γ are as follows. As shown in Figure 7b, at zero EF, d x 2 −y 2 and part of d xz and d yz are occupied and form states on the Fermi level. As a positive EF is applied, those states are shifted to higher energy and become unoccupied. This modification may reduce the SOCs of u, xy|L z I |o, x 2 − y 2 and u(o), xz|L z I |o(u), yz . Due to this reduction, the MCAE may reduce significantly. For the negative EF is vice versa.
We confirmed the change of occupation number from the partial density of states (PDOS) under EF. Figure 12 shows the PDOSs of minority-spin states on Fe 1 for η Fe = −5.9% at zero field (shade area), positive EF (solid red line), and negative EF (solid blue line). Under the positive EF, there is a clear change of states near the Fermi level. However, under the negative EF, only a small change is observed. This may be the reason for a small change of MCAE under negative EF. Similar to this work, the enhancements in MCAE and γ were also observed in metal/oxide interfaces as strain effect [45,46].

Conclusions
We systematically calculated the MA energy and VCMA coefficient in the Fe/MgO interface with the Cr buffer layer. We introduced two kinds of perturbation, namely, the alloying and strain effects. In the calculation of MA energy, we included both contributions from the MCAE originating from SOC and the SMAE originating from spin dipole-dipole interaction. Using the GCFT scheme, we evaluated the atom-and k-resolved MCAEs. We found the MCAE values from the GCFT scheme are well reproduced by those from the TE scheme. Both structure-I and structure-II showed the perpendicular MCAE. As the alloying is introduced, the MCAE increased by a factor of two. After combining with the SMAE, the MA energy in structure-I becomes negative while in structure-II it remains positive. The latter value has the same sign (perpendicular anisotropy) as the experimental measurement with the same structure. From the atom-resolved data, the enhancement of MCAE mainly originates from the Fe interface with MgO (Fe 1 ). Moreover, the atomic k-resolved MCAE showed that the enhancement is a result of enhancements and reductions at the positive and negative MCAE contributions located at Γ-Y and M-X lines in 2DBZ, respectively. As the strain is introduced in structure-II, the MCAE increased with increasing the compressive strain and becomes the maximum value of 2.2 mJ/m 2 . For the VCMA coefficients, as the compressive strain increases, the sign changes to negative and the absolute value becomes enhanced to 170 fJ/Vm. The atomic k-resolved VCMA showed that the enhancement appears at the certain region/specific regions in 2DBZ which correspond to a large amount of certain d-orbital states near the Fermi level. We believe that our present results can give shed light on the design of material for Me-RAM.

Appendix A. Atomic k-Resolved MCAE of Structure-I
In the GCFT scheme, we calculated the atomic k-resolved MCAE for each Fe atom. The results for structure-I are shown in Figure A1. In the figure, the bare data of no symmetrization and asymmetrized component are also included as well as the symmetrized component discussed in Section 3. The set of eigenvalues determined by the GCFT scheme contains the effect of odd parity in k-space. After summing over the k-space (2DBZ), eventually, such effect does not contribute to MCAE. However, as shown at Fe 1 in Figure A1, the odd parity part (asymmetrized component) is dominant over a wide region of 2DBZ in the bare data. Such contribution may be discussed with the Rashba effect at the atomic layer of Fe 1 and the effective electric field imposed on the layer.

Appendix B. Atom k-Resolved MCAE of Structure-II
The results of atomic k-resolved MCAE for the structure-II of η Fe = −5.9% are shown in Figure A2.