Intersublattice magnetocrystalline anisotropy using a realistic tight-binding method based on maximally localized Wannier functions

TL;DR

This paper develops a realistic tight-binding approach based on maximally localized Wannier functions to analyze intersublattice magnetocrystalline anisotropy in transition metal compounds, revealing complex intersublattice effects.

Contribution

The study introduces a tight-binding method that accurately evaluates two-ion MAE contributions, highlighting the significance of intersublattice effects beyond single-ion models.

Findings

Inter-sublattice MAE contributions vary in sign among compounds.

Inter-sublattice MAE can be significant, affecting overall magnetic anisotropy.

The method aligns with first-principles results and captures Fermi surface features.

Abstract

Using a realistic tight-binding Hamiltonian based on maximally localized Wannier functions, we investigate the two-ion magnetocrystalline anisotropy energy (MAE) in $L1_0$ transition metal compounds. MAE contributions from throughout the Brillouin zone are obtained using magnetic force theorem calculations with and without perturbation theory. The results from either method agree with one another, and both reflect features of the Fermi surface. The intra-sublattice and inter-sublattice contributions to MAE are evaluated using a Green's function method. We find that the sign of the inter-sublattice contribution varies among compounds, and that its amplitude may be significant, suggesting MAE can not be resolved accurately in a single-ion manner. The results are further validated by scaling spin-orbit-coupling strength in density functional theory. Overall, this realistic tight-binding method provide an effective approach to evaluate and analyze MAE while retaining the accuracy of corresponding first-principles methods.

Figures

Tables

Table 1. Table 1: Spin-orbit coupling constants ξ i ↑ ↓ superscript subscript 𝜉 𝑖 ↑ absent ↓ \xi_{i}^{\uparrow\downarrow} (meV) in various compounds calculated in FP-LMTO. On the 3 ​ d 3 𝑑 3d sites (Mn, Fe, and Co), ξ i ↑ ↑ / ξ i ↑ ↓ ≈ ξ i ↑ ↓ / ξ i ↓ ↓ superscript subscript 𝜉 𝑖 ↑ absent ↑ superscript subscript 𝜉 𝑖 ↑ absent ↓ superscript subscript 𝜉 𝑖 ↑ absent ↓ superscript subscript 𝜉 𝑖 ↓ absent ↓ \xi_{i}^{\uparrow\uparrow}/\xi_{i}^{\uparrow\downarrow}\approx\xi_{i}^{\uparrow\downarrow}/\xi_{i}^{\downarrow\downarrow} varies between 1.07 and 1.13 for the d 𝑑 d channel. SOC energies (isotropic) E so superscript 𝐸 so E^{\text{so}} (meV) calculated in TB are also listed.

Compound	1st element		2nd element		$E^{\text{so}}$
	${\xi }_p$	${\xi }_d$	${\xi }_p$	${\xi }_d$
FePt	197.2	55.0	2626.4	574.9	193.8
CoPt	190.5	72.1	2793.1	580.8	210.6
FePd	168.9	54.7	898.0	200.8	22.4
FeNi	232.9	55.9	218.2	91.6	11.7
MnGa	193.1	41.0	209.6	84.1	5.3
MnAl	203.0	41.3	31.3	0.4	2.7

Table 2. Table 2: MAE ( μ 𝜇 \mu eV/f.u.) in L ​ 1 0 𝐿 subscript 1 0 L1_{0} systems calculated in our tight-binding framework with (denoted as TB-PT) and without (denoted as TB) perturbation approach. MAE values calculated using FP-LMTO with both PBE (denoted as FP) and BH (denoted as FP-LDA) exchange-correlation functionals, and using vasp with PBE functional, are also listed for comparison. The SOC constants used in TB were obtained from FP using PBE functionals.

Compound	TB	FP	TB-PT	VASP	FP-LDA
FePt	2495	2556	2692	2656	2746
CoPt	836	788	1098	858	1156
FePd	175	164	186	194	173
FeNi	66	68	65	80	68
MnGa	399	381	415	431	421

Table 3. Table 3: Sublattice-resolved MAE (normalized to 1) in L ​ 1 0 𝐿 subscript 1 0 L1_{0} systems calculated using perturbation theory in TB. For each compound A ​ B 𝐴 𝐵 AB , MAE is resolved into intrasublattice contributions, K A ​ - ​ A subscript 𝐾 𝐴 - 𝐴 K_{A\text{-}A} and K B ​ - ​ B subscript 𝐾 𝐵 - 𝐵 K_{B\text{-}B} , and intersublattice contribution K A ​ - ​ B subscript 𝐾 𝐴 - 𝐵 K_{A\text{-}B} . We also define the contribution from individual sublattice A 𝐴 A as K A = ( K A ​ - ​ A + K A ​ - ​ B / 2 ) subscript 𝐾 𝐴 subscript 𝐾 𝐴 - 𝐴 subscript 𝐾 𝐴 - 𝐵 2 K_{A}=(K_{A\text{-}A}+K_{A\text{-}B}/2) . The hypothetical FeCo structure with c / a = 1.1 𝑐 𝑎 1.1 c/a=1.1 is also included.

$AB$	$K_{A\text{-}A}$	$K_{A\text{-}B}$	$K_{B\text{-}B}$	$K_A$	$K_B$
FePt	0.11	-0.55	1.44	-0.16	1.16
CoPt	0.17	-0.77	1.59	-0.21	1.21
FePd	1.80	-2.58	1.78	0.51	0.49
FeNi	4.51	-5.88	2.37	1.57	-0.57
MnGa	0.66	0.24	0.10	0.78	0.22
MnAl	0.98	0.02	0.00	0.99	0.01
FeCo	0.10	0.63	0.28	0.41	0.59

Table 4. Table 4: Site-and-orbital-resolved charge q 𝑞 q and magnetic moment m 𝑚 m ( μ B subscript 𝜇 B \mu_{\text{B}} /atom) in L ​ 1 0 𝐿 subscript 1 0 L1_{0} compounds calculated in TB and vasp . The moments of s 𝑠 s -orbitals are negligible and not shown. Spin-orbit is not included in calculation.

Compound	Method	$1_{st}$ element							$2_{nd}$ element							Total
		$q_s$	$q_p$	$m_p$	$q_d$	$m_d$	$q_1$	$m_1$	$q_s$	$q_p$	$m_p$	$q_d$	$m_d$	$q_2$	$m_2$	$q$	$m$
FePt	VASP	0.41	0.45	-0.01	6.12	2.92	6.98	2.92	0.61	0.50	-0.05	7.73	0.43	8.84	0.36	15.82	3.28
	TB	0.83	0.81	-0.04	6.38	2.95	8.03	2.95	0.98	0.81	-0.09	8.18	0.45	9.97	0.34	18.00	3.29
CoPt	VASP	0.41	0.44	-0.01	7.20	1.91	8.05	1.89	0.61	0.50	-0.04	7.72	0.46	8.82	0.41	16.88	2.30
	TB	0.83	0.80	-0.04	7.48	1.92	9.11	1.87	0.96	0.78	-0.06	8.14	0.47	9.89	0.40	19.00	2.26
FePd	VASP	0.42	0.41	-0.01	6.09	2.97	6.92	2.96	0.37	0.31	-0.04	7.96	0.42	8.63	0.36	15.55	3.33
	TB	0.81	0.67	-0.04	6.35	3.05	7.83	3.00	0.85	0.70	-0.09	8.63	0.42	10.17	0.29	18.00	3.29
FeNi	VASP	0.47	0.49	-0.02	6.17	2.69	7.12	2.67	0.51	0.51	-0.07	8.31	0.72	9.33	0.62	16.45	3.29
	TB	0.80	0.70	-0.05	6.38	2.74	7.89	2.67	0.85	0.73	-0.09	8.53	0.71	10.11	0.57	18.00	3.24
MnGa	VASP	0.29	0.31	0.02	5.13	2.49	5.73	2.52	1.08	1.31	-0.11	0.16	0.02	2.56	-0.14	8.29	2.39
	TB	0.55	0.47	0.02	5.45	2.60	6.46	2.64	1.36	1.89	-0.15	0.28	0.03	3.54	-0.15	10.00	2.49
MnAl	VASP	0.33	0.34	0.02	5.12	2.37	5.78	2.41	0.40	0.42	-0.04	0.00	0.00	0.82	-0.06	6.60	2.35
	TB	0.56	0.47	0.03	5.49	2.49	6.52	2.55	1.21	1.87	-0.14	0.41	0.03	3.48	-0.17	10.00	2.39

Equations

$$V_{\text{so}}=\xi \mathbf{L}\cdot \mathbf{S}=\frac{{\hslash }^2}{2M^2{c}^2}\frac{1}{r}\frac{\mathrm{d}V}{\mathrm{d}r}\mathbf{L}\cdot \mathbf{S},$$

$$K=\sum _{k,b}\left({ϵ}_{k,b}^{110}{f}_{k,b}^{110}-{ϵ}_{k,b}^{001}{f}_{k,b}^{001}\right),$$

$$K=\sum _{k,b}\left(({ϵ}_{k,b}^{110}-{ϵ}_{\text{F}}^0)f_{k,b}^{110}-({ϵ}_{k,b}^{001}-{ϵ}_{\text{F}}^0)f_{k,b}^{001}\right),$$

$$E^{\text{so}}=-\frac{1}{2\pi }\mathrm{Im}{\int }_{-\infty }^{{ϵ}_F}\mathrm{d}ϵ\mathrm{Tr}\left[\tilde{G}(ϵ)V_{\text{so}}\right],$$

$$E^{\text{so}}$$

$$G(ϵ)={\left(ϵ-H\right)}^{-1}.$$

$$E^{\text{so}}(\hat{\text{n}})=\sum _{il;j{l}^{\prime },\sigma {\sigma }^{\prime }}{E}_{il,j{l}^{\prime }}^{\sigma ,{\sigma }^{\prime }}(\hat{\text{n}}),$$

$$E_{il,j{l}^{\prime }}^{\sigma {\sigma }^{\prime }}(\hat{\text{n}})=-\frac{1}{2\pi }\mathrm{Im}{\int }_{-\infty }^{{ϵ}_{\text{F}}}\mathrm{d}ϵ\int \mathrm{d}\mathbf{k}\mathrm{Tr}{E}_{il,j{l}^{\prime }}^{\sigma {\sigma }^{\prime }}(\mathbf{k},ϵ;\hat{\text{n}}),$$

$$E_{il,j{l}^{\prime }}^{\sigma {\sigma }^{\prime }}(\mathbf{k},ϵ;\hat{\text{n}})=G_{il,j{l}^{\prime }}^{\sigma }(\mathbf{k},ϵ)V_{j{l}^{\prime }}^{\sigma {\sigma }^{\prime }}(\hat{\text{n}})G_{j{l}^{\prime },il}^{{\sigma }^{\prime }}(\mathbf{k},ϵ)V_{il}^{{\sigma }^{\prime }\sigma }(\hat{\text{n}}).$$

$$K=\sum _{ij,\sigma {\sigma }^{\prime }}{K}_{ij}^{\sigma {\sigma }^{\prime }}=\sum _{ij,\sigma {\sigma }^{\prime }}{E}_{ij}^{\sigma {\sigma }^{\prime }}({\hat{\text{n}}}_{110})-E_{ij}^{\sigma {\sigma }^{\prime }}({\hat{\text{n}}}_{001}).$$

$$\begin{array}{rl}2U& =V({\hat{\text{n}}}_{110})+V({\hat{\text{n}}}_{001}),\\ 2A& =V({\hat{\text{n}}}_{110})-V({\hat{\text{n}}}_{001}).\end{array}$$

$$K=\sum _{ij,\sigma {\sigma }^{\prime }}{\tilde{K}}_{ij}^{\sigma {\sigma }^{\prime }}$$

$${\tilde{K}}_{ij}^{\sigma {\sigma }^{\prime }}=$$

$$\mathrm{Tr}\left[G_{ij}^{\sigma }(\mathbf{k},ϵ)U_j^{\sigma {\sigma }^{\prime }}{G}_{ji}^{{\sigma }^{\prime }}(\mathbf{k},ϵ)A_i^{{\sigma }^{\prime }\sigma }\right].$$

$${\tilde{K}}_{ij}^{\uparrow \downarrow }+{\tilde{K}}_{ij}^{\downarrow \uparrow }=K_{ij}^{\uparrow \downarrow }+K_{ij}^{\downarrow \uparrow }.$$

$$V_{\text{so}}(\mathbf{\lambda })=\sum _i{\lambda }_i{\xi }_i\mathbf{L}\cdot \mathbf{S},$$

$$K(\mathbf{\lambda })$$

$${\alpha }_{ij}=K_{ij}=\sum _{\sigma {\sigma }^{\prime }}{K}_{ij}^{\sigma {\sigma }^{\prime }}=\sum _{\sigma {\sigma }^{\prime }}{\tilde{K}}_{ij}^{\sigma {\sigma }^{\prime }}.$$

$$\mathbf{L}\cdot \mathbf{S}=\frac{1}{2}({\mathbf{J}}^2-{\mathbf{L}}^2-{\mathbf{S}}^2)=\frac{{\hslash }^2}{2}\left(\begin{array}{cc}{L}_z& L_{-}\\ L_{+}& -L_z\end{array}\right).$$

$$⟨l,m|L_z|l,m⟩$$

$$⟨l,m-1|L_{-}|l,m⟩$$

$$⟨l,m+1|L_{+}|l,m⟩$$

$$Y_{\ell m}$$

$${\mathbf{O}}^{(\mathbb{R})}$$

$${\bf U}_{\mathbb{C}\leftarrow\mathbb{R}}=\frac{1}{\sqrt{2}}\left(\begin{array}[c]{ccccccc}i&0&0&0&0&0&1\\ 0&i&0&0&0&1&0\\ 0&0&i&0&1&0&0\\ 0&0&0&\sqrt{2}&0&0&0\\ 0&0&i&0&-1&0&0\\ 0&-i&0&0&0&1&0\\ i&0&0&0&0&0&-1\\ \end{array}\right).$$

$$H_{\text{so}}(\hat{\mathbf{n}})=\frac{\xi }{2}U(\theta ,\phi )(\mathbf{L}\cdot \mathbf{S})U^{†}(\theta ,\phi ),$$

$$U\left(\theta,\varphi\right)=\left(\begin{array}[]{cc}e^{i\frac{\phi}{2}}\cos\left(\frac{\theta}{2}\right)&e^{-i\frac{\phi}{2}}\sin\left(\frac{\theta}{2}\right)\\ -e^{i\frac{\phi}{2}}\sin\left(\frac{\theta}{2}\right)&e^{-i\frac{\phi}{2}}\cos\left(\frac{\theta}{2}\right)\\ \end{array}\right),$$

$$\begin{gathered}H_{\text{so}}({\bm{\hat{\textnormal{\bf n}}}}_{110})\\ =\frac{\xi}{2}\left(\begin{array}[]{cc}\frac{\sqrt{2}}{2}\left({L_{x}+L_{y}}\right)&-L_{z}+\frac{i\sqrt{2}}{2}\left(L_{x}-L_{y}\right)\\ -L_{z}-\frac{i\sqrt{2}}{2}\left(L_{x}-L_{y}\right)&-\frac{\sqrt{2}}{2}\left({L_{x}+L_{y}}\right)\\ \end{array}\right).\end{gathered}$$

Full text

I Introduction

Magnetocrystalline anisotropy (MA) arises from the interplay between spin-orbit coupling (SOC) and crystal-field effects and is one of the most fundamental intrinsic magnetic properties van Vleck (1937). Materials with high MA have been used in many applications including permanent magnets McCallum et al. (2014) and magnetic recording media. One of the latest examples is the recently realized magnetic two-dimensional (2D) van der Waals (vdW) class of materials, in which MA is required to stabilize the long-range magnetic ordering down to atomically thin dimensions. These materials can be exploited as platforms for true 2D magnetism and for innovative applications such as energy-efficient, ultra-compact, spin-based electronics. In general, it is of great interest to evaluate and resolve MA, and to unravel the underlying mechanisms in a given system. Ultimately, such understanding will guide the control and tuning of MA, accelerating the development of new materials and their applications Rau et al. (2014).

Density functional theory (DFT) has proven to be a valuable tool to investigate and predict MA energy (MAE) in various systems. The relativistic effects of valence electrons are often treated using various approximations to reduce computational complexity and cost. Instead of directly solving the four-component Dirac equation self-consistently, one usually treats SOC as perturbation and starts first with the two-component scalar-relativistic (SR) Hamiltonian Koelling and Harmon (1977), omitting SOC but including all other relativistic effects such as mass-velocity and Darwin terms. SOC can be added directly into the SR Hamiltonian or included in a subsequent step using the basis (often a subset of it) of SR wave functions (second variation) Li et al. (1990); Shick et al. (1997). Because the charge- and spin-density variations caused by SOC vanish to first order in the SOC strength Wang et al. (1996a), the magnetic force theorem (MFT) Mackintosh and Andersen (1980); Weinert et al. (1985); Daalderop et al. (1990) is often applied on top of the second-variation method to calculate MAE as the difference of one-electron band energies. Finally, perturbation theory (PT) Yosida et al. (1965); Abate and Asdente (1965); Takayama et al. (1976); Bruno (1989); Cinal et al. (1994); van der Laan (1998); Ke and van Schilfgaarde (2015); Jesche et al. (2015) is also often used with the MFT to compute and analyze the change of band energies due to SOC. Overall, depending on the system size and approximations used, MAE computation can be quite demanding, because of the enlarged dimension of Hamiltonian and the reduced symmetry due to SOC, and the denser $k$ mesh needed for high accuracy.

Empirical or semi-realistic tight-binding (TB) methods Bruno (1989) were widely used to study MAE long before MAE became accessible to more sophisticated DFT methods Daalderop et al. (1990). Pioneering work Bruno (1989); Cinal et al. (1994) using TB provided a fundamental understanding of MAE in various systems. However, empirical TB Hamiltonians are generally hard to parametrize and often have insufficient accuracy to describe band structure, usually limiting TB to obtaining qualitative results in systems with large MAE. The recently developed maximally localized Wannier functions (MLWFs) method Marzari and Vanderbilt (1997); Marzari et al. (2012); Souza et al. (2001) has been widely used to effectively construct TB Hamiltonians to compute many properties such as Fermi surfaces, Berry curvature, and transport. With a smaller basis, it can describe an isolated set of bands and/or entangled bands in a given energy window. Indeed, this method is also very suitable for MAE calculations, considering that MAE is, after all, a ground-state quantity, determined by the occupied states. However, due to the minuteness of MAE, it is not clear how accurately the realistic TB can describe MAE in systems such as transition metal bulk compounds. Here, we demonstrate that the realistic TB framework based on the MLWFs method can produce accurate MAE in comparison to DFT, thereby providing an efficient framework to compute and analyze MAE.

Various decomposition schemes have been used to resolve MAE into $k$ space, atomic sites, orbital, and spin channels, providing insight and guidance on tuning MAE. The MFT enables resolution of MAE into individual bands on each $k$ point in reciprocal space, allowing a band-structure-origin analysis of MAE. Site-resolved MAE is often calculated using methods such as evaluation of the on-site SOC energy Antropov et al. (2014), second-order PT, scaling the SOC strength, or others Wang et al. (1996b); Åberg et al. (2015). For example, $R_{2}$Co17 (with $R$ = Y or Ce) compounds have very small MAE and ab initio analysis found that one particular Co sublattice, the so-called dumbbell sites, has a large negative contribution to uniaxial anisotropy Ke et al. (2016). Thus, proper substituents that preferentially occupy the dumbbell sites and eliminate the negative contributions will significantly improve the uniaxial anisotropy, as observed in experiments. Another interesting example is the oscillating MAE behavior in diluted nitridometalates Li2[(Li${}_{1-x}T_{x}$)N], with $T$ = Mn, Fe, Co, or Ni. The MAE can be solely attributed to individual $T$ atoms and the band-filling effect on MAE can be quantitatively described in a single-ion MAE model Ke and van Schilfgaarde (2015); Jesche et al. (2015).

Site-resolved MAE values, together with exchange parameters, can also be used as inputs for subsequent large-scale atomic spin simulations to calculate temperature-dependent magnetic properties. When interfacing ab initio results with atomic spin simulation, MAE is usually treated in the single-ion model in the atomic spin Hamiltonian. However, it has been found that a two-ion anisotropy model is needed to properly describe properties such as temperature-dependent MAE in some systems Mryasov et al. (2005). Thus, it is of great interest to resolve MAE into sites and, in particular, beyond the single-ion model.

In this work, we investigate two-ion MAE in various $L1_{0}$ compounds, one of the most widely studied systems, with MAE values ranging from several tens of $\mu\mathrm{e}\mathrm{V}$ to a couple of $\mathrm{meV}$ per formula unit Daalderop et al. (1991); Solovyev et al. (1995); Kota and Sakuma (2014). Our approach is based on second-order PT using a Green’s function method implemented within the realistic TB framework. We demonstrate that this approach achieves accuracy similar to DFT and provides a highly efficient means to compute and analyze the two-ion anisotropy in transition metal systems.

II Computational Details

We first construct the real-space scalar-relativistic TB Hamiltonian using the MLWFs method. The corresponding Green’s function is also constructed for use in the PT approach. MAE is calculated using the MFT in TB, with and without PT, referred to hereafter as PT and MFT, respectively, for simplicity. DFT methods, including both vasp Kresse and Hafner (1993); Kresse and Furthmüller (1996) and an all-electron full-potential LMTO (FP-LMTO) method Methfessel et al. (2000), are used to calculate MAE and compare with TB. To compare with the two-ion MAE values obtained using PT in TB, we also calculate intersublattice MAE contribution by scaling the SOC strength in vasp. All DFT calculations are carried out within the generalized gradient approximation (GGA) using the functional of Perdew, Burke, and Ernzerhof (PBE) Perdew et al. (1996) unless local density approximation (LDA) von Barth and Hedin (1972) is specified.

II.1 TB Hamiltonian and SOC

The MLWFs are constructed through a postprocessing procedure Marzari and Vanderbilt (1997); Souza et al. (2001); Marzari et al. (2012) using the output of a self-consistent scalar-relativistic vasp calculation. For each $L1_{0}$ compound, 18 MLWFs corresponding to $s$-, $p$-, and $d$-type orbitals for each of the two atoms in the unit cell were generated using wannier90 Mostofi et al. (2014). The spread functional for entangled energy bands is minimized by a two-step procedure Souza et al. (2001). An outside energy window with a larger number of bands was selected to ensure good description of the band structure of the “frozen” inner energy window, spanning from the bottom of the valence band to a few eV above the Fermi level. A real-space Hamiltonian $H({\bf R})$ with dimensions 18$\times$18 is constructed to accurately represent the band structures in this specified “frozen” energy window. Then, $H({\bf k})$ is obtained by Fourier transformation. The energy bands are recalculated within TB to ensure that DFT bands can be accurately reproduced before further MAE calculations.

The TB Hamiltonian is represented in a basis of orthonormalized atomic functions $|i,l,m,\sigma\rangle$, where $i$ labels atomic sites, $l,m$ angular and magnetic quantum numbers (in cubic harmonics), and $\sigma$ the spin. The SOC part of the Hamiltonian, which can be directly added into $H$ or included using PT, can be written as

[TABLE]

where ${\bf L}\cdot{\bf S}$ depends explicitly on the direction of spin quantization axis (details can be found in Appendix A). The radial part of $V_{\text{so}}$, the SOC constants $\xi_{i,l}^{\sigma\sigma^{\prime}}$, are calculated using FP-LMTO. For simplicity, we ignore the energy and spin dependence of $\xi$. Furthermore, the occupation numbers of the Pt-$p$ orbitals, which have a large SOC constant, are overestimated in TB in comparison to DFT. We renormalize the Pt-$p$ occupation numbers based on the DFT value when adding SOC into the Hamiltonian.

II.2 MAE

Turning on SOC lowers the system energy. Here, we refer to this energy change as the SOC energy $E^{\text{so}}$, which depends on the spin direction. For uniaxial geometry, MAE can be defined as $K=E^{\text{so}}_{110}-E^{\text{so}}_{001}$, with $E^{\text{so}}_{110}$ and $E^{\text{so}}_{001}$ indicating the SOC energies along the spin directions $[110]$ and $[001]$, respectively. We use $[110]$ as the reference direction for the basal plane. A positive $K$ value indicates the system has uniaxial anisotropy with the easiest spin direction being out of plane.

II.2.1 Magnetic force theorem

After SOC is added into the TB Hamiltonian, tetrahedron integration with Blöchl correction Blöchl et al. (1994) is used to determine the Fermi level $\epsilon_{\text{F}}$, band weights, and band sums. MAE is calculated as

[TABLE]

where, $b$ is the band index, $k$ refers to the wave vector in the first Brillouin zone (BZ), and $f_{k,b}$ is the corresponding band occupancy. In order to resolve the MAE into $k$ space, Eq. (2) needs to be modified. A grand-canonical ensemble version Barreteau et al. (2016) is used:

[TABLE]

where $\epsilon^{0}_{\text{F}}$ is the Fermi level calculated without SOC. The total MAE value does not depend on the reference energy because the total number of valence electrons $\sum f_{k,b}$ is conserved when spin is along different directions; however, the $k$-resolved MAE does depend on the choice of reference energy Subkow and Fähnle (2009); Åberg et al. (2015). We use Eq. (3) to properly decompose MAE because $k$ resolution of MAE based on Eq. (2) will only reflect the change of band occupancy $\Delta f_{k,i}$. The Eq. (2) result will be dominated by the relaxation effect of the Fermi surface due to turning on SOC with magnetization along various directions, and MAE contributions will only be significant near the Fermi surface.

II.2.2 Perturbation theory

Using second-order PT, we can express orbital moment, SOC energy, and their anisotropies in terms of the susceptibility Ke and van Schilfgaarde (2015) calculated using the unperturbed band structure. The SOC energy $E^{\text{so}}$ due to the spin-orbit interaction $V_{\text{so}}$ can be written as

[TABLE]

where $\widetilde{G}(\epsilon)$, the full Green’s function, includes SOC and can be constructed from the non-perturbed Green’s function ${G}(\epsilon)$. Using second-order PT (here we consider only systems with a uniaxial geometry), the SOC energy can be written as

[TABLE]

where the Green’s function is constructed using

[TABLE]

A complex contour integration on an elliptical path Zeller et al. (1982) is used for the integration. By exploiting the fact that $G$ is spin diagonal and $V$ is ($i,l$) diagonal, the SOC energy can be written as

[TABLE]

where $E_{il,jl^{\prime}}^{\sigma,\sigma^{\prime}}({\bm{\hat{\textnormal{\bf n}}}})$ is the contribution from the sublattice-orbital-spin pair ($il\sigma,jl^{\prime}\sigma^{\prime}$). We have

[TABLE]

and

[TABLE]

$V=\xi_{i,l}^{\sigma\sigma^{\prime}}(\epsilon)\mathbf{L}{\cdot}\mathbf{S}({\bm{\hat{\textnormal{\bf n}}}})$ couples states within the same $l$ channel at the same site. Here, for simplicity, we treat SOC strength as a constant $\xi_{il}$ for each $l$ channel at site $i$, ignoring its energy and spin dependence. $\mathbf{L}{\cdot}\mathbf{S}$ can be written as a function of magnetization direction ${\bm{\hat{\textnormal{\bf n}}}}$. The MAE can be written as

[TABLE]

We define the isotropic and anisotropic parts of $V({\bm{\hat{\textnormal{\bf n}}}})$ as $U$ and $A$, respectively:

[TABLE]

Then, Eq. (10) can also be written as

[TABLE]

with

[TABLE]

Here, we have $\widetilde{K}_{ij}^{\sigma\sigma^{\prime}}={K}_{ij}^{\sigma\sigma^{\prime}}$ when $\sigma=\sigma^{\prime}$. Unlike ${K}_{ij}^{\uparrow\downarrow}={K}_{ij}^{\downarrow\uparrow}$, we have $\widetilde{K}_{ij}^{\uparrow\downarrow}\neq\widetilde{K}_{ij}^{\downarrow\uparrow}$, however,

[TABLE]

According to Eq. (13), the strength of the intersublattice MAE contribution $K_{ij}$ depends on the SOC strength of both sublattices $\xi_{i}$ and $\xi_{j}$ (contained in $V_{\text{so}}$ or $U$ and $A$) and on the intersublattice Green’s function $G_{ij}$. The element types of the sublattices determine the SOC strength $\xi$ while $G_{ij}$ is relevant to the hopping or hybridization between two sublattices and depends on the detail of electronic structure. If there is no hybridization between the two sublattices or they are coupled only through $s$ orbitals (i.e., the angular parts of $U$ and $A$ vanish) elements, the corresponding intersublattice contribution becomes negligible.

II.2.3 Scaling SOC strength

Instead of using Eq. (13), the intrasublattice and intersublattice MAE contributions can also be obtained by scaling the SOC strength $\xi_{i}$ on each site $i$ by a factor $\lambda_{i}$, and by fitting the MAE as a function of scaling vector $\bm{\lambda}=[\lambda_{1},\lambda_{2},...,\lambda_{n}]^{\intercal}$:

[TABLE]

[TABLE]

Comparing Eq. (16) to Eqs. (10), (12), and (13), the coefficients $\alpha_{ij}$ are nothing but the corresponding terms containing $\xi_{i}\xi_{j}$ in $K$. Thus, we have

[TABLE]

The SOC-scaling procedure is often used in DFT, probably due to its rather straightforward implementation. Obviously, the scaling procedure can be generalized from sites to orbitals to obtain contributions from individual orbitals.

II.3 Crystal structure of $L1_{0}$ compounds

We chose to focus on $L1_{0}$ magnetic compounds because they are one of the most widely studied systems Daalderop et al. (1991); Solovyev et al. (1995); Kota and Sakuma (2014). Their simple CuAu-type crystal structure is shown in Fig. 1. The primitive cell is body-centered tetragonal (bct) and contains one formula unit (f.u.) while the conventional cell is face-centered tetragonal (fct) and contains two f.u. For all $L1_{0}$ magnetic compounds that we study in this work, experimental lattice parameters have been used. The $c/a$ ratio values (with respect to the bct primitive cell) are in the range of 1.28–1.414. Considering that hypothetical bct-FeCo structures with different $c/a$ ratios have received significant attention in the past few years, we also investigated a hypothetical bct-FeCo structure with a $c/a$ ratio of 1.1.

III Results

III.1 Electronic structure and magnetic moment

For all compounds we investigated in this work, the scalar-relativistic band structures recalculated in TB are essentially in perfect agreement with those obtained from DFT, within the specified energy window. In both TB and DFT, the tetrahedron integration with Blöchl correction method Blöchl et al. (1994) is used to determine the Fermi level. Here, we use FePt as an example to illustrate our implementation of SOC in the realistic TB framework. The new Fermi level obtained in TB with charge neutrality deviates from the original DFT Fermi level by less than $0.01\text{\,}\mathrm{eV}$. As shown in Fig. 2, without SOC, the band structures are nearly identical between TB bands and the all-electron full potential bands calculated using FP-LMTO. The vasp bands (not shown) are essentially the same also. The SOC constants $\xi_{il}$ calculated using FP-LMTO are used as input parameters to construct relativistic TB bands.

Directly using the SOC constants calculated in FP-LMTO, the resulting band structure’s SOC splitting is overestimated just above and below the Fermi level at $R$ and $M$ points, respectively. By investigating the eigenvectors of the non-SOC band in FP-LMTO and TB, we found that TB overestimates the eigenvector component of the Pt-$p$ states, which have a very large SOC constant of about $2.6\text{\,}\mathrm{e}\mathrm{V}$. (The site-and-orbital-resolved charges and moments are listed in Appendix B.) This is because $p$ orbitals are much more extended than $d$ orbitals and we do not include the interstitial sites into the projection, and as a result, the interstitial components also fold into the atomic sites. Thus, we simply renormalize the occupation numbers (or, equivalently, renormalize the SOC constants) of the Pt-$p$ channel using the ratio of atomic $p$ charge between vasp and TB. This adjustment improves agreement between DFT bands and TB bands. Correspondingly, the MAE also improves as we discuss later. In comparison to a previously reported empirical TB method Zemen et al. (2014) that used Slater-Koster parametrization, our TB method shows significantly better agreement with DFT.

SOC parameters $\xi_{il}^{\uparrow\downarrow}$ calculated in FP-LMTO in various compounds are summarized in Table 1. Generally, $\xi$ integration is only significant near a nucleus where electrons move fast. As a result, for a given element in various compounds, $\xi_{l}$ barely changes, especially for the $d$ channel, enabling transferability. The $d$ orbitals of $3d$ elements are more spin polarized than those of $5d$ elements, and the ratio of $\xi_{i}^{\uparrow\uparrow}/\xi_{i}^{\uparrow\downarrow}\approx\xi_{i}^{\uparrow\downarrow}/\xi_{i}^{\downarrow\downarrow}$ of $3d$ elements varies between 1.07 and 1.13 in various $L1_{0}$ compounds. For simplicity, we use the value of $\xi_{il}^{\uparrow\downarrow}$ for all spin channels.

The SOC energy $E^{\text{so}}$ is proportional to $\xi^{2}$ within second-order PT. Among the $L1_{0}$ compounds we studied, CoPt has the largest $\xi$ values as well as the largest $E^{\text{so}}$. The $E^{\text{so}}$ values are much larger than MAE, indicating that the isotropic part of $E^{\text{so}}$ is much larger than the anisotropic part. In comparison to other compounds, MnGa and MnAl have a rather high anisotropic/isotropic ratio, suggesting that they have a more “efficient” band structure, in the sense that the Fermi level is close to the bandfilling position that gives the largest MAE value, as we show in subsection III.2. On the other hand, it is challenging to analyze the relationship between MAE and band structure for FePt and CoPt as each MAE value is only a small fraction of $E^{\text{so}}$.

III.2 MAE and band-filling effect

As shown in Table 2, the MAE values calculated using TB agree well with DFT, especially with those obtained by the all-electron FP-LMTO method. They are also comparable to previous DFT calculations using various methods Kota and Sakuma (2014); Antropov et al. (2014). We found that MAE values calculated using PT generally agree very well with the MFT results. The largest difference in MAE is in CoPt, in which PT gives a MAE 20% larger than the MFT result. Interestingly, the MAE of CoPt also strongly depends on the exchange-correlation functionals used. LDA increases the MAE values by 30%. Likely, this is because CoPt has a large SOC and its MAE depends on the detailed, subtle band features near the Fermi level. We will study CoPt in more detail in a later section.

Figure 3 shows the band-filling dependence of MAE in FePt, CoPt, FePd, and FeNi. The oscillation Heine et al. (1984) of MAE can be associated with the local susceptibility Ke and van Schilfgaarde (2015); Inoue et al. (2015). For the entire bandfilling range, MAE values calculated using DFT and TB agree well. On the other hand, although PT gives a good description of MAE at the actual electron filling for each compound, it disagrees with the MFT result in certain bandfilling ranges in FePt and CoPt. Specifically, the disagreement is pronounced from four to eight electrons in FePt and from four to ten in CoPt. A similar finding was reported by a previous LMTO-ASA study Kota and Sakuma (2014). In comparing results for $3d$, $4d$, and $5d$ compounds, second-order PT is best suited for the lower SOC strength found in the lighter compounds, and that’s where we obtained the best agreement between MFT and PT. In general, we found that we could not improve agreement by using a denser $k$ mesh.

Although MnGa does not contain heavier $4d$ and $5d$ elements, its MAE is larger than FePd. As shown in Fig. 4, the Fermi level is located close to the filling with maximum MAE. PT generally agrees well with MFT. We also consider the renormalization of the $p$ (for both Ga and Mn) charge when including SOC, and it slightly decreases the MAE. All TB results are within $\pm 10\%$ of the FP-LMTO results. While MnGa MAE is about $10\text{\,}\mathrm{\char 37\relax}$ of the SOC energy, the ratio is much higher than in FePt and CoPt.

III.3 Reciprocal-space resolved MAE and its correlation with Fermi surface

As shown in Table 2, MAE values calculated using second-order TB with either MFT or PT generally agree well with each other and with the corresponding DFT calculations. To further validate the applicability of PT, we compare the $k$-resolved MAE using both methods for a more stringent test. Remarkably, for all $L1_{0}$ compounds we study, the two methods produce very similar results, further suggesting that the overall effect of Fermi surface relaxation is non-significant and that second-order PT is valid. CoPt has the largest difference of MAE values between MFT and PT among the compounds we study. Figure 5 shows the $k_{z}$ dependence of MAE contributions calculated in FePd and CoPt. Although the $k_{z}$-resolved MAE in CoPt shows a larger difference, the two methods still generally agree with each other. The larger difference is likely due to its rather large SOC for the Pt atom and more complex Fermi surface, which results in a larger Fermi surface relaxation effect. Similar to FePd, the $k_{z}$-resolved MAE calculated using MFT and PT in other $L1_{0}$ compounds (not shown) are also nearly identical. And, the $Z$-$R$-$A$ ($k_{z}=0.5$) plane has the largest positive contribution to MAE.

The $(k_{x},k_{y})$-resolved MAE at $k_{z}=0$ in CoPt calculated using both MFT and PT are shown in Fig. 6, which also includes the corresponding non-SOC Fermi contour. Figure 6(a) shows the resolved MAE calculated using MFT. Obviously, the four-fold symmetry is broken due to SOC. The twofold or mirror symmetry along the [110] direction is simply an artifact of our choice of [110] as the reference in-plane spin direction. The Fermi contour plot, as shown in Fig. 6(b), is calculated by integrating the electron density at the Fermi level with a width of $0.02\text{\,}\mathrm{eV}$. To better compare with the Fermi surface and $k$-resolved MAE, we also symmetrized $k$-space contributions with symmetry operations that are compatible with SOC when the spin quantization axis is along the $z$ direction. The symmetrized $k$-resolved MAE calculated in MFT and PT are shown in Figs. 6(c) and 6(d), respectively.

As shown in Figs. 6(b) and 6(d), the correlation between the $(k_{x},k_{y})$-resolved MAE and the Fermi contour is apparent and two features stand out. First, large changes of MAE contributions occur at the Fermi contour. This is because MAE contributions are obtained by integrating Eq. (13) up to $\epsilon_{\text{F}}$. When the $k$ path crosses the Fermi contour, some bands become occupied or unoccupied, and their MAE contributions appear or disappear. Second, from the point of view of PT, the strongest contributions to MAE are from those virtual transitions between the unoccupied and occupied states near the Fermi level that are coupled by SOC Ke and van Schilfgaarde (2015); Daalderop et al. (1994). As a result, the largest contributions are located at and near the degenerate states across the Fermi contour, as expected from PT, and as shown in Fig. 6(d).

As shown in Figs. 6(c) and 6(d), the overall $(k_{x},k_{y})$-resolved MAE calculated using the two methods share great similarity. The largest differences exist around the Fermi contour, where the relaxation effect is large. The large contributions (near the degenerate states) observed in PT disappear in MFT. Indeed, the MFT results show more complex features than the non-SOC Fermi contour, which corresponds to the Fermi surfaces when SOC is on and the magnetization direction being along in- and out-of-plane directions. The SOC-induced lifting of band degeneracy, especially near the $\epsilon_{\text{F}}$, is often discussed to explain MAE in various systems Wang et al. (1993); Daalderop et al. (1994); Moos et al. (1996); Ravindran et al. (2001).

Thus, we demonstrate that MFT and PT give very similar results for not only total MAE but also for $k$-resolved MAE. To achieve this, it is important to use the reference energy $\epsilon^{0}_{\text{F}}$ as in Eq. (3). Otherwise, if one just resolves MAE using Eq. (2), a very different dependence can be obtained, and $k$-resolved MAE will only manifest the Fermi surface, near which the change of band occupancy is significant. MAE has often been resolved into $k$ space along certain line paths between high-symmetry points. Not surprisingly, using Eq. (2) to resolve MAE will result in spikes at points where bands cross the Fermi level. Further resolution of MAE into atomic sites by projecting eigenvectors of each $k$ point may produce unphysical results Subkow and Fähnle (2009); Åberg et al. (2015); Šipr et al. (2014).

III.4 Two-ion MA: Intersublattice contribution

Site-resolved MAE values, together with exchange parameters, can be used to construct an atomic spin Hamiltonian for subsequent Monte Carlo or spin-dynamics simulations to calculate the temperature dependence of magnetic properties. Methods such as evaluating the anisotropy of on-site SOC energy Antropov et al. (2014), which is a local quantity, have been used to resolve the MAE contribution from each individual sublattice. Here, we use PT in TB to resolve MAE into sublattices and validate the decomposition using the SOC-strength-scaling approach in vasp.

As discussed above, PT can well describe the MAE in these systems. To quantify the single-ion and two-ion contributions of MAE, we first use PT within TB to resolve MAE into intrasublattices and intersublattice contribution. Results are summarized in Table 3. Interestingly, all intrasublattice contributions are positive in all elements except for the $s$-like Al site in MnAl, where it vanishes. The sign of the intersublattice contribution varies and its amplitude is generally comparable to or even larger than that of either individual intrasublattice. For FePt and CoPt, the major contributions are from the Pt sites. The intersublattice contributions are negative for FePt, CoPt, FePd, and FeNi. Especially in FeNi, the amplitude of the negative intersublattice contribution is larger than each individual intrasublattice contribution. In contrast, the intersublattice contribution is positive in MnGa. An even larger positive intersublattice contribution is found in hypothetical FeCo with $c/a=1.1$.

To validate TB results, we also investigate the intrasublattice and intersublattice MAE contributions by scaling the SOC strength in vasp. Figure 7 shows the normalized MAE as a function of the SOC-scaling factors (between 0 and 1) for $L1_{0}$ materials using the second-variation method in vasp. For all compounds, the sign and relative amplitude of intrasublattice and intersublattice contributions agrees well with TB-PT results. Furthermore, we fit MAE as a function of SOC-scaling factors (with $0.9<\lambda_{i}<1.1$) using Eq. (16) and further confirm that the second-order terms agree very well with TB results listed in Table 3. The fourth-order terms are generally small especially for $4d$ and $3d$ compounds. Owing to stronger SOC, FePt and CoPt have larger and negative fourth-order contributions: $\sim$8% of total MAE. A previous study also found a small and negative high-order contribution to MAE in FePt Ayaz Khan et al. (2016). The good agreement between TB and vasp further validates the accuracy of the PT approach for those systems. Unlike the scaling procedure, the PT approach resolves all contributions in a single calculation. Thus, for the same analysis, once the TB Hamiltonian is constructed, the TB-PT approach is orders of magnitude faster than the SOC-scaling approach in DFT.

We further investigate FeNi, in which the intersublattice MAE dominates. Along each direction, the intrasublattice SOC energy is large and dominated by the $\lambda^{2}$ term while the intersublattice term is rather small. However, the majority of intrasublattice terms cancel out between the two directions, while the intersublattice term does not. Hence, the intersublattice term becomes dominant in the anisotropy of SOC. In other words, the intrasublattice terms in SOC are large but more isotropic, while the intersublattice terms are smaller but more anisotropic with respect to the magnetization quantization direction.

A large intersublattice MAE contribution may suggest the need to go beyond the single-ion MAE model when interfacing ab initio methods with atomic spin simulation. To simulate temperature-dependent MAE or other magnetic properties, exchange coupling is often included over the first one or two nearest-neighbor shells while MAE is often included using the single-ion MAE term $k_{i}(s^{z}_{i})^{2}$. For systems with strong intersublattice contribution, one may also need to include two-ion terms such as $k_{ij}s^{z}_{i}s^{z}_{j}$ into the atomic spin Hamiltonian.

IV Conclusions

Using $L1_{0}$ systems as a test case, we demonstrate that the ab initio TB framework, constructed using the maximally localized Wannier functions method, can be used to efficiently and accurately compute and resolve MAE in transition metal systems. With the magnetic force theorem, TB quantitatively reproduces DFT results over the full band-filling range from the bottom of valence band to a few eV above the Fermi level. When calculating $k$-resolved MAE in TB, the magnetic force theorem and perturbation theory results agree with one another, and both yield MAE contour maps that are consistent with the Fermi surface. We also resolve MAE into intrasublattice and intersublattice contributions using perturbation theory in TB and a scaled spin-orbit strength procedure in DFT. The results using these two methods are in excellent agreement. We found that the sign of the intersublattice contribution differs among compounds, and its amplitude may be comparable to or even larger than the intrasublattice contributions, suggesting the need to go beyond the single-ion MAE model. Depending on the system size, once the TB Hamiltonian is constructed, it can speed up the calculation by orders of magnitude, providing an efficient, accurate, and high-resolution method to calculate MAE. We expect that it can be applied to more complex compounds and structures to compute and analyze MAE. Finally, this realistic TB method can also be interfaced with ab initio methods beyond DFT, such as the much more expensive self-consistent $GW$ methods van Schilfgaarde et al. (2006); Kotani et al. (2007), to greatly accelerate the calculations and analysis of MAE or other SOC-related properties using those methods.

Acknowledgments

I thank B. Harmon for helpful discussions. This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, and Early Career Research Program. The initial development of the software was supported by the Laboratory Directed Research and Development Program of Ames Laboratory. Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No. DE-AC02-07CH11358.

Appendix A Spin-orbit coupling operators in real spherical harmonics representation

[TABLE]

In the presentation of complex spherical harmonics $Y_{\ell}^{m}$, the non-vanished matrix elements of $L_{z}$, $L_{+}$, and $L_{-}$ are

[TABLE]

wannier90 uses the real spherical harmonics $Y_{\ell m}$, also known as tesseral spherical harmonics, which can be written in terms of the complex spherical $Y_{\ell}^{m}$ as

[TABLE]

The angular momentum matrices in the real-spherical-harmonics representation, ${\bf O}^{(\mathbb{R})}$, can be obtained by directly evaluating the angular momentum operator on real spherical functions in Eq. (20), or transforming ${\bf O}^{(\mathbb{C})}$, the corresponding operator matrix from complex representation:

[TABLE]

From Eq. (20), the transfer matrix ${\bf U}_{\mathbb{C}\leftarrow\mathbb{R}}$ can be written as

[TABLE]

For $p$ and $d$ orbitals, only the corresponding subblocks of Eq. (22) are needed.

Similarly, for quantization along other directions, the matrix can be rotated by using Wigner matrix

[TABLE]

and

[TABLE]

where $\theta$ and $\varphi$ are the angles of the direction of magnetization when the unit vector is defined by $\hat{\bf n}=(\sin\theta\cos\varphi,\sin\theta\sin\varphi,\cos\theta)$. When the spin quantization axis is along [110] direction, the SOC Hamiltonian can be written as

[TABLE]

Appendix B Charge and moment calculated in TB and VASP

Table 4 lists the site-resolved charge and magnetic moments in $L1_{0}$ compounds calculated in TB and vasp.