{'fraction_non_alphanumeric': 0.03956695450444645, 'fraction_numerical': 0.010826137388838768, 'mean_word_length': 4.622705314009662, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The generalized gradient approximation kernel in time-dependent density functional theory

A complete understanding of a material requires both knowledge of the excited states as well as of the ground state. In particular, the low energy excitations are of utmost importance while studying the electronic, magnetic, dynamical, and thermodynamical properties of the material. Time-Dependent Density Functional Theory (TDDFT), within the linear regime, is a successful \textit{ab-initio} method to access the electronic charge and spin excitations. However, it requires an approximation to the exchange-correlation (XC) kernel which encapsulates the effect of electron-electron interactions in the many-body system. In this work we derive and implement the spin-polarized XC kernel for semi-local approximations such as the adiabatic Generalized Gradient Approximation (AGGA). This kernel has a quadratic dependence on the wavevector, {\bf q}, of the perturbation, however the impact of this on the electron energy loss spectra (EELS) is small. Although the GGA functional is good in predicting structural properties, it generality overestimates the exchange spin-splitting. This leads to higher magnon energies, as compared to both ALDA and experiment. In addition, interaction with the Stoner spin-flip continuum is enhanced by AGGA, which strongly suppresses the intensity of spin-waves.


I. INTRODUCTION
Recent developments in the field of laser-induced spindynamics have opened up the world of femtomagnetism 1 , whereby the spin degree-of-freedom is controlled using ultrafast laser pulses 2 . As the name suggests, femtomagnetism concerns charge and spin dynamics on the femtosecond (= 10 −15 s) time scale, corresponding to energies in the meV range. Electronic excitations in this energy regime can be clasified as either localized single-particle like transitions, e.g. Stoner spin-flips, or those exhibiting a collective nature, such as excitons or magnons 3 , which are spread over many atomic sites. These collective excitations with large wavelength occur at relatively lower energies as compared to the singleparticle excitations and hence are important in studying the material's properties. To exploit the vast potential femtomagnetism offers, it is vital that we are able to accurately describe these collective excitations, in order to understand, and ultimately control, them.
Theoretical studies of charge and spin excitations can be performed using either simple models like the Heisenberg model 4 or Landau-Lifshitz-Gilbert equation 5 , or computationally more demanding, parameter-free, abinitio methods. In contrast to ab-initio methods, model based approaches are limited by their lack of generality, as they are usually tailored to study only specific problems and cannot be used universally for all materials.
Time-dependent density functional theory (TDDFT) [6][7][8] is an ab-initio method which can predict the excited state properties of a material. Since its theoretical foundation in 1984 6 , it has been successfully applied to study excited state properties of a wide range of materials 9 . Compared to other ab-initio methods, such as many-body perturbation theory, TDDFT provides a similar level of accuracy, but at far less computational cost.
The real time evolution of electronic charge and spin densities is calculated using TDDFT by solving the single-particle Kohn-Sham (KS) equations. The effects of electron-electron interactions come into this noninteracting KS system via an effective potential, the so called Hartree exchange-correlation (XC) potential. Although TDDFT is an exact theory for treating systems under the influence of strong time-dependent external potentials [10][11][12][13][14][15][16][17] , it is most commonly applied within the weak perturbation limit. When working in this linear regime, one requires the functional derivative of the XC potential, the so called XC kernel. In a practical TDDFT calculation, an approximation to the XC potential and kernel is required.
There are many different flavors of XC energy functionals in ground-state DFT, which can be divided into the local density approximation (LDA), generalized gradient approximations (GGAs), meta-GGAs, hybrids and Fock-like approximations, comprising the so called Jacob's ladder of approximations, where the level of accuracy increases as we climb from LDA to hybrids. The performance of these approximations in static ground-state DFT has been well-studied, however much less is known about their behavior in TDDFT (when combined with the adiabatic approximation). Most of the research so far concerns the optical absorption spectra, and, in particular, the failure of simple XC kernels to predict bound excitons. From these studies, we know the importance of describing the long wavelength limit of the XC kernel correctly in order to obtain reasonable exciton binding energies, leading to a number of new approximations [18][19][20][21][22][23][24][25] . However for magnetic excitations, only the ALDA XC kernel has been properly studied, e.g. for calculations of the magnon spectra [26][27][28][29][30] , where for many cases it overestimates magnon energies as compared to experiments.
In the static DFT case, it is well known that going from LDA to GGA improves many ground state properties 31 . Hence, in the work presented here, we climb up to the next rung of Jacobs's Ladder within TDDFT and ask if including gradient corrections to the XC kernel improves the charge and spin excitation spectra. The paper is organized as follows: Section II gives the basic equations of TDDFT and how they may be used to calculate excitation energies via the linear-response susceptibilities. We will also derive the adiabatic GGA (AGGA) XC kernel for non-collinear spin systems in this section. In Section III, we apply the AGGA kernel to first study the Electron Energy Loss Spectra (EELS) for medium-(diamond) and large-(LiF) bandgap insulators. A comparison is made with experiments and the ALDA kernel. Then a more comprehensive study is made for magnetic excitations of simple bulk ferromagnetic systems Fe, Co and Ni, and Heusler materials. Again we compare with experimental results, as well as previous theoretical works, all of which used the ALDA XC kernel. Finally, in Section IV, we give some concluding remarks on the performance of AGGA.

II. THEORETICAL FORMULATION
When the external perturbation is small, the response of the system to this stimulus is studied through a suitable response function. In general, the response of a system to external stimulus can be expanded in a Taylor series with respect to the perturbation. The coefficients of this expansion are the response functions which have useful information embedded in them, such as the optical absorption spectra, Pockels effect, optical rectification, second harmonic generation, Kerr effect etc. In this paper we primarily focus on the first order response functions and particularly on the charge-charge response, (δρ/δv ext ), and the spin-spin response, (δm/δB ext ), where ρ, v ext , m, B ext correspond to the charge density, electric scalar potential, magnetization density and magnetic field, respectively. For noncollinear systems, the full response is by: , and has a 4 × 4 structure as shown in Fig. (1).
The non-interacting KS linear response functions can easily be calculated in terms of the KS spinors, φ(r), using first-order perturbation theory: where f p , f q denote the occupation number of the p th , q th band, respectively, and σ µ are the fourdimensional counterparts of Pauli spin matrices.
TDDFT relates this non-interacting response function of the KS system to that of the interacting system via a Dyson-like equation: where f µν H (r, r ) = δ µ0 δ ν0 v(r, r ) is the Hartree kernel and v(r, r ) = 1/|r − r | is the Coulomb potential, ω corresponds to frequency and f µν XC (r, r , ω) is the XC kernel, which is the Fourier transform of: is the combined XC potential for the scalar, v XC , and magnetic, B XC , fields.
Since it is convenient to work in reciprocal space for periodic systems, all quantities are represented as matrices via Fourier transformation. The transformed interacting response has the form: Here G, G are the reciprocal lattice vectors and q is the wave vector of the perturbation.
Conventionally, the excitations are studied in the decoupled limit where the off-diagonal terms of Fig. (1) are set to zero (i.e. δm/δv ext = δρ/δB ext =0). This allows us to separate the dielectric response and magnetic response.
Experimental observables may then be extracted from the response functions, for example, the inverse dielectric function is The imaginary part of −1 gives the EELS whereas the imaginary part of gives the absorption spectrum. Likewise, the neutron scattering cross-section is proportional to the transverse magnetic response 32 where Here, the transverse terms comprise χ −+ (= 2χ xx + 2iχ xy ) and χ +− (= 2χ xx − 2iχ xy ). The term χ zz does not contribute to the spin-flip excitations, rather it's the transverse terms of the magnetic susceptibility which give rise to the Stoner and magnon excitations.

A. GGA Kernel
Interactions between the electrons are encapsulated in the XC potential. Knowledge of the exact form of v XC , and hence f XC , would lead to the exact solution of all many-body problems. However, the exact form of v XC is not known and approximations are required for all practical calculations. Here we will derive the XC kernel for GGA functionals within the adiabatic approximation (AA), which is semi-local in space and local in time. For the spin unpolarized case, the XC energy functional, E XC , depends not only on the density, n(r), but also on its gradient, ∇n(r), at each point r in space. The XC potential and kernel can be obtained from first and second order functional derivatives, respectively, of E XC with respect to density i.e., where E XC [ρ] = e XC (ρ, ∇ρ)( r)d 3 r for the GGA functional and e XC is the XC energy density.
The variation of the XC energy is defined by: Taylor expanding the energy density up to first order gives leading to the expansion of energy functional, and Carrying out integration by parts of the second term gives us: Comparing Eq. (14) with Eq. (10), we find the XC potential as: where σ = ∇ρ · ∇ρ is often used in practice. Variation of this potential to first order will give the kernel: Integrating these terms individually by introducing a delta function gives us the kernel for GGA functional [33][34][35] Repeating the above derivation for the spin polarized case (see Appendix A) gives us two equations comprising the symmetric terms f αα XC (r, r ), f ββ XC (r, r ) and the asymmetric terms f αβ XC (r, r ), f βα XC (r, r ) of the XC kernel matrix, where α and β label the up and down spins, respectively, which is extended to non-collinear systems via the Kübler method 36 , see Eqs. (28)-(30).

III. COMPUTATIONAL DETAILS
All calculations are performed using the all-electron full-potential linearized augmented plane wave electronic structure code ELK 43 with PW91 (LDA) 44 and PBE (GGA) 45 functionals. For diamond and LiF, a fcc crystal structure with lattice spacing of 3.5599Å and 4.0259 A, respectively, is used. A dense k-point grid is required to obtain the response functions, hence the BZ is sampled on a k-point grid of 25 × 25 × 25 for both. The method to obtain response functions is a two-step procedure, firstly a ground-state calculation is done to obtain the converged density and potentials. The scissor operator has been used to correct the optical band-gap by 1.306 eV and 5.06 eV for diamond and LiF, respectively. Then the EELS spectra of LiF and diamond are obtained from the LDA and GGA kernels using the corrected band gaps.
The magnon spectra are highly sensitive to a number of parameters, hence convergence has to be checked with respect to k-point grid and the number of G-vectors. In these calculations we have used a 40×40×40 k-point grid. The response functions are expanded in G-space with the length of G-vector up to 6 Bohr −1 . The lattice constant used for Co 2 MnSi is 5.640Å and for NiMnSb is 5.897 A. The spectra are corrected by choosing convergence parameters so that it best satisfies Goldstone's theorem.

A. Semiconductor Spectra
It is well known that the q-dependent behavior of the XC kernel is of vital importance for predicting the optical response of materials. For example, the XC kernel must go as 1/q 2 in order to capture excitonic effects 7,22,25,46-48 in the long-wavelength limit (q → 0). However, the first  rung on Jacob's ladder, the ALDA, does not display any q-dependence, owing to the local approximation for the XC energy. This explains why ALDA does not yield excitonic peaks 7 . The second rung consists of semi-local functionals which include information not just about the density but also its gradients. In this case, it has been shown that the AGGA kernel shows q 2 behaviour 49 . Hence we explore if this has an impact on the spectra.
In Fig. 2(a) we plot the EELS for (i) LiF, which is a large bandgap material with a bound exciton, and (ii) diamond which is a medium bandgap material with excitonic effects appearing as a shift in the spectral weight towards lower energies. For LiF, we can see that AGGA shifts the peak energies for q = 0.24ΓX, 0.48ΓX towards the correct experimental values and also makes the peaks more pronounced, as compared to ALDA. Both ALDA and AGGA fail to capture the excitonic peak at 13 eV as neither has the correct 1/q 2 behavior in the longwavelength limit. Outside the first BZ (q = 1.52ΓX), ALDA and AGGA exhibit similar behaviour. For diamond, neither AGGA nor ALDA captures the shift in spectral weight as can be seen in Fig. 2(b). In fact there is little difference between the results obtained using the two approximations.
To conclude this section, the q-dependence of the AGGA kernel slightly improves upon the ALDA results, although neither captures the excitonic effects.

B. Spin Excitation Spectra
For a computational method to be fully ab-initio, given the atomic composition, it should first predict the equilibrium geometry. Having found the minimum energy crystal structure from ground-state calculations, we can then calculate the excited state properties, all without reference to experimental data. Only those methods which follow this prescription can be considered fully predictive.
We begin by reviewing ground-state DFT calculations with experimental and optimized lattice parameters. The results are summarized in Table I. From this we conclude that (i) GGA is very good in reproducing the structures of materials whereas (ii) LDA is better in predicting the magnetic moments. Within TDDFT the magnon spectra of a system can be calculated from the transverse response function χ −+ (q, ω), which is found using the Dyson-like equation, Eq. (3). The excitations in χ −+ (q, ω) originate from two sources, (1) renormalized poles of the KS response χ 0 corresponding to the Stoner continuum of single-particle spin-flips and (2) additional peaks created by the XC kernel corresponding to spin-wave excitations.
To find the magnon dispersion, we calculate Im{χ −+ (q, ω)} for each q value, extract the magnon peak position and then plot these as a function of q. This is shown in Fig. (3) for nickel, cobalt, and iron along with experimental measurements and ALDA calculations reported in the literature [26][27][28][29]39,53 . We first note that the ALDA results reported in this paper are consistent with previously reported results. In the following, we will discuss each material individually before commenting on the general behavior of the AGGA XC kernel.
For Ni, Fig. 3(a), both ALDA and AGGA show quantitatively the same behaviour from the BZ center to |q| = 0.4. As we move further away from the zone center, the AGGA kernel tends to deviate from ALDA until it becomes ≈ 80 meV higher in energy at the zone boundary. For Co, Fig. 3(b), the experimental results are well captured by both ALDA and AGGA calculations. For Fe, Fig. 3(c), in contrast to ALDA which reproduces the experimental values, the AGGA dispersion overestimates the magnon energies. Beyond the zone center, the transverse response function obtained using AGGA becomes too broad to assign a single energy to the excitation peaks.
The strength and width of the peak in Im{χ −+ (q, ω)} is related to the scattering amplitude and lifetime of the magnon, respectively. To visualize how these properties change throughout the BZ, we can make a 2D contour plot of Im{χ −+ (q, ω)}. These are shown for both ALDA and AGGA in Fig. (4) for Ni, Co and Fe.
Beginning again with Ni, we see that the peaks in  Table I.
Im{χ −+ (q, ω)} obtained by using ALDA (Fig. 4(a)) are stronger in intensity and better resolved than AGGA ( Fig. 4(b)). There exists a high probability of observing a magnon at the BZ boundary with ALDA whereas it is suppressed significantly beyond q = 0.5ΓX with AGGA. We also observe a strong suppression in the magnon intensity between |q| = 0.1 and |q| = 0. shows two branches, both correctly predict an abrupt change in the magnon dispersion around this value of q. This is due to the Stoner spin-flip transitions having comparable energy to the magnon energy causing strong interference and intensity suppression at these values of q (as can be seen later in Fig. 7). In contrast to Ni, the experimental dispersion for fcc Co obtained by Balashov et al. 41 along [100] does not show any optical branches. Both ALDA and AGGA behave the same and show good agreement with experiment, with AGGA being slightly lower in energy. Observing the full transverse response function over the whole BZ (Figs. 4(c) and 4(d)) we again see a reduction in the peak strength and suppression of the magnons by AGGA as compared to ALDA. Qualitatively, AGGA also reproduces the jump in magnon energy witnessed in experiments around |q| = 0.6, although at a higher value of |q|(= 0.8).
For Fe, we see significant broadening in the AGGA (Fig. 4(f)) transverse response for |q| > 0.5, to such an extent that it becomes impossible to assign a peak position. This may explain the experimental observation that magnons were not observed for |q| > 0.46. At this value of q, we observed, with AGGA, strong suppression of the magnon due to interaction with the Stoner continuum.
We now test the predictive power of LDA and GGA by comparing their behavior for Ni, Co, and Fe at the experimental and optimized lattice parameters (see Fig. 5). We find that the AGGA magnon spectra are more sensitive to the lattice parameters than the ALDA. In most cases, the AGGA results at the corresponding GGA parameter (a 0 ) theo are lower in energy and closer to the experimental results than at (a 0 ) exp , although still overestimated. Suprisingly, for the case of Co, the GGA result is better than the LDA result. .
Next we investigate Heusler (or half-Heusler) materials Co 2 MnSi and NiMnSb, which due to their geometry of interlocking magnetic fcc lattices, can (in principle) have multiple magnon branches. In Fig. 6 the magnon spectra of Co 2 MnSi and NiMnSb are plotted along with the experimental and ALDA calculations. For Co 2 MnSi (see Fig. 6(a)) both, an acoustic branch and an optical branch are observed. An increase in the energies of acoustic branch is noted when compared with earlier ALDA results 54 , based on loss function. However, energies of the optical branch with AGGA are within the same range as reported with ALDA. One interesting aspect is that the upper branch, reported with ALDA in literature, was not seen with ALDA in this work, rather it is seen with AGGA. For NiMnSb (Fig. 6(b)), both ALDA and AGGA predict only an acoustic branch, as is also the case experimentally. This is likely due to Ni not possessing a strong local moment as most of the total moment is localized on the Mn atoms. In this case AGGA severely overestimates the magnon energies.
Finally, we offer the underlying explanation as to why AGGA tends to overestimate the magnon energies. The role of the XC kernel is to transform the excitation structure of Im{χ −+ 0 } into the true response. From the Stoner single-particle excitations, contained in Im{χ −+ 0 }, it must create the magnon peak. At q = 0, the gap in Im{χ −+ 0 } is related to the exchange splitting between spin-up and spin-down states. This splitting dictates the position of the Stoner continuum across the BZ. In Figs. 7(a) and (b), we plot Im{χ −+ 0 } for LDA and GGA, where we observe that the Stoner gap has increased by approximately 60 meV for Nickel. This increment stems from the fact that GGA increases the exchange splitting in Ni by 59.9 meV compared to LDA, leading to the shift in Stoner continuum towards high energies. Similar behavior was also observed for other materials, e.g for Fe there is ≈ 150 meV increase, and even a 50% increase of the LDA Stoner gap in χ −+ 0 for the half-metal NiMnSb. This increase in the Stoner excitation energies translates into an increase in the magnon energies. Given that LDA is well-known to overestimate the exchange splitting, the further enhancement on going from LDA to GGA leads to large overestimation of the magnon energies.

V. CONCLUDING REMARKS
We have studied the charge and spin excitation spectra using the gradient dependent AGGA XC kernel within the linear response regime of TDDFT. The calculated EELS for LiF and diamond show that the AGGA kernel performs slightly better than the ALDA kernel, although, as would be expected, neither captures excitonic effects. For magnon dispersions, AGGA generally overestimates the magnon energies. This is due to the fact that the GGA XC functional overestimates the exchange splitting. Furthermore, the intensity of the peaks is greatly suppressed in the spectra obtained by the AGGA XC kernel due to interaction of spin-waves with the Stoner continuum. This suppression is also observed in experiments, suggesting AGGA might provide better qualitative understanding than ALDA. Heusler materials consisting of multiple magnetic sublattices were also studied where it was found that AGGA is better at resolving higher-energy optical magnon branches. When experimental reference data for the system geometry and lattice parameters is unavailable, we found that using GGA consistently, i.e. for both the ground-state (i.e. lattice parameter as well as density) as well as to calculate the response function, gave better results than using the experimental parameters.
In this work, we principally investigated the spinspin response of collinear ferromagnetic systems, which greatly simplified the XC kernel. However, the AGGA XC kernel derived here is valid for all systems, and, in particular, has interesting terms for the the spin-charge response and for non-collinear systems. This will be explored in future work.
It is important to test the performance of adiabatic functionals in TDDFT as their behavior can be quite different from the ground-state case. Only by implementing, assessing, and understanding this behavior can we gain insight into the relevant features necessary for accu-rate XC kernels, which can guide us in developing new approximations in TDDFT.