Spin wave dispersion of 3d ferromagnets based on QSGW calculations

TL;DR

This paper uses QSGW calculations to accurately determine spin wave dispersions in 3d ferromagnets, improving agreement with experimental data over traditional methods.

Contribution

It introduces a QSGW-based approach for calculating spin wave dispersions, providing more accurate spin stiffness constants for ferromagnetic materials.

Findings

QSGW yields better spin stiffness constants than LDA.

Spin wave dispersions in Co match experimental data.

Discrepancies in FeCo suggest need for further research.

Abstract

We calculate transverse spin susceptibility in the linear response method based on the ground states determined in the quasi-particle self-consistent $GW$ (QSGW) method. Then we extract spin wave (SW) dispersions from the susceptibility. We treat bcc Fe, hcp Co, fcc Ni, and B2-type FeCo. Because of the better description of the independent-particle picture in QSGW, calculated spin stiffness constants for Fe, Co, and Ni give much better agreement with experiments in QSGW than that in the local density approximation (LDA), where the stiffness for Ni in LDA is two times bigger than the experiment. For Co, both acoustic and optical branches of SWs agree with the experiment. As for FeCo, we have some discrrepancy between the spin stiffness in QSGW and that in the experiment. We may need further theoretical and experimental investigations on the discrepancy.

Tables

Table 1. Table 1: t 2 ​ g subscript 𝑡 2 𝑔 t_{2g} level of Fe at Γ Γ \Gamma for the minority spin, and that at N for the majority spin. These are in LDA, in QSGW, in addition to the experimental data by ARPES ARPES_Fe . Energy is relative to E Fermi subscript 𝐸 Fermi E_{\rm Fermi} .

	band energy [eV]
	LDA	QSGW	Expt. ARPES_Fe
$\Gamma$(Minority)	-0.32	-0.11	-0.19
N (Majority)	-0.74	-0.68	-0.57

Table 2. Table 2: Calculated stiffness constant D 𝐷 D for Fe, Ni, Co and FeCo. The results by other groups are shown together; the LR LR_Friedrich , with the LF Pajda , and with the time-dependent DFT (TDDFT) Niesert (on average). In addition we show inelastic neutron scattering data Fe_expt1 ; Fe_expt2 ; Ni_expt2 ; Ni_expt3 ; hcpCo_Shirane ; hcpCo_Shirane ; feco_Lowde .

	$D$ [meV$\cdot$Å2]
Material	LR (LDA)	LR (QSGW)	Expt.	LR (GGA) LR_Friedrich	LF Pajda	TDDFT Niesert
bcc Fe	155	222	230 (RT) Fe_expt1	248	250	189
			280 (4.2 K) Fe_expt2
fcc Ni	873	449	433 Ni_expt2		756	1097
			555 Ni_expt3
hcp Co [100]	565	486	478 hcpCo_Perring
hcp Co [001]	752	532	410 hcpCo_Perring
			510 hcpCo_Shirane
B2 FeCo	407	307	450-500 feco_Lowde

Equations

$${\mathcal{H}}_0(\mathbf{r},{\mathbf{r}}^{\prime })=-\frac{{\nabla }^2}{2}+V_{\mathrm{ext}}+V_{\mathrm{H}}+V_{\mathrm{xc}}(\mathbf{r},{\mathbf{r}}^{\prime }),$$

$$\mathcal{H}(\mathbf{r},{\mathbf{r}}^{\prime };\omega )=-\frac{{\nabla }^2}{2}+V_{\mathrm{ext}}+V_{\mathrm{H}}+\Sigma (\mathbf{r},{\mathbf{r}}^{\prime };\omega ).$$

$$\displaystyle V^{\rm xc}=\frac{1}{2}\sum_{ij}\ket{\psi_{i}}\Bigl{\{}{\rm Re}\bigl{[}\Sigma(\varepsilon_{i})\bigr{]}_{ij}+{\rm Re}[\Sigma(\varepsilon_{j})]_{ij}\Bigr{\}}\bra{\psi_{j}},$$

$$R^{+-}(1,2)=\frac{\delta ⟨{\hat{\sigma }}^{+}(1)⟩}{\delta B^{+}(2)},$$

$$⟨{\hat{\sigma }}^{+}(1)⟩=-i\sum _{\alpha ,\beta }{\sigma }_{\beta \alpha }^{+}{G}_{\alpha \beta }(1,1^{+})(\alpha ,\beta \in \{\uparrow ,\downarrow \}),$$

$$R_{\uparrow \downarrow }^{(4)}(12,34)=K_{\uparrow \downarrow }(12,34)$$

$$+\int \int K_{\uparrow \downarrow }(12,56)W(5^{+},6)R_{\uparrow \downarrow }(56,34)d5d6,$$

$$-K_{\uparrow \downarrow }(12,34)=-i{G}_{\uparrow }^0(1,3)G_{\downarrow }^0(4,2^{+}),$$

$$-K_{\uparrow \downarrow }({\mathbf{r}}_1,{\mathbf{r}}_2;{\mathbf{r}}_3,{\mathbf{r}}_4;\omega )$$

$$=\sum _{\mathbf{k},n}^{\mathrm{occ}}\sum _{{\mathbf{k}}^{\prime },n^{\prime }}^{\mathrm{unocc}}\frac{{\Psi }_{\mathbf{k}n\downarrow }^{\ast }({\mathbf{r}}_2){\Psi }_{\mathbf{k}n\downarrow }({\mathbf{r}}_4){\Psi }_{{\mathbf{k}}^{\prime }{n}^{\prime }\uparrow }({\mathbf{r}}_1){\Psi }_{{\mathbf{k}}^{\prime }{n}^{\prime }\uparrow }^{\ast }({\mathbf{r}}_3)}{\omega -({\epsilon }_{{\mathbf{k}}^{\prime }{n}^{\prime }\uparrow }-{\epsilon }_{\mathbf{k}n\downarrow })+i\delta }$$

$$+\sum _{\mathbf{k},n}^{\mathrm{unocc}}\sum _{{\mathbf{k}}^{\prime },n^{\prime }}^{\mathrm{occ}}\frac{{\Psi }_{\mathbf{k}n\downarrow }^{\ast }({\mathbf{r}}_2){\Psi }_{\mathbf{k}n\downarrow }({\mathbf{r}}_4){\Psi }_{{\mathbf{k}}^{\prime }{n}^{\prime }\uparrow }({\mathbf{r}}_1){\Psi }_{{\mathbf{k}}^{\prime }{n}^{\prime }\uparrow }^{\ast }({\mathbf{r}}_3)}{-\omega -({\epsilon }_{\mathbf{k}n\downarrow }-{\epsilon }_{{\mathbf{k}}^{\prime }{n}^{\prime }\uparrow })+i\delta },$$

$${\Psi }_{\mathbf{k}n}(\mathbf{r})$$

$$w_{\mathbf{R}i}^{\mathbf{k}}(\mathbf{r})=\frac{1}{\sqrt{N}}\sum _{\mathbf{T}}{w}_{\mathbf{R}i}(\mathbf{r}-\mathbf{R}-\mathbf{T})\exp (i\mathbf{k}\cdot \mathbf{T}),$$

$$\sum _{{\mathbf{R}}^{\prime }j}{H}_{\mathbf{R}i{\mathbf{R}}^{\prime }j}^{\mathbf{k}}{a}_{{\mathbf{R}}^{\prime }j}^{\mathbf{k}n}$$

$$-K_{\mathbf{R}ij,{\mathbf{R}}^{\prime }kl}^{\uparrow \downarrow }(\mathbf{q},\omega )$$

$$=\frac{1}{N}\sum _{\mathbf{k}}^{\mathrm{BZ}}\sum _n^{\mathrm{occ}}\sum _{n^{\prime }}^{\mathrm{unocc}}\frac{a_{\mathbf{R}j\beta }^{\mathbf{k}n\ast }{a}_{{\mathbf{R}}^{\prime }l\beta }^{\mathbf{k}n}{a}_{\mathbf{R}i\alpha }^{\mathbf{k}+\mathbf{q}{n}^{\prime }}{a}_{{\mathbf{R}}^{\prime }k\alpha }^{\mathbf{k}+\mathbf{q}{n}^{\prime }\ast }}{\omega -({\epsilon }_{\mathbf{q}+\mathbf{k}{n}^{\prime }\uparrow }-{\epsilon }_{\mathbf{k}n\downarrow })+i\delta }$$

$$+\frac{1}{N}\sum _{\mathbf{k}}^{\mathrm{BZ}}\sum _n^{\mathrm{unocc}}\sum _{n^{\prime }}^{\mathrm{occ}}\frac{a_{\mathbf{R}j\beta }^{\mathbf{k}n\ast }{a}_{{\mathbf{R}}^{\prime }l\beta }^{\mathbf{k}n}{a}_{\mathbf{R}i\alpha }^{\mathbf{k}+\mathbf{q}{n}^{\prime }}{a}_{{\mathbf{R}}^{\prime }k\alpha }^{\mathbf{k}+\mathbf{q}{n}^{\prime }\ast }}{-\omega -({\epsilon }_{\mathbf{k}n\downarrow }-{\epsilon }_{\mathbf{q}+\mathbf{k}{n}^{\prime }\uparrow })+i\delta }.$$

Full text

I Introduction

Spin wave (SW) is one of the important factors to control magnetic properties of material. SW is excited at considerably low temperature compared to room temperature (RT), and its energy range typically lies in a few hundred meV. When one magnetic moment tilted from the parallel spin configuration, the exchange interaction triggers the SW propagation throughout the material as collective excitation. We can observe SWs in bulk materials by inelastic neutron scattering experiment, e.g., in bcc Fe Fe_expt1 , fcc Ni Ni_expt2 , and even half-metals like perovskite $\rm{{La_{0.7}Sr_{0.3}MnO_{3}}}$ lsmo_expt . In addition to collective excitation, another magnetic excitation like spin-flip excitation is called Stoner excitation, whose excitation energy is related to the exchange splitting ${\Delta}E_{\rm x}$. We can experimentally observe Stoner excitation by the high energy experiment such as spin-polarized electron energy loss spectroscopy (SPEELS) Stoner_Vollmer . High energy SWs are strongly damped because of the hybridization with the Stoner excitation.

Let us explain how we determine the spin stiffness $D$ experimentally. From the macroscopic point of view, the Bloch’s $T^{\frac{3}{2}}$ rule Bloch_rule in the temperature dependence of magnetization M(T) is derived from the SW theory. For the wave vector ${\bf q}\sim 0$, the SW dispersion ${\omega({\bf q})}$ behaves as ${\omega({\bf q})}=D{\bf q}^{2}$. Since this behavior of ${\omega({\bf q})}$ results in the $T^{\frac{3}{2}}$ rule in low temperature, we can determine $D$ by analyzing the temperature dependence of magnetization Fe_expt2 .

We mainly have three methods to calculate ${\omega({\bf q})}$ in the first-principles methods. The first one is the Lichtenstein formula (LF) Lichtenstein . Assuming the Heisenberg model, we calculate exchange interaction ${J_{ij}}$ or its Fourier transform ${J({\bf q})}$ based on the magnetic force theorem MFT . Here $i,j$ are for site indices. Then ${\omega({\bf q})}$ is calculated from ${J({\bf q})}$. In Ref. Lichtenstein, , they calculated ${J_{ij}}$ up to the second nearest neighbors, resulting in $D$, which are in good agreement with experiments for Fe and Ni. Later, Pajda et al. investigated the convergence of $D$ for a range of neighbors and found that converged $D$ are in good agreement with experiments for Fe but overestimated for Ni Pajda .

The second one is the frozen magnon method (FMM) FMM_Halilov , which assumes the Heisenberg model as in LF. In FMM, we employ adiabatic approximation; namely, we neglect motions of the magnetic moment compared to electron motions. Then we calculate ${J({\bf q})}$ from the constraint spin-spiral configurations with the fixed magnitude of the magnetic moment. Once we get ${J({\bf q})}$, we solve the eigenvalue problem for deriving ${\omega({\bf q})}$. This method works well for bcc Fe FMM_Halilov ; Grotheer . Note that we can not describe the decay of collective SWs (Stoner damping) in both of these two methods.

The third one is the linear response (LR) method for transverse spin susceptibility $R^{+-}({\bf q},\omega)$ LR_Gross . The LR method directly gives ${\omega({\bf q})}$ in the reciprocal space. Cooke et al. first introduced the LR method for calculating $R^{+-}({\bf q},\omega)$, and they discussed Stoner damping in SWs in bcc Fe and fcc Ni LR_Cooke . Savrasov treated spin fluctuations based on the many-body perturbation theory and reproduced the experimental ${\omega({\bf q})}$ LR_Savrasov . Karlsson and Aryasetiawan also calculated $R^{+-}({\bf q},\omega)$ based on the Green function method LR_Karlsson . From a view of computational efficiency, Şaşıoǧlu et al. proposed a LR method with maximally-localized Wannier function (MLWF) LR_Sasioglu . In the method, we decrease to the second power of the number of a Wannier basis set and we can decrease the calculation cost. With this efficient method, they can use fine ${\bf q}$ mesh for calculating $R^{+-}({\bf q},\omega)$.

These three methods mainly have been applied to the ground states given in the local density approximation (LDA). However, the ground state given in LDA is not necessarily good enough. For example, Sponza et al. shows that $3d$-bandwidth and ${\Delta}E_{\rm x}$ in LDA are not good enough to calculate ${\omega({\bf q})}$ Sponza . In antiferromagnetic transition metal oxides such as NiO and MnO, the calculated ${\omega({\bf q})}$ does not agree with the experiment due to too small ${\Delta}E_{\rm x}$ and too small bandgap TMO_QSGW . Serious disagreement is also found in the ${\omega({\bf q})}$ in $\rm{{La_{0.7}Sr_{0.3}MnO_{3}}}$, for which LDA fails to reproduce the half-metallic electronic structure of that compound LSMO_QSGW . It is possible to start from the ground states of LDA+$U$; however, we sometimes have difficulty in determination of $U$. It may suggest a limitation of LDA+$U$ itself.

To overcome such limitations in LDA, Kotani et al. calculated ${\omega({\bf q})}$ for strongly-correlated materials in an LR method for the ground states determined in the quasi-particle self-consistent $GW$ (QSGW) method TMO_QSGW ; LSMO_QSGW . Then we see reasonable agreement with experiments for NiO and MnO because QSGW gives good descriptions of the band quantities such as ${\Delta}E_{\rm x}$ and bandgaps Deguchi_2016 . We expect such good agreement for wide-range of materials. However, Kotani’s LR method used in Refs. TMO_QSGW ; LSMO_QSGW is too simple to apply a wide range of materials. Thus we implemented the efficient LR method to calculate $R^{+-}({\bf q},\omega)$ based on the MLWF given by Şaşıoǧlu et al. LR_Sasioglu in QSGW calculation package $ecalj$ compiled by Kotani et al. ecalj . We demonstrate how the method works for typical ferromagnets such as bcc Fe, fcc Ni, hcp Co, and B2 FeCo (CsCl structure) and we discuss the difference between LDA and QSGW. Except for FeCo, the SWs in QSGW agree with experiments. We find some discrepancies for FeCo.

II Computational methods

II.1 quasiparticle self-consistent GW (QSGW)

Until now, varieties of $GW$ calculations based on the Hedin’s GW approximation hedin_gw1 ; hedin_gw2 have been performed since it is introduced to the first-principles calculations by Hyberstein and Louie HybersteinLouie1986 . Most of the $GW$ calculations are so-called one-shot $GW$. Starting from $G^{0}$ for the one-body Hamiltonian in LDA $\mathcal{H}^{\rm LDA}_{0}$, we calculate corrections to the eigenvalues of $\mathcal{H}^{\rm LDA}_{0}$ to reproduce quasiparticle energies. In the one-shot $GW$, the self-energy for the corrections is given as $\Sigma(1,2)=iG^{0}(1,2)W(1^{+},2)$, where we use notation $1\equiv({\bf r}_{1},t_{1})$. The screened Coulomb interaction $W(1^{+},2)$ is calculated as $W=(1-vP)^{-1}v$ from the bare Coulomb interaction $v$ and the polarization function $P=-iG^{0}\times G^{0}$. The one-shot $GW$ has a shortcoming since the one-shot $GW$ is just a perturbation on top of $\mathcal{H}^{\rm LDA}_{0}$.

To overcome the shortcoming of the one-shot $GW$, we utilize QSGW method qsgw1 ; qsgw2 ; qsgw3 implemented in $ecalj$ package ecalj . Let us summarize QSGW method. At first, recall the above $GW$ procedure which can be applicable to any static one-body Hamiltonian $\mathcal{H}_{0}({\bf r},{\bf r}^{\prime})$ as

[TABLE]

where we have the external potential $V_{\rm ext}$, the Hartree potential $V_{\rm H}$, and the non-local exchange-correlation potential $V_{\rm xc}({\bf r},{\bf r}^{\prime})$. With $\Sigma(1,2)=iG^{0}(1,2)W(1^{+},2)$ where $G^{0}=1/(\omega-\mathcal{H}_{0})$, we have the energy-dependent one-body Hamiltonian $\mathcal{H}({\bf r},{\bf r}^{\prime};\omega)$ as

[TABLE]

That is, $GW$ approximation gives a procedure $\mathcal{H}_{0}\rightarrow\mathcal{H}$. QSGW requires “quasiparticle self-consistency”, that is, minimization of the difference between $\mathcal{H}_{0}$ and $\mathcal{H}$. The minimization gives the procedure $\mathcal{H}\rightarrow\mathcal{H}_{0}$, replacing the $\omega$-dependent $\Sigma$ in Eq. (2) with the static non-local exchange-correlation potential $V^{\rm xc}$ as

[TABLE]

where eigenvalues $\varepsilon_{i}$ and eigenfunctions $\psi_{i}$ are those of $\mathcal{H}_{0}$. This defines a procedure to give a new $\mathcal{H}_{0}$, $\mathcal{H}\rightarrow\mathcal{H}_{0}$. Thus we finally have a ’quasiparticle self-consistency’ cycle $\mathcal{H}_{0}\rightarrow\mathcal{H}\rightarrow\mathcal{H}_{0}\rightarrow\mathcal{H}\rightarrow\cdots$ (or $G^{0}\rightarrow G\rightarrow G^{0}\rightarrow\cdots$) until converged.

II.2 Dynamical magnetic susceptibility

In LR, we follow the procedure given in Ref. LR_Friedrich ; LR_Sasioglu . Here we treat the transverse spin susceptibility $R^{+-}(1,2)$, which describes the response of the expectation value of a spin density operator $\hat{\sigma}^{+}(1)$ to the the external magnetic field $B^{-}(2)$ as,

[TABLE]

where $1=({\bf r}_{1},t_{1})$. See Eq. (20) in Ref. LR_Friedrich, . Here the expectation value of $\hat{\sigma}^{+}(1)$ is given as

[TABLE]

where $G(1,1^{+})$ is the single-particle Green function from 1 to $1^{+}$. For our calculation below, it is convenient to consider four-points representation $R^{(4)}_{{\uparrow\downarrow}}(12,34)$. The trace of matrix $R^{(4)}_{{\uparrow\downarrow}}(11,33)$ leads to two-point representation $R^{+-}(1,2)$.

In order to obtain $R^{(4)}_{{\uparrow\downarrow}}(12,34)$, we solve the Bethe-Salpeter equation where we use the static screened Coulomb interaction $W(1^{+},2)$ which is $\propto\delta(t_{1}-t_{2})$. It is

[TABLE]

where $K_{{\uparrow\downarrow}}(12,34)$ is the non-interacting two-particle (particle-hole with opposite spin) propagator given as

[TABLE]

where we consider $t_{1}=t_{2}$ and $t_{3}=t_{4}$, $i.e.$, $K_{{\uparrow\downarrow}}({\bf r}_{1},{\bf r}_{2};{\bf r}_{3},{\bf r}_{4};t_{1}-t_{3})$. The Fourier transform is from $t_{1}-t_{3}$ to $\omega$. We symbolically solve Eq. (6) to be $R=K+KWK+KWKWK+\cdots=K(1-WK)^{-1}$, where the numerator $K$ describes the Stoner excitations, whereas zeros of the denominator $(1-WK)$ gives the collective excitation.

This $K_{{\uparrow\downarrow}}$ is given as

[TABLE]

where ${\bf k},{\bf k}^{\prime}$ are in the first Brillouin zone, $n(n^{\prime})$ is the band index summed over occupied (unoccupied) states, $\varepsilon_{{\bf k}n\uparrow}$ ($\varepsilon_{{\bf k}n\downarrow}$) is the $n$th majority (minority) band energy at ${\bf k}$, and $\Psi$ is the eigenfunction of $\mathcal{H}_{0}$.

As mentioned in Ref. LR_Sasioglu , in order to satisfy the Goldstone theorem $\omega({\bf q})\rightarrow 0$ $({\bf q}\rightarrow 0)$, we need to introduce a factor $\eta$ for $R=K(1-\eta WK)^{-1}$. In principle, the Goldstone theorem should be automatically satisfied with the LR method since we expect that the LR method evaluates the second derivative of the total energy of the ground states. However, our LR is not formulated to reproduce the second derivative exactly; furthermore, QSGW is not formulated to minimize the total energy. This simple scaling by introducing $\eta$ is a quick remedy to satisfy the theorem; their deviations from unity show the size of vertex corrections, which should be added to the interaction $W$. The calculated $\eta$ of LDA (QSGW) are 1.15 (1.19), 1.41 (1.87), 1.26 (1.33), and 1.05 (0.87) for Fe, Ni, Co, and FeCo, respectively. These $\eta$ are in good agreement with previous calculations 1.28, 1.5, and 1.33 for Fe LR_Friedrich , Ni LR_Sasioglu , and FeCo LR_Friedrich . The deviations are not small enough. We may need to treat the vertex correction accurately in order to override the ambiguity due to this quick remedy in the future.

II.3 Wannier representation

Based on Refs. MLWF1 ; MLWF2 , we generate MLWFs from eigenfunctions of LDA or QSGW. Once we generate MLWFs, we can obtain the Wannier representation of $R^{{\uparrow\downarrow}}$ as follow.

In the Wannier basis, we expand eigenfunctions as

[TABLE]

where $a_{{\bf R}i}^{{\bf k}n}$ is the expansion coefficient, ${\bf R}$ is atomic position in a primitive cell, $i$ is the Wannier orbital ($e.g.$ $i=3d_{xy}$) of each atom on ${\bf R}$. $w_{{\bf R}{i}}^{{\bf k}}({\bf r})$ is represented as a complete set of orthogonal basis {$w_{{\bf R}i}({\bf r})$},

[TABLE]

where ${\bf T}$ is the lattice translation vector and $N$ is the normalization constant satisfying the Born von Karman boundary condition. By using the orthogonality, the eigenvalue equations $\mathcal{H}\Psi_{{\bf k}n}({\bf r})={\varepsilon_{{\bf k}n}}\Psi_{{\bf k}n}({\bf r})$ can be rewritten with this Wannier representation,

[TABLE]

where the Hamiltonian matrix with Wannier basis $H_{{\bf R}i{\bf R}^{\prime}j}^{{\bf k}}$ is the Fourier transform of $H_{{\bf R}i{\bf R}^{\prime}j}^{{\bf T}-{\bf T}^{\prime}}\equiv\bra{w_{{\bf R}i}({\bf r}-{\bf R}-{\bf T})}\mathcal{H}\ket{w_{{\bf R}^{\prime}j}({\bf r}-{\bf R}^{\prime}-{\bf T}^{\prime})}$.

Substituting Eqs. (9) and (10) to Eq. (8) and using Fourier transform of real-space, we will obtain the time-ordered linear response function for a non-interacting system represented in a restricted Hilbert space,

[TABLE]

We calculate the imaginary part of $-K^{{\uparrow\downarrow}}_{{\bf R}ij,{\bf R}^{\prime}kl}({\bf q},\omega)$ by a tetrahedron method and obtain its real part by the Hilbert transform. The matrix element of $R^{{\uparrow\downarrow}}_{{\bf R}ij,{\bf R^{\prime}}kl}$ is calculated through $R=K(1-\eta WK)^{-1}$, where $W$ is calculated in the random phase approximation (RPA) in the product basis technique developed in Ref. Kotani_2001 .

II.4 Calculation details

All of the calculation procedures above are implemented in the first-principles package $ecalj$ ecalj ; Deguchi_2016 . The $ecalj$ is based on the linearized augmented plane-wave and muffin-tin orbital (MTO) method (PMT method), which combines augmented plane wave (APW) and MTO basis sets. We also generate MLWFs in $ecalj$. We perform LDA and QSGW calculations for band structures with $20\times 20\times 20$ and $16\times 16\times 16$ $k$-point mesh respectively. We consider 9 MLWFs ($spd$) for the $3d$ elemental materials (Fe and Ni) and 18 MLWFs for hcp Co and binary FeCo. In the calculations of $K^{{\uparrow\downarrow}}$, we use $48\times 48\times 48$ $q$-point mesh for the 3$d$ elemental material and $24\times 24\times 24$ for binary FeCo. We use static and onsite $W$, $i.e.$, we take $W_{ijkl}(\omega)=W_{{\bf R}ij,{\bf R}kl}(\omega=0)$. We use experimental lattice parameters, $a=2.867$ Å, $a=3.524$ Å, $a=2.850$ Å for Fe, Ni, and FeCo, respectively. For hcp Co, we use $a=2.507$ Å and $c=4.070$ Å.

III Results and discussion

III.1 bcc Fe

Figs. 1(a), (b), and (c) show the majority and minority band structures and the partial density of states in QSGW for Fe, while Figs. 1(d), (e), and (f) in LDA as well. Calculated total magnetic moments in LDA and QSGW are both 2.22 $\mu_{\rm B}$ for Fe, in agreement with the experimental value 2.22 $\mu_{\rm B}$ mmom_Danan , in contrast to 2.93 $\mu_{\rm B}$ in the fully self-consistent $GW$ method scGW . Our results are consistent with Ref. Sponza by Sponza et al. The superposed Wannier band structures in Eq. (11) by broken lines are entirely on the original band structures by bold grey lines. Size of colored circles show the weights of each MLWF. In Table 1, we show the $t_{2g}$ of minority spin at $\Gamma$ and that of majority spin at N in LDA and QSGW. QSGW gives better agreement with the angle-resolved photoemission spectroscopy (ARPES) data ARPES_Fe . The $3d$-bandwidth in QSGW is a little smaller than that in LDA. Except for this difference, the overall shapes of the majority and the minority bands are similar in both LDA and QSGW.

Fig. 2(a) shows $-{\rm Im}[K^{+-}({\bf q}=0,\omega)]$ in LDA and in QSGW, where $K^{+-}({\bf q},\omega)$ means the trace of the matrix $K^{{\uparrow\downarrow}}$ given as $K^{+-}({\bf q},\omega)=\sum_{{\bf R},i,j}K^{{\uparrow\downarrow}}_{{\bf R}ii,{\bf R}jj}({\bf q},\omega)$. We use a little different definition from Refs. LR_Sasioglu, ; LR_Friedrich, ; LR_Friedrich2018, , thus it is not meaningful to compare absolute value of $K^{+-}({\bf q},\omega)$ with their results. As shown in the figure, QSGW gives smaller ${\Delta}E_{\rm x}$ and 3$d$-bandwidth, which is consistent with results by Sponza et al. Roughly speaking, the shape of $-{\rm Im}[K^{+-}({\bf q}=0,\omega)]$ agree with the shape of density of states (DOS) of majority spin. The peak around 2 eV originates from the $t_{2g}^{\uparrow}$-$t_{2g}^{\downarrow}$ and $e_{g}^{\uparrow}$-$e_{g}^{\downarrow}$ transition, i.e., vertical transitions to the unoccupied minority states above the Fermi energy $E_{\rm Fermi}$ from the occupied majority states just below the $E_{\rm Fermi}$ in Fig. 1. The second peak around 4 eV is stemmed from another $e_{g}^{\uparrow}$-$e_{g}^{\downarrow}$ transition to $E_{\rm Fermi}+2$ eV in minority states from $E_{\rm Fermi}-2$ eV in majority states.

We see two features in the difference between LDA and QSGW in $-{\rm Im}[K^{+-}({\bf q}=0,\omega)]$ shown in Fig. 2(a). One is that the width of the peak around $2$ eV in QSGW is wider than that in LDA. The difference of DOS in LDA and QSGW can not explain this fact; it can be due to the difference of eigenfunctions. The peak becomes wider in QSGW, probably because of the general tendency of QSGW that it makes a more significant difference between occupied $3d$ states and unoccupied $3d$ states. The former is more localized, and the latter more extended in comparison with the case in LDA. The other is the width due to the $3d$ band; corresponding to the width of $3d$ band shown in the inset of Fig. 2(a), we see narrower width in $-{\rm Im}[K^{+-}({\bf q}=0,\omega)]$ in QSGW.

Figs. 2(b) and (c) show the Stoner excitation spectrum $-{\rm Im}[K^{+-}({\bf q},\omega)]$ in LDA and QSGW. Our LDA results give good agreement with Fig. 6 in Ref. LR_Friedrich2018, . We see red triangle-like strong intensity around $\Gamma$, especially in LDA. The center of peak moves up as a function of ${\bf q}$. This is because shifted $q$ from $\Gamma$ requires corresponding energy shift to trace the peak of $-{\rm Im}[K^{+-}({\bf q},\omega)]$ as a function of $\omega$. This is explained in Fig. 7 of Ref. LR_Friedrich2018, .

Fig. 3 shows ${\rm Im}[R^{+-}({\bf q},\omega)]$ in LDA (a) and in QSGW (b), where $R^{+-}({\bf q},\omega)$ means the trace of the matrix $R^{{\uparrow\downarrow}}$ given as $R^{+-}({\bf q},\omega)=\sum_{{\bf R},i,j}R^{{\uparrow\downarrow}}_{{\bf R}ii,{\bf R}jj}({\bf q},\omega)$. We superpose experimental data Fe_expt1 ; expt_Loong on it. We also superpose the SW dispersion calculated with the LF Pajda in LDA, and that with FMM in LDA FMM_Halilov . These are not only in (a) but also in (b) as a guide of eye. As shown in Fig. 3, the peak broadening due to the Stoner damping can be seen even below 100 meV because bcc Fe is a weak ferromagnet, whose majority and minority $3d$ have relatively large DOS at $E_{\rm Fermi}$ as shown in the inset of Fig. 2(a). This results in relatively large low-energy Stoner excitations. It means that SWs are getting to be hybridized well with Stoner excitation immediately after departing from $\Gamma$. The strong damping around H is also seen in the previous calculation combining the the generalized gradient approximation (GGA) and the MLWF approach with 6 MLWFS ($sd$) LR_Friedrich . Our LDA calculation indicates Kohn anomalies in $\Gamma$-H, H-N, and $\Gamma$-N, which are also found in the other calculations Pajda ; FMM_Halilov ; Grotheer . We checked calculations with denser q-point mesh (60$\times$60$\times$60) and confirmed the strong anomaly at 2/3 along $\Gamma$-N in LDA, and especially in QSGW. Ref. LR_Friedrich2018 explains how such anomalies can be traced back to the band structures, although they have not given explicit analysis. Real metals such as Fe can have complicated band structures, resulting in too complicated Fermi-surface-nestings like phenomena to be analyzed. Thus, we also have not yet got into such analysis. We are somehow skeptical whether it is worth to do or not.

In Table 2, we summarize calculated results of stiffness constant $D$, with another LR result based on the GGA LR_Friedrich , and with that of the time-dependent density functional theory (TDDFT) Niesert . To obtain $D$, we fit the calculated SW dispersion by quadratic functions. For the fitting, we just take peaks for small $q$ as $|{\bf q}|<0.20(\frac{2\pi}{a})$ where little Stoner damping occurs. Details for Fe and Ni are in supplements Supplemental . LDA gives $D=155$ meV$\cdot$Å2, which is a little smaller than experiments $D=230$, $280$ meV$\cdot$Å2 Fe_expt1 ; Fe_expt2 . On the other hand, QSGW gives $D=222$ meV$\cdot$Å2 in much better agreement with the experimental values. Note that we see a contradiction between our LR (LDA) and the other two previous calculations, the LR (GGA) and the LF. Our values $D=155$ meV$\cdot$Å2 is too low in comparison with the other data 248, 250 meV$\cdot$Å2, although the smaller difference from $D=189$ meV$\cdot$Å2 in TDDFT. However, we currently have no definite idea to resolve the discrepancy from these previous works.

III.2 fcc Ni

The calculated magnetic moment for Ni in LDA is in agreement with the experiment, 0.62 $\mu_{\rm B}$ mmom_Danan . On the other hand, QSGW gives 0.80 $\mu_{\rm B}$. Sponza et al. Sponza indicates that this is reasonable because we have not taken into account the longitudinal quantum spin fluctuation. In LDA, we may have accidentally had a good agreement because of too small exchange splitting cancels the fact that calculations do not include the fluctuation.

Fig. 4(a) shows the $-{\rm Im}[K^{+-}({\bf q}=0,\omega)]$ in Ni. Peaks at 0.7 eV and 0.8 eV in LDA and QSGW are the Stoner gaps, corresponding to the difference of peaks between majority and minority spins in DOS shown in its inset. ${\Delta}E_{\rm x}$ given in LDA and QSGW are about two times larger than 0.3 eV, which is the value obtained by ARPES at $L_{3}$ point exchange_Ni . Sponza et al. Sponza indicates that the overestimation is due to the missing of spin fluctuations. Figs. 4(b) and (c) show $-{\rm Im}[K^{+-}({\bf q},\omega)]$ in LDA and QSGW. Our LDA results give good agreement with Fig. 6 of Ref. LR_Friedrich2018, . We see that strong intensity around $\Gamma$ get broadened as a function of ${\bf q}$ as in the case of homogeneous electron gas shown in Fig. 5 of Ref. LR_Friedrich2018, . In QSGW, ${\bf q}$-dependence of $-{\rm Im}[K^{+-}({\bf q},\omega)]$ looks slightly weakened around $\Gamma$, probably because of the reflection of flattened (weak $q$-dependent) 3$d$ band.

In Fig. 5 (a), we show ${\rm Im}[R^{+-}({\bf q},\omega)]$ in LDA. We can identify the SW dispersion in the whole BZ in contrast to the case of Fe in Fig. 3. Our SW dispersion in LDA is consistent with a previous LR calculation by Savrasov LR_Savrasov and a TDDFT calculation by Niesert Niesert . As superposed in Fig. 5, results with FMM FMM_Halilov and with the LF Pajda give a little lower ${\omega({\bf q})}$. Let us compare QSGW result shown in Fig. 5(b) with (a), where we can use black lines as a guide of eye. ${\omega({\bf q})}$ curvature around $\Gamma$ is smaller in QSGW. In fact, Table 2 shows that QSGW gives very smaller $D=449$ meV$\cdot$Å2 around $\Gamma$ than $D=873$ meV$\cdot$Å2 in LDA. This is in agreement with the experimental values $D=433$, $555$ meV$\cdot$Å2 Ni_expt2 ; Ni_expt3 . This is the reflection of weak ${\bf q}$-dependence of $-{\rm Im}[K^{+-}({\bf q},\omega)]$ around $\Gamma$ in the previous paragraph. Along $\Gamma$-L, QSGW successfully trace an experiment Ni_expt4 even up to the half of the BZ boundary. Although (b) may be taken as a simple elongation of (a) at a glance, it is not true if we take the behavior around $\Gamma$ into account. In Ref. LR_Karlsson , Karlsson and Aryasetiawan gives good agreement with the SW dispersion along [100] by adjusting the ${\Delta}E_{\rm x}$ of Ni. However, such a procedure may give a simple shrinkage. Thus the physical mechanism in QSGW is very different from their method even though both our QSGW and their method reproduce the experimental $D$.

III.3 hcp Co

Fig. 6(a) shows the $-{\rm Im}[K^{+-}({\bf q}=0,\omega)]$ in Co and Figs. 6(b) and (c) show $-{\rm Im}[K^{+-}({\bf q},\omega)]$ in LDA and QSGW. The calculated magnetic moments per Co atom is 1.67 $\mu_{\rm B}$ in LDA, 1.76 $\mu_{\rm B}$ in QSGW. These are a little larger than the experiment 1.58 $\mu_{\rm B}$ hcpCo_moment . It is reasonable in the sense that the QSGW value relative to experiment is 1.76 $\mu_{\rm B}$/1.58 $\mu_{\rm B}$, in between 2.22 $\mu_{\rm B}$ /2.22 $\mu_{\rm B}$ (Fe) and 0.80 $\mu_{\rm B}$ /0.62 $\mu_{\rm B}$ (Ni). Let us compare peaks of $3d$ shown in insets with those for Fe and Ni (Figs. 2 and 4). In QSGW, $3d$ bands are narrower than LDA in both of majority, and minority spins in Co and Ni, in contrast to the case of Fe where little narrowing of DOS in the minority spins. It is probably because the bcc structure has more hybridization with $sp$ bands than fcc and hcp. In Co, the largest peaks of $3d$ are pushed down by QSGW relative to LDA, with keeping the exchange splitting. Thus changes of $-{\rm Im}[K^{+-}({\bf q}=0,\omega)]$ from QSGW to LDA are similar in Fe and Co. As we already noted in Sec.III.1, we admit several universal tendencies of QSGW relative to LDA, however, such changes of DOS and $-{\rm Im}[K^{+-}({\bf q}=0,\omega)]$ are hardly predicted without calculations in practice.

In Fig. 7(a), we show ${\rm Im}[R^{+-}({\bf q},\omega)]$ in LDA together with plots of the SW dispersion given by the FMM FMM_Halilov (black broken lines) and by the LF Pajda (black lines). In these plots, two branches appear because of two atoms per primitive cell. The LF traces peaks of our ${\rm Im}[R^{+-}({\bf q},\omega)]$ very well especially along $\Gamma$-A-K-H-A. At M around, the black lines are slightly lower than the peak of ${\rm Im}[R^{+-}({\bf q},\omega)]$ seen at $\sim$800 meV. Near $\Gamma$, ${\rm Im}[R^{+-}({\bf q},\omega)]$ shows no optical branch. Experimental data shown by oval circles hcpCo_Perring ; hcpCo_Shirane are a little lover than the plots and peaks of ${\rm Im}[R^{+-}({\bf q},\omega)]$.

In contrast, we have an impressive agreement with the experiment in QSGW. As seen in Fig. 7(b), oval circles are on the peak of ${\rm Im}[R^{+-}({\bf q},\omega)]$ in QSGW. The calculated $D$ shown in Table 2 in QSGW are 486 meV$\cdot$Å2 along [100], and 532 meV$\cdot$Å2 along [001]. These give much better agreements with experiments, consistent with the agreement in Fig. 7(b). This agreement of the SW energy is probably originated from narrower $3d$ band in QSGW, resulting weaker ${\bf q}$-dependence of $-{\rm Im}[K^{+-}({\bf q},\omega)]$, rather than LDA.

III.4 B2 FeCo

We treat B2 FeCo in the CsCl structure. Calculated magnetic moments per cell are 4.44 $\mu_{\rm B}$ in LDA, 4.80 $\mu_{\rm B}$ in QSGW. The latter is close to experiment 4.70 $\mu_{\rm B}$ feco_Goldman . It is consistent with other compounds TMO_QSGW ; LSMO_QSGW where QSGW give agreements with experiments as for magnetic moments when LDA gives underestimation. Alternatively, we may take FeCo as a case between Fe and Co. Since QSGW/experiment = 2.22 $\mu_{\rm B}$/2.22 $\mu_{\rm B}$ for Fe, = 1.76 $\mu_{\rm B}$/1.58 $\mu_{\rm B}$ for Co, we may say that slight overestimation 4.80 $\mu_{\rm B}$/4.70 $\mu_{\rm B}$ is reasonable.

Fig. 8(a) shows $-{\rm Im}[K^{+-}({\bf q}=0,\omega)]$ in LDA and QSGW. In its inset, ${\Delta}E_{\rm x}$ is $\sim$ 2.8 eV in QSGW while $\sim$2.2 eV in LDA. The difference results in the difference of peaks in $-{\rm Im}[K^{+-}({\bf q}=0,\omega)]$. Figs. 8(b) and (c) show $-{\rm Im}[K^{+-}({\bf q},\omega)]$ in LDA and QSGW, although we see no specific features worth to be mentioned.

Fig. 9 shows ${\rm Im}[R^{+-}({\bf q},\omega)]$ in (a) LDA and in (b) QSGW, together with the previous SW calculation in the FMM Grotheer . ${\rm Im}[R^{+-}({\bf q},\omega)]$ in LDA shows the lower peaks of ${\omega({\bf q})}$ than FMM. ${\rm Im}[R^{+-}({\bf q},\omega)]$ in LDA gives $D=407$ meV$\cdot$Å2 is a little smaller than 500 meV$\cdot$Å2 by Grotheer Grotheer . The optical branch is weakened as in the case of Fe. Weak peak around $\sim$ 600 meV are close to ${\omega({\bf q})}$ in FMM.

In QSGW, there is lower ${\omega({\bf q})}$ in the whole BZ as in the case of Co. Table 2 shows that $D=307$ meV$\cdot$Å2 in QSGW is much smaller than the experiment 450-500 meV$\cdot$Å2 by inelastic neutron scattering feco_Lowde . Considering success on Fe, Ni, and Co, this FeCo was the case that we could expect a good agreement with experiments. We have not yet found a reason why QSGW gives such discrepancy from the experiment.

IV Summary

In order to calculate SW dispersion in QSGW, we have implemented an effective numerical method for calculating $R^{+-}({\bf q},\omega)$ in a package $ecalj$. This is in the linear response formulation based on the maximally localized Wannier functions as given in Ref. LR_Sasioglu, .

Then we apply the method to Fe, Ni, Co, and FeCo. We compare peak of ${\rm Im}[R^{+-}({\bf q},\omega)]$ with inelastic neutron scattering data and with the spin stiffness $D$. For Fe, Ni, and Co, QSGW gives much better agreements with the experiment rather than LDA does. Notably, too large $D$ of Ni in LDA is reduced by half, resulting in a good agreement with the experiment. We see similar agreement for Co in comparison with the neutron scattering data. For FeCo, we have not yet understood why $D$ in QSGW disagree with the experiment.

Such good agreements are owing to the reliable description of the electronic structure in QSGW. QSGW gives a good description of $3d$-bandwidth, ${\Delta}E_{\rm x}$ and magnetic moments, except the case of Ni where we have a too large magnetic moment. Our method developed here is promising in the sense that it covers wide range of materials from metals treated here to transition-metal oxides where LDA can be hardly applicable.

Acknowledgements.

This work was partly supported by the Building of Consortia for the Development of Human Resources in Science and Technology project, implemented by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan. This work was partly supported by JST CREST Grant number JPMJCR1812 and by JSPS KAKENHI Grant Number JP18H05212. T. Kotani thanks to supporting by JSPS KAKENHI Grant Number 17K05499. We also thank the computing time provided by Research Institute for Information Technology (Kyushu University). We want to thank T. Fukazawa for giving us useful comments.