Time-Dependent Magnons from First Principles

We propose an efficient and non-perturbative scheme to compute magnetic excitations for extended systems employing the framework of time-dependent density functional theory. Within our approach, we drive the system out of equilibrium using an ultrashort magnetic kick perpendicular to the ground-state magnetization of the material. The dynamical properties of the system are obtained by propagating the time-dependent Kohn–Sham equations in real time, and the analysis of the time-dependent magnetization reveals the transverse magnetic excitation spectrum of the magnet. We illustrate the performance of the method by computing the magnetization dynamics, obtained from a real-time propagation, for iron, cobalt, and nickel and compare them to known results obtained using the linear-response formulation of time-dependent density functional theory. Moreover, we point out that our time-dependent approach is not limited to the linear-response regime, and we present the first results for nonlinear magnetic excitations from first principles in iron.


INTRODUCTION
Collective magnetic excitations play a central role in our understanding of the stability of magnetic materials; for example, they determine the Curie (Neél) temperature of (anti-) ferromagnets. 1 In recent years, the experimental confirmation of intrinsic two-dimensional magnets 2−7 has attracted a lot of attention, since it adds magnetism to the toolbox of van der Waals heterostructures, which bears the promise for creating efficient novel devices with tailored electronic, optical, and magnetic properties. 8 To understand the coupled spin and charge dynamics of extended electronic systems, numerical simulations can provide insight into the involved physical processes and complement experimental investigation. Broadly speaking, there are two approaches to study the spin magnetization dynamics: (1) A phenomenological approach, where the spin degrees of freedom are encoded in a model Hamiltonian, such as the Heisenberg model, and the magnetization dynamics is computed, for example, by solving the classical Landau−Lifshitz−Gilbert equation. 9 (2) The first-principles approach, based on the solution of the microscopic Pauli−Schrodinger equation for the electrons, where magnetism emerges due to the intrinsic magnetic moment of the interacting electrons. In our work, we follow the latter route; however, the direct solution of the (time-dependent) Schrodinger equation for interacting electrons is not possible in practice, at least for extended systems. We employ instead the exact reformulation of quantum mechanics based on density functional theory (DFT). 10−12 Since its inception, DFT has evolved into the most widespread approach to numerically study extended systems from first principles. Its generalization to time-dependent or out-ofequilibrium phenomena, the so-called time-dependent DFT (TD-DFT), 13−15 is currently the only viable and reliable route to study the dynamics of large-scale quantum mechanical systems from first principles despite the intrinsic limitations introduced by approximating interaction effects beyond the classical Hartree interaction. Specifically, the linear-response formulation of TD-DFT has been used to study not only optical properties but also magnetic excitations such as magnons, which are collective fluctuations of the spin magnetization. 16−23 Over the last couple of years, first studies investigating magnetization dynamics from first principles in real time have emerged. For example, the ultrafast demagnetization due to intense laser pulses 24 has been successfully investigated using TD-DFT approaches. 25 −29 In the present work, we are interested in modeling transverse magnetic excitations, specifically magnons, which are long-wavelength collective excitations of magnetic materials with a typical energy of tens to hundreds of milli-electronvolts. We propose a timedependent alternative to the widely used linear-response formulation, in spirit closely related to the work of Bertsch et al. for optical excitations, 30 which allows us to go beyond the linear regime and address the response of the system to arbitrary spatially and time-dependent electromagnetic fields. To induce magnetization dynamics in the system, we employ a "transverse magnetic kick". The subsequent time evolution of the spin magnetization, governed by the time-dependent Schrodinger equation, is then analyzed to obtain the density of states (DOS) of transverse magnetic excitations. We show that for small magnetic kicks we recover the results from the linearresponse theory. However, our approach does not rely on the assumption of small perturbations and we can analyze magnetic excitations beyond the linear response. Surprisingly, we find that the frequency of magnons in the ferromagnetic metal iron gets detuned to higher energies, in contrast to the prediction of simple models for the magnons based on the Heisenberg model.
The outline of the paper is as follows: In Section 2, we introduce our prescription to induce transverse magnetic excitations in a real-time TD-DFT framework. Results for the simple transition metals iron, cobalt, and nickel are presented in Section 3, and we provide our conclusions and outlook in Section 4. Furthermore, there are two appendices detailing more technical aspects of our approach. In Appendix A, we provide a more in-depth discussion of the TD-DFT linearresponse formalism, and in Appendix B, we discuss the socalled generalized Bloch theorem, which allows for the efficient simulation of long-wavelength magnons using the chemical unit cell.

REAL-TIME MAGNONS WITHIN TD-DFT
We consider a Hamiltonian of the form where T m 1 2 2 = Π̂is the kinetic energy operator and n(r) and σ(r) are the local operators yielding the electronic density, n(r, t) = ⟨Φ(t)|n(r)|Φ(t)⟩, and the spin magnetization, m(r, t) = ⟨Φ(t)|σ(r)|Φ(t)⟩, which couple to the scalar potential v(r, t) and the magnetic field B(r, t), respectively. The momentum operator Π̂contains the vector potential and spin−orbit coupling, 31−34 but in the absence of an external laser pulse and ignoring the spin−orbit coupling, it is simply given by Π̂= −iℏ∇. Finally, Ŵrepresents the interaction between the electrons, a Coulombic density−density interaction. Within TD-DFT, the time-dependent densities, here the electronic density n(r, t) and the spin magnetization m(r, t), are determined by solving a time-dependent Schrodinger equation for effectively noninteracting electrons, the so-called Kohn− Sham (KS) equation respectively. The Hartree potential v H is the classical potential due to the electronic density n. The so-called exchange-correlation (xc) potential and magnetic field, v xc and xc , are functionals of the density and spin magnetization and need to be approximated in practice. The xc terms take electron−electron interactions beyond the classical Hartree interaction into account. Within DFT, a magnetic material is characterized by an effective field, s x c = , that is, the effective magnetic field is purely internal (due to electron−electron scattering). Put differently, the xc magnetic field is crucial to stabilize the magnetic order.
A common approach to study transverse magnetic excitations is to employ the linear-response formulation of TD-DFT. In a nutshell, one computes the transverse magnetic susceptibility χ +− , which determines the transverse spin magnetization response induced by a weak transverse magnetic field, that is In eq 3, we introduced m ± = m x ± im y and i x y 1 2

= [ ± ]
± , which diagonalizes the transverse magnetic response for systems with a collinear ground-state magnetization. 35 The transverse magnetic excitation spectrum is encoded in the susceptibility χ +− . It is worth noting that the transverse magnetic excitation provides complementary information to the longitudinal magnetic response. 36 The longitudinal magnetic response of finite systems has been studied, for example, using a dipolar kick in refs 37, 38. We provide a more detailed discussion of the linear response within TD-DFT in Appendix A.
Next, we describe an alternative approach for the investigation of transverse magnetic excitation based on the time propagation of the KS equations of TD-DFT (2 In the previous equation, we denote all explicitly timeindependent contributions to the Hamiltonian by Ĥs. To induce transverse magnetic excitations in the system, we add an explicit time-dependent contribution, Δ̂(t), to the Hamiltonian based on the building block which corresponds to a time-dependent Zeeman coupling due to a magnetic field characterized by the wave vector q and frequency ω, that is The operators σ̂±(r) = σx(r) ± iσŷ(r) are the local spin-raising and -lowering operators, respectively. Taking a superposition of the building blocks (5) over all frequencies, we arrive at where we damped high frequencies with a characteristic time τ > 0. Averaging over the frequencies leads to the Poisson kernel representation of the delta function in time, δ τ (t) = τ/ π(t 2 + τ 2 ), with a pulse width ∝ τ. In the limit of vanishing pulse width (τ → 0) and infinite field strength ( ) → ∞ ⊥ , we arrive at The angle θ is defined as the finite limit of / τ ℏ ⊥ when τ →0 and → ∞ ⊥ . One can interpret this somewhat artificial "transverse magnetic kick" as exciting all transverse magnetic modes with a wave vector q and an arbitrary frequency. Working in the limit of an infinitesimally sharp kick also allows us to compute the effect of this transverse magnetic kick on the KS states analytically. The propagator from the initial time t 0 = 0 to the time just after the transverse magnetic kick, t 0 + = 0 + , is simply given by where we have used the integral only over an infinitesimal time interval and therefore only the contribution proportional to the δ-function survives. In eq 9a, we give the definition of the infinitesimal propagator representing the transverse magnetic kick, and in eq 9b, we show its explicit form acting on twocomponent spinor wave functions in real space. It is instructive to look at the spin magnetization just after the kick where m 0 (r) = ⟨Ψ(t = 0)|σ|Ψ(t = 0)⟩, the ground-state spin magnetization (oriented here along the z-axis), and R q (r) is a spatially varying rotation matrix characterized by the wave vector q Equation 11 corresponds to a rotation around a spatially dependent rotation axis u = [cos(q·r),sin(q·r), 0] T by a constant angle θ. A pictorial presentation of the effect of the transverse magnetic kick on the spin magnetization is given in Figure 1.
After the transverse magnetic kick, the system is no longer in an eigenstate of Ĥs, and a subsequent propagation with a fieldfree Hamiltonian, that is, a Hamiltonian with no additional (explicitly) time-dependent external field or potential, leads to nontrivial dynamics of the observables. In this work, we are explicitly interested in the dynamics of the transverse spin magnetization encoded in m ± (r). For a small opening angle θ, we expect to recover the linear-response features contained in The main numerical challenge for such a calculation is that typical magnon frequencies are in the range of tens to hundreds of milli-electronvolts while we are interested in wave vectors within the first Brillouin zone. To resolve frequencies on the order of milli-electronvolts, we need to propagate the system up to a picosecond. Moreover, if we are interested in a wave vector, which is a fraction 1/m of the reciprocal lattice vector for the ground-state magnetic unit cell of the magnetic material, a we have to construct an m × 1 × 1 supercell (m repetitions along the wave vector q), which can host the magnetic structure with a periodicity given by q. However, there is a way to circumvent the construction of supercells, which goes by the name of the generalized Bloch theorem (GBT). The GBT has been introduced first by Sandradskii for the calculation of ground-state spin waves. 39 In Appendix B, we show that the GBT can also be used for the time propagation of the states induced by the magnetic kick given in eqs 9a and 9b.

RESULTS FOR SIMPLE MAGNETIC TRANSITION METALS
We have implemented the transverse magnetic kick, described in Section 2, and the GBT, detailed in Appendix B, in Octopus, 40 which is a real-time, real-space implementation of TD-DFT. The initial kick, defined in eq 9b, is straightforward to implement in a real-space code because it is simply a local rotation of two-component Pauli spinors. The GBT is more involved since it corresponds to modified boundary conditions on the spinors (cf. Figure 7 in Appendix B). In principle, the GBT allows for arbitrary wave vectors q for the initial kick or, Figure 1. Schematic representation of the effect of the initial transverse magnetic kick on the ground-state spin magnetization of a solid. The initial spin magnetization is shown as hollow black arrows (ferromagnetic "atomic" moments). The transverse magnetic kick rotates the spin by an angle θ around a rotation axis in the x−y plane. The rotation axis itself (represented by green, dashed lines with black arrowheads) is rotating in the plane perpendicular to the ground-state spin magnetization with the wavelength given by 2π/|q|. The resulting spins, depicted as filled red arrows, form a spin wave.
Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article equivalently, the induced magnetic excitation. In practice, however, we found that the propagation turns out to be unstable if q is not commensurate with the k-point grid. A transverse magnetic kick with the wave vector q mixes spin-up and spin-down components of a spinor, which differ by a wave vector of q. Only if the wave vectors differing by q are representable by the k-point grid, this mixing of the spin-up and spin-down components is numerically accurate. This means that small wave vectors require a very dense k-point grid. A way around this numerical challenge could be to use two k-point meshes, one for spin-up components and the other for spin-down components, which are shifted by q. Subsequently, we show results for the well-studied reference magnetic transition metals iron (Fe), cobalt (Co), and nickel (Ni), frequently used benchmarks for experimental and theoretical approaches to investigate magnetic excitations. In all of the reported results, we employed the "twisted" boundary conditions to benefit from the GBT, and the calculations are performed in primitive cells of the materials unless stated differently. For each transition metal, we induce transverse magnetic excitations at a fixed wave vector q applying the kick described above with an opening angle of θ = 2 × 10 −2 . This value of the opening angle was determined by decreasing θ until a clear spectral feature, such as a magnon peak at small wave vectors, did not change its position further, indicating that the linear-response regime is reached. We record the time evolution of the spin magnetization, specifically where the spin magnetization is computed from the KS states [cf. eq 2c]. Next, we compute the Fourier transform of m ± (q, t), that is where the time signal is artificially damped, described by a rate T Γ = α , such that it decays to a numerically small value at the final time T of the propagation interval. Hence, we can formally take T → ∞ for the upper integration limit in eq 13. Note that a finite value of α fixes the value of the broadening, which in our studies was chosen to be α = 5. The m ± (q, ω) values are then proportional to the transverse spin susceptibilities χ ±∓ (q, ω; G = G′ = 0), provided the initial kick is not too strong. Below, we show heat maps of the imaginary part of m + (q, ω) as a function of the wave vector q along high-symmetry lines (x-axis) and the energy (y-axis), which correspond to the density of states (DOS) of magnetic excitations. All results are obtained using the adiabatic localdensity approximation (LDA) for the effective potentials, allowing us to compare our results to prior studies. Unless specified otherwise, we propagated the system after the initial magnetic kick for T = 18 000 au ≈ 435 fs. The spectral resolution, Δω, can be estimated by requiring that half a period is resolved within the propagation time, T, which yields Δω = ℏπ/T ≈ 5 meV.
3.1. Nickel. We compute Ni using an FCC unit cell with a theoretical lattice constant of a = 3.436 Å, 23 a 16 × 16 × 16 kpoint mesh with four shifts, and a norm-conserving Hamann− Schluẗer−Chiang−Vanderbilt (HSCV) pseudopotential. 45,46 The real-space grid is sampled along the primitive lattice vectors of the cell using a grid spacing of 0.27 Bohr. This specific k-point grid allows us to compute 16 distinct wave vectors along q x X = Γ . In Figure 2, we show a heat map of the imaginary part of (the negative of) eq 13. For small wave vectors, x < 0.25, we find a pronounced peak in the DOS at low energies, which corresponds to the collective magnon mode. At larger wave vectors, this peak gets broadened due to the vicinity of the Stoner continuum, that is, the continuum of single spin-flip excitations (cf. right panel of Figure 2). From the lowest three wave vectors, we can fit a quadratic dispersion relation, ω(q) = Dq 2 , shown as the black parabola. Toward the Brillouin zone boundary, the strength of the mode decreases and its energy settles around 400 meV. The value for the fitted spin stiffness 47 (D = 890 meV Å 2 ) and the energy toward the Brillouin zone boundary is in good agreement with the available linear-response studies using the adiabatic LDA. 18,19,[21][22][23]48 3.2. Cobalt. For simulating the magnetization dynamics of Co, we employed an FCC unit cell with a theoretical lattice constant of a = 3.429 Å 23 and, similar to the calculations for Ni, a 16 × 16 × 16 k-point mesh with four shifts, a grid spacing of 0.27 Bohr, and an HSCV pseudopotential. Figure 3 depicts the DOS of transverse magnetic excitations for Co. At small wave vectors, we find a clear magnon peak, which fades out  41 The discrepancy between the INS results and the TD-DFT results within the adiabatic LDA is sometimes attributed to the overestimated spin gap within LDA. 42,43 Yellow stars are results obtained from inelastic scanning tunneling spectroscopy (ISTS). 44 Right panel: Density of states for several fixed wave vectors in arbitrary units.

Journal of Chemical Theory and Computation
pubs.acs.org/JCTC Article around halfway through the Brillouin zone along X Γ . The value of the spin stiffness estimated from the lowest three wave vectors (D = 325 meV Å 2 ) is in good agreement with linearresponse studies using the adiabatic LDA 18,23,48 and, in contrast to Ni, close to experimental values obtained from inelastic neutron scattering. 49 3.3. Iron. In contrast to our study on Ni and Co, we did not use the elementary (chemical) unit cell for Fe but instead used a cubic supercell containing two Fe atoms. This allowed us to compute the transverse magnetic response at q H 1 2 = Γ without using the GBT, hence numerically verifying the validity of the GBT. We found no sizable difference between the supercell calculation and the GBT, thus validating our implementation of the GBT. For our simulation, we used the experimental lattice constant of a = 2.867 Å, a 16 × 16 × 16 kpoint grid, a grid spacing of 0.27 Bohr, and a PseudoDojo 51 pseudopotential Fe, which we found to be softer for Fe than the HSCV pseudopotential, used for Ni and Co. The use of the specified k-point mesh enabled the computation of 16 wave vectors along H Γ and eight wave vectors along N Γ . No shifts are applied to the k-point mesh, which yields the same density of wave vectors as in our studies of Ni and Co. We find that only for small wave vectors around Γ, that is, roughly up to H 1 4 Γ and N 1 3 Γ , respectively, we get a clear peak corresponding to the collective magnon mode. For larger wave vectors, the magnon peak gets broadened by the Stoner continuum (see Figure 4). Again, as for Ni and Co, our results for Fe are in good agreement with the available linear-response calculations using the adiabatic LDA within TD-DFT. 18,19,[21][22][23]48 The main advantage of our approach to compute transverse magnetic excitations is that we are not limited to the linearresponse regime. A simple way to go beyond the linear regime is to increase the opening angle of the magnetization spiral formed due to the initial kick (cf. Figure 1). Since we found that the positions of relevant spectral features are already converged for a propagation time of T = 9000 au ≈ 220 fs in the linear regime, all following results are obtained with this shorter propagation time. In Figure 5 = Γ , where we have a clear magnon peak in the linear regime, we find that the magnon becomes stiffer for wider opening angles (stronger perturbations). This is quite interesting because a simple, classical Heisenberg model with nearest-neighbor (ferromagnetic) couplings between the spins predicts a softening of the magnon frequency for wider opening angles. Even more surprising is the behavior of the DOS at q H 1 2 = Γ , where we do not have a clear magnon peak in the linear regime due to the damping by the single spin-flip excitations (Stoner continuum). We find that the spectral signal sharpens as we increase the opening angle such that around θ = 0.1, we find a clear single peak, similar to the magnon peaks at smaller wave vectors in the linear regime. This suggests that the Stoner continuum can no longer efficiently damp the collective magnon modes. Increasing the opening angle further, we observe a stiffening of the peak similar to the behavior of the magnon peak at smaller wave vectors.

CONCLUSIONS AND OUTLOOK
In our work, we proposed a novel approach to study magnetic excitations from first principles. The basic idea is to induce magnetization dynamics by an ultrashort kick with a transverse magnetic field. We implemented two approaches in the realspace real-time TD-DFT code Octopus: One is based on the use of supercells that are big enough to contain the wave vector of interest. The other uses the flexibility of the real-space grid to implement the "twisted" boundary conditions of the GBT, which allows for the computation of magnetic excitation of any fixed wave vector using the primitive cell. For the second approach, we found the sampling of the Brillouin zone to be important.
As a consistency check, we performed calculations in the linear-response regime, where we can compare against known results and found an excellent agreement with previously published results. In the linear-response regime, our approach has several features that render it different from the usual way of computing the linear response. b First of all, it scales linearly with the number of occupied states and does not require unoccupied states, in contrast to the construction of the linearresponse function, which usually requires a summation over unoccupied (virtual) states. This becomes increasingly important for simulating large-scale systems and is therefore extremely relevant in case the GBT cannot be applied, for example, when spin−orbit interactions are included. Second, since we are solving the Schrodinger equation, and not the  Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article Dyson equation, for the linear-response function, we only need an approximation for the exchange-correlation potentials and not for the exchange-correlation kernel. Formally, the exchange-correlation kernel is given by the functional derivative of the exchange-correlation potentials with respect to the densities, but in practice, it might be rather cumbersome to implement the required kernel. For example, this plays a role when using orbital functionals, for example, meta-or hyper-GGAs, 53 or mixed approaches like hybrid functionals and LDA +U 54 because for these approaches, which may be referred to as generalized KS schemes, 55,56 it is not obvious how to construct the proper kernel for the linear-response Dyson equation. Finally, even if the exchange-correlation potentials and the corresponding exchange-correlation kernels are available, the ground-state KS equation, which yields the KS response function, and the Dyson equation from which the interacting response function emerges are usually represented in different numerical basis sets. This can lead to a (numerical) violation of the Goldstone theorem, which states that the magnon frequency vanishes for q → 0 in the absence of spin− orbit coupling, while it is by construction satisfied in our approach. Most importantly, our approach is not limited to the linearresponse regime. Already, our proof-of-principles studies for iron revealed interesting trends as transverse magnetic modes are excited more violently. In the region of the energy− momentum space, where a clear magnon peak is present already in the linear regime, we found a stiffening of the magnon frequency. While general nonlinearities can either soften or stiffen the frequency of a harmonic oscillator, the classical Heisenberg model with the nearest-neighbor (ferromagnetic) coupling, for example, predicts a softening of the magnon frequency when the magnon is excited more strongly in contrast to our results. Similarly, magnon−magnon interactions tend to soften the magnon frequency. However, our study on iron, a prototypical ferromagnetic metal, shows the opposite behavior. Even more striking are our findings in the region of the energy−momentum space where collective magnetic excitations are strongly damped by single-electron spin-flip excitations in the linear regime: A stronger excitation seems to "un-burrow" the collective mode, since a clear peak emerges in the density of states. This will be further analyzed in future works.
We stress that our approach is not limited to (anti-) ferromagnetic materials. For example, the machinery can readily be applied to systems exhibiting a nontrivial magnetic order in the ground state, such as spin spirals or skyrmions. However, starting from an already noncollinear ground state implies that transverse and longitudinal magnetic excitations are mixed and so the interpretation of the results might be more difficult.
Finally, using our approach, it is straightforward to include spin−orbit coupling or to combine the magnetic kick with other external drivers, such as periodic lasers (Floquet engineering 57 ) or ultrashort intense laser pulses (e.g., pumpprobe spectroscopy). This opens up exciting possibilities to study the interplay of electronic and magnetic excitations in and out of equilibrium within the first-principles framework.

MAGNONS
In this appendix, we provide a concise summary of the linearresponse formulation of TD-DFT. We emphasize that all new results presented in this work are obtained from the time propagation presented in the main body of this paper and not within the linear-response approach discussed in this appendix. The motivation for this appendix is twofold: First of all, we would like to give a short recap of the linear-response formulation of magnetic excitations to provide the minimal background for readers not deeply familiar with magnetic excitations, for example, the definition of magnons within the linear response formalism. Second, we provide the linearresponse formulas using our "notational" conventions to allow for a direct comparison to other approaches formulated in the linear response regime.
In the linear regime, the response of the density and the spin magnetization to a perturbation by a scalar potential and magnetic field is given by i k j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z (14) In eq 14, we have employed a compressed notation where we omit the dependence on (r, t) of the induced density and spin magnetization and the dependence on (r′, t′) of the perturbing scalar potential and magnetic field. Note that we have included the coupling constant μ B in the definition of the perturbing magnetic field, B B μ = , which means that the "magnetic field" has the same units as the scalar potential v, that is, both have the units of energy. The response matrix χ ̲ depends on, both, (r, t) and (r′,t′), and the * denotes a contraction (integral) over (r′,t′). Within TD-DFT, the response function of the interacting system χ ̲ is computed from the response function of the effective KS electrons χ ̲ s and the so-called (Hartree-)exchange-correlation kernel f̲ Hxc by solving the following equation The KS response is determined from the eigenstates and eigenvalues of the static KS equation. A common reformulation of the linear response eq 14 is to introduce the following linear combination of the densities and potentials In terms of these densities and potentials, the KS Hamiltonian takes the form i k j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z The Hartree-exchange-correlation kernel can couple the different channels. In fact, already the Hartree contribution leads to a coupling between the ↑↑ and ↓↓ channels, but it vanishes for the +− and −+ channels, since it stems from a density−density interaction. In this work, we focus on the adiabatic local-density approximation (ALDA) for the exchange-correlation contribution to the effective potentials.
In the linear-response formulation of TD-DFT, this leads to an exchange-correlation kernel f̲ xc , which is local in space and time, that is  ) where f̲ xc unif is the static (ω → 0) long-wavelength (q → 0) limit of the exchange-correlation kernel for a uniform electron gas with the density n and spin magnetization m evaluated at the local ground-state density and spin magnetization, n 0 (r) and m 0 (r), respectively, of the system. Like the Hartree contribution, it is block-diagonal in the "up"−"down" sector, but, in contrast to the Hartree contribution, it has a nonvanishing diagonal in the +− and −+ channels. This contribution is crucial for the description of magnons, which are collective excitations of the transverse spin magnetization: without f xc , there would be no magnons, much like there would be no plasmons without the Hartree contribution. In Figure 6, we show a sketch of the imaginary part of χ +− , which represents the density of states of magnetic excitations perpendicular to the ground-state magnetization.
■ APPENDIX B: TIME-DEPENDENT GENERALIZED

BLOCH THEOREM
The generalized Bloch theorem (GBT) has been introduced already some time ago by Sandratskii 39 to study magnetic systems exhibiting a spin spiral configuration in the ground state. Moreover, this approach can also be used to compute a constrained ground state, that is, the minimal energy state compatible with a magnetization forming a spin spiral, which can be used to compute adiabatic or frozen magnons. To the best of our knowledge, the GBT has not been applied within a time-dependent framework. Hence, within this appendix, we answer the question regarding the conditions under which the application of the GBT is possible within a time-dependent approach.
Let us start by recalling why the Bloch theorem is applicable for the propagation of solids: If the KS Hamiltonian for all times commutes with translations by crystal lattice vectors, the propagator commutes itself with these translations. Hence, if the propagation is started in an eigenstate of the translation operator, for example, a Bloch state, it remains an eigenstate of the translation operator. Now, let us look at how the states after the transverse magnetic kick, defined in eq 9a and 9b, behave under translations by a lattice vector. We assume that the system is collinear in the ground state, which means that we ignore the spin−orbit coupling or a noncollinear spin magnetization due to frustration (for example, in a twodimensional hexagonal lattice with anti-ferromagnetic nearest- Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article neighbor coupling). This implies that the states before the magnetic kick are Bloch states, which are either pure "spin-up" or pure "spin-down". After the initial kick, the states are given as follows The "partial" spin magnetizations m α/β 0 are determined exclusively by and therefore are lattice periodic