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 
N-terminated surface of III-V nitrides from first-principles calculations. PMA ranges from 24.1 meV/u.c. in Fe/BN to 53.7 meV/u.c. in Fe/InN. Symmetry-protected degeneracy between x2 − y2 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 3d 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 Hso = λL ⋅ S. In the atomic limit of a single Fe, six Fe 3d valence electrons could ideally have a total spin of Sz = 2 and angular momentum of Lz = 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, x2 − y2, and 3z2 − r2 orbitals as the eigenstates. All of the new orbitals have zero Lz 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 3d magnetic element, Δ ~ 1 eV and λ ~ 0.03 eV (10). It is thus not surprising to achieve 1-meV PMA in most 3d 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 
N-terminated surface of III-V nitrides XN, where X = B, Al, Ga, and In (Fig. 1A). First-order perturbation of SOC is exactly the mechanism responsible, and the atomic energy limit is approached.
(A) Surface structures of one-monolayer (1ML) Fe deposited on 
N-terminated surface of III-V nitride XN (X = B, Al, Ga, and In) substrate in top, side, and perspective views, respectively. In top view, three dash-dotted lines indicate three mirror planes of the C3v point group on Fe. (B) Crystal field diagram for Fe(3d)’s spin majority (spin up) and minority (spin down) channels. The spin minority channel is further split by SOC. EF refers to the Fermi level. (C) The relative total energy per unit cell (u.c.) of Fe(1ML)/III-V nitride thin film as a function of 
, where θ denotes the angle between magnetization orientation and the z direction. Relative energy at 
gives the value of PMA. The dotted curves are their fit to 
.
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 
is shown in Fig. 1, where θ is the angle between the magnetization and z direction. At small θ, a linear relation between magnetoanisotropy energy and 
is observed. It is consistent with well-adopted descriptions of PMA in most thin film systems (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 
is present, and high-order anisotropy 
contributes to PMA as well. Therefore, we fit the magnetoanisotropy energy to 
(19, 28). PMA, the energy difference between perpendicular magnetization and in-plane magnetization, is thus given by Ku1 + Ku2. Fitting parameters for Fe/BN, Fe/AlN, Fe/GaN, and Fe/InN are listed in Table 1. No significant changes of Ku1 are found, whereas the fourth-order term grows significantly from BN to InN. Except for Fe/BN, the contribution of 
to PMA is considerable and is even about twice the second-order contribution in Fe/InN. The resulting PMA values range from 24.1 to 53.7 meV/u.c. They are more than one order of magnitude larger than the PMA of Fe/MgO, which is around 1 to 2 meV/u.c. (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.
The values of Ku1, Ku2, total PMA in units of millielectron volt per unit cell and millijoule per square meter, spin moments (ms), orbital moments (ml), and occupation number 
of the 
-dominated state in spin minority are listed.
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(3d) 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 x2 − y2 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.4e− 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.
(A) The positive and (B) negative part of the charge difference between the total charge density of Fe(1ML)/GaN and the sum of charge densities of a suspended 1ML Fe and pure GaN supercell. Charge deficiency in xz/yz-like orbitals shown in (A) is transferred to x2 − y2/xy-like orbitals shown in (B). (C to F) Orbital-resolved projected density of state (PDOS) of 3z2 − r2, xz/yz, and x2 − y2/xy orbitals, respectively, in the absence of SOC. The bandwidth Δ of x2 − y2/xy is labeled in (F). (E) PDOS for Fe(3d) orbitals with SOC included. Positive and negative values of PDOS refer to spin majority and spin minority channels, respectively, and the Fermi level is set to zero. (G) Occupation number of each Fe(3d) orbital and (H) PMA as a function of the Hubbard U when SOC is included.
Bulk refers to the bulk Fe and GaN; surface refers to the clean GaN 
N-terminated surface; interface refers to 1ML Fe on the GaN substrate. The last column gives the results with SOC included.
To explain the charge transfer from xz/yz to x2 − y2/xy orbitals and thereby identify valence states of Fe cations, we investigated the crystal field and orbital-resolved PDOS of Fe(3d) 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 3z2 − r2 are far above the Fermi surface and almost unoccupied. Double-degenerate orbitals, labeled as e orbitals, predominated by x2 − y2/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 C3v point group of the crystal field around each Fe cation. In reality, an overlap between xz/yz and x2 − y2/xy orbitals is present but small. According to the density matrix of Fe(3d) orbitals and the corresponding occupation number, e states are mixed with 3% xz/yz orbitals, and e′ orbitals contain 3% x2 − y2/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 x2 − y2/xy is 0.095, and that of xz/yz is 0.564. Therefore, x2 − y2/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 x2 − y2 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 
or equivalently
(1)
In terms of spherical harmonics, this state is quite close to 
, which has a nonzero expectation value of the SOC energy. Similarly, the higher splitting state is close to 
but hybrid with e′ states. Occupation numbers of those three states are 0.042, 0.033, and 0.032, respectively; that is, almost empty. The 3z2 − r2 (that is, 
)−dominated state is insensitive to SOC, and its occupation number is 0.293. Therefore, the net orbital magnetic moment on Fe(3d) along the z direction is 1.54 μB, consistent with Lz = 2 due to the splitting into 
and 
in 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 Fe2+ cation, with the fully filled state in the spin majority and only one electron occupied on 
in the spin minority. This result is demonstrated by Fig. 1B. With the occupation and orbital components derived, one can estimate the SOC energy ΔE of 
by Δ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(3d) is on the scale of 20 meV, two orders of magnitude smaller than the bandwidth (~ 2.2 eV) of x2 − y2/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); 
, where U − J = 4.0 eV is the U value chosen for our first-principles calculations and 
denotes the occupation number of orbital m in spin channel σ. The energy of the 
-dominated state is − 1.62 eV and that of the 
-dominated state is 1.98 eV, leading to a total splitting of 3.6 eV. It is consistent with the splitting in PDOS (Fig. 2D) and well exceeds the bandwidth of the original x2 − y2/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 
- and 
-dominated states is reduced when 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 x2 − y2/xy orbitals. Eventually, when the Hubbard U is zero, almost equal populations of 
and 
states are recovered, and the magnitude of PMA is only 1.88 meV/u.c., entering the regime dominated by the second-order correction of SOC. These results also confirm that the value 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
) surface of BN, AlN, and InN share the same behavior. The PDOS of Fe(3d) x2 − y2/xy orbitals without SOC is shown in Fig. 3 (A to C). Again, the x2 − y2 and xy orbitals are degenerate because of the C3v 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 x2 − y2/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 x2 − y2/xy orbitals into 
and 
.
(A to C) PDOS of x2 − y2/xy orbitals without SOC. The bandwidth in the spin minority channel is indicated by double arrows. (D to F) PDOS of Fe(3d) with SOC. In (D), the shadow region gives the projection onto x2 − y2/xy orbitals where the magnitude of splitting is indicated by the double arrow.
It is confirmed by the PDOS of Fe(3d) 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 
-dominated state, which is 0.724 for Fe/BN, smaller than 0.854 for AlN, 0.904 for GaN, and 0.930 for InN. Lifting of the degeneracy between x2 − y2 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 x2 − y2/xy orbitals and therefore results in a nearly full splitting between 
- and 
-dominated states after SOC is turned on. A PMA of 53.7 meV/u.c. for Fe/InN is consequently the largest and almost hits the atomic limit of anisotropy energy of Fe. Giant PMA for all four systems follows the first-order perturbation scheme of SOC. The magnitude of PMA in units of millijoule per square meter is also listed in Table 1. In these units, PMA in Fe/BN is no longer the smallest because of small unit cell size in BN.
We further investigated the thickness dependence of PMA by using the Fe/GaN 
system as an example. Slab supercells with 2ML and 3ML Fe cations on top of GaN were built following the hexagonal closed packing along wurtzite GaN 
direction. As a result, PMA values for Fe(2ML)/GaN and Fe(3ML)/GaN are 37 meV/u.c. (58.4 mJ/m2) and 21 meV/u.c. (33.2 mJ/m2), respectively, which are still large enough to be considered in the first-order SOC perturbation regime.
Fe thin films grown on the (
) N-terminated surface of III-V nitrides XN can be formed experimentally by the adlayer-enhanced lateral diffusion method (33, 34). Because the X-terminated surface is found to be the most stable structure, at least for 1 × 1 (
) surface of GaN and InN (35, 36), 1 ML of Fe is grown on the X-terminated 
surface first. Then, the N atoms from N2 or NH3 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, 
surface of III-V nitrides provides a crystal field of C3v symmetry. Using noncollinear spin-polarized first-principles calculations, we discovered a giant PMA in the 1ML Fe thin film on this 
N-terminated surface of III-V nitride substrate. PMA ranged from 24.1 meV/u.c. in BN to 53.7 meV/u.c. in InN substrate. They are exceedingly large compared to existing PMA thin films and approach to the atomic limit of SOC energy of Fe. Electronic structure calculations and ligand field analysis show that each Fe cation has a net orbital angular momentum Lz = 2 that originated from the symmetry-guaranteed splitting of x2 − y2 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(3d) 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 
surface. Dangling bonds of the X-terminated (0001) surface at bottom of the substrate were saturated by pseudohydrogen atoms. Vacuum layers in the supercells were more than 14 Å thick, and the lattice constant in the plane was fixed to be 2.554, 3.113, 3.183, and 3.456 Å for BN, AlN, GaN, and InN, respectively, as determined from the GGA-PBE result of the wurtzite phase of III-V nitrides. Positions of X and N atoms were borrowed from their bulk values. Positions of Fe atoms were optimized under collinear magnetic calculations without SOC until the force on each Fe atom was less than 1 meV/Å. It was confirmed that the N top was the lowest energy site of Fe atoms, as shown in Fig. 1A.
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
[Fe(1ML)/GaN] without SOC included by GGA + 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(3d) 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.
DOS (gray dashed line) and PDOS (solid blue line) of Fe(1ML)/GaN without SOC derived by (A) GGA-PBE, (B) GGA + U at U = 4.0 eV and (C) HSE06, respectively. Consistency between (B) and (C) shows that GGA + U at U = 4.0 eV can well describe Fe(3d) orbitals in this system.
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
- 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).
