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 (Néel) 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–Schrödinger 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) Schrödinger 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-of-equilibrium 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 time-dependent 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 Schrödinger 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 linear-response 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 . 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 linear-response formalism, and in Appendix B, we discuss the so-called 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 formwhere is the kinetic energy operator and n̂() and σ̂() are the local operators yielding the electronic density, n(, t) = ⟨Φ(t)|n̂()|Φ(t)⟩, and the spin magnetization, (, t) = ⟨Φ(t)|()|Φ(t)⟩, which couple to the scalar potential v(, t) and the magnetic field (, 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(, t) and the spin magnetization (, t), are determined by solving a time-dependent Schrödinger equation for effectively noninteracting electrons, the so-called Kohn–Sham (KS) equationHere, Φs, represents the KS Pauli spinors labeled by a generic index n and f is the corresponding occupation number specifying the initial state. The effective potentials vs and are related to the external potential, v, and magnetic field, , of the interacting Hamiltonian (1) by vs = v + vH[n] + vxc[n, ] and , respectively. The Hartree potential vH is the classical potential due to the electronic density n. The so-called exchange-correlation (xc) potential and magnetic field, vxc and , 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, , 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 isIn eq , we introduced m± = m ± im and , 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). We can rewrite the time-dependent Schrödinger equation for the KS electrons asIn the previous equation, we denote all explicitly time-independent 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 blockwhich corresponds to a time-dependent Zeeman coupling due to a magnetic field characterized by the wave vector and frequency ω, that isThe operators σ̂±() = σ̂() ± iσ̂() are the local spin-raising and -lowering operators, respectively. Taking a superposition of the building blocks (5) over all frequencies, we arrive atwhere 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) = τ/π(t2 + τ2), with a pulse width ∝ τ. In the limit of vanishing pulse width (τ → 0) and infinite field strength , we arrive atThe 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 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 t0 = 0 to the time just after the transverse magnetic kick, t0+ = 0+, is simply given bywhere 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 two-component spinor wave functions in real space. It is instructive to look at the spin magnetization just after the kickwhere 0() = ⟨Ψ(t = 0)||Ψ(t = 0)⟩, the ground-state spin magnetization (oriented here along the z-axis), and () is a spatially varying rotation matrix characterized by the wave vector Equation corresponds to a rotation around a spatially dependent rotation axis = [cos(·),sin(·), 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 .

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π/||. The resulting spins, depicted as filled red arrows, form a spin wave.

After the transverse magnetic kick, the system is no longer in an eigenstate of Ĥs, and a subsequent propagation with a field-free 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±(). 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 ), which can host the magnetic structure with a periodicity given by . 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 and 9b.

Results for Simple Magnetic Transition Metals

We have implemented the transverse magnetic kick, described in Section , 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 , 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 in Appendix B). In principle, the GBT allows for arbitrary wave vectors for the initial kick or, equivalently, the induced magnetic excitation. In practice, however, we found that the propagation turns out to be unstable if is not commensurate with the -point grid. A transverse magnetic kick with the wave vector mixes spin-up and spin-down components of a spinor, which differ by a wave vector of . Only if the wave vectors differing by are representable by the -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 -point grid. A way around this numerical challenge could be to use two -point meshes, one for spin-up components and the other for spin-down components, which are shifted by .

Figure 7

Schematic representation of the boundary conditions used for the (a) Bloch theorem and (b) generalized Bloch theorem. We implemented these “twisted” boundary conditions for which an extra phase factor is added to the boundary points. This extra phase depends on the label (α or β) and the component of the Pauli spinor, as explained in the main text.

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 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, specificallywhere the spin magnetization is computed from the KS states [cf. eq ]. Next, we compute the Fourier transform of m±(, t), that iswhere the time signal is artificially damped, described by a rate , 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 . Note that a finite value of α fixes the value of the broadening, which in our studies was chosen to be α = 5. The m±(, ω) values are then proportional to the transverse spin susceptibilities χ±∓(, ω; = ′ = 0), provided the initial kick is not too strong. Below, we show heat maps of the imaginary part of m+(, ω) as a function of the wave vector 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 local-density 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.

Nickel

We compute Ni using an FCC unit cell with a theoretical lattice constant of a = 3.436 Å,[23] a 16 × 16 × 16 -point mesh with four shifts, and a norm-conserving Hamann–Schlüter–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 -point grid allows us to compute 16 distinct wave vectors along . In Figure , we show a heat map of the imaginary part of (the negative of) eq . 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 ). From the lowest three wave vectors, we can fit a quadratic dispersion relation, ω(q) = Dq2, 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−23,48]

Figure 2

Left panel: density of states of magnetic excitations for nickel. Green stars represent experimental results from inelastic neutron scattering (INS).[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.

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 -point mesh with four shifts, a grid spacing of 0.27 Bohr, and an HSCV pseudopotential. Figure depicts the DOS of transverse magnetic excitations for Co. At small wave vectors, we find a clear magnon peak, which fades out around halfway through the Brillouin zone along . The value of the spin stiffness estimated from the lowest three wave vectors (D = 325 meV Å2) is in good agreement with linear-response studies using the adiabatic LDA[18,23,48] and, in contrast to Ni, close to experimental values obtained from inelastic neutron scattering.[49]

Figure 3

Density of states of magnetic excitations for cobalt. The black line depicts a quadratic fit, ω(q) = Dq2, of the first three magnon peaks. Yellow stars are results estimated from inelastic scanning tunneling spectroscopy (ISTS) measurements.[44]

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 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 -point 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 -point mesh enabled the computation of 16 wave vectors along and eight wave vectors along . No shifts are applied to the -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 and , 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 ). 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−23,48]

Figure 4

Density of states of magnetic excitations for iron. We show the DOS along two high-symmetry lines, that is, along and . Experimental results (isotropic) from inelastic neutron scattering are shown as green stars.[50] For each high-symmetry line, we fitted a quadratic magnon dispersion relation. The anisotropy is slightly reduced compared to the linear-response TD-DFT presented in ref (48).

The main advantage of our approach to compute transverse magnetic excitations is that we are not limited to the linear-response 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 ). 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 , we show the DOS for the fixed wave vectors and , for various opening angles θ. For a wave vector of , 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 , 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.

Figure 5

Left panel: dependence of the DOS on the opening angle for a wave vector for iron. Right panel: density of states at wave vector for various opening angles θ.

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 real-space 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 linear-response 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 Schrödinger equation, and not the 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 → 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 linear-response 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., pump-probe 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.