# An implementation of the density functional perturbation theory
in the PAW framework

Xiaoqiang Liu  International Center for Quantum Materials, School of Physics,
Peking University, Beijing 100871, China  Yihao Lin  International Center for
Quantum Materials, School of Physics, Peking University, Beijing 100871, China
Ji Feng<EMAIL_ADDRESS>International Center for Quantum Materials, School
of Physics, Peking University, Beijing 100871, China  Hefei National
Laboratory, Hefei 230088, China

###### Abstract

Quantifying materials’ dynamical responses to external electromagnetic fields
is central to understanding their physical properties. Here we present an
implementation of the density functional perturbation theory for the
computation of linear susceptibilities using the projector augmented-wave
method. The Sternheimer equations are solved self-consistently through a
nested iterative procedure to compute the first-order wavefunctions, from
which the linear susceptibilities are obtained. As a demonstration, we compute
the spin wave spectral functions of two magnetic metals. The computed magnon
spectra for half-metallic CrO2 and a Heusler intermetallic Cu2MnAl show
gapless Goldstone modes when spin rotation symmetry is preserved and display
reasonable agreement with available experimental data. The Landau damping is
computed to be small in CrO2, but significant in Cu2MnAl producing an
asymmetric Lorentzian spectral lineshape. The access to linear
susceptibilities as well as first-order wavefunctions offers a range of novel
possibilities in quantitative understanding of materials’ electronic
properties from ab initio methods.

##  I Introduction

A microscopic understanding of electrical and magnetic characteristics of
materials plays a key role in condensed matter physics, furnishing a unified
insight into a wide range of phenomena. Indeed, the generalized density
response functions (a.k.a. susceptibilities) of a many-electron system to
external electromagnetic fields, [1] in a broad sense encompasses the
information of the usual longitudinal dielectric function and magnetic
permeability, as well as the cross terms for electromagnetic coupling. The
properties of collective excitations, e.g. plasmon and magnon, can also be
procured from the susceptibilities. Therefore, computing the susceptibilities,
which involve both charge and spin degrees of freedom of electrons, is
essential to a full characterization of electronic properties. Though
relatively straightforward for a non-interacting system, computing the
susceptibilities for an interacting many-electron system is a nontrivial task
due to interaction effects.

The Kohn-Sham density-functional theory (DFT) [2] is by far the most widely
employed ground state electronic structure method for materials and molecules.
By mapping the ground state energy to a non-interacting system described by
Kohn-Sham equations, the Kohn-Sham ansatz enables the variational formulation
of the static responses of a many-electron system. The time-dependent (td) DFT
[3] is subsequently developed, in which the electrodynamics is described by td
Kohn-Sham equations. In these theories, the Hartree and exchange-correlation
potentials are functionals of density, and treated as self-consistent fields.
The self-consistent first-order perturbation in td DFT leads to the density
functional perturbation theory (DFPT), [4] which is then the machinery for
linear response calculations leading directly to the full response functions.

DFPT-based methods have been successfully applied to calculate the dielectric
function [5, 6] and phonon dispersion, [7, 8] with the results in quantitative
agreement with experimental observations. In the computations of dielectric
function, the first-order wavefunctions with respect to $\bm{k}$ are solved,
while deformation potential from frozen atomic displacements is used as the
external field to compute phonon dispersions. Dynamical responses to external
electromagnetic fields from DFPT, which account for the screening of both
charge and spin, have attracted considerable interest recently. [9, 10, 11,
12, 13, 14, 15, 16] In addition to quantifying electrical and magnetic
properties of materials, the linear susceptibilities computed from DFPT also
find applications in the $GW$ approximations. [17, 18, 19]

Approaches to the DFPT can be broadly grouped into two categories: a Dyson-
like equation is solved in the first, [10, 11, 16] whereas the Sternheimer
equations are solved in the second. [9, 12, 13, 14, 15] The Dyson-like
equation approach starts with response functions computed for the Kohn-Sham
ground state. Though formally transparent and amenable to various iterative
techniques, the Dyson-like equation approach suffers from two shortcomings. It
requires a large number of unoccupied states and huge planewave bases for
adequate convergence. [16] More serious is the subtle basis set
incompatibility between the Kohn-Sham states and the DFPT process, which gives
rise to an artifactual spin excitation gap in systems with spin rotation
symmetry. The latter problem can only be partially mended with delicate
engineering of the interaction kernel. [11, 20] In the second category, the
first-order wavefunctions are procured (often iteratively) by solving the
Sternheimer equations, from which the density is updated with charge mixing
iterations and the response functions are computed upon convergence. In this
case, one is forced to deal with wavefunctions and various pseudopotentials
with all the technicalities, [21, 22, 23, 24] in addition to the nested
iterative procedure. Since the full Kohn-Sham response functions are never
required, this method is free of the burden of summation over a huge number of
empty states. In addition, the first-order wavefunctions and densities
computed are actually bonus, which can be useful for computing a variety of
properties.

A particularly popular planewave-based approach to DFT is based on the
projector augmented-wave (PAW) method. [23, 24] Combining the formal
simplicity of pseudopotentials and the versatility of the linearized
augmented-planewave method, the PAW method offers both efficiency and accuracy
to Kohn-Sham DFT calculations on extended solids, and a wide range of
capabilities in various implementations. [25, 26, 27] Despite its popularity,
DFPT in the PAW framework has remained to be developed, which is accomplished
in this work by solving the td Sternheimer equations to compute the linear
susceptibilities of crystalline materials accounting for both charge and spin
degrees of freedom. The paper is organized as follows. In Sec. II.1 the
general theory of DFPT is introduced, from the viewpoint of the dressed spin
in td external electromagnetic fields. The screening built on the notion of
dressed spin leads to the Sternheimer equations in the frequency and momentum
domain and explicit formula for the linear susceptibilities. In Sec. II.2 the
PAW method is reviewed, based on which the formulation of DFPT in the PAW
framework is described, along with a few implementation details. As a first
calibration our implementation, the spin wave spectral functions are extracted
from the computed linear susceptibilities. Two examples are presented: half-
metallic CrO2 (Sec. III.1) shows a clean spin wave spectrum and minimal Landau
damping, and a full Heusler intermetallic Cu2MnAl (Sec. III.2) shows
significant Landau damping in the spin excitations that can be quantified with
a simple asymmetric Lorentzian lineshape. Lastly, a summary is provided with
an eye on room for development from algorithm and physics points of view.

##  II Theory and implementation

###  II.1 Density functional perturbation theory

The td DFT offers an efficient description of the dynamics of an interacting
many-electron system in the presence of external fields. [3] As a self-
consistent perturbation theory of td DFT, DFPT is introduced in this section,
wherein the Sternheimer equations are specialized to crystalline systems under
a monochromatic, periodic external electromagnetic field.

In td DFT, to account for both charge and spin degrees of freedom, the
generalized density is given by

$\displaystyle\rho(\bm{r},t)$
$\displaystyle=\sum_{n}\theta_{n}\operatorname{tr}\\{\psi_{n}^{\dagger}(\bm{r},t)\sigma\psi_{n}(\bm{r},t)\\}$
(1) $\displaystyle=(\rho_{0},\bm{m})=(\rho_{0},\rho_{1},\rho_{2},\rho_{3})$

where $\theta_{n}$ is the occupancy of the spinor single-particle state
$\psi_{n}$, and the four-vector spin
$\sigma=(\sigma_{0},\sigma_{1},\sigma_{2},\sigma_{3})$ ($\sigma_{0}$ is the
identity matrix and $\sigma_{\alpha}$ with $\alpha=1,2,3$ are the Pauli
matrices). The atomic units 111$e=\hbar=m_{e}=1$ are adopted in this paper, so
$\rho_{0}$ is the total charge density, and $\bm{m}$ is magnetization density.
The dynamics of $\psi_{n}$ is prescribed by the td Kohn-Sham equations [3]

$\mathrm{i}\partial_{t}|\psi_{n}(t)\rangle=[H+\delta
H(t)]|\psi_{n}(t)\rangle.$ (2)

The ground state Kohn-Sham Hamiltonian [2] in Eq. (2) is
$H=-\frac{1}{2}\nabla^{2}+v[\rho^{(0)}](\bm{r})$ where the self-consistent
potential is a functional of the ground state density $\rho^{(0)}$, composed
of ionic, Hartree and exchange-correlation (xc) potentials, namely,
$v^{\text{i}}$, $v^{\text{H}}$ and $v^{\text{xc}}$. Though the xc potential
$v^{\text{xc}}$ as a functional of density is in principle nonlocal in space
and time, [29, 30, 31, 32, 33] the commonly adopted local and adiabatic
approximation (ALDA) is assumed in this work. [34, 35, 36]

The first-order Hamiltonian $\delta H$ comprises two contributions. The first
arises from the coupling of the four-vector spin $\sigma$ with external fields

$v^{\text{ext}}(\bm{r},t)=-B_{\alpha}(\bm{r},t)\sigma_{\alpha},$ (3)

where the four-vector electromagnetic field
$B(\bm{r},t)=(-\phi,\frac{1}{2}\bm{B})$. 222Here, $\bm{B}$ is the Zeeman
field. We ignore the nonlocal coupling between the magnetic field and
orbitals, which requires the current density functional theory and is beyond
the scope of the present work. The indices $\alpha,\beta=0,1,2,3$ are
implicitly summed over when repeated, but we will keep other summations
explicit. The self-consistent inclusion of the density dependence in the
Hartree and xc potentials means that $\delta H$ also includes a second
contribution from induced density $\delta\rho(\bm{r},t)$ that screens
$v^{\text{ext}}$. In the adiabatic linear response theory, this can be
formulated in terms of a _dressed spin_ $\tau_{\alpha}$

$\tau_{\alpha}(\bm{r},\bm{r}^{\prime},t-t^{\prime})=-\frac{\delta
H(\bm{r},t)}{\delta
B_{\alpha}(\bm{r}^{\prime},t^{\prime})}=\sigma_{\alpha}\delta(\bm{r}-\bm{r}^{\prime})\delta(t-t^{\prime})\\\
-\sigma_{\gamma}\int
f_{\gamma\beta}(\bm{r},\bm{r}^{\prime\prime})\chi_{\beta\alpha}(\bm{r}^{\prime\prime},\bm{r}^{\prime},t-t^{\prime})\text{d}\bm{r}^{\prime\prime}.$
(4)

Here $\chi_{\alpha\beta}(\bm{r},\bm{r}^{\prime},t)$ is the linear
susceptibility that we pursue in this work, defined via

$\delta\rho_{\alpha}(\bm{r},t)=\int\\!\chi_{\alpha\beta}(\bm{r},\bm{r}^{\prime},t-t^{\prime})B_{\beta}(\bm{r}^{\prime},t^{\prime})\text{d}\bm{r}^{\prime}\text{d}t^{\prime}.$
(5)

The interaction kernel has two components,
$f_{\alpha\beta}=f^{\text{H}}_{\alpha\beta}+f^{\text{xc}}_{\alpha\beta}$,
namely, the Hartree and xc kernels in the ALDA

$f_{\alpha\beta}(\bm{r},\bm{r}^{\prime})=\frac{\delta_{\alpha 0}\delta_{\beta
0}}{|\bm{r}\\!-\\!\bm{r}^{\prime}|}+\frac{1}{2}\delta(\bm{r}-\bm{r}^{\prime})\operatorname{tr}\left[\sigma_{\alpha}\frac{\partial
v^{\text{xc}}}{\partial\rho_{\beta}}\right].$ (6)

In terms of the dressed spin, the first-order Hamiltonian in unscreened
external fields is written

$\delta
H(\bm{r},t)=-\int\\!\tau_{\alpha}(\bm{r},\bm{r}^{\prime},t-t^{\prime})B_{\alpha}(\bm{r}^{\prime},t^{\prime})\text{d}\bm{r}^{\prime}\text{d}t^{\prime}.$
(7)

Now with the first-order Hamiltonian, the first-order wavefunctions can be
obtained by solving the Sternheimer equations [38, 39, 9, 8, 40, 12]

$(\mathrm{i}\partial_{t}-H)|\psi^{(1)}_{n}(t)\rangle=\delta
H(t)|\psi_{n}^{(0)}(t)\rangle,$ (8)

in which $\psi_{n}^{(\ell)}(t)$ is the $\ell$th-order wavefunction, with
$|\psi_{n}^{(0)}(t)\rangle=e^{-\text{i}\varepsilon_{n}t}|\psi_{n}\rangle$.
With the first-order wavefunctions, we will be able to compute the induced
density via the variation of Eq. (1), from which the first-order Hamiltonian
will be updated.

Figure 1: Flow chart for a nested-loop DFPT calculation in the Sternheimer
equation approach. The outer loop depicted here is the charge mixing. The
inner loop is incurred in solving the Sternheimer equations in each step of
the outer iteration.

With the above introduction, a DFPT calculation can then be performed in a
nested iterative process depicted in Fig. 1. In the initializing step, $\delta
H$ is constructed from the external electromagnetic fields, whence the
Sternheimer equations are solved in the inner iteration. This produces a set
of first-order wavefunctions, from which the induced density is calculated. A
charge mixing strategy [41, 42] is employed to revise the induced density,
with the linear susceptibility and dressed spin computed subsequently. Then a
new $\delta H$ is constructed from the updated spin to enter a second round of
the outer iteration. The above process is repeated till convergence.

Now we will explain how the Sternheimer equations are solved in a crystalline
solid. For electrons in a crystal, the initial Kohn-Sham states are Bloch
functions such that
$|\psi_{n}\rangle\mapsto|\psi_{n\bm{k}}\rangle=e^{\text{i}\bm{k}\cdot
r}|u_{n\bm{k}}\rangle,$ where $n$ is the band index and $\bm{k}$ is the
crystal momentum, and $u_{n\bm{k}}$ is the cell-periodic part of the Bloch
function. The system is subject to spatially periodic and monochromatic
external electromagnetic fields

$B_{\alpha}(\bm{r},t)=\mathcal{B}_{\alpha}e^{\mathrm{i}(\bm{p}\cdot\bm{r}-\omega
t)}+\mathrm{c.c.},$ (9)

where $\bm{p}=\bm{q}+\bm{g}$ with $\bm{g}$ being a reciprocal lattice vector
for $\bm{q}$ to in the first Brillouin zone. In this case, the following
expansion is useful

$\zeta(\bm{r},t)=\sum_{l}e^{\mathrm{i}(\bm{l}\cdot\bm{r}-\nu
t)}\zeta(\bm{r},l)$ (10)

for $\zeta=\delta H,\delta\rho_{\alpha}$, in which $\zeta(\bm{r},l)$ is a
cell-periodic function with

$l\equiv(\nu,\bm{l})=\pm(\omega,\bm{q}).$

Then the first-order Hamiltonian in the $l$ channel is

$\delta
H(\bm{r},l)=-e^{-\mathrm{i}\bm{l}\cdot\bm{r}}\tau_{\alpha}(\bm{r},\operatorname{sign}(\nu\omega)\bm{p},\nu)\mathcal{B}_{\alpha}(l),$
(11)

where $\mathcal{B}_{\alpha}(\omega,\bm{q})=\mathcal{B}_{\alpha}$,
$\mathcal{B}_{\alpha}(-\omega,-\bm{q})=(\mathcal{B}_{\alpha})^{*}$, and

$\tau_{\alpha}(\bm{r},\bm{p}^{\prime},\nu)=\int\tau_{\alpha}(\bm{r},\bm{r}^{\prime},t)e^{\mathrm{i}(\bm{p}^{\prime}\cdot\bm{r}^{\prime}+\nu
t)}\text{d}\bm{r}^{\prime}\text{d}t.$ (12)

The first-order wavefunction is expanded as

$\displaystyle|\psi_{n\bm{k}}^{(1)}(t)\rangle=\int\frac{\text{d}\nu}{2\pi}e^{-\mathrm{i}(\nu+\varepsilon_{n\bm{k}})t}|\psi^{(1)}_{n\bm{k}}(\nu)\rangle.$
(13)

For the fields in Eq. (9), $|\psi^{(1)}_{n\bm{k}}(\nu)\rangle$ is nonzero only
when $\nu=\pm\omega$ and can be written as

$\displaystyle|\psi_{n\bm{k}}^{(1)}(\nu)\rangle=e^{\mathrm{i}(\bm{k}+\bm{l})\cdot\bm{r}}|u^{(1)}_{n\bm{k}}(l)\rangle.$
(14)

Then the Sternheimer equations become

$\left(\nu+\mathrm{i}\eta+\varepsilon_{n\bm{k}}-H_{\bm{k}+\bm{l}}\right)|u_{n\bm{k}}^{(1)}(l)\rangle=\delta
H(\bm{r},l)|u_{n\bm{k}}\rangle,$ (15)

in which
$H_{\bm{k}}=e^{-\mathrm{i}\bm{k}\cdot\bm{r}}He^{\mathrm{i}\bm{k}\cdot\bm{r}}$
is the Bloch Hamiltonian. A positive infinitesimal $\eta$ is introduced on the
left-hand side as a convergence factor to embody the causality structure of
the linear susceptibility. In actual calculations, $\eta$ takes finite values
to ensure convergence with finite $\bm{k}$-mesh especially for metals and
bestows finite broadening on spectral peaks.

From the first-order wavefunctions, we compute the first-order induced density
as

$\delta\rho_{\alpha}(\bm{r},l)=\sum_{n\bm{k}}\theta_{n\bm{k}}\operatorname{tr}\left\\{u_{n\bm{k}}^{\dagger}(\bm{r})\sigma_{\alpha}u^{(1)}_{n\bm{k}}(\bm{r},l)+{u^{(1)}_{n\bm{k}}}^{\dagger}(\bm{r},-l)\sigma_{\alpha}u_{n\bm{k}}(\bm{r})\right\\}.$
(16)

The linear susceptibility can be extracted from the Fourier components of
$\delta\rho_{\alpha}(\bm{r},l)$

$\delta\rho_{\alpha}(\bm{g}^{\prime},l)=\chi_{\alpha\beta}(\bm{g}^{\prime},\operatorname{sign}(\nu\omega)\bm{g},l)\mathcal{B}_{\beta}(l),$

whence

$\displaystyle\chi_{\alpha\beta}(\bm{g}^{\prime},\bm{g}^{\prime\prime},l)=$
$\displaystyle\int
e^{-\mathrm{i}(\bm{g}^{\prime}+\bm{l})\cdot\bm{r}^{\prime}}\chi_{\alpha\beta}(\bm{r}^{\prime},\bm{r}^{\prime\prime},\nu)e^{\mathrm{i}(\bm{g}^{\prime\prime}+\bm{l})\cdot\bm{r}^{\prime\prime}}\text{d}\bm{r}^{\prime}\text{d}\bm{r}^{\prime\prime}.$
(17)

It is worth mentioning that the Sternheimer equations in (15) admit the
following formal solutions in the $l$ channel

$|u^{(1)}_{n\bm{k}}(l)\rangle=-\mathcal{B}_{\alpha}(l)\sum_{n^{\prime}}\frac{|u_{n^{\prime}\bm{k}+\bm{l}}\rangle\tau_{\alpha}(\bm{k},l)_{n^{\prime}n}}{\nu+\varepsilon_{n\bm{k}}\\!-\\!\varepsilon_{n^{\prime}\bm{k}+\bm{l}}+\mathrm{i}\eta},$
(18)

in which the dressed spin matrix element is

$\tau_{\alpha}(\bm{k},l)_{n^{\prime}n}=\int
u^{\dagger}_{n^{\prime}\bm{k}+\bm{l}}(\bm{r})e^{-\mathrm{i}\bm{l}\cdot\bm{r}}\tau_{\alpha}(\bm{r},\operatorname{sign}(\nu\omega)\bm{p},\nu)u_{n\bm{k}}(\bm{r})\text{d}\bm{r}.$
(19)

Because it requires the summation over a large number of empty states, this
formal solution is not used in practice. 333As will be explained later, in our
implementation the exact solution is used to quickly correct the contributions
of occupied and a small number of low-lying particle states in the iterative
solution of the Sternheimer equations. The screened susceptibility then has
the same expression as the (bare) Kohn-Sham susceptibility, except that the
(unscreened) external fields are now coupled to the dressed spin; that is,

$\chi_{\alpha\beta}(\bm{g}^{\prime},\operatorname{sign}(\nu\omega)\bm{g},l)=-\sum_{nn^{\prime}\bm{k}}(\theta_{n\bm{k}}-\theta_{n^{\prime}\bm{k}+\bm{l}})\frac{\langle
u_{n\bm{k}}|e^{-\mathrm{i}\bm{g}^{\prime}\cdot\bm{r}}\sigma_{\alpha}|u_{n^{\prime}\bm{k}+\bm{q}}\rangle\tau_{\beta}(\bm{k},l)_{n^{\prime}n}}{\nu+\varepsilon_{n\bm{k}}-\varepsilon_{n^{\prime}\bm{k}+\bm{l}}+\mathrm{i}\eta}.$
(20)

In arriving at the last expression, we have used the fact that
$\tau_{\alpha}(\bm{k},l)_{n^{\prime}n}=\tau_{\alpha}(\bm{k}+\bm{l},-l)^{*}_{nn^{\prime}}$.

###  II.2 DFPT with PAW method

The previous Subsection presents a sketch of the DFPT for periodic systems
without recourse to computational details. In practical calculations, however,
various technologies have been developed such that one can perform DFT (and
therefore DFPT) calculations on valence electrons only for crystalline
materials using planewave bases with the aid of pseudopotentials, to yield
satisfactory accuracy with high efficiencies. It is well known that the norm-
conserving pseudopotential [21] requires large planewave bases particularly
for localized orbitals in transition elements, while the application of
ultrasoft pseudopotential [22] is partly limited by the rather laborious
construction. In contrast, the PAW method, combining the pseudopotential and
linearized augmented-plane-wave methods, is free from the above difficulties
and has been used widely. Thus, performing DFPT within the PAW framework [23,
24] for inhomogeneous and td electromagnetic fields is evidently useful though
notably nontrivial. DFPT with the PAW method has been implemented within the
Vienna ab initio simulation package (VASP) [25] for atomic displacements in
the static and long-wavelength limit to calculate zone-center phonon energies.
The extension to inhomogeneous and td electromagnetic fields is accomplished
in this work based on VASP 5.4.4, and a few implementation details warrant
further clarification. Here, we will briefly review the PAW method and
describe how it is used in our DFPT calculations. Although the formalism for
the ground state quantities is identical to those in the literature and
notationally notorious, we feel compelled to provide some of these details,
particularly in view of the time- and position-dependent external fields
involved in our implementation.

The PAW method is based on a linear transformation between the all-electron
(AE) Hilbert space orthogonal to core states and pseudo (PS) Hilbert space.
The AE and PS wavefunctions are related by

$|\psi_{n\bm{k}}\rangle=T|\tilde{\psi}_{n\bm{k}}\rangle,$ (21)

with the linear operator defined as

$T=1+\sum_{i}(|\phi_{i}\rangle-|\tilde{\phi}_{i}\rangle)\langle\tilde{p}_{i}|.$
(22)

The index $i$ is a shorthand encapsulating the atomic site located at
$\bm{R}_{i}$ as well as the quantum numbers $(nlm)$ of the local orbitals and
spin. $\phi$, $\tilde{\phi}$ and $\tilde{p}$ are AE partial waves, PS partial
waves and projector functions, respectively, and should all be understood as
spinors. In order to perform the DFPT calculations on a crystalline solid in
the electromagnetic fields in Eq. (9), the key is to find the cell-periodic
part of each PAW-pseudized quantity, especially the nonlocal ones.

Upon application of the time-independent transformation $T$ to the first-order
wavefunctions the pseudized Sternheimer equations read

$(\mathrm{i}\partial_{t}S-\tilde{H})|\tilde{\psi}^{(1)}_{n\bm{k}}(t)\rangle=\delta\tilde{H}(t)|\tilde{\psi}_{n\bm{k}}^{(0)}(t)\rangle,$
(23)

with $S=T^{\dagger}T$, $\tilde{H}=T^{\dagger}HT$ and
$\delta\tilde{H}(t)=T^{\dagger}\delta{H}(t)T$.

According to the implementation of PAW method in VASP, [24] we have

$\displaystyle S$ $\displaystyle=1+\sum_{ij}\left|\tilde{p}_{i}\right\rangle
q_{ij}\left\langle\tilde{p}_{j}\right|,$ (24) $\displaystyle\tilde{H}$
$\displaystyle=-\frac{1}{2}\Delta+\tilde{v}^{\mathrm{eff}}+\sum_{ij}\left|\tilde{p}_{i}\right\rangle
D_{ij}\left\langle\tilde{p}_{j}\right|,$

in which the nonlocal potential $D_{ij}=\hat{D}_{ij}+\tilde{D}_{ij}$ with
$\tilde{D}_{ij}=D^{1}_{ij}-\tilde{D}^{1}_{ij}$. The quantities $q_{ij}$,
$D^{1}_{ij}$ and $\tilde{D}^{1}_{ij}$ are defined in reference. [24]
Apparently, the local potential $\tilde{v}^{\mathrm{eff}}(\bm{r})$ is a
functional of pseudo density $\tilde{n}(\bm{r})$ and compensation density
$\hat{n}(\bm{r})$, while $\tilde{D}_{ij}$ is a function of density matrix
$\varrho$, i.e.

$\displaystyle\tilde{v}^{\mathrm{eff}}=$
$\displaystyle\tilde{v}^{\mathrm{eff}}[\tilde{n}+\hat{n}],$ (25)
$\displaystyle\hat{D}_{ij}=$
$\displaystyle\int\frac{1}{2}\operatorname{tr}\\{\sigma_{\alpha}\tilde{v}^{\mathrm{eff}}(\bm{r})\\}Q_{ij}^{\alpha}(\bm{r})\text{d}\bm{r},$
$\displaystyle\tilde{D}_{ij}=$ $\displaystyle\tilde{D}_{ij}(\varrho).$

We have hidden the functional dependence on the pseudized core densities in
$\tilde{v}^{\mathrm{eff}}(\bm{r})$, which are kept frozen during the DFPT
calculations. $\tilde{n}(\bm{r})$, $\hat{n}(\bm{r})$ and $\varrho_{ij}$ are
given by, respectively,

$\displaystyle\tilde{n}(\bm{r})=$
$\displaystyle\sum_{n\bm{k}}\theta_{n\bm{k}}\operatorname{tr}\\{\tilde{\psi}_{n\bm{k}}^{\dagger}(\bm{r})\sigma\tilde{\psi}_{n\bm{k}}(\bm{r})\\},$
(26) $\displaystyle\hat{n}(\bm{r})=$
$\displaystyle\sum_{i,j}\varrho_{ij}Q_{ij}(\bm{r}),$
$\displaystyle\varrho_{ij}=$
$\displaystyle\sum_{n\bm{k}}\theta_{n\bm{k}}\langle\tilde{\psi}_{n\bm{k}}|\tilde{p}_{i}\rangle\langle\tilde{p}_{j}|\tilde{\psi}_{n\bm{k}}\rangle.$

Here
$Q_{ij}(\bm{r})=\sum_{L}\operatorname{tr}\\{\sigma\hat{Q}_{ij}^{L}(\bm{r})\\}$
with $\hat{Q}_{ij}^{L}(\bm{r})$ defined to construct the compensation density
$\hat{n}(\bm{r})$ in the reference. [24] It should be noticed that for
$\varrho_{ij}$, only elements with $\bm{R}_{i}=\bm{R}_{j}$ are useful in our
calculations, while $D_{ij}$ are nonzero only when $\bm{R}_{i}=\bm{R}_{j}$.

Now we derive the expression for the first-order Hamiltonian
$\delta\tilde{H}(t)$. For the fields in Eq. (9), the first-order local
densities and potentials follow the same expansion as in Eq. (10), while the
first-order density matrix $\delta\varrho_{ij}$ and nonlocal potential $\delta
D_{ij}$ can be expanded as

$\delta\zeta_{ij}(t)=\sum_{l}e^{\mathrm{i}(\bm{l}\cdot\bm{R}_{i}-\nu
t)}\delta\zeta_{ij}(l).$ (27)

Here, the factor $e^{\mathrm{i}\bm{l}\cdot\bm{R}_{i}}$ in the $l$ channel is
introduced such that $\delta\zeta_{ij}(l)$ is cell-periodic, i.e.,
$\delta\zeta_{ij}(l)=\delta\zeta_{i^{\prime}j^{\prime}}(l)$ if the positions
of atomic site for $i,j$ and $i^{\prime},j^{\prime}$ differ by a lattice
vector.

The contribution of external electromagnetic fields in $\delta\tilde{H}(t)$
can be calculated directly

$\displaystyle\delta\tilde{H}^{\text{ext}}(\bm{r},t)=$
$\displaystyle\sum_{l}e^{-\mathrm{i}\nu
t}\left[e^{\mathrm{i}\bm{l}\cdot\bm{r}}v^{\text{ext}}(\bm{r},l)\right.+\sum_{ij}\left.e^{\mathrm{i}\bm{l}\cdot\bm{R}_{i}}\left|\tilde{p}_{i}\right\rangle
D^{\text{ext}}_{ij}(l)\left\langle\tilde{p}_{j}\right|\right],$ (28)

with

$\displaystyle v^{\text{ext}}(\bm{r},l)=$
$\displaystyle-\mathcal{B}_{\alpha}(l)\sigma_{\alpha}e^{\operatorname{sign}(\nu\omega)\mathrm{i}\bm{g}\cdot\bm{r}},$
(29) $\displaystyle D^{\text{ext}}_{ij}(l)=$
$\displaystyle\langle\phi_{i}|e^{\mathrm{i}\bm{l}\cdot(\bm{r}-\bm{R}_{i})}v^{\text{ext}}(\bm{r},l)|\phi_{j}\rangle-\langle\tilde{\phi}_{i}|e^{\mathrm{i}\bm{l}\cdot(\bm{r}-\bm{R}_{i})}v^{\text{ext}}(\bm{r},l)|\tilde{\phi}_{j}\rangle.$

The contribution to $\delta\tilde{H}(t)$ from the induced densities has a
similar expression as in Eq. (28) and can be calculated via an explicit finite
difference, as $\tilde{H}$ is a functional of $\tilde{n}$, $\hat{n}$ and
$\rho_{ij}$. The first-order densities are found to be

$\displaystyle\delta\tilde{n}(\bm{r},l)=$
$\displaystyle\sum_{n\bm{k}}\theta_{n\bm{k}}\operatorname{tr}\\{{\tilde{u}_{n\bm{k}}}^{\dagger}(\bm{r})\sigma\tilde{u}^{(1)}_{n\bm{k}}(\bm{r},l)+\tilde{u}^{(1),\dagger}_{n\bm{k}}(\bm{r},-l)\sigma\tilde{u}_{n\bm{k}}(\bm{r})\\},$
(30) $\displaystyle\delta\hat{n}(\bm{r},l)=$
$\displaystyle\sum_{i,j}e^{\mathrm{i}\bm{l}(\bm{R}_{i}-\bm{r})}\delta\varrho_{ij}(l)Q_{ij}(\bm{r}),$
$\displaystyle\delta\varrho_{ij}(l)=$
$\displaystyle\sum_{n\bm{k}}\theta_{n\bm{k}}[\langle\tilde{u}_{n\bm{k}}|\tilde{p}_{i\bm{k}}\rangle\langle\tilde{p}_{j\bm{k}+\bm{l}}|\tilde{u}_{n\bm{k}}^{(1)}(l)\rangle+\langle\tilde{u}_{n\bm{k}}^{(1)}(-l)|\tilde{p}_{i\bm{k}-\bm{l}}\rangle\langle\tilde{p}_{j\bm{k}}|\tilde{u}_{n\bm{k}}\rangle],$

where we define
$|\tilde{p}_{i\bm{k}}\rangle=e^{-\mathrm{i}\bm{k}\cdot(\bm{r}-\bm{R}_{i})}|\tilde{p}_{i}\rangle$.
To linear order in external fields, the first order effective local potential
is given by

$\delta\tilde{v}^{\mathrm{eff}}(l)\approx\tilde{v}^{\mathrm{eff}}[\tilde{n}+\hat{n}+\delta\tilde{n}(l)+\delta\hat{n}(l)]-\tilde{v}^{\mathrm{eff}}[\tilde{n}+\hat{n}].$
(31)

Similarly, the first-order nonlocal potentials can be approximated as

$\displaystyle\delta\hat{D}_{ij}(l)\approx$ $\displaystyle\int
e^{\mathrm{i}\bm{l}\cdot(\bm{r}-\bm{R}_{i})}\frac{1}{2}\operatorname{tr}\\{\sigma_{\alpha}\delta\tilde{v}^{\mathrm{eff}}(\bm{r},l)\\}Q_{ij}^{\alpha}(\bm{r})\text{d}\bm{r},$
(32) $\displaystyle\delta\tilde{D}_{ij}(l)\approx$
$\displaystyle\tilde{D}_{ij}(\varrho+\delta\varrho(l))-\tilde{D}_{ij}(\varrho).$

Though introduced as forward differences in Eq.(31) and (32), these quantities
are evaluated using 4th-order centered finite differences, with a step length
of a thousandth the density variables.

With the above results, the final Sternheimer equations become

$\left(\nu
S_{\bm{k}+\bm{l}}+\varepsilon_{n\bm{k}}S_{\bm{k}+\bm{l}}-\tilde{H}_{\bm{k}+\bm{l}}\right)|\tilde{u}_{n\bm{k}}^{(1)}(l)\rangle=\delta\tilde{H}_{\bm{k}}(l)|\tilde{u}_{n\bm{k}}\rangle,$
(33)

with

$\displaystyle S_{\bm{k}+\bm{l}}$
$\displaystyle=1+\sum_{ij}\left|\tilde{p}_{i\bm{k}+\bm{l}}\right\rangle
q_{ij}\left\langle\tilde{p}_{j\bm{k}+\bm{l}}\right|,$ (34)
$\displaystyle\tilde{H}_{\bm{k}+\bm{l}}$
$\displaystyle=-\frac{1}{2}\Delta_{\bm{k}+\bm{l}}+\tilde{v}^{\mathrm{eff}}+\sum_{ij}\left|\tilde{p}_{i\bm{k}+\bm{l}}\right\rangle
D_{ij}\left\langle\tilde{p}_{j\bm{k}+\bm{l}}\right|,$
$\displaystyle\delta\tilde{H}_{\bm{k}}(l)$
$\displaystyle=v^{\text{ext}}(l)+\delta\tilde{v}^{\mathrm{eff}}(l)+\sum_{ij}\left|\tilde{p}_{i\bm{k}+\bm{l}}\right\rangle[D^{\text{ext}}_{ij}(l)+\delta\hat{D}_{ij}(l)+\delta\tilde{D}_{ij}(l)]\left\langle\tilde{p}_{j\bm{k}}\right|.$

Here $\tilde{u}_{n\bm{k}}$ and $\tilde{u}_{n\bm{k}}^{(1)}(l)$ are the cell-
periodic parts of corresponding pseudo wavefunctions, respectively.

The pseudized Sternheimer equations in $\pm(\omega,\bm{q})$ channels are
solved separately in each iteration using a variant of residual minimization
method with a direct inversion in the iterative subspace (RMM-DIIS), [44, 45]
which is already implemented in VASP 5.4.4. The Löwdin perturbation theory is
also performed to correct the first-order wavefunctions in the subspace of
occupied states and low-lying excitations to speed up convergence

$|\tilde{u}^{(1)}_{n\bm{k}}(l)\rangle\rightarrow|\tilde{u}^{(1)}_{n\bm{k}}(l)\rangle-\sum_{n^{\prime}}|\tilde{u}_{n^{\prime}\bm{k}+\bm{l}}\rangle\langle\tilde{u}_{n^{\prime}\bm{k}+\bm{l}}|S_{\bm{k}+\bm{l}}|\tilde{u}^{(1)}_{n\bm{k}}(l)\rangle+\sum_{n^{\prime}}\frac{|\tilde{u}_{n^{\prime}\bm{k}+\bm{l}}\rangle\langle\tilde{u}_{n^{\prime}\bm{k}+\bm{l}}|\delta\tilde{H}|\tilde{u}_{n\bm{k}}\rangle}{\nu+\varepsilon_{n\bm{k}}-\varepsilon_{n^{\prime}\bm{k}+\bm{l}}}.$
(35)

In the last equation above, the summations on $n^{\prime}$ run over the
occupied bands plus a few empty bands.

Apparently, solving Eq. (33) requires a $\bm{k}$-grid supplemented by two
additional grids shifted by $\pm\bm{q}$ when $\bm{q}$ itself is not on the
$\bm{k}$-grid. Doing so, however, not only increases the computational burden
but also, more seriously, obliterates the exact cancellation of the
contribution of the occupied manifold to the density change due to a
$\bm{k}$-grid discretization error. The latter can be easily avoided by
employing a pair of grids with a $\bm{q}$ shift, which also reduces the
calculation partly. Then Eq. (33) is solved on the $\bm{k}$-grid in $+\bm{q}$
channel, and on the $\bm{k}+\bm{q}$ grid in $-\bm{q}$ channel. It is observed
that in this dual grid setup, the above cancellation is well preserved.

The xc potentials in ALDA are functionals of real-valued densities. Thus,
calculating $\delta\tilde{v}^{\mathrm{eff}}(l)$ like in Eq. (31) requires
caution as $\delta\tilde{n}(l)$ and $\delta\hat{n}(l)$ are usually complex. In
fact, the real and imaginary part of $\delta\tilde{v}^{\mathrm{eff}}(l)$ are
calculated separately,

$\mathfrak{F}\delta\tilde{v}^{\mathrm{eff}}(l)\approx\tilde{v}^{\mathrm{eff}}[\tilde{n}+\hat{n}+\mathfrak{F}\delta\tilde{n}(l)+\mathfrak{F}\delta\hat{n}(l)]-\tilde{v}^{\mathrm{eff}}[\tilde{n}+\hat{n}]$
(36)

where $\mathfrak{F}=\operatorname{Re},\operatorname{Im}$ takes the real or
imaginary part, respectively. In the case of nonlocal potential,
$\delta\tilde{D}_{ij}(l)$ and $\delta\varrho_{ij}(l)$ are first decomposed
into two independent Hermitian matrices (i.e. Hermitian part and anti-
Hermitian part multiplied by $-\mathrm{i}$), and then finite differenced
separately in an analogous fashion.

Symmetry reduction is also performed in our implementation, where the
summation over $\bm{k}$ points in Eq. (30) is restricted to the symmetry-
irreducible part of the Brillouin zone. The symmetry group here is the
subgroup of the magnetic group of the studied crystal in which the external
electromagnetic fields in Eq. (9) is invariant.

##  III Application to spin-wave spectrum calculation

Our implementation enables computing the linear susceptibilities
$\chi_{\alpha\beta}$ with the self-consistent inclusion of the interaction
kernel. Directly inverting $\chi_{\alpha\beta}$ yields the dielectric tensor,
which is composed of the usual charge sector $\epsilon_{00}$, spin sector
$\epsilon_{\alpha\beta}$, and the spin-charge sector $\epsilon_{0\beta}$, each
embodying unique physics. Computing $\chi_{\alpha\beta}$ then can have diverse
applications in evaluating materials properties pertaining to both charge and
spin fluctuations, or in subsequent many-body calculations beyond the Kohn-
Sham mean fields. One immediate application that has received considerable
attention is the calculation of spin-wave excitation. [9, 10, 11, 12, 13, 14,
15, 16] According to the fluctuation-dissipation theorem, the spin-spin
correlation function, directly accessible by various spin-sensitive inelastic
scattering probes, [46, 47, 48] is related to the imaginary part of the linear
susceptibilities,

$S_{\scriptscriptstyle+-}(\bm{p},\omega)=\frac{\operatorname{Im}\chi_{\scriptscriptstyle+-}(\bm{g},\bm{g},\omega,\bm{q})}{1-e^{-\hbar\omega/k_{\mathrm{B}}T}}.$
(37)

[b]  Authors (Year) Potential Basis set Software Savrasov (1998) [9] Full
potential LMTO1 LMTO Magnons Cao et al. (2018) [12] NCPP2, USPP3 Planewave QE4
Gorni et al. (2018) [13] NCPP Planewave QE Tancogne-Dejean et al. (2020) [14]
NCPP Real space Octopus Singh et al. (2020) [15] Full potential LAPW5 Elk
1Linear muffin-tin orbital. 2Norm-conserving pseudopotential. 3Ultrasoft
pseudopotential. 4Quantum ESPRESSO. 5Linearised augmented-plane wave.

Table 1: Reported implementations of DFPT for spin wave spectra calculations
by solving the Sternheimer equations.

Although for magnetic systems dominated by local moments the magnon can be
described effectively by localized spin models, this method is subject to
debate when delocalization sets in, and ultimately of questionable validity
for itinerant magnetism. In these latter cases, which include a wide range of
magnetic materials, the DFPT route becomes invaluable for computing the spin
wave spectra ab initio. To our knowledge, there have been just a handful of
works devoted to implementing DFPT scheme for this purpose by solving the
Sternheimer equations, as summarized in Table 1. In these efforts,
implementations are limited to the full-potential, [9, 15] or norm-conserving
and ultrasoft pseudopotentials. [12, 13, 14]

In this section, we present an initial application of our implementation of
DFPT in the PAW framework to the calculations of spin-wave spectra for a
couple of magnetic materials. For ferromagnets with the spin polarized along
$z$ direction, transverse magnetic field can be applied by choosing
$\mathcal{B}=(0,1,-\mathrm{i},0)$ in Eq. (9), from which
$\chi_{\scriptscriptstyle+-}$ is calculated directly. In these cases, the only
remaining symmetry operation keeping the crystal and the transverse magnetic
field unchanged is identity transformation. Thus, there is no room for
symmetry reduction. For notationaly convenience, we define
$\chi_{\scriptscriptstyle+-}(\bm{p},\omega)\equiv\chi_{\scriptscriptstyle+-}(\bm{g},\bm{g},\omega,\bm{q})$
as $\bm{p}=\bm{q}+\bm{g}$.

###  III.1 Half-metallic chromium dioxide

As shown in the inset in Fig. 2, chromium dioxide, CrO2, is a ferromagnetic
half-metallic oxide with a rutile crystal structure, where each Cr atom is
situated at the center of an octahedral cage formed by oxygen atoms. [49, 50]
Widely used as a magnetic recording material, CrO2 also has various potential
applications in spintronics and magnetoelectronics [51, 52] due to its half-
metallic properties.

Figure 2:  Spin-resolved densities of states of CrO2. The spin-flip gap is
found to be about 310 meV. Inset: the crystal structure of CrO2, highlighting
the CrO6 octahedra.

The experimental lattice parameters $a=b=4.4218~{}$Å and $c=2.9182~{}$Å [50]
are used in our calculations. The planewave energy cutoff is set to be $500$
eV and a $13\times 13\times 20$ $\Gamma$-centered mesh of $\bm{k}$-points is
used. The spin-resolved density of states of CrO2 is computed from the
collinear spin-polarized calculation and shown in Fig. 2, where the half-
metallic band structure is clearly seen. The magnetic moment of Cr is found to
be 2 $\mu_{\mathrm{B}}$. We then turn to the noncollinear calculations, and
compute the transverse spin susceptibility
$\chi_{\scriptscriptstyle+-}(\bm{p},\omega)$ along [100] and [001] directions
for $\omega\leq 400$ meV on 10 meV intervals. The broadening parameter $\eta$
introduced in Eq. (15) is set to be 50 meV in the calculations in this
subsection.

Figure 3: (a) Im$\chi_{\scriptscriptstyle+-}(\bm{p},\omega)$ of CrO2 along
[100] and [001] directions. The black vertical line indicates the center of
the first Brillouin zone. (b) The folded magnon energy dispersion
$\omega_{\bm{q}}$ of CrO2 in the first Brillouin zone, extracted from
$\operatorname{Im}\chi_{\scriptscriptstyle+-}$ shown in (a). The inset shows
the quadratic fits to $\omega_{\bm{q}}$ at small $q$ along [100] and [001]
directions, respectively, for calculations without Hubbard $U$ correction and
with $U^{\mathrm{eff}}=2.1$ eV. The squares are the calculated data.

Fig. 3(a) shows the computed Im$\chi_{\scriptscriptstyle+-}(\bm{p},\omega)$,
without spin-orbit interaction, along two $\bm{p}$ paths. In general,
$\chi_{\scriptscriptstyle+-}(\bm{p},\omega)$ is not periodic in $\bm{p}$. Then
the branch in the first Brillouin Zone is composed of acoustic magnon modes,
while the branch in the second Brillouin zone optical modes. The profile of
the magnon peak at a given $\bm{p}$ is nearly perfect Lorentzian over the
entire energy range. The extracted half width at half maximum $\eta_{\bm{p}}$
are almost a constant and equal to the artificial broadening parameter $\eta$,
indicating that the Landau damping in CrO2 is negligible. This is expected
given that half-metallic CrO2 has a large spin-flip gap around 310 meV, as
shown in Fig. 2.

The location of the maxima of magnon peaks, $\omega_{\bm{p}}$, are then
recorded and folded to the first Brillouin zone ($\omega_{\bm{q}}$). As shown
in Fig. 3(b), we find one acoustic magnon branch and one optical branch,
consistent with the fact that the unit cell in CrO2 contains two magnetic Cr
atoms. There is no magnon gap at the Brillouin zone boundaries, which is
result of the $n$ glide symmetry. The energy of long-wavelength acoustic
magnons is quadratic in $q$ with a gapless Goldstone mode,
$\omega_{\bm{q}}=D_{\parallel}(q_{x}^{2}+q_{y}^{2})+D_{z}q_{z}^{2}$, as
expected for a ferromagnet with spin rotation symmetry. The spin stiffness
coefficients are found to be $D_{\parallel}=82$ meV$\cdot$Å2 along [100] and
$D_{z}=92$ meV$\cdot$Å2 along [001] directions, respectively. From this, we
can estimate the average spin stiffness coefficient to be 85 meV$\cdot$Å2,
which is close to the experimental measured result ($\sim 112.5$ meV$\cdot$Å2
[53]).

For comparison, additional electron correlation is included statically within
the LSDA+U formalism [54] in both the ground state calculation and subsequent
DFPT calculations, with $U^{\mathrm{eff}}=2.1$ eV. [55] As shown in the inset
of Fig. 3(b), the gapless Goldstone mode is obtained again, accompanied by an
average spin stiffness coefficient $D=391$ meV$\cdot$Å2, almost five times as
large as the one without Hubbard $U$ correction. The energy of magnon in CrO2
seems to be highly overestimated in LSDA+U calculations.

As a further test, we examine the Goldstone gap as a result of breaking the
spin rotation symmetry by introducing spin-orbit interaction. The atomic spin-
orbit interaction for Cr is fairly weak (on the order of tens of meV). Since
the gap in the Goldstone magnon is second order in the spin-orbit coupling, it
is small for CrO2. In order to visualize the effect of spin-orbit coupling, we
introduce a parameter $\lambda$ to artificially tune its strength (or speed of
light), as in $H=H^{0}+\lambda H^{\mathrm{soc}},$ where $\lambda=1$
corresponds to the actual strength of spin-orbit coupling in CrO2. The
calculated Im$\chi_{\scriptscriptstyle+-}(p=0,\omega)$ as a function of
$\omega$ for different $\lambda$ values are shown in Fig. 4(a). Apparently
there is blue shift of the Goldstone mode with increasing $\lambda$,
indicating the emergence of a Goldstone gap. The Goldstone gap indeed shows a
quadratic dependence on $\lambda$, as demonstrated by the gap-vs-$\lambda^{2}$
plot in Fig. 4(b). The extrapolated Goldstone gap in CrO2 is about 0.1 meV.

Figure 4: (a) Im$\chi_{\scriptscriptstyle+-}(p=0,\omega)$ as a function of
$\omega$ for different $\lambda$ values. $\lambda=1$ corresponds to the actual
strength of spin-orbit coupling in CrO2. The squares are the calculated data
and the solid lines are the fits with Lorentzian line shape. (b) The Goldstone
gap as a function of $\lambda^{2}$ in CrO2. The solid line is the linear fit.

###  III.2 Heusler intermetallic Cu2MnAl

The ternary intermetallic Cu2MnAl is a Mn-based full-Heusler alloy with the
$L2_{1}$ structure type (see inset in Fig. 5). The experimental lattice
parameter for the conventional cubic cell (space group $Fm\bar{3}m$) is
$a=5.968$ Å. [56] Cu2MnAl is ferromagnetic below the relatively high Curie
temperature (603 K). [56] Apart from being regarded as a prototype for
understanding the electronic correlations in Heusler intermetallics, [57]
Cu2MnAl is also being used as a neutron polarizer and monochromator material.
[58, 59]

Figure 5: Spin-resolved densities of states of Cu2MnAl. Inset: the crystal
structure of full Heusler Cu2MnAl, showing a conventional cubic unit cell for
the $L2_{1}$ structure.

A planewave energy cutoff of $350$ eV and a $15\times 15\times 15$
$\Gamma$-centered $\bm{k}$-grid are used in our calculations. Spin-orbit
coupling is not included. The spin-resolved densities of states of Cu2MnAl
computed from the collinear spin-polarized calculation confirms the
ferromagnetism of Cu2MnAl, as shown in Fig. 5. The magnetic moment is carried
primarily by Mn atoms and computed to be 3.4 $\mu_{\mathrm{B}}$/Mn. The
transverse spin susceptibility $\chi_{\scriptscriptstyle+-}(\bm{p},\omega)$
along [100], [110] and [111] crystallographic directions is then computed in
noncollinear calculations for $\omega\leq 300$ meV on 5 meV intervals, with a
broadening parameter $\eta$ of 50 meV.

Fig. 6(a) shows the computed magnon spectral function
Im$\chi_{\scriptscriptstyle+-}(\bm{p},\omega)$ along the three principal
directions. The acoustic magnon branch is seen clearly only at low energies
near the Brillouin zone center. The spectral peaks of these low-energy modes
can be adequately fitted with the Lorentzian lineshape as in the CrO2 case.
The dispersion of the long-wavelength modes is quadratic and isotropic, as
demonstrated in the inset of Fig. 7 (a). A spin stiffness coefficient $D=268$
meV$\cdot$Å2 is procured from the quadratic fit, which is about 1.5 times
larger than the experimental value of 175 meV$\cdot$Å2. [60]

Figure 6: (a) Im$\chi_{\scriptscriptstyle+-}(\bm{p},\omega)$ of Cu2MnAl along
[100], [110] and [111] directions. Reciprocal lattice vectors of conventional
cubic cell are adopted here. (b,c)
Im$\chi_{\scriptscriptstyle+-}(\bm{p},\omega)$ as a function of $\omega$ for
$\bm{p}=(1,0,0)$, $(1,1,0)$ and $(1,1,1)$. The squares are the calculated
data. The solid lines are the fits with (b) symmetric or (c) asymmetric
Lorentzian function.

Notably, at higher energies and near Brillouin zone boundaries, the magnon
peaks become fuzzier and broader, attesting to substantial Landau damping in
this material. In stark contrast to the CrO2 case with almost no Landau
damping, the coupling to the Stoner continuum in Cu2MnAl bestows a magnon peak
at a given $\bm{p}$ an asymmetric profile that defies a Lorentzian fit, as
shown in Fig. 6(b). The stronger Landau damping in this system is consistent
with the absence of the spin-flip gap as shown in Fig.5. Viewing the coupling
to the Stoner continuum as a Fano-type resonance, we superimpose a linear
function on the Lorentzian to describe the asymmetric line shape, as

$A(\bm{p},\omega)=\frac{a_{\bm{p}}\eta_{\bm{p}}}{(\omega-\omega_{\bm{p}})^{2}+\eta_{\bm{p}}^{2}}+\xi_{\bm{p}}(\omega-\omega_{\bm{p}})$
(38)

with $\omega_{\bm{p}},a_{\bm{p}},\eta_{\bm{p}}$ and $\xi_{\bm{p}}$ as fitting
parameters.

As it turns out, this simple modification leads to satisfactory fitting for
the entire spectrum, as evidenced in Fig.6(c). The extracted magnon dispersion
is then shown in Fig. 7(a), which coincides with the calculated results of
Buczek et al. [10, 61] and agrees well with the experimental observations
along [100] direction. [60] Along the [110] and [111] directions where the
Landau damping seems more pronounced, our computed dispersion shows
significant discrepancy from the experimental one. For low-energy modes,
$\eta_{\bm{p}}$ is dominated by the artificial broadening parameter $\eta$ and
the asymmetry is small, as shown in Fig. 7(b). With increasing energies,
however, the broadening quickly exceeds $\eta$ and the asymmetry becomes
pronounced, especially near Brillouin zone boundaries, both providing
quantitative characterization of the Landau damping.

Figure 7: (a) The magnon energy dispersion $\omega_{\bm{p}}$ of Cu2MnAl along
[100], [110] and [111] directions, extracted from the asymmetric Lorentzian
fits. $\triangle$ is the experimental data of Tajima et al. [60] Inset is the
quadratic fit to the long-wavelength acoustic magnons. (b) Broadening
$\eta_{\bm{p}}$ and asymmetry $\xi_{\bm{p}}$ of the magnon peaks in Cu2MnAl
along [100], [110] and [111] directions. The error bars show 95% confidence
bounds on the fitting parameters from the fits.

##  IV Summary and outlook

In conclusion, we report an implementation of the DFPT method in the PAW
framework, which is capable of computing the full linear susceptibilities of
real materials. A nested iterative procedure is employed to self-consistently
solve the Sternheimer equations, to procure linear susceptibilities along with
the first-order wavefunctions and densities in monochromatic and periodic
external electromagnetic fields. The time cost of each DFPT calculation (given
an external field direction, momentum and frequency) is comparable to that of
a corresponding Kohn-Sham DFT calculation.

As a demonstration, we compute the spin wave spectra for CrO2 and Cu2MnAl.
Gapless magnon dispersion is obtained for both materials from the calculations
without spin-orbit coupling. The spin stiffness coefficient extracted from the
quadratic fit is in agreement with experimental value for CrO2 but 1.5 times
larger for Cu2MnAl. The Landau damping in CrO2 is insignificant due to its
half-metallic nature, while in Cu2MnAl is remarkable at high energies and can
be quantified with a simple asymmetric Lorentzian fit. LSDA+U method as well
as the effect of spin-orbit coupling are examined in CrO2, from which the
former highly overestimates the magnon energy, while the latter gives rise to
a Goldstone gap quadratic in spin-orbit coupling strength $\lambda$.

There is clearly room for future developments, to make the current
implementation more efficient and versatile. From an algorithm viewpoint, the
occupied subspace is not projected out in Sternheimer equations in the current
implementation. As the contribution to the first-order wavefunctions from the
occupied states does not contribute to the first-order densities, projecting
out the occupied subspace [8] can potentially improve the efficiency and
stability of the nested iteration. As an additional benefit of the projection,
it also renders the principle integrals explicit and amenable to analytic
techniques, which can further reduce the number of $\bm{k}$-points required
and improve efficiencies. Alternative iterative techniques should be tested in
general, for both inner and outer loops, especially in conjunction with the
projection.

From a physics viewpoint, a few tasks are on immediate agenda and new
possibilities are clearly on the horizon, beyond the initial demonstrations
presented herein. For the spin-wave spectral functions, it will be valuable to
compare the computed spectra with experimental results for more materials. A
particularly interesting comparison can be made between the dispersion
relations obtained ab initio from our DFPT implementation and those from
Heisenberg models parametrized from constrained DFT energies on the basis of
the magnetic force theorem. [62, 63, 64] Such comparisons should be examined
in detail for materials in the localized and the itinerant limits, as well as
for the continuum falling in between. Further systematic studies for the
gradient correction (as in generalized gradient approximations) and for the
Hubbard correction in LSDA+U method can reveal the effect of correlation on
the spin-wave spectra. As the first-order wavefunctions are also produced in
our code, it is also tempting to evaluate other physical properties, related
to density and current responses, such as the magnetoelectric coupling and
related transport coefficients. A particular connection may be made by
observing that

$W=f+f\chi f$ (39)

is the screened kernel, which now includes the charge, spin and cross
screening effects. This will enable analyzing the many-electron effects in
magnetic materials with strong spin-orbit coupling, and potentially evaluating
novel bound states from the screened charge/spin interactions.

###### Acknowledgements.

We acknowledge the financial support from the National Natural Science
Foundation of China (Grant No. 11725415), the National Key R&D Program of
China (Grant Nos. 2018YFA0305601 and 2021YFA1400100), and the Innovation
Program for Quantum Science and Technology (Grant No. 2021ZD0302600).

## References

  * Kubo [1957] R. Kubo, Statistical-mechanical theory of irreversible processes. I. general theory and simple applications to magnetic and conduction problems, J. Phys. Soc. Japan 12, 570 (1957).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. 140, A1133 (1965).
  * Runge and Gross [1984] E. Runge and E. K. U. Gross, Density-functional theory for time-dependent systems, Phys. Rev. Lett. 52, 997 (1984).
  * Gross and Kohn [1985] E. K. U. Gross and W. Kohn, Local density-functional theory of frequency-dependent linear response, Phys. Rev. Lett. 55, 2850 (1985).
  * Hybertsen and Louie [1987] M. S. Hybertsen and S. G. Louie, Ab initio static dielectric matrices from the density-functional approach. I. formulation and application to semiconductors and insulators, Phys. Rev. B 35, 5585 (1987).
  * Gajdoš  _et al._ [2006] M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller, and F. Bechstedt, Linear optical properties in the projector-augmented wave methodology, Phys. Rev. B 73, 045112 (2006).
  * Giannozzi  _et al._ [1991] P. Giannozzi, S. de Gironcoli, P. Pavone, and S. Baroni, Ab initio calculation of phonon dispersions in semiconductors, Phys. Rev. B 43, 7231 (1991).
  * Baroni  _et al._ [2001] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. 73, 515 (2001).
  * Savrasov [1998] S. Y. Savrasov, Linear response calculations of spin fluctuations, Phys. Rev. Lett. 81, 2570 (1998).
  * Buczek  _et al._ [2009] P. Buczek, A. Ernst, P. Bruno, and L. M. Sandratskii, Energies and lifetimes of magnons in complex ferromagnets: A first-principle study of heusler alloys, Phys. Rev. Lett. 102, 247206 (2009).
  * Rousseau  _et al._ [2012] B. Rousseau, A. Eiguren, and A. Bergara, Efficient computation of magnon dispersions within time-dependent density functional theory using maximally localized wannier functions, Phys. Rev. B 85, 054305 (2012).
  * Cao  _et al._ [2018] K. Cao, H. Lambert, P. G. Radaelli, and F. Giustino, Ab initio calculation of spin fluctuation spectra using time-dependent density functional perturbation theory, plane waves, and pseudopotentials, Phys. Rev. B 97, 024420 (2018).
  * Gorni  _et al._ [2018] T. Gorni, I. Timrov, and S. Baroni, Spin dynamics from time-dependent density functional perturbation theory, Eur. Phys. J. B 91, 249 (2018).
  * Tancogne-Dejean  _et al._ [2020] N. Tancogne-Dejean, F. G. Eich, and A. Rubio, Time-Dependent Magnons from First Principles, J. Chem. Theory Comput. 16, 1007 (2020).
  * Singh  _et al._ [2020] N. Singh, P. Elliott, J. K. Dewhurst, E. K. U. Gross, and S. Sharma, Ab-initio real-time magnon dynamics in ferromagnetic and ferrimagnetic systems, Phys. Status Solidi B 257, 1900654 (2020).
  * Skovhus and Olsen [2021] T. Skovhus and T. Olsen, Dynamic transverse magnetic susceptibility in the projector augmented-wave method: Application to Fe, Ni, and Co, Phys. Rev. B 103, 245110 (2021).
  * Hybertsen and Louie [1986] M. S. Hybertsen and S. G. Louie, Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies, Phys. Rev. B 34, 5390 (1986).
  * Shishkin and Kresse [2006] M. Shishkin and G. Kresse, Implementation and performance of the frequency-dependent $GW$ method within the paw framework, Phys. Rev. B 74, 035101 (2006).
  * Shishkin and Kresse [2007] M. Shishkin and G. Kresse, Self-consistent $GW$ calculations for semiconductors and insulators, Phys. Rev. B 75, 235102 (2007).
  * Skovhus  _et al._ [2022] T. Skovhus, T. Olsen, and H. M. Rønnow, Influence of static correlation on the magnon dynamics of an itinerant ferromagnet with competing exchange interactions: First-principles study of MnBi, Phys. Rev. Materials 6, 054402 (2022).
  * Hamann  _et al._ [1979] D. R. Hamann, M. Schlüter, and C. Chiang, Norm-conserving pseudopotentials, Phys. Rev. Lett. 43, 1494 (1979).
  * Vanderbilt [1990] D. Vanderbilt, Soft self-consistent pseudopotentials in a generalized eigenvalue formalism, Phys. Rev. B 41, 7892 (1990).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, _et al._ , Advanced capabilities for materials modelling with Quantum ESPRESSO, J. Phys.: Condens. Matter 29, 465901 (2017).
  * Gonze  _et al._ [2020] X. Gonze, B. Amadon, G. Antonius, F. Arnardi, L. Baguet, J.-M. Beuken, J. Bieder, F. Bottin, J. Bouchet, E. Bousquet, _et al._ , The Abinit project: Impact, environment and recent developments, Comput. Phys. Commun. 248, 107042 (2020).
  * Note [1] $e=\hbar=m_{e}=1$.
  * Vignale and Kohn [1996] G. Vignale and W. Kohn, Current-dependent exchange-correlation potential for dynamical linear response theory, Phys. Rev. Lett. 77, 2037 (1996).
  * Vignale  _et al._ [1997] G. Vignale, C. A. Ullrich, and S. Conti, Time-dependent density functional theory beyond the adiabatic local density approximation, Phys. Rev. Lett. 79, 4878 (1997).
  * Kootstra  _et al._ [2000] F. Kootstra, P. L. de Boeij, and J. G. Snijders, Application of time-dependent density-functional theory to the dielectric function of various nonmetallic crystals, Phys. Rev. B 62, 7071 (2000).
  * van Faassen  _et al._ [2002] M. van Faassen, P. L. de Boeij, R. van Leeuwen, J. A. Berger, and J. G. Snijders, Ultranonlocality in time-dependent current-density-functional theory: Application to conjugated polymers, Phys. Rev. Lett. 88, 186401 (2002).
  * Marques  _et al._ [2006] M. A. Marques, C. A. Ullrich, F. Nogueira, A. Rubio, K. Burke, and E. K. U. Gross,  _Time-dependent density functional theory_ (Springer, 2006).
  * Ando [1977] T. Ando, Inter-subband optical absorption in space-charge layers on semiconductor surfaces, Z. Phys. B 26, 263 (1977).
  * Zangwill and Soven [1980a] A. Zangwill and P. Soven, Density-functional approach to local-field effects in finite systems: Photoabsorption in the rare gases, Phys. Rev. A 21, 1561 (1980a).
  * Zangwill and Soven [1980b] A. Zangwill and P. Soven, Resonant photoemission in Barium and Cerium, Phys. Rev. Lett. 45, 204 (1980b).
  * Note [2] Here, $\bm{B}$ is the Zeeman field. We ignore the nonlocal coupling between the magnetic field and orbitals, which requires the current density functional theory and is beyond the scope of the present work.
  * Sternheimer [1954] R. M. Sternheimer, Electronic polarizabilities of ions from the hartree-fock wave functions, Phys. Rev. 96, 951 (1954).
  * de Gironcoli [1995] S. de Gironcoli, Lattice dynamics of metals from density-functional perturbation theory, Phys. Rev. B 51, 6773 (1995).
  * Giustino  _et al._ [2010] F. Giustino, M. L. Cohen, and S. G. Louie, GW method with the self-consistent sternheimer equation, Phys. Rev. B 81, 115105 (2010).
  * Broyden [1965] C. G. Broyden, A class of methods for solving nonlinear simultaneous equations, Math. Comput. 19, 577 (1965).
  * Johnson [1988] D. D. Johnson, Modified Broyden’s method for accelerating convergence in self-consistent calculations, Phys. Rev. B 38, 12807 (1988).
  * Note [3] As will be explained later, in our implementation the exact solution is used to quickly correct the contributions of occupied and a small number of low-lying particle states in the iterative solution of the Sternheimer equations.
  * Pulay [1980] P. Pulay, Convergence acceleration of iterative sequences. the case of scf iteration, Chem. Phys. Lett. 73, 393 (1980).
  * Wood and Zunger [1985] D. M. Wood and A. Zunger, A new method for diagonalising large matrices, J. Phys. A: Math. Gen. 18, 1343 (1985).
  * Mook and Nicklow [1973] H. A. Mook and R. M. Nicklow, Neutron scattering investigation of the magnetic excitations in iron, Phys. Rev. B 7, 336 (1973).
  * Braicovich  _et al._ [2009] L. Braicovich, L. J. P. Ament, V. Bisogni, F. Forte, C. Aruta, G. Balestrino, N. B. Brookes, G. M. De Luca, P. G. Medaglia, F. M. Granozio, M. Radovic, M. Salluzzo, J. van den Brink, and G. Ghiringhelli, Dispersion of magnetic excitations in the cuprate ${\mathrm{La}}_{2}{\mathrm{CuO}}_{4}$ and ${\mathrm{CaCuO}}_{2}$ compounds measured using resonant x-ray scattering, Phys. Rev. Lett. 102, 167401 (2009).
  * Ament  _et al._ [2011] L. J. P. Ament, M. van Veenendaal, T. P. Devereaux, J. P. Hill, and J. van den Brink, Resonant inelastic x-ray scattering studies of elementary excitations, Rev. Mod. Phys. 83, 705 (2011).
  * Schwarz [1986] K. Schwarz, CrO2 predicted as a half-metallic ferromagnet, J. Phys. F 16, L211 (1986).
  * Cloud  _et al._ [1962] W. H. Cloud, D. Schreiber, and K. Babcock, X-ray and magnetic studies of CrO2 single crystals, J. Appl. Phys. 33, 1193 (1962).
  * Singh  _et al._ [2016] A. Singh, C. Jansen, K. Lahabi, and J. Aarts, High-quality ${\mathrm{CrO}}_{2}$ nanowires for dissipation-less spintronics, Phys. Rev. X 6, 041012 (2016).
  * Rabe  _et al._ [2000] M. Rabe, J. Dreßen, D. Dahmen, J. Pommer, H. Stahl, U. Rüdiger, G. Güntherodt, S. Senz, and D. Hesse, Preparation and characterization of thin ferromagnetic CrO2 films for applications in magnetoelectronics, J. Magn. Magn. Mater. 211, 314 (2000).
  * Watts [2002] S. M. Watts,  _Magnetotransport in half-metallic ferromagnetic oxides_ , Ph.D. thesis, Florida State University (2002).
  * Dudarev  _et al._ [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study, Phys. Rev. B 57, 1505 (1998).
  * Korotin  _et al._ [1998] M. A. Korotin, V. I. Anisimov, D. I. Khomskii, and G. A. Sawatzky, ${\mathrm{CrO}}_{2}$: A self-doped double exchange ferromagnet, Phys. Rev. Lett. 80, 4305 (1998).
  * Buschow  _et al._ [1983] K. Buschow, P. van Engen, and R. Jongebreur, Magneto-optical properties of metallic ferromagnetic materials, J. Magn. Magn. Mater. 38, 1 (1983).
  * Weber  _et al._ [2015] J. A. Weber, A. Bauer, P. Böni, H. Ceeh, S. B. Dugdale, D. Ernsting, W. Kreuzpaintner, M. Leitner, C. Pfleiderer, and C. Hugenschmidt, Spin-resolved fermi surface of the localized ferromagnetic heusler compound ${\mathrm{Cu}}_{2}\mathrm{MnAl}$ measured with spin-polarized positron annihilation, Phys. Rev. Lett. 115, 206404 (2015).
  * Delapalme  _et al._ [1971] A. Delapalme, J. Schweizer, G. Couderchon, and R. Perrier de la Bathie, Étude de l’alliage de heusler (Cu2MnAl) comme monochromateur de neutrons polarisés, Nucl. Instrum. Methods 95, 589 (1971).
  * Neubauer  _et al._ [2012] A. Neubauer, F. Jonietz, M. Meven, R. Georgii, G. Brandl, G. Behr, P. Böni, and C. Pfleiderer, Optical floating zone growth of high-quality Cu2MnAl single crystals, Nucl. Instrum. Methods Phys. Res. Sect. A 688, 66 (2012).
  * Tajima  _et al._ [1977] K. Tajima, Y. Ishikawa, P. J. Webster, M. W. Stringfellow, D. Tocchetti, and K. R. A. Zeabeck, Spin waves in a heusler alloy Cu2MnAl, J. Phys. Soc. Japan 43, 483 (1977).
  * Buczek [2009] P. A. Buczek,  _Spin dynamics of complex itinerant magnets_ , Ph.D. thesis, Universitäts- und Landesbibliothek Sachsen-Anhalt (2009).
  * Liechtenstein  _et al._ [1984] A. I. Liechtenstein, M. I. Katsnelson, and V. A. Gubanov, Exchange interactions and spin-wave stiffness in ferromagnetic metals, J. Phys. F 14, L125 (1984).
  * Liechtenstein  _et al._ [1987] A. Liechtenstein, M. Katsnelson, V. Antropov, and V. Gubanov, Local spin density functional approach to the theory of exchange interactions in ferromagnetic metals and alloys, J. Magn. Magn. Mater. 67, 65 (1987).
  * Bruno [2003] P. Bruno, Exchange interaction parameters and adiabatic spin-wave spectra of ferromagnets: A “renormalized magnetic force theorem”, Phys. Rev. Lett. 90, 087205 (2003).