Validity of perturbative methods to treat the spin–orbit interaction: application to magnetocrystalline anisotropy

Our reference ab initio MAE value is obtained by substracting the total energies E_tot between two fully-relativistic self-consistent calculations, which include SOI, for two different orientations of the magnetization,

$$\begin{eqnarray}&&\mathrm{MAE}={E}_{{\rm{tot}}}^{x}-{E}_{{\rm{tot}}}^{z},\end{eqnarray} \tag{ 2 }$$

where the spins are aligned along the OX and OZ directions shown in the L1₀ unit cell model of figure 1. The main shortcoming of this method is its computational cost, since equation (2) implies substracting two large numbers, which requires demanding convergence criteria. In fact, the obtention of well-converged MAEs from equation (2) is crucial in this work (see details in the next section), since they are used as a benchmark for the approximations described in the next subsections.

A scalar relativistic ground state (GS) is converged in a spin-polarized calculation without SOI. The so-obtained charge density is used to initialize a fully-relativistic calculation (i.e. non-collinear) by turning it into a block-diagonal charge density matrix. Then, the spin–orbit H_SO term calculated for a given magnetization axis is added to the scalar-relativistic hamiltonian H₀ and new eigenvalues are calculated by diagonalization without further self-consistent cycles. We denote the resulting total energy change ${\rm{\Delta }}{E}_{{\rm{tot}}}^{x,z}$ and the corresponding charge density change Δρ^x,z. The MAE is approximated as the difference in the band energies, ${E}_{{\rm{band}}}^{x,z}$, between the two orientations of the magnetization, bearing in mind that the Fermi levels are in general different for the two orientations, as they are computed independently,

$$\begin{eqnarray}&&\mathrm{MAE}\simeq {\rm{\Delta }}{E}_{{\rm{tot}}}^{x}-{\rm{\Delta }}{E}_{{\rm{tot}}}^{z}\simeq {E}_{{\rm{band}}}^{x}-{E}_{{\rm{band}}}^{z}=\displaystyle \sum _{k}^{{N}_{k}}\displaystyle \sum _{n}^{{N}_{b}}[{f}^{x}({\epsilon }_{{kn}}^{x}){\epsilon }_{{kn}}^{x}-{f}^{z}({\epsilon }_{{kn}}^{z}){\epsilon }_{{kn}}^{z}].\end{eqnarray} \tag{ 3 }$$

Here, the sum runs over one-electron eigenvalues ${\epsilon }_{{kn}}^{x,z}$, calculated with the spins aligned along the OX and OZ directions, respectively, and integrated over the entire first Brillouin Zone (1BZ). N_b bands are considered (index n) and a discrete grid of N_k points (index k) is used to sample the 1BZ. f^x,z are the Fermi–Dirac distribution functions, which depend on the magnetization axes through the Fermi energy, while the finite electronic temperature kT acts as a smearing parameter. The approximation is based on the fact that ${\rm{\Delta }}{E}_{{\rm{tot}}}^{x}$ and ${\rm{\Delta }}{E}_{{\rm{tot}}}^{z}$ are correct to order ${({\rm{\Delta }}{\rho }^{x})}^{2}$ and ${({\rm{\Delta }}{\rho }^{z})}^{2}$, respectively. The method is thus sometimes called 'second variation' [10] or 'force theorem' [9, 11]. Equation (3) is correct to order Δρ^x,z while the ${({\rm{\Delta }}{\rho }^{x,z})}^{2}$-order corrections have a small effect, since there are cancellations from the two magnetization directions [9]. If one further assumes that the self-consistency cycles introduce negligible modifications in the charge density matrix and in the exchange and correlation potential, then equations (2) and (3) should provide very similar MAE values, although with a considerable reduction in the computational cost in the latter case.

A widely used alternative approximation treats the SOI as a second-order perturbation (2PT) to the many-body GS, $| {{\rm{\Psi }}}^{(0)}\rangle$. The general expression for the 2PT energy correction is

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{\rm{SO}}}^{(2)}=\displaystyle \sum _{i\ne 0}\displaystyle \frac{\langle {{\rm{\Psi }}}^{(0)}| {H}_{{\rm{SO}}}| {{\rm{\Psi }}}^{(i)}\rangle \langle {{\rm{\Psi }}}^{(i)}| {H}_{{\rm{SO}}}| {{\rm{\Psi }}}^{(0)}\rangle }{{E}_{0}-{E}_{i}},\end{eqnarray} \tag{ 4 }$$

where the sum runs over excited states. In a many-body language the GS $| {{\rm{\Psi }}}^{(0)}\rangle$ is formed as a Slater determinant by occupation of the lowest-lying one-electron Kohn–Sham eigenstates up to the Fermi level. Each excited state $| {{\rm{\Psi }}}^{(i)}\rangle$ is then constructed by creating electron–hole (e–h) pairs using the unoccupied Kohn–Sham eigenstates. Thus, the ${E}_{0}\mbox{--}{E}_{i}$ term in the denominator is simply the energy difference between the occupied and unoccupied eigenvalues associated to the particular e−h excitation [13, 18–20]. The perturbative expansion may then be written in terms of the GS Kohn–Sham eigenstates $| {kn}\sigma \rangle$ as follows:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{\rm{SO}}}^{(2)}=\displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{k}^{{N}_{k}}\displaystyle \sum _{n,n^{\prime} }^{{N}_{b}}\displaystyle \sum _{\sigma ,\sigma ^{\prime} }\displaystyle \frac{f({\epsilon }_{{kn}\sigma })[1-f({\epsilon }_{{kn}^{\prime} \sigma ^{\prime} })]}{{\epsilon }_{{kn}\sigma }-{\epsilon }_{{kn}^{\prime} \sigma ^{\prime} }}\times \langle {kn}\sigma | {H}_{{\rm{SO}}}| {kn}^{\prime} \sigma ^{\prime} \rangle \langle {kn}^{\prime} \sigma ^{\prime} | {H}_{{\rm{SO}}}| {kn}\sigma \rangle ,\end{eqnarray} \tag{ 5 }$$

where $\sigma ,\sigma ^{\prime}$ stand for the spin indices.

The second-order formula given by equation (5) is applicable only to a non-degenerate GS. A degenerate GS in an extended metallic system happens when there is a band crossing at a certain k-point precisely at the Fermi level and this band pair is coupled by H_SO (i.e. the corresponding matrix element is not zero). This can happen eventually, and in this situation the eigenstate pair should be treated separately by first-order degenerate state perturbation theory. However, in a calculation with a large N_k the contribution to ${\rm{\Delta }}{E}_{{\rm{SO}}}^{(2)}$ of these exactly degenerate states would be negligible, since only a handful of band crossings are expected at the Fermi level, and they contribute with a factor of order ξ/N_k [37]. A sufficiently fine k-grid can map the spin–orbit band splitting effect nearby the crossing, so that equation (5) can be safely used.

In practice, it is convenient to express equation (5) in a basis whose elements have well defined lm quantum numbers (spherical harmonics Y_lm). A natural choice are AOs, which already constitute the basis set in a number of DFT-based packages, leading to Bloch Kohn–Sham eigenstates of the form:

$$\begin{eqnarray}&&| {kn}\sigma \rangle =\displaystyle \frac{1}{\sqrt{{N}_{k}}}\displaystyle \sum _{{\bf{R}}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{R}}}\displaystyle \sum _{\alpha }^{{N}_{\alpha }}{c}_{\alpha }^{n\sigma }(k)| \alpha ({\bf{R}}),\sigma \rangle ,\end{eqnarray} \tag{ 6 }$$

where R runs over the lattice vectors, $| \alpha ({\bf{R}})\rangle$ denotes an AO located in unit cell R and N_α is the total number of AOs in the basis set. Inserting equation (6) into 5 yields:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{\rm{S}}{\rm{O}}}^{(2)}=\displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{k}^{{N}_{k}}\displaystyle \sum _{\sigma {\sigma }^{{\rm{{\prime} }}}}\displaystyle \sum _{n,{n}^{{\rm{{\prime} }}}}^{{N}_{\alpha }}\displaystyle \frac{f({\epsilon }_{kn\sigma })[1-f({\epsilon }_{{kn}^{{\rm{{\prime} }}}{\sigma }^{{\rm{{\prime} }}}})]}{{\epsilon }_{kn\sigma }-{\epsilon }_{{kn}^{{\rm{{\prime} }}}{\sigma }^{{\rm{{\prime} }}}}}\displaystyle \sum _{\alpha \beta ,{\alpha }^{{\rm{{\prime} }}}{\beta }^{{\rm{{\prime} }}}}^{{N}_{\alpha }} \langle \alpha | {H}_{{\rm{S}}{\rm{O}}}^{\sigma {\sigma }^{{\rm{{\prime} }}}}(k)| \beta \rangle \langle {\alpha }^{{\rm{{\prime} }}}| {H}_{{\rm{S}}{\rm{O}}}^{{\sigma }^{{\rm{{\prime} }}}\sigma }(k)| {\beta }^{{\rm{{\prime} }}} \rangle {n}_{\alpha \beta }^{{nn}^{{\rm{{\prime} }}}\sigma {\sigma }^{{\rm{{\prime} }}}}(k){n}_{{\alpha }^{{\rm{{\prime} }}}{\beta }^{{\rm{{\prime} }}}}^{{n}^{{\rm{{\prime} }}}n{\sigma }^{{\rm{{\prime} }}}\sigma }(k),\end{eqnarray} \tag{ 7 }$$

where the k-space H_SO(k) matrix elements are given by:

$$\begin{eqnarray}&&\langle \alpha | {H}_{{\rm{SO}}}^{\sigma \sigma ^{\prime} }(k)| \beta \rangle =\displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{{\bf{R}}}{{\rm{e}}}^{{\rm{i}}{\bf{k}}\cdot {\bf{R}}}\langle \alpha (0),\sigma | {H}_{{\rm{SO}}}| \beta ({\bf{R}}),\sigma ^{\prime} \rangle \end{eqnarray} \tag{ 8 }$$

and the generalized projected charges by:

$$\begin{eqnarray}&&{n}_{\alpha \beta }^{{nn}^{\prime} \sigma \sigma ^{\prime} }(k)={\left({c}_{\alpha }^{n\sigma }(k)\right)}^{* }{c}_{\beta }^{n^{\prime} \sigma ^{\prime} }(k).\end{eqnarray} \tag{ 9 }$$

It is usual to further assume the so-called on-site approximation, whereby the l_i · s_i operators only mix states within the same l-shell of a given atom contained in the origin unit cell. Equation (8) then becomes:

$$\begin{eqnarray}&&\langle \alpha | {H}_{{\rm{SO}}}^{\sigma \sigma ^{\prime} }(k)| \beta \rangle \approx \langle \alpha (0),\sigma | {H}_{{\rm{SO}}}| \beta (0),\sigma ^{\prime} \rangle \approx {\xi }_{\alpha ,l}\langle \alpha {lm}\sigma | {\bf{l}}\cdot {\bf{s}}| \beta l^{\prime} m^{\prime} \sigma ^{\prime} \rangle {\delta }_{\alpha \beta }{\delta }_{{ll}^{\prime} }\end{eqnarray} \tag{ 10 }$$

where indices $\alpha ,\beta$ in the last term now refer to the principal quantum number in a given atom (in a multi-ζ scheme, also the particular ζ [38]), and lm stand for the orbital and magnetic quantum numbers of each AO. ξ_α,l is the SOI strength for this l-shell resulting from the integration of the radial part in the $\langle \alpha {lm}\sigma | {H}_{{\rm{SO}}}| \alpha {lm}^{\prime} \sigma ^{\prime} \rangle$ matrix elements (independent of ${mm}^{\prime}$ and $\sigma \sigma ^{\prime}$). The angular part of these matrix elements, $\langle \alpha {lm}\sigma | {\bf{l}}\cdot {\bf{s}}| \alpha {lm}^{\prime} \sigma ^{\prime} \rangle$, take simple analytical expressions and tabulated formulas as a function of the spherical harmonics Y_lm involved can be found, for example, in [39, 40]. The on-site approximation simplifies considerably the 2PT formula:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{\rm{SO}}}^{(2)}=\displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{k}^{{N}_{k}}\displaystyle \sum _{\sigma \sigma ^{\prime} }\displaystyle \sum _{n,n^{\prime} }^{{N}_{\alpha }}\displaystyle \frac{f({\epsilon }_{{kn}\sigma })[1-f({\epsilon }_{{kn}^{\prime} \sigma ^{\prime} })]}{{\epsilon }_{{kn}\sigma }-{\epsilon }_{{kn}^{\prime} \sigma ^{\prime} }}{A}^{{nn}^{\prime} \sigma \sigma ^{\prime} }(k){A}^{n^{\prime} n\sigma ^{\prime} \sigma }(k),\end{eqnarray} \tag{ 11 }$$

where we have defined:

$$\begin{eqnarray}&&{A}^{{nn}^{\prime} \sigma \sigma ^{\prime} }(k)=\displaystyle \sum _{\alpha l}{\xi }_{\alpha l}\displaystyle \sum _{{mm}^{\prime} }^{2l+1}\langle \alpha {lm}\sigma | {\bf{l}}\cdot {\bf{s}}| \alpha {lm}^{\prime} \sigma ^{\prime} \rangle {n}_{\alpha {lm}\alpha {lm}^{\prime} }^{{nn}^{\prime} \sigma \sigma ^{\prime} }(k).\end{eqnarray} \tag{ 12 }$$

We observe that the 2PT equations (7) and (11) are perturbative in the SOI strength, which appears explicitly in the form of the parameters ${\xi }_{\alpha l}$ in this equation.

Next, we consider the case when the unperturbed spin-polarized calculation is realized employing a plane-wave basis set, as many DFT codes do. Instead of using a Bloch-function representation of the Kohn–Sham eigenstates, we use a set of N_w maximally localized Wannier functions (MLWF) as formulated in [41]. MLWFs are constructed to yield the exact eigenvalues as the ab initio calculation. Usually, a previous disentanglement procedure is carried out, whereby a handful of relevant bands within an energy window are isolated from the rest [42]. For the systems under study, we focus on the bands that originate from the d-valence electrons of both metal atoms (allowing also for some degree of s–d hybridization) and belong to a window of about 10 eV below and 5 eV above the Fermi energy. Afterwards, it is straightforward to obtain new Bloch functions on a k-grid as dense as desired by interpolation [43]. This procedure allows us to estimate the 2PT MAE using as input solely a scalar-relativistic first-principles calculation and does not require a highly dense 1BZ k-sampling.

The jth MLWF localized at the unit cell R that results from band disentanglement and wannierization for states with spin σ is:

$$\begin{eqnarray}&&| {w}_{j}^{\sigma }({\bf{R}})\rangle =\displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{k}^{{N}_{k}}{{\rm{e}}}^{-{\rm{i}}{\bf{k}}\cdot {\bf{R}}}\displaystyle \sum _{n}^{{N}_{b}}{Q}_{j}^{n\sigma }(k)| {kn}\sigma \rangle ,\end{eqnarray} \tag{ 13 }$$

where $| {kn}\sigma \rangle$ stands for the Kohn–Sham eigenstates already interpolated in the dense k-grid [43] and ${Q}_{j}^{n\sigma }({\bf{k}})$ are the coefficients that relate the Wannier and Bloch functions.

Typically, atomic-like wavefunctions, formed by a radial function and spherical harmonics Y_lm to describe the angular component, are used to initialize the wannierization procedure. Here, we use d-orbitals centred at the atomic sites and a few s-waves at interstitial positions. The purpose of the latter functions is to facilitate the wannierization and will not take part in the 2PT MAE calculation. We fix l = 2 and drop this index in the following. Thus, the j label accounts for the α-th atom in the cell R and the m = 0, ±1, ±2 quantum number. If the deviation of the MLWFs from actual atomic wavefunctions is small, we can approximate the matrix elements $\langle {w}_{\alpha m}^{\sigma }({\bf{R}})| {H}_{{\rm{SO}}}| {w}_{\alpha m^{\prime} }^{\sigma ^{\prime} }({\bf{R}})\rangle$ by the ones in the AO representation $\langle \alpha m\sigma | {\bf{l}}\cdot {\bf{s}}| \alpha m^{\prime} \sigma ^{\prime} \rangle$ and take advantage of their simple analytic expressions [39, 40]. Note that, in general, the resulting MLWFs do not keep a well-defined orbital character because they have to account for both the intra- and inter-AO hybridization present in the system [44]. Nevertheless, a suitable choice of axes in the systems under study allows to obtain MLWFs that keep the AO character and justify this approximation. Substituting equation (13) in (5) and using this approach, the second-order energy correction associated to the SOI is

$$\begin{eqnarray}&&{\rm{\Delta }}{E}_{{\rm{S}}{\rm{O}}}^{(2)}=\displaystyle \sum _{{\alpha }_{1}{\alpha }_{2}}\displaystyle \sum _{\sigma {\sigma }^{{\rm{{\prime} }}}}\displaystyle \sum _{{m}_{1}{m}_{2}}\displaystyle \sum _{{m}_{1}^{{\rm{{\prime} }}}{m}_{2}^{{\rm{{\prime} }}}}{\xi }_{{\alpha }_{1}}{\xi }_{{\alpha }_{2}}{F}_{{m}_{1}{m}_{2}{m}_{1}^{{\rm{{\prime} }}}{m}_{2}^{{\rm{{\prime} }}}}^{\sigma {\sigma }^{{\rm{{\prime} }}}{\alpha }_{1}{\alpha }_{2}} \langle {\alpha }_{1}{m}_{1}^{{\rm{{\prime} }}}{\sigma }^{{\rm{{\prime} }}}| {\bf{l}}\cdot {\bf{s}}| {\alpha }_{1}{m}_{1}\sigma \rangle \langle {\alpha }_{2}{m}_{2}^{{\rm{{\prime} }}}\sigma | {\bf{l}}\cdot {\bf{s}}| {\alpha }_{2}{m}_{2}{\sigma }^{{\rm{{\prime} }}} \rangle ,\end{eqnarray} \tag{ 14 }$$

where the m_i indexes run over the 5 d-orbitals of each atom α_i in the unit cell. The F coefficients contain the details of the basis change from Kohn–Sham states to MLWF states and, implicitly, they allow us to use a dense k-grid by interpolation:

$$\begin{eqnarray}&&{F}_{{m}_{1}{m}_{2}{m}_{1}^{{\rm{{\prime} }}}{m}_{2}^{{\rm{{\prime} }}}}^{\sigma {\sigma }^{{\rm{{\prime} }}}{\alpha }_{1}{\alpha }_{2}}=\displaystyle \frac{1}{{N}_{k}}\displaystyle \sum _{k}^{{N}_{k}}\displaystyle \sum _{{nn}^{{\rm{{\prime} }}}}^{{N}_{w}}\displaystyle \frac{f({\epsilon }_{kn\sigma })[1-f({\epsilon }_{{kn}^{{\rm{{\prime} }}}{\sigma }^{{\rm{{\prime} }}}})]}{{\epsilon }_{kn\sigma }-{\epsilon }_{{kn}^{{\rm{{\prime} }}}{\sigma }^{{\rm{{\prime} }}}}}{Q}_{{\alpha }_{1}{m}_{1}^{{\rm{{\prime} }}}}^{{n}^{{\rm{{\prime} }}}{\sigma }^{{\rm{{\prime} }}}}(k){\left({Q}_{{\alpha }_{1}{m}_{1}}^{n\sigma }(k)\right)}^{\star }{Q}_{{\alpha }_{2}{m}_{2}^{{\rm{{\prime} }}}}^{n\sigma }(k){\left({Q}_{{\alpha }_{2}{m}_{2}}^{{n}^{{\rm{{\prime} }}}{\sigma }^{{\rm{{\prime} }}}}(k)\right)}^{\star }.\end{eqnarray} \tag{ 15 }$$

Equation (14) is similar to the 2PT MAE formula for extended systems developed by van der Laan [20], and here we generalize the result to the case of more than one atom per unit cell. It is also straigthforward to show that the on-site approximation equation (11) reduces to the above expression after replacing the Bloch coefficients ${c}_{\alpha }^{n\sigma }(k)$ in equation (9) by the ${Q}_{j}^{n\sigma }(k)$ ones and restricting the summation over the angular momentum numbers to the $l=l^{\prime} =2$ case.

Note that equation (14) (or equation (7)) is a 'four-legged' expression in the sense that is contributed by two different e–h pairs. There are other 2PT formulations based on the use of a localized basis set of orbitals [13, 15, 18, 19]. However, it is worth to mention that our formulation of the MAE does not neglect spin-flip contributions, unlike the one proposed by Bruno [18]. Formulas that neglect wavefunction phase effects in the Kohn–Sham-to-local-basis projection have also been proposed [13, 15]. By doing so, the 'four-legged' equation (14) becomes only 'two-legged' and equation (15) can be written in simpler terms, namely the projected densities of states on the local orbitals.

Table 4 contains the collection of the MAE values of the L1₀ alloys calculated with the methodologies presented in section 3. All the approaches provide the same behaviour in the MCA of each alloy, albeit there are some quantitative differences in the corresponding MAE values, which will be discussed below. The easy magnetization axis is OZ in the five studied systems. Overall, MAE values are smaller than 0.5 meV except for FePt, a well-known prototype of large MCA, which shows a MAE of nearly 3 meV in good agreement with other ab initio values available in the literature [46–48, 59].

The first important observation is that the MAE of Fe-based L1₀ alloys is not directly correlated with the SO strength of the alloying metal. Indeed, it is the electronic structure what governs the MCA of these alloys, overruling the effect of the SO strength. We find larger MAE values for FeCo and FeCu than for FePd and FeAu in spite of ${\xi }_{\mathrm{Pd},\mathrm{Au}}$ being larger than ${\xi }_{\mathrm{Co},\mathrm{Cu}}$ (see table 3). In the case of FeCo it is known that a large MAE can be achieved by a strain on the cell along OZ. For c/a = 1.2, the geometry chosen for this work, a maximum is obtained because degenerate states coupled by H_SO lie on the Fermi level [30]. As a matter of fact, in this respect, FeCo is not an isolated case [15, 60] .

The FSC values obtained in the VASP and SG calculations are in good agreement for FeCo, FeCu, and FePd, where discrepancies of 0.07 meV or smaller are found. For FePt and FeAu larger deviations of around 0.3 meV exist. The reason for this discrepancy is difficult to identify. Small quantitative differences in the band structures provided by both codes would be unimportant for most physical properties but, unfortunately, they become significant for the estimation of the MAEs.

However, the MAE values predicted by the FT method (see table 4) are, in general, very close to the FSC values, with typical deviations well below 0.1 meV. It is striking to find such a good agreement even for the heaviest metal systems, where the charge density change is expected to be larger. There are only a few exceptions, namely the VASP calculation for FePt, FeCo and FeAu and the SG one for FeAu for which, nonetheless, differences remain smaller than 0.2 meV. For the formers, we assign the discrepancy to calculation details rather than to the limitation of using the non-self-consistent charge density. Among the technical details, the smearing method used for the Fermi level determination is crucial, since the FT approach is sensitive to the small Fermi energy shift when magnetization occurs in one or other direction. In the FeAu case with SG, where full k-convergence is achieved, we may consider a 0.12 meV deviation as the upper limit to the accuracy of the FT, probably due to the larger ξ_Au SOI strength. Notwithstanding, the H_SO term is fully accounted for by this approach, although not self-consistently.

Table 5 shows the MAE values for the thin films. The deviation of the FT values from the FSC ones is about 0.1 meV for FeCo and 0.3 meV for FePt bilayers. From these results, it is clear that the dimensionality reduction can undermine the performance of FT. This has been observed in Fe and Co adatoms on Rh(111), Pd(111) [12] and Pt(111) [61], too, with agreement between FT and FSC values largely depending on the adsorption site and atomic species [12]. We recall that under the FT the electron density is not allowed to relax when SOI is included. In the FSC calculation we expect these relaxations to be significant in the low-dimensional cases while in bulk 3D systems, with the concomitant symmetry constraints, they are strongly reduced. Therefore, we find the FT approach more reliable for the bulk (table 4) than for thin films (table 5).

Next, we address the performance of the 2PT approximation to the MAE, first focussing on the SG MAEs in table 4: MAE values are in almost perfect agreement with FSC and FT results, the ony exception is the FeAu alloy with a 0.2 meV difference. Thus, the same behaviour between the FT and 2PT methods indicates that their accuracy is very similar despite the neglect of high-order H_SO terms in the latter, both providing excellent results at a much lower computational cost compared to the reference FSC calculations. The 2PT MAE values obtained with the wannierized bands from the VASP calculation yield worse agreement than the other methods becoming more pronounced the heavier the metal in the Fe alloy (see table 4). Specifically, the easy axis for FePd is switched to OX while the MAE of FePt is considerably underestimated and that of FeAu overestimated. Although this method allows the use of very dense k-point grids by interpolation, the quality of the wannierization procedure is crucial for numerical accuracy. Nonetheless, the main factor that undermines the final MAE values is the approximation in the H_SO matrix elements: the assumption that the MLWFs correspond to AOs with well-defined lm quantum numbers and, to a lesser extent, the on-site approximation (equation (10)) and the neglect of sp-orbital contributions. All in all equation (14) is a coarse approximation. Despite the apparent resemblance of MLWFs with atomic wavefunctions by visual inspection (see figure 5 in the supplementary information), the deformations of these 'd-orbitals' are significant, due to the fact that Fe-alloys have strongly-hybridized bands. For FeCo and FeCu, nevertheless, the MAE behaviour in the energy windows analysed is similar to the SG ones (not shown). The performance of the 2PT method for the thin films (table 5) is not as good as in the bulk case. As mentioned above, in a low-dimensional scenario we expect the spin–orbit effects on the electronic structure to be stronger. Thus, due to its perturbative nature, it is more difficult for the 2PT method to capture these effects fully.

Now, using the information provided by the L1₀ alloys we address the important issue of how close to the Fermi level are the states determining the MAE values, the above mentioned usual assumption. The 2PT-derived MAE is determined by eigenstates around the Fermi level within an energy range given by the SOI strength ξ and the ${({\epsilon }_{{kn}\sigma }-{\epsilon }_{{kn}^{\prime} \sigma ^{\prime} })}^{-1}$ factors. Indeed, equation (7) gives less weight to the e−h pairs coupled by H_SO matrix elements that lie farther in energy (the contributions of the energy prefactors are shown in figure 6 of the Supplementary Information). Intuitively, deep states in the valence band might be regarded as negligible for the MAE. However, a third factor needs to be considered: the avaible number of states. To analyse the competition of these three effects, we have calculated the 2PT MAE allowing only initial (final) states within an energy window below (above) the Fermi energy when evaluating equation (7). The results are shown in figure 2 (left) for the SG case, while those with VASP are qualitatively the same (not shown). When imposing an energy window on the occupied states, a plateau in the MAE value is not reached until the window is 5–7 eV wide, depending on the system. These energy ranges comprise the d-band widths below the Fermi level (see the densities of states (DOS) projected on the d-orbitals in figure 2 (right). With heavier atoms, deeper states contribute non-trivially to the MAE, even reversing its sign, as it is the case of FeAu. Therefore, the final value of the MAE of each system depends on its electronic structure details, since the dispersion and binding energies of the individual band pairs coupled by H_SO largely vary from one system to another. The effect of constraining the accessible empty states in figure 2 is less dramatic and reveals a common behaviour in all the systems. With a window above the Fermi level, the MAE plateau is reached at ∼2.5 eV. As shown in the DOS plots, this energy range corresponds to the empty spin-minority states of Fe in all the alloys and other empty metal d-states also contribute, although to a lesser extent, if available (Co and Pt).

Zoom In Zoom Out Reset image size

Therefore, we conclude that, as a general rule, the high DOS counterbalances the decay of the ${({\epsilon }_{{kn}\sigma }-{\epsilon }_{{kn}^{\prime} \sigma ^{\prime} })}^{-1}$ factors and dominates over the ξ prefactors. We can see these trends in the systems under study. In brief, figure 2 shows that states distant from the Fermi level by as much as several eV (that is, spanning the whole d-band of the alloy) cannot be neglected in a 2PT calculation, not even in the case of atoms with weak SOI. We note also that, because of the DOS effect, the observed MAE for FeAu is not the largest, despite of the heavy Au atoms. From this figure, we conclude that the accessibility of empty Fe spin-minority states is a common feature that allows for sizeable MCAs in the Fe-based alloys, while the details of the occupied d-bands of each case determine the final MAE value.

Another appealing aspect of the 2PT formulation is that it allows to split the MAE into contributions arising from pairs of atoms in the metallic alloy: Fe–Fe, Fe–Me and Me–Me. This is straighforward when using MLWFs and equation (14) (see figures 7(a)–(e) in the supplementary information), while if equation (7) is used, the decomposition can be realized by restricting $(\alpha ,\alpha ^{\prime} )$ to a given pair of atoms and performing the summation over all the other $(\beta ,\beta ^{\prime} )$ AO contributions. The result is shown in figures 3(a)–(e). As expected, Fe and Co have a balanced weight on the final FeCo MAE, while Cu hardly contributes to that of FeCu. In the rest of alloys the decompositions show the contribution we could expect from the knowledge of the electronic structure, at least qualitatively. The Pt–Pt and Fe–Pt terms are the main contributors in the FePt alloy, because ξ_Pt is an order of magnitude larger than ξ_Fe and because the d electrons of both species are strongly hybridized. We find, in agreement to [14], that the large contribution of the SOI of Pt atoms to the perpendicular anisotropy (OZ easy axis) is partially balanced by the effect of the Fe–Pt hybrid bands, which favour in-plane anisotropy. In FeAu, the Au–Au term is also very large due to the magnitude of ξ_Au, but now we see that the contribution of Fe–Au terms is weaker, since there is little hybridization with Au-d electrons.

Zoom In Zoom Out Reset image size

Finally, we compare our perturbative results with the widely used 2PT equation by Bruno for bands of d-orbital character [18], which assumes that all the accesible holes are minority spin states. Therefore, it neglects e–h excitations that involve spin-flip, which leads to ${\rm{\Delta }}{E}^{(2)}\propto \xi \langle {\bf{L}}\rangle$. Table 3 shows the MAE values obtained in this approximation with the DFT-calculated AO moment values μ_L for each magnetization direction. Although the density of majority-spin states above the Fermi level is marginal (shown in figure 2), the results of the Bruno formula differ significantly from the other methods, and in some cases it does not even predict the correct easy magnetization axis. In addition to the breakdown of the approximation itself, we can ascribe the discrepancies to the tendency of DFT to underestimate orbital moment values.

Figures 3(f)–(j) shows the contributions of the e–h spin pairs to the MAE of equation (7) (the equivalent result with the MWLF approach is shown in the supplementary information figures 7(f)–(j)). We observe that the spin-flip terms have, indeed, a non-negligible contribution in all cases. In FeCo and FeCu the spin-down final (empty) states clearly dominate the total MAE value, since the contribution from the spin-up final sates is almost negligible due to their very low DOS (see figure 2). However, the interpretation of the e–h spin pair contributions for the rest of alloys is less trivial; in the case of FePt, in particular, the up-up term surprisingly accounts for most of the MAE. A detailed energy-resolved analysis of this counterintuitive result (not shown) reveals that, despite the spin-down final states still yield the largest contributions in absolute value for all alloys, they tend to cancel or even provide net negative values (FePd, FePt, and FeAu) when integrated around the Fermi level, whereas the transitions involving spin-up final states increase in magnitude as the Me atom becomes heavier and tend to take positive values.

Figure 4 shows the MAE as a function of the number of electrons N_e calculated for the L1₀ cases within the FT and the 2PT approaches. Here, we have followed a rigid band approach, in which the Fermi level is recalculated for each N_e employing the eigenvalues and eigenstates of the unperturbed neutral calculation. Hence, the plots represent the MAE behaviour of each alloy under conditions of doping or application of a bias voltage, which are common working conditions in magnetic devices, albeit, due to the rigid band approach, only small deviations from the neutral situation are physically meaningful. Nevertheless, this approach provides a deep understanding of the physical origin of the MCA and yields valuable information about the accuracy of the MAE values, as we show below. As expected from the disccusion in the previous section, there are only small differences between the FT curves calculated with SG and VASP.

Zoom In Zoom Out Reset image size

Switching of the easy magnetization axis occurs several times as a function of band filling in all cases. Interestingly, a MAE reversal already takes place by removal or addition of one or two electrons per unit cell. Furthermore, the MAE undergoes drastic changes in magnitude, even attaining values which are one order of magnitude larger than the neutral ones, specially in the ${N}_{e}\approx 10$ region where the d-bands show the highest density of states.

With a methodological aim in mind, the study of the MAE in the non-neutral case allows us to study the validity the 2PT perturbative approach with greater confidence than in the neutral case, since now we can compare a MAE curve with a rich structure instead of just a single value. The SG 2PT curves (dashed lines in the left-hand panel of figure 4) are in very good agreement with the FT curves for FeCo, FeCu, and FePd, while large discrepancies appear when heavier atoms are present. This is particularly evident in FeAu for N_e = 5 − 10, which corresponds to the spectrum region where the d-bands of Au lie. Considering the strong dependence of the MAE on the electron band structure details discussed in the previous section and the fact that the agreement with FT is not homogeneous as N_e changes, both facts suggest that 2PT loses its validity for elements with strong SOI. Thus, the coincidence in the neutral-case FePt and FeAu MAEs could be considered to be fortuitous, in the sense that the coincidence is a consequence of the specific band structure of the alloy, as it is the case of FePt (also observed in [62]) and FeAu⁵ .

This restriction of the 2PT validity to light atoms is not a serious disadvantage. This approach facilitates the MAE evaluation with fine k-point grids and narrow smearing widths at a lower computational cost, since it requires a single DFT calculation without SOI. We recall that a weak MCA does not necessarily follow from a small ξ, as we see in the systems under study. In extended systems like the current ones, band dispersion governs the MCA.

Both FT and 2PT describe SOI with a perturbative treatment of either the charge density changes or the ξ parameter, respectively. At this point, it is important to note another fundamental difference between the FT and 2PT formalisms. In the former method, the eigenvalues change with the magnetization axis and some degeneracies may be broken. In the latter method, the unperturbed band structure is not explicitely modified. Instead, e–h pair excitations of the GS $| {{\rm{\Psi }}}^{(0)}\rangle$ are induced by H_SO, with probabilities given by their matrix elements. In other words, 2PT is a many-body approach, while FT is a one-electron approach. Thus, if we take the limit of single ions and uniaxial anisotropy, the 2PT approach described here tends to the widely-used formalism of the spin hamiltonian for non-degenerate states ${H}_{{\rm{ion}}}={{DS}}_{z}^{2}+E({S}_{x}^{2}-{S}_{y}^{2})$ [63, 64]. This magnetic anisotropy model is widely used in the study of single-atom magnetism and spin excitations [21–23].

Visual evidence of the inherent difference between FT and 2PT is presented in figure 5. This figure shows, for the case of FeCo and the VASP-Wannier calculation, the MAE densities as a function of N_e in the reciprocal space along a few high-symmetry directions inside the 1BZ. To guide the eye, the bands without SOI have been transformed from energies to the corresponding filling N_e and superimposed on the MAE density graph. As expected, for the FT case (top panel) the regions of non-zero MAE are localized close to the unperturbed bands, and take positive or negative values depending on the relative H_SO induced shifts in the eigenvalues between the OZ and OX magnetization directions. The map also reveals abrupt switching of the MAE densities close to several band-crossings. For example, this happens near the Γ point and N_e ≃ 9, where a pair of majority spin bands is split by a few meV for magnetization along OZ. Since the splitting is nearly symmetric in energy about the non-perturbed bands, the contribution to MAE is negative as the bottom band is filled and changes sign when the upper one starts to become occupied. When both bands are completely filled the MAE vanishes. In the 2PT approach the MAE density (lower panel) takes a very different aspect since the bandstructure is not modified and, as shown by figure 2, e–h pairs with an energy difference as large as a few eV can contribute to MAE at a given N_e (see also supplementary information figure 6). Thus, the map presents broad plateaus with a non-negligible MAE density in areas devoid of bands, while sign changes are always localized precisely at the bands, since they take place when the combined e–h distribution functions (term $f({\epsilon }_{{kn}\sigma })[1-f({\epsilon }_{{kn}^{\prime} \sigma ^{\prime} })]$ in equation (5)) also changes sign as the band is crossed. Nevertheless, for FeCo the two approaches yield a similar overall description of the MAE(N_e) behaviour, as seen in the k-integrated curve of figure 4 (right), in spite of treating bandstructures in a different way by construction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

With the 2PT calculation that uses MLWFs poorer agreement is found in the profile of the MAE(N_e) curves, due to the limitations of this methodology. Still, the qualitative behaviour is reproduced in all alloys except FeAu, where similar deviations as for the SG case are found (see right-hand panel of figure 4). Therefore, it could be used under suitable conditions to make predictions on the MAE behaviour as a function of doping or bias voltage at a low computational cost from a DFT calculation which does not include SOI. Those conditions are (i) weak SOI strength and (ii) resemblance between the MLWFs and Y_2m spherical harmonics. When the first condition is met, the MAE obtained in the on-site approximation (equation (10)) performs as accurately as FT. This has been checked with SG calculations (see Supplementary Information figure 8), where the drawback of point (ii) is not present. Interestingly, the effect of inter-atomic d-orbital hybridization on the MAE is well captured by this methodology via the F factors in equation (15), while the crude approximation done for the SOI matrix elements seems less important, as it can be deduced from the good agreement with the FT curves in the FeCo and FeCu panels of figure 4 (right).

Finally, we analyse the effect of low dimensionality in the non-neutral case. Figure 6 shows the MAE dependence on the number of electrons per unit cell for the four studied bilayers. As discussed above in section 4.1, the impact of SOI in the two-dimensional bandstructure is stronger, which is manifested in the large magnitude of the oscillations in the MAE curves of figure 6. Although FT is a less reliable method for the thin films, we observe the same main trend as in the bulk results of figure 4, namely, that the agreement is better for FeCo, since the SOI strength is inside the perturbative regime.