## Abstract

Large perpendicular magnetic anisotropy (PMA) in transition metal thin films provides a pathway for enabling the intriguing physics of nanomagnetism and developing broad spintronics applications. After decades of searches for promising materials, the energy scale of PMA of transition metal thin films, unfortunately, remains only about 1 meV. This limitation has become a major bottleneck in the development of ultradense storage and memory devices. We discovered unprecedented PMA in Fe thin-film growth on the *x*^{2} − *y*^{2} and *xy* orbitals and its lift by the spin-orbit coupling play a dominant role. As a consequence, PMA in Fe/III-V nitride thin films is dominated by first-order perturbation of the spin-orbit coupling, instead of second-order in conventional transition metal/oxide thin films. This game-changing scenario would also open a new field of magnetism on transition metal/nitride interfaces.

## INTRODUCTION

Magnetic anisotropy is a relativistic effect originating from spin-orbit coupling (SOC). Perpendicular magnetic anisotropy (PMA) in magnetic thin films has led to rich physics and become a key driving force in the development of magnetic random-access memory (MRAM) devices (*1*–*3*). Establishment of PMA in nanostructures and nanopatterned magnetic multilayers paves a new avenue toward nanomagnetism, in which fascinating physics such as spin Hall switching and skyrmions is blooming (*4*–*7*). Because the strength of the SOC is a quartic function of the atomic number, it is not surprising to have large magnetic anisotropy in heavy metals such as rare earth materials, which are commonly used for permanent magnets (*8*, *9*). It is, however, a challenge to enable large anisotropy, especially PMA, in commonly used 3*d* transition metals such as Fe thin films.

Actually, the strength of PMA is determined by the energy correction from the SOC, which couples the orbital angular momentum **L** to the spin momentum **S** via *H*_{so} = λ**L** ⋅ **S**. In the atomic limit of a single Fe, six Fe 3d valence electrons could ideally have a total spin of *S*_{z} = 2 and angular momentum of *L*_{z} = 2.The atomic limit of the SOC energy λ**L** ⋅ **S** is thus 75 meV, given by the SOC coefficient λ ≈ 19 meV (*10*). However, in all existing discussions of Fe-based thin films on MgO substrates, such as Fe/MgO-based (*11*–*15*) and CoFeB/MgO-based systems (*16*–*18*), the size of PMA is only 1 meV, far below the atomic limit. This leaves a vast window to escalate PMA in transition metal thin films unexploited.

Under crystal field, five *d* orbitals are superposed and form *xy*, *yz*, *xz*, *x*^{2} − *y*^{2}, and 3*z*^{2} − *r*^{2} orbitals as the eigenstates. All of the new orbitals have zero *L*_{z} because of the time reversal symmetry. If these orbitals are nondegenerate, the first-order energy correction from SOC vanishes, leaving the second-order perturbation as the dominant contribution (*11*, *19*). This is the scenario in most thin film systems (*11*–*18*, *20*). The energy scale of PMA is then λ^{2}/Δ, where Δ is the bandwidth of the state crossing the Fermi level. For a typical 3*d* magnetic element, Δ ~ 1 eV and λ ~ 0.03 eV (*10*). It is thus not surprising to achieve 1-meV PMA in most 3*d* magnet thin films.

To escalate the PMA, one thus would like to find a regime in which the first-order perturbation of SOC is dominant. In this regime, PMA is proportional to λ instead and has the chance to approach the atomic limit of SOC energy. It occurs when partially filled degenerate orbitals exist around the Fermi level. A successful example has already been demonstrated in a single adatom (*21*–*24*) or dimer (*25*, *26*) deposited on specific substrates. However, once a thin film is formed, PMA in these systems is greatly reduced and brought back to the second-order perturbation scheme.

Here, we report giant PMA in Fe ultrathin films grown on the wurtzite

## RESULTS

As a central result of this work, total energies with different magnetization directions of 1 monolayer (1ML) Fe on BN, AlN, GaN, and InN, respectively, were obtained. The relative magnetoanisotropy energy as a function of *z* direction. At small θ, a linear relation between magnetoanisotropy energy and *3*, *27*). The lowest energy lies at θ = 0 so that the uniaxial anisotropy along the *z* direction is clearly identified. When θ approaches to angle π/2, a clear deviation from *19*, *28*). PMA, the energy difference between perpendicular magnetization and in-plane magnetization, is thus given by *K*_{u1} + *K*_{u2}. Fitting parameters for Fe/BN, Fe/AlN, Fe/GaN, and Fe/InN are listed in Table 1. No significant changes of *K*_{u1} are found, whereas the fourth-order term grows significantly from BN to InN. Except for Fe/BN, the contribution of *11*–*15*). In particular, PMA of Fe(1ML)/InN is 53.7 meV/u.c., approaching the atomic limit of the SOC energy of 75 meV for an isolated Fe. The PMA values for the other three are also on the same order as that limit.

## DISCUSSION

The giant PMA in Fe/III-V nitrides is considerably beyond the energy scale of the second-order perturbation of SOC. To understand the origin of this giant PMA, we studied the electronic structure of Fe(3*d*) orbitals. Without loss of generality, the Fe(1ML)/GaN system was analyzed in detail below. Figure 2 (A and B) displays the difference between the total charge density of our Fe(1ML)/GaN system and the sum of charge densities of a suspended 1ML Fe and a pure GaN supercell. Electron density is reduced in blue contours, whereas it is increased in yellow ones. Thus, charge transfer occurs from the blue contour to the yellow contour during the formation of the Fe-GaN interface. The yellow contour indicates the formation of strongly polarized Fe–N bonds and the enhancement of in-plane *x*^{2} − *y*^{2} and *xy* orbitals. From blue contours, a significant reduction of Fe’s itinerant and *xz*/*yz* electrons was witnessed. The reduction of itinerant electrons is reasonable because Fe electrons saturate the dangling bonds from N atoms on the N-terminated surface, so that Fe atoms lose electrons and become cations. The ionic behaviors of Fe are doubly confirmed by the Bader charge (*29*–*31*) results (Table 2), of which the difference corresponds to the charge increasing/decreasing on one atom. About 0.4*e*^{−} electrons per Fe atom is transferred to N atoms on the interface. These interfacial N atoms thus have almost the same number of valence electrons as that in bulk GaN. In addition, there is no additional interatomic charge transfer when SOC is included, shown on the last column in Table 2.

To explain the charge transfer from *xz*/*yz* to *x*^{2} − *y*^{2}/*xy* orbitals and thereby identify valence states of Fe cations, we investigated the crystal field and orbital-resolved PDOS of Fe(3*d*) without SOC first. As shown in Fig. 2 (C to F), all five *d* orbitals in the spin majority channel are far below the Fermi level and fully filled. Rich physics is present in the spin minority channel. Three orbitals predominated by *xz*, *yz*, and 3*z*^{2} − *r*^{2} are far above the Fermi surface and almost unoccupied. Double-degenerate orbitals, labeled as *e* orbitals, predominated by *x*^{2} − *y*^{2}/*xy* orbitals are low-lying (Fig. 2F). They are crossing the Fermi level and thus partially occupied. In Fig. 2D, the twofold degeneracy of the *xz*/*yz*-predominated orbitals, labeled as *e*′ orbitals, is explicitly shown. Double degeneracies of *e* and *e*′ orbitals are protected by the two-dimensional irreducible representation *E* of the *C*_{3v} point group of the crystal field around each Fe cation. In reality, an overlap between *xz*/*yz* and *x*^{2} − *y*^{2}/*xy* orbitals is present but small. According to the density matrix of Fe(3*d*) orbitals and the corresponding occupation number, *e* states are mixed with 3% *xz*/*yz* orbitals, and *e*′ orbitals contain 3% *x*^{2} − *y*^{2}/*xy* components.

In spin minority channels, *e* states are almost half-filled, with an occupation number of 0.457, whereas occupation number of *e*′ states is only 0.035. As a comparison, in suspended 1ML Fe, the occupation number of *x*^{2} − *y*^{2}/*xy* is 0.095, and that of *xz*/*yz* is 0.564. Therefore, *x*^{2} − *y*^{2}/*xy* orbitals have escalated occupation once Fe is deposited on GaN. It is consistent with the charge density contours discussed earlier.

Once SOC is included, one can expect the lift of degeneracy between *x*^{2} − *y*^{2} and *xy* orbitals due to nonzero off-diagonal matrix elements. It is confirmed by the PDOS shown in Fig. 2E, where almost no states in the spin minority channel are found near the Fermi level. Degenerate *e* states are split, and a large splitting of about 3.0 eV is present. According to the density matrix, the occupation number of the lower splitting state is 0.904, and the corresponding eigenstate is

In terms of spherical harmonics, this state is quite close to *e*′ states. Occupation numbers of those three states are 0.042, 0.033, and 0.032, respectively; that is, almost empty. The 3*z*^{2} − *r*^{2} (that is, *d*) along the *z* direction is 1.54 μ_{B}, consistent with *L*_{z} = 2 due to the splitting into *e* states near the Fermi level. The spin moment in the unit cell from this self-consistent calculation is 3.84 μ_{B}. It is quite consistent with *S* = 2, the high-spin configuration of Fe^{2+} cation, with the fully filled state in the spin majority and only one electron occupied on *E* of *E* = 0.904 × λ(2 × 2α + 1 × 2β) ≈ 32.9 meV, where λ ≈ 19 meV is the SOC coefficient (*10*). It matches well with the final PMA of 32.5 meV for Fe(1ML)/GaN. Therefore, PMA in Fe/GaN thin films is dominated by the first-order perturbation of SOC.

According to the discussion above, large band splitting and the partial occupation in consequence are the precursor of large PMA. However, one should note that the SOC of Fe(3*d*) is on the scale of 20 meV, two orders of magnitude smaller than the bandwidth (~ 2.2 eV) of *x*^{2} − *y*^{2}/*xy* orbitals, labeled as Δ in Fig. 2F. SOC alone can hardly generate such a large splitting of the entire band. This is resonated by the fact that, in simple non–self-consistent calculations, the PMA contributed by SOC alone is only 3.0 meV, a typical value in the second-order correction scheme. Because the spin splitting changes slightly after SOC is included and no structural reconstruction driven by SOC happens, SOC cannot be the major driving force of the band splitting. The only interaction on the scale of electron volts under investigation is the on-site electron-electron correlation interaction described by the Hubbard *U*. It thus suggests that the correlation interaction between occupied electrons must play a significant role on large PMA here.

The energy contribution from Hubbard *U* can be given by the single-particle expression under the Dudarev formation of L(S)DA + *U* (*32*); *U* − *J* = 4.0 eV is the *U* value chosen for our first-principles calculations and *m* in spin channel σ. The energy of the *x*^{2} − *y*^{2}/*xy* orbitals. On the other hand, no band splitting takes place if SOC is turned off because both *e* and *e*′ states receive the same energy shift from the Hubbard *U* and acquire the same occupation number due to the degeneracy. Therefore, the band splitting is triggered by SOC but amplified by the Hubbard *U*.

To further confirm this conclusion, we performed the SOC-included self-consistent calculations with multiple values of the Hubbard *U*. As shown in Fig. 2G, the splitting between *U* decreases. Figure 2H shows that PMA keeps a high value when *U* = 2.0 to 5.0 eV and drops significantly when the Hubbard *U* = 1.0 eV, becoming smaller than the bandwidth of the original *x*^{2} − *y*^{2}/*xy* orbitals. Eventually, when the Hubbard *U* is zero, almost equal populations of *U* = 4.0 eV is reasonable for our calculations because PMA is insensitive to *U* when *U* is higher than 2.0 eV. Thus, the Hubbard *U* plays a key role on the SOC-driven band splitting and the consequential large PMA in this thin film system.

Electronic structures of 1ML Fe on N-terminated (000*d*) *x*^{2} − *y*^{2}/*xy* orbitals without SOC is shown in Fig. 3 (A to C). Again, the *x*^{2} − *y*^{2} and *xy* orbitals are degenerate because of the *C*_{3}*v* symmetry. They are fully filled in the spin majority channel and partially filled in the spin minority one. Although valence states of Fe cations in these structures are similar to that of Fe(1ML)/GaN, the bandwidth shows strong dependence on the lattice constant of the substrate. It is ~5.0, 3.0, and 1.5 eV for Fe/BN, Fe/AlN, and Fe/InN, respectively. Therefore, Fe/BN with the smallest lattice constant has the largest bandwidth due to the large overlap of *x*^{2} − *y*^{2}/*xy* orbitals lying in the plane. A bandwidth of 5.0 eV there exceeds the magnitude of the Hubbard *U* so that the combination of SOC + *U* can hardly induce a large splitting of *x*^{2} − *y*^{2}/*xy* orbitals into

It is confirmed by the PDOS of Fe(3*d*) orbitals with SOC, as shown in Fig. 3 (D to F). Relevant magnetic and orbital properties for each Fe(1ML)/III-V nitride are listed in Table 1. Considerable splitting around 3 eV at Fermi level of spin minority is found in Fe/AlN, Fe/GaN, and Fe/InN but is reduced to ~1.5 eV in Fe/BN, as shown in Fig. 3E. According to Table 1, such small splitting in Fe/BN is consistent with the occupation number of the *x*^{2} − *y*^{2} and *xy* orbitals gives a net orbital magnetic moment along the *z* direction. It is 0.91 μ_{B} for Fe/BN, which is much smaller than 1.44 μ_{B} for Fe/AlN, 1.54 μ_{B} for Fe/GaN, and 1.51 μ_{B} for Fe/InN. Spin moment for Fe/BN is 3.56 μ_{B}. It is also the smallest among all III-V nitrides. As a result, Fe/BN has the lowest PMA, at 24.1 meV/u.c., of all four materials under investigation. Still, it is an order of magnitude larger than any other transition metal thin films ever reported. On the other side, the large lattice constant of InN leads to small energy dispersion of the degenerate *x*^{2} − *y*^{2}/*xy* orbitals and therefore results in a nearly full splitting between

We further investigated the thickness dependence of PMA by using the Fe/GaN ^{2}) and 21 meV/u.c. (33.2 mJ/m^{2}), respectively, which are still large enough to be considered in the first-order SOC perturbation regime.

Fe thin films grown on the (*33*, *34*). Because the X-terminated surface is found to be the most stable structure, at least for 1 × 1 (*35*, *36*), 1 ML of Fe is grown on the X-terminated _{2} or NH_{3} plasma are deposited into the space between an Fe adlayer and the top X layer and gradually diffuse laterally under the Fe adlayer. Finally, a whole ML of N can be formed between Fe and X layers, resulting in a good surface morphology of Fe(1ML)/III-V nitrides.

## CONCLUSION

In summary, *C*_{3v} symmetry. Using noncollinear spin-polarized first-principles calculations, we discovered a giant PMA in the 1ML Fe thin film on this *L*_{z} = 2 that originated from the symmetry-guaranteed splitting of *x*^{2} − *y*^{2} and *xy* orbitals. The on-site correlation interaction amplifies the splitting so that PMA is dominated by first-order perturbation, which is linearly proportional to the strength of SOC of Fe. Thickness dependence shows that PMA keeps a large value in multiple Fe layers deposited on GaN. It eases the experimental realization of our theoretical prediction.

In the rapidly developing technology of MRAM, the lack of large PMA becomes a bottleneck in downsizing the binary bits. Giant PMA discovered here suggests that a 2.0 nm × 2.0 nm flake of Fe(1ML)/III-V nitride has a total uniaxial magnetic anisotropy energy of about 1.2 eV, reaching the criteria for 10-year data retention at room temperature (*3*). Therefore, giant PMA in this thin film can ultimately lead to nanomagnetism and promote revolutionary ultrahigh storage density in the future. Furthermore, large anisotropy energy could lead to large coercivity. Fe/III-V nitride could lead to a new type of permanent magnet without rare earth element potentially.

## METHODS

The calculations were carried out in the framework of the noncollinear spin-polarized first-principles calculations with the projector-augmented wave pseudopotential (*37*) implemented in the Vienna ab initio simulation package (*38*). We used the generalized gradient approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) formation (*39*) plus Hubbard *U* (GGA + *U*) (*32*) with *U* = 4.0 eV on Fe(3*d*) orbitals.

To build the slab supercell, we used four X-N (X = B, Al, Ga, and In) principal layers as the substrate and deposited 1 to 3 ML of Fe on the N-terminated

Charge density of the SOC-free ground state was used as the initial state. Self-consistent total energy calculations were used to derive the noncollinear calculation with SOC included. Γ-centered 25 × 25 × 1 K-point meshes in the two-dimensional Brillouin zone were used with an energy cutoff of 600 eV for the plane-wave expansion. The accuracy of the total energy is thus guaranteed to be better than 0.1 meV/u.c.

To obtain a reasonable *U* value for GGA + *U* calculations, we compared density of state (DOS) and PDOS of 1ML Fe on GaN*U* at *U* = 4.0 eV with those by GGA and the Heyd, Scuseria, and Ernzerhof hybrid functional (HSE06) (*40*), respectively. HSE06 well describes the occupied states in magnetic systems by hybridizing the exact Fock exchange energy with GGA. To integrate over occupied states, we used tetrahedron smearing in GGA and GGA + *U* and the Gaussian smearing with a σ value of 0.05 eV in HSE06. According to the PDOS shown in Fig. 4, spin majority channels of Fe(3*d*) for GGA, GGA + *U* at *U* = 4.0 eV, and HSE06 were mainly located at −4.0 to −1.0 eV, −7.0 to −4.0 eV, and −7.5 to −4.5 eV, respectively. Meanwhile, three peaks of the spin minority channel were located at −0.4, 1.6, and 2.7 eV in HSE06 and at −0.3, 1.3, and 2.3 eV in GGA + *U*, quite consistent with each other. As a comparison, all spin minority states were hybridized near the Fermi level from pure GGA calculations, which were completely different from the HSE06 result. Therefore, we used GGA + *U* (*U* = 4.0 eV) for all calculations except for special notes.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under award no. DE-SC0016424 and used the Extreme Science and Engineering Discovery Environment under grant no. TG-PHY170023.

**Author contributions:**J.-X.Y. and J.Z. conceived the project together. J.-X.Y. carried out the first-principles calculations. Both authors discussed the results and wrote the manuscript.

**Competing interests:**J.-X.Y. and J.Z. are inventors on a U.S. provisional patent application related to this work filed by the University of New Hampshire (application no. 62/589,399; filed 21 November 2017). The authors declare that they have no other competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper. Additional data related to this paper may be requested from the authors.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).