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

A second-order perturbation (2PT) approach to the spin–orbit interaction (SOI) is implemented within a density-functional theory framework. Its performance is examined by applying it to the calculation of the magnetocrystalline anisotropy energies (MAE) of benchmark systems, and its efficiency and accuracy are compared with the popular force theorem method. The case studies are tetragonal FeMe alloys (Me=Co, Cu, Pd, Pt, Au), as well as FeMe (Me=Co, Pt) bilayers with (111) and (100) symmetry, which cover a wide range of SOI strength and electronic band structures. The 2PT approach is found to provide a very accurate description for 3d and 4d metals and, moreover, this methodology is robust enough to predict easy axis switching under doping conditions. In all cases, the details of the bandstructure, including states far from the Fermi level, are responsible for the finally observed MAE value, sometimes overruling the effect of the SOI strength. From a technical point of view, it is confirmed that accuracy in the MAE calculations is subject to the accuracy of the Fermi level determination.


Introduction
In a model system of interacting magnetic moments various contributions can be identified that lead to the observation of the magnetic anisotropy, i.e. the existence of a preferential magnetization direction in the system. These terms are the classical dipole-dipole interaction, resulting in the so-called shape anisotropy, and quantum-mechanical ones with origin in the electronic spin-orbit interaction (SOI). These latter include the anisotropic exchanges, both the symmetric and antisymmetric (known as Dzyaloshinskii-Moriya interaction [1,2]), and the magnetocrystalline anisotropy (MCA). Here, this fundamental property (MCA) is used as a probe for approximate methods to treat the SOI, but our methods are potentially valid to treat other physical properties determined by SOI, like anisotropic g-tensors or anisotropic exchange interaction [3], and Elliot-Yafet theory of spin relaxation [4].
The historical development of non-volatile memories based on the property of magnetorresistivity has been closely related to the ability of balancing out those competing contributions. The key achievement was the perpendicular magnetic anisotropy in thin film heterostructures, since the out-of-plane MCA at the interface tends to dominate the in-plane shape anisotropy. Indeed, the efficient spin orientation control in the ferromagnetic (FM) electrodes of magnetic tunnelling junctions is still a technological challenge. Additionally, the use of external magnetic fields for this purpose is evolving towards voltage control of spintronic devices [5,6] and spin transfer torques induced by spin-polarized currents [7]. These advances motivate efficient, robust and accurate modelling of the physics of the MCA for different materials/interfaces under external stimuli, such as external fields or strain.
The interplay of MCA and dimensionality is considered as a promising route for the development of spintronics. Furthermore, MCA plays a central role in the magnetic properties of two-dimensional materials, since it can prevent thermal fluctuations from destroying long-range magnetic order [8]. In this context, density functional theory (DFT) calculations that include SOI constitute a powerful theoretical tool to characterize the MCA. Nonetheless, the small magnitude of the associated magnetocrystalline anisotropy energy (MAE) imposes stringent convergence to the DFT calculations at the expense of high computational demands. In practice, fully relativistic DFT calculations that include SOI are carried out by first computing self-consistently the Hamiltonian of the system including the scalar relativistic term and next adding the spin-orbit contribution. The latter may be treated self-consistently (FSC) or non-self-consistently [9,10] within the so-called force theorem (FT) [11] to obtain the MAE. FT is often considered to produce a good estimate of the FSC MAE value for every system, at least for bulk crystals, although it has been shown to be less accurate in the case of low dimensional systems [12]. In any case, this comparison between FT and SCF calculations can only be done for relatively simple systems, because the FSC tends to be computationally too demanding. Alternatively, Green's function methods that treat SOI at the level of second-order perturbation theory (2PT) have been formulated in the literature under different flavours [13][14][15][16][17]. These approaches rely on the knowledge of the spin-orbit effects on atomic orbitals (AOs) [18][19][20], which can be quasi-analytically accounted. However, to our knowledge, the versatility of those perturbative methods has not been exploited within a DFT environment. Considering the high computational expense involved in DFT calculations that include SOI under the aforementioned techniques, it would be timely to explore the 2PT route.
In this work, we address approximate methods to treat spin-orbit effects in DFT calculations and establish their efficiency, accuracy, and applicability range using MAE values of several Fe-based alloys as benchmarks. Overall, we tackle fundamental concepts related to the magnetic anisotropy, such as the relationship between the one-electron wavefunctions used in ab initio or tight binding schemes with the multiplets (spin states) used in spin hamiltonian formulations [21][22][23]. More specifically, our aim is two-fold: (i) set the technical grounds for implementing a 2PT approach in regular DFT codes, with a detailed description of the handling of two different families of localized basis sets. Thereby, the 2PT many-body expression for the MAE is evaluated using oneelectron operators and wavefunctions. (ii) Taking the FSC MAEs of the model systems as reference values, study the accuracy of the FT and 2PT approximations and their dependence on the physics of the problem. This requires to take into account careful convergence criteria.
As case studies, we choose tetragonal transition-metal alloys FeMe (Me=Co, Cu, Pd, Pt, and Au) with L1 0 structure (see figure 1) and we also consider thin FeMe (Me=Co, Pt) bilayers with (111) and (100) orientations. With this choice we can cover a wide range of two orders of magnitude in SOI strengths by using elements of different atomic weigths and also we analyse differerent hybrid bands characters (d-d for FeCo, FePd, and FePt, and d-s for FeCu and FeAu). Additionally, the MCA of this type of alloys is well characterized and, thus, provides us with a reliable quantitative reference to compare the approximate methods considered. In fact, for solids with cubic symmetry, MCA is an effect of fourth order in the SOI strength but a tetragonal distortion can enable a second-order MCA. This idea is behind the materials used in (or proposed for) some of the aforementioned devices. For example, multiferroics allow for strain-mediated magneto-electric coupling [24,25] and, more recently, tetragonal Heusler alloys have attracted attention for combining high MCA and half-metallicity Figure 1. Structure of the L1 0 tetragonal unit cell and cartesian axes that provide a good match between the MLWFs and the Y 2m orbital functions. [26,27]. The family of L1 0 alloys considered in this work is also versatile, since their structure and magnetic properties can be tailored by varying the stoichiometry, as it is the case of the strongly magnetostrictive 'galfenol' (Fe 1−x Ga x ) [28,29], Fe 1−x Co x [30][31][32], and Cu-Ni films [33]. Several theoretical studies have shown that a bias voltage could significantly affect the MCA [34,35] and, in fact, an electric-field-induced MCA switching has been realized in Fe 30 Co 70 alloy films [36].
Our first principles calculations show that the second-order perturbative (2PT) approach to treat spin orbit interaction is very reliable for the lighter (3d and 4d) transition metals but it breaks down for heavier 5d elements. In general, the MAE values calculated using the force theorem (FT) agree reasonably well with the fully self-consistent calculated values also for the heavier elements, but minor quantitative discrepancies between one and other are apparent in the two-dimensional thin bilayers. The different performance of these two approximate methods is explained by their construction: while 2PT is perturbative in the spin-orbit strength, FT is perturbative in the charge density. The general character of these approaches suggests that they can be applied to any system featuring dispersive bands. Regarding the nature of the MCA in the alloys under study, we find its physical origin in the availability of empty Fe electron states, although the whole valence bandstructure contributes to its final magnitude, this latter calculated with tenths of meV precision.
The paper is organized as follows: section 2 describes the three methods used here to compute the MAE, namely FSC (section 2.1), FT (section 2.2) and, in more detail, two different implementations of a 2PT formula in DFT codes (section 2.3). Details of the DFT calculation parameters for the Fe-based alloys are presented in section 3. The results for the charge neutral and non-neutral cases are shown in sections 4.1 and 4.2, respectively. Finally, conclusions are drawn in section 5.

Theoretical background
The SOI Hamiltonian is generally written as a sum over one-electron operators: where l i and s i are the orbital and spin momentum operators, respectively, acting on the ith electron in the system and ξ i accounts for the SOI strength. In practice, as most of the relevant electronic and magnetic properties of solids derived from SOI originate from valence electrons, only outer-shell and semi-core electrons are considered in our first principles calculations. Since ξ i is proportional to the radial derivative of the potential, it increases with the atomic number. Furthermore, it is often a good approximation to take the same value for all the electrons within the same l-shell.
In the next subsections we describe the three methods considered in this work to evaluate the MAEs from first-principles, all of them including SOI at different levels of approximation. In particular, we will examine under which conditions a second-order perturbative approach, where ξ i acts as the perturbation constant, breaks down.

Self-consistent MAE (FSC)
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, where the spins are aligned along the OX and OZ directions shown in the L1 0 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.

Non-self-consistent MAE based in the FT
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 blockdiagonal charge density matrix. Then, the spin-orbit H SO term calculated for a given magnetization axis is added to the scalar-relativistic hamiltonian H 0 and new eigenvalues are calculated by diagonalization without further self-consistent cycles. We denote the resulting total energy change E x z tot , D and the corresponding charge density change Δρ x, z . The MAE is approximated as the difference in the band energies, E x z band , , 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, Here, the sum runs over one-electron eigenvalues 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 E x tot D and E z tot D are correct to order ) , 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 x z , 2 r D ( ) -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.

Second-order perturbative MAE (2PT)
A widely used alternative approximation treats the SOI as a second-order perturbation (2PT) to the many-body GS, 0 Y ñ | ( ) . The general expression for the 2PT energy correction is where the sum runs over excited states. In a many-body language the GS 0 Y ñ | ( ) 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 i Y ñ | ( ) is then constructed by creating electron-hole (e-h) pairs using the unoccupied Kohn-Sham eigenstates. Thus, the E E i 0 -term in the denominator is simply the energy difference between the occupied and unoccupied eigenvalues associated to the particular e−h excitation [13,[18][19][20]. The perturbative expansion may then be written in terms of the GS Kohn-Sham eigenstates knsñ | as follows: where , s s¢ 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 E SO 2 D ( ) 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: where R runs over the lattice vectors, R a ñ | ( ) 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: where the k-space H SO (k) matrix elements are given by: and the generalized projected charges by: 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: where indices , a b 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 lm H lm SO a s a s á ¢ ¢ñ | | matrix elements (independent of mm¢ and ss¢). The angular part of these matrix elements, lm lm l s a s a s á ¢ ¢ñ | · | , 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: where we have defined: We observe that the 2PT equations (7) and (11) are perturbative in the SOI strength, which appears explicitly in the form of the parameters l x a in this equation.
Next, we consider the case when the unperturbed spin-polarized calculation is realized employing a planewave 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: where knsñ | stands for the Kohn-Sham eigenstates already interpolated in the dense k-grid [43] and Q k j ns ( ) 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 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 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 kgrid by interpolation:  (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 k n a s ( ) in equation (9) by the Q k j ns ( ) ones and restricting the summation over the angular momentum numbers to the l l 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.

Computational details
In all the procedures described in this section, to prevent the MAE values from being biased by the structural parameters, we have kept the bulk L1 0 lattice constants a, c fixed in all calculations (see figure 1 and table 1). The model for FeCo has the Fe bcc unit cell volume and c/a=1.2 to maximise the anisotropy, as suggested by the literature on strain effects [30]. For FeCu, we keep the Cu-Cu as in bulk fcc Cu and c/a=1.34 [45]. Finally, we use published lattice constants for FePd, FeAu [46] and FePt [47]. The interplanar distance of the FeCo and FePt bilayer models with (100) structure has been set to d c 2 = . For the (111) films, the Co and Pt bulk lattice constants have been used and the interplanar distance has been optimized (see table 2).
Two DFT codes have been used in this work: SIESTA-Green (SG) [38,48] and VASP [49,50]. The former uses multi-ζ non-orthogonal strictly localized numerical AOs as basis set and replaces core electrons by pseudopotentials, while the latter uses plane-waves and projector-augmented wave potentials to describe the ion cores [51]. Both codes feature SOI implementations [48,52] that allow to obtain the MAE values directly from equation (2) or the SOI-corrected eigenvalues of equation (3). By working with both codes we ensure that the conclusions about the MCA are not biased by the basis set type. The exchange and correlation functional used in all calculations is that of Perdew, Burke and Ernzerhof [53].
In the VASP calculations the p semi-core states are included as valence electrons. We have set the plane-wave energy cut-off to 400 eV in all the systems, except for 420 eV in FeAu, and suitable k-point Monkhorst-Packgrids [54] according to each lattice dimensions and calculation type ((24×24×20) for FeCo and (24×24×18) for the other L1 0 units cells, and (48×48×1) for the thin film calculations). The tetrahedron method with Blöchl corrections is used to obtain the Fermi level [55]. For the SG calculations we have adopted a double-ζ polarized scheme to generate the AO basis set using a confinement energy of 100 meV although, as Table 1. Unit cell lattice parameters (see figure 1) of the considered model bulk systems and calculation parameters used in the wannierization. n s is the number of s-wave-like Wannier functions, introduced in addition to the d-orbital-like ones, placed at interstitial sites of the structure. [w 0 , w 1 ] are the frozen windows used for disentanglement with respect to the Fermi energy. Two intervals are shown for FeAu that correspond to different windows for the spin-majority and spin-minority bands, respectively. k 1 is the grid used in the reference DFT calculation and k 2 the interpolated grid used to evaluate the MAE in the 2PT approximation. p min is the minimum projection factor p m a s (equation (16)) found for each system. opposed to VASP, no p semi-core states are considered. Pseudo-core corrections are included for all the atoms involved, while a very fine real space mesh is employed by setting the mesh cut-off to 4000 Ry.
In the SG case, SOI matrices are calculated under the fully-relativistic pseudo-potential (FR-PP) method described elsewhere [48]. This approach goes beyond the on-site approximation [56] (equation (10)) as it takes into account intra-and inter-atomic interactions between different l-shells. Although the off-site terms tend to be small, test calculations show that neglecting them can induce errors in the MAE of around 0.2 meV or even larger (around 1 meV) in particular cases. Nevertheless, the on-site approximation allows to extract the SOI strengths ξ αl , which we provide in table 3 for the d orbitals.
Even when working with the FSC method for SOI, the calculation parameters must be carefully tested. The Fermi energy smearing is a decisive technical factor for the MAE of some of the systems. This issue has been addressed with both codes for the FSC and FT methods. We find satisfactory convergence when the tetrahedron method is used for the FSC MAE with VASP [55]. By doing so, we avoid kT-dependence in the resulting MAE values (we have checked that the total energy extrapolation to kT=0 gives, in general, good agreement with the results of the tetrahedron method). With SG, since the use of finer k-point grids can be afforded at a reasonable computational and memory cost, a high convergence could be systematically achieved in the k-grid. Convergency values below 0.01 meV are obtained with kT values entering the Fermi-Dirac distribution function as low as 1 meV. In the FT calculations we find that the smearing function, whether Fermi-Dirac or Methfessel-Paxton of a given order [57], has a non-negligible influence, but it becomes less important for sufficiently fine kgrids and small kT. A detailed convergence analysis for all phases can be found in the Supplementary Information sections 1, 2 and figures 1-4 (stacks.iop.org/NJP/21/073054/mmedia).
For the more elaborate 2PT methodology, we have implemented equation (7) for the SG calculations and the semi-analytical form of equation (14) for VASP in combination with Wannier90 [41,58]. The k-grids needed to obtain a faithful representation of the electronic structure with MLWFs can be less dense than the ones typically needed to obtain the MAEs. The latter, also used in the explicit evaluation of equations (14) and (15), can be chosen as dense as desired by MLWF interpolation [43]. Besides, due to the strongly hybridized d-bands in the L1 0 Fe-alloys, prior to wannierization it is useful to perform a disentanglement of the bands [42] within an energy window that contains the relevant states. A numerical drawback in the whole procedure is that the quality of the wannierization is system-dependent. Therefore, for each alloy we have chosen suitable settings, shown in table 1, to meet the usual sanity requirements of a wannierization, namely, little overlap between MLWFs (at least between the functions that emulate the d-orbitals), bandstructure reproducibility, and spatial localization. The same strategy has been adopted for the four thin films (see table 2).
It is worth mentioning that the five atomic d-orbital wavefunction geometries, i.e. the representations of the angular functions Y r m 2 (ˆ) in cartesian coordinates, depend on the convention taken for the OX, OY, OZ axes in each code. Obviously, the electronic structure that results from hybridization of the atomic wavefunctions must be independent of those conventions. However, if we align the crystallographic directions and the cartesian axes of VASP and Wannier90 as shown in figure 1, we can represent the bonding states as overlapping Y r R   which range between zero and one. As shown in tables 1 and 2, despite the dispersive character of the bandstructures, in this work we find projections above 0.80 after wannierization. Table 4 contains the collection of the MAE values of the L1 0 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][47][48]59].

Neutral systems
The first important observation is that the MAE of Fe-based L1 0 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 Pd,Au x being larger than Co,Cu x (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 kconvergence 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.  (7), i.e. a SG calculation) from the atom-pair terms (panels (a)-(e)) and the created e-h spin-pair terms (panels (f)-(j)). Here 'dn' stands for spin down in the graph. Each colour corresponds to a FeMe system. 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  (7), with a Fermi-Dirac smearing of kT = 10 meV is used here. Right: VASP calculations, i.e. the dashed line is obtained from equation (14), with a Fermi-Dirac smearing of kT=50 meV. Table 4. MAE values in meV for the bulk neutral systems. Labels SG and V indicate SIESTA-Green and VASP calculations, respectively, with the cut-off energy and k-point samplings discussed in the text. In the 2PT SG (VASP) calculations a Fermi-Dirac smearing with kT=10 meV (kT=50 meV) has been used. (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 Figure 5. k-resolved MAE as a function of the number of valence electrons N e along high-symmetry directions inside the first Brillouin zone for the FeCo alloy. The top and bottom panels show the FT and 2PT approaches, respectively, obtained from equation (14) (VASP calculation) with a Fermi-Dirac smearing of kT=50 meV. The correspondence between the Wannier-interpolated energy eigenvalues without SOI and the number of valence electrons is shown as black (majority spin) and green (minority spin) bandstructures. The horizontal line indicates charge neutrality. In the colour scale, red and blue regions of the spectrum account for a contribution to anisotropy easy axis along OZ or OX, respectively. Table 5. MAE values in meV for the neutral bilayer models. Labels SG and V indicate SIESTA-Green and VASP calculations, respectively, with the paramters discussed in the text. As in the bulk structures, the 2PT SG (VASP) calculations used a Fermi-Dirac smearing with kT=10 meV (kT=50 meV). 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 0 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 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).
Therefore, we conclude that, as a general rule, the high DOS counterbalances the decay of the ) 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)   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.
Finally, we compare our perturbative results with the widely used 2PT equation by Bruno for bands of dorbital 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 E L 2 x D µ á ñ ( ) . 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 0 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.

Non-neutral systems
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 10 e » 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 5 . 5 In the FeAu alloy, the Au-d states are mostly confined in a band between −7 and −4 eV below the Fermi energy (see figure 2). On the one hand, those states are subject to strong couplings by SOI, since ξ Au =615 meV. On the other hand, because of the kn ) factors, those states have a weaker effect on the MAE for N e values close to charge neutrality (N e =19) than for smaller N e values. E.g. a band filling N e =10 corresponds to a downward shift of the Fermi level of −3.7 eV, close to the Au-d states. Therefore, fast sharp oscillations are observed in the MAE curve at N e =5−10, while smooth behaviour and apparent agreement with FT exists at N 12 e > . Importantly, this does not mean that the Au-d states have a negligible contribution, as evidenced by the absence of a plateau in the occupied states curve of the FeAu panel of figure 2, which represents the N e =19 case. For the FePt alloy, since the Pt-d band is less localized in energies, the disagreement between 2PT and FT is visible throughout the MAE(N e ) curve. 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 0 Y ñ | ( ) 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 ) [63,64]. This magnetic anisotropy model is widely used in the study of single-atom magnetism and spin excitations [21][22][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 f 1 kn kn (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.
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.

Conclusions
In this work, we have included the SOI in DFT calculations at different levels of approximation to obtain the magnetocrystalline anisotropy energies (MAE). As reference, we use fully-relativistic fully-self-consistent (FSC) calculations. We have examined the accuracy of a many-body second-order perturbation (2PT) treatment of the SOI on the scalar-relativistic GS and we have found that it is as good as the well established FT approach. In both cases (2PT and FT), we have confirmed that an accurate Fermi level determination is crucial to obtain well converged MAE values.
As case studies, we have considered several FeMe tetragonal ferromagnetic alloys with L1 0 structure, as well as two FeMe (Me=Co, Pt) bilayers with (111) and (100) symmetry. By choosing Me=Co, Cu, Pd, Pt, and Au, we cover the scenarios of s-d and d-d hybrid band effects and a range of atomic SOI strength parameters ξ. We find that the 2PT approximation describes accurately the MAEs of FeCo, FeCu, (ξ∼0.05 eV) and satisfactorily that of FePd (ξ∼0.1 eV), but fails for the alloys containing 5d metals (ξ∼0.5 eV), while FT to provides reliable MAE values in general and some minor quantitative discrepancies with respect to FSC values for the bilayers. The difference in the performance of the two approximations has the following origin: while 2PT is perturbative in ξ, FT is perturbative in the charge density changes by SOI.
The MAE in this family of alloys is determined not only by the empty spin-minority Fe states, which lie about 2 eV above the Fermi level, but also by the whole valence band, which lies several eV below the Fermi level. Thus, the MCA is determined by electronic states that lie from the Fermi level much further than the SOI strength parameter. The details of the bandstructure are, in essence, responsible for the final MAE value.
Finally, the 2PT approximation is sound enough to account for the anisotropy behaviour of the weak SOI alloys under deviations from charge neutrality. We show that magnetic switching can be induced by addition or removal of electrons. This effect could be tuned by, for example, doping or strain, to find an scenario under a minimal bias voltage. These properties have ample applications in magnetoelectric and magnetostrictive devices. The 2PT approximation to the SOI, when valid, can be used to study large numbers of these cases efficiently, as it uses as sole input a scalar-relativistic DFT calculation.