Ab‐Initio Real‐Time Magnon Dynamics in Ferromagnetic and Ferrimagnetic Systems

Magnonics—an emerging field of physics—is based on the collective excitations of ordered spins called spin waves. These low‐energy excitations carry pure spin currents, paving the way for future technological devices working at low energies and on ultrafast timescales. The traditional ab‐initio approach to predict these spin‐wave energies is based on linear‐response time‐dependent density functional theory (LR‐TDDFT) in the momentum and frequency regime. Herein, the simulation of magnon dynamics using real‐time time‐dependent density functional theory is demonstrated, thus extending the domain of ab‐initio magnonic studies. Unlike LR‐TDDFT, this enables us to observe atom‐resolved dynamics of individual magnon modes and, using a supercell approach, the dynamics of several magnon modes can be observed simultaneously. The energies of these magnon modes are concurrent with those found using LR‐TDDFT. Next, the complex dynamics of the superposition of magnon modes is studied, before finally studying the element‐resolved modes in multisublattice magnetic systems.

DOI: 10.1002/pssb.201900654 Magnonics-an emerging field of physics-is based on the collective excitations of ordered spins called spin waves. These low-energy excitations carry pure spin currents, paving the way for future technological devices working at low energies and on ultrafast timescales. The traditional ab-initio approach to predict these spin-wave energies is based on linear-response time-dependent density functional theory (LR-TDDFT) in the momentum and frequency regime. Herein, the simulation of magnon dynamics using real-time time-dependent density functional theory is demonstrated, thus extending the domain of ab-initio magnonic studies. Unlike LR-TDDFT, this enables us to observe atom-resolved dynamics of individual magnon modes and, using a supercell approach, the dynamics of several magnon modes can be observed simultaneously. The energies of these magnon modes are concurrent with those found using LR-TDDFT. Next, the complex dynamics of the superposition of magnon modes is studied, before finally studying the element-resolved modes in multisublattice magnetic systems.
via ultrafast laser pulses. [24][25][26] Thus, we extend RT-TDDFT to calculate magnon dynamics in real time, resulting in a method that is both predictive and that captures linear as well as nonlinear dynamics.
This article is arranged as follows: we begin by studying a simple ferromagnet [face centered cubic (FCC) Fe] and discuss the individual modes that may be excited in a 4-atom supercell. This demonstrates the insight into magnon modes that can be gained from their real-space real-time visualizations. With this tool in hand, we then study the atom-resolved dynamics in a multisublattice magnetic material (Co 50 Ni 50 ) and explore the coupled and decoupled modes. Finally, we go beyond ferromagnetic coupling and explore the magnon modes of a ferrimagnetic system (Mn 3 Ga), where the magnetic moment on one Mn atom is antiparallel to the other two Mn atoms in each unit cell.

Magnons in RT-TDDFT
Time-dependent density functional theory (TDDFT) [27][28][29][30] is an ab initio method for solving the dynamics of many-electron systems. It does so via an exact mapping from the interacting system to a noninteracting system, known as the Kohn-Sham (KS) system, defined such that it yields the same density dynamics. For noncollinear spin systems, TDDFT may be extended such that both the exact charge density and magnetization density vector field dynamics are reproduced by the KS system. In this case, the TDKS equation reads where ϕ j ðr, tÞ are two-component Pauli spinor TDKS orbitals, A ext ðtÞ is the external laser field, written as a purely timedependent vector potential, σ are Pauli matrices, v S ðr, tÞ ¼ v ext ðrÞ þ v H ðr, tÞ þ v XC ðr, tÞ is the KS effective scalar potential, and B S ðr, tÞ ¼ B ext ðr, tÞ þ B XC ðr, tÞ is the KS effective magnetic field. The external scalar potential, v ext ðrÞ, includes the electronnuclei interaction, whereas B ext ðr, tÞ is an external magnetic field which interacts with the electronic spins via the Zeeman interaction. The Hartree potential, v H ðr, tÞ is the classical electrostatic interaction. Finally, we have the XC potentials, the scalar v XC ðr, tÞ, and the XC magnetic field, B XC ðr, tÞ, which require approximation. In the adiabatic approximation, these may be calculated using a DFT XC energy functional.
In LR-TDDFT, perturbation theory is used to derive the analytic response functions of the densities to small perturbations of the external potentials. A Dyson-like equation relates the response of the KS system to that of the interacting system. Information regarding the frequency and lifetimes of excitations can be extracted from these response functions. [9,11,15,31] In particular, the magnon quasiparticle excitations can be found from the interacting spin-spin response function, which relates the change in transverse magnetization to applied magnetic fields.
To simulate magnon modes using RT-TDDFT simulations, a two-step procedure is followed. First, a supercell is constructed such that it contains the magnon modes we wish to study. This can be understood from the behavior of the transverse magnetic moments of a simplified magnon mode m x,y ðR α , tÞ $ e iðq⋅R α ÀωtÞ (2) where q is the wavevector of the magnon mode, R α are atomic positions, and ω is the magnon frequency. In a supercell, only those modes where q is commensurate with the lattice vectors R i are included, i.e., only when q ⋅ R i ¼ 2πn. For example, to study the magnon mode with q ¼ X where X ¼ ð0, π=a, 0Þ is the X point of the BZ of a cubic lattice, we construct a supercell with lattice vectors R 1 ¼ að1, 0, 0Þ, R 2 ¼ að0, 2, 0Þ, and R 3 ¼ að0, 0, 1Þ. This can then be generalized for arbitrary unit cells. In this article, we will work with FCC or similar lattices, where X ¼ ð0, 0, 2π=aÞ. As we will see, one advantage the supercell approach has over other possible methods, such as utilizing the time-dependent generalized Bloch theorem, is that several modes can be studied together. When using the generalized Bloch theorem, a q vector must be specified as input to the calculation which constrains the magnetization density to spatially oscillate with phase ϕ p ¼ q ⋅ r, thus preventing other magnon modes to exist. The supercell approach removes this restriction on the spatial behavior of the magnon modes. However, the disadvantage is that low q modes require large supercells which make them computationally demanding. Once a suitable supercell is chosen, a ground-state DFT calculation is first converged, with the spin oriented along the z-axis, for example. We then add small perturbative transverse (x and y) magnetic fields on each atom and perform a single additional ground-state iteration. Alternatively, a ground-state calculation with fixed transverse moments via constraining magnetic fields could be used or more correctly a spatially varying magnetic field could be applied for a short time period. These approaches would allow the z-component of the spin to decrease, allowing nonlinear magnon affects to be studied. The purpose of these applied magnetic fields is to mimic the effect of an external perturbation that excites the magnon modes and thus prepares the system in a state that includes magnon excitations. By carefully choosing the amplitude and direction of these applied fields, we can excite individual magnon modes exclusively.
This initial state is then propagated in time using Equation (1). Nonlinear effects due to applied external fields, such as laser pulses, can be included during these propagations. In this work, we focus on the field-free magnon dynamics; however, interaction with an applied pulse was studied by Singh et al. [26] From the TDKS orbitals, we can calculate the time-varying magnetization density and visualize the magnon in real space and real time. For simplicity, the magnetization density is integrated in a sphere centered on each atom to obtain time-varying atomic moments. If we Fourier transform m x ðtÞ and m y ðtÞ to find the power spectrum, we can recover the magnon frequency.

Computational Details
The KS orbitals are time propagated [32] with a time step of 1.209 attoseconds using the ELK electronic structure code. [33] An 8 Â 8 Â 8 k-point grid is used, and the adiabatic local density approximation (ALDA) functional is used to approximate the XC potential. As in the study by Singh et al., [34] among others, [11,14,15] ALDA can overestimate magnon frequencies, especially in the case of Ni. This is due to the overestimation of the exchange splitting in the ground state. The lattice constants used for constructing the supercell for a) iron and Co 50 Ni 50 are a ¼ 7.2902 a.u. . [35] The length of G þ k vectors is set to 2.917 a.u., and the length of G vectors used for expanding the interstitial density and potentials is 12 a.u. To reduce computational costs, the typical propagation time is 50 À 100 fs; however, for low-energy modes this can be increased if required.
To easily visualize interatomic magnon modes, the time-dependent magnetization density may be integrated within the muffin-tin sphere centered on each atom to define an atomic moment. The Fourier transform power spectrum of these time-dependent atomic moments is then used to identify the magnon frequencies where f ðtÞ is a third-order polynomial damping function [36] included to reduce high-frequency noise introduced by the finite period of simulation. However, this function will also artificially increase the width of any peaks if the simulation time is not sufficiently long. The width of these peaks is related to the magnon lifetime; thus, longer times are needed to accurately extract these features. For more complex magnon modes, the Fourier transform of magnetization density may be performed, and the spatial form as a function of frequency may be examined. At higher energies, the transverse oscillations of the Stoner excitations will also be seen.

Iron
We begin with the study of FCC Fe, where the unit cell is extended along the c-axis to obtain a supercell consisting of four atoms. We expect to see four independent modes as there are 4 q-vectors commensurate with this supercell. Below we excite each mode individually by tailoring the initial perturbation. The modes found in our simulations are shown in the cartoon in Figure 1.
The arrangement of spin for all the four atoms after the initial perturbation can be represented in the xy-plane using ", ! for þy, þx and #, ← for Ày, Àx directions, respectively. As time evolves, these spins will precess in the counterclockwise direction (note that in these simulations, the moment points in the negative z-direction). We will now discuss each mode in more detail.
3.1.1. Goldstone mode: " " " " The initial state required to excite the Goldstone mode is all spins pointing in the same direction; this corresponds to a wavevector q ¼ Γ ¼ 0. This is created by perturbing each atom with the same magnetic field. Propagating in time, we do not see any precession of the spins (Figure 2a) as these are excitations of zero energy/frequency. This can also be seen from the Fourier transform of y-moments in Figure 2b. The m y of all the atoms stays fixed to its initial value (0.044 μ B ) which is the moment induced by the applied magnetic field. The Goldstone mode is simply a tilting of the ground-state magnetization along a new direction, this therefore costs no energy when the system has magnetic isotropy. This mode can also be understood as a consequence of Goldstone's theorem which states that a zero-energy mode must exist when a continuous symmetry is broken. For systems with magnetic anisotropy, where the spins align in a preferential www.advancedsciencenews.com www.pss-b.com direction, this mode will no longer be a zero-energy mode, as tilting of spins off the easy axis will cost a small amount of energy. This may also be understood as a consequence of the violation of Larmor's theorem due to the combination of spin-orbit coupling and the crystal field. [37,38] 3.1.2. Optical mode: " # " # Next an initial state for studying the optical mode is prepared by applying magnetic fields in alternating y and Ày directions on the neighboring atoms. The time propagation of orbitals is performed and transverse moments as a function of time are obtained, as shown in Figure 3a. Here the neighboring atoms oscillate 180 out of phase with each other, but the first and third atoms, and second and fourth atoms behave identically. The wavevector for this excitation is q ¼ X, where X corresponds to the zone boundary of the reciprocal cell of the FCC lattice. As shown in Figure 3a, at each instant in time, m x and m y of each atom are 90 out of phase with each other, and they precess counterclockwise in agreement with the simple picture of magnon oscillation of Equation (2). The frequency of this mode is 110 meV, as can be obtained by performing a Fourier transform of x or y moments, see Figure 3b. Also in Figure 3a, we can also see modulation of the amplitude and frequency of the mode, this is due to interaction with the Stoner continuum of spin-flip excitations.
We label the wavevector of this mode þQ, which is equal to half the vector connecting the Γ and X points of the BZ. In this case, the neighboring atoms have a 90 phase with each other. This implies that the y moment of the first atom and x moment of the second atom will be in phase with each other during time propagation, as shown in Figure 4a. These oscillations correspond to a frequency of 76 meV, as shown in Figure 4b. As expected, this mode has a lower frequency compared with the optical mode and less interaction with the Stoner continuum.
3.1.4. ÀQ ¼ ÀΓX/2 mode: " ← # ! Finally, for FCC Fe, the dynamics of the ÀQ mode is shown in Figure 5a. Similar to þQ mode, the neighboring atoms are now À90 (or 270 ) out of phase with each other. The þQ and ÀQ modes are degenerate modes and oscillate with the same frequency of 76 meV (see Figure 5b). For systems with a strong Dzyaloshinskii-Moriya interaction, the degeneracy of these two modes would be lifted.
To summarize this first section, we have demonstrated for a simple elemental ferromagnet that RT-TDDFT may be used to calculate magnon dynamics. The four modes seen in our simulations, and shown in Figure 1, behave as expected according to Equation (2). Before moving to more sophisticated materials, we  cross-check the frequencies found in our RT-TDDFT calculations against those calculated using LR-TDDFT. By definition, RT-TDDFT in the linear regime should agree with LR-TDDFT. In Figure 6, we plot the magnon dispersion obtained with LR-TDDFT along the ΓX direction along with those obtained using RT-TDDFT, shown as the two black dots in Figure 6. As shown, there is good agreement between the two methods, validating our real-time calculations. We often find that LR-TDDFT calculations of magnon frequencies are difficult to converge due to the numerical difficulty in capturing the pole structure of the response function and the sum-over-states nature of the LR equations, hence the numerical noise shown in Figure 6. Now with the machinery built to interpret the time-varying moments in RT-TDDFT calculations, we can observe the element-specific behavior in multimagnetic atomic systems.

CoNi
A more complex magnetic material is Co 50 Ni 50 , which is a multisublattice ferromagnet. In the ground state, the Co atom has a magnetic moment of 1.87 μ B , whereas the Ni atom is 0.71 μ B . We again construct a 4-atom supercell that allows the same wavevectors q ¼ Γ, 1=2ΓX, X as previously studied (this is due to the L1 0 primitive cell simply being an extension of the FCC lattice).
In such a magnetic system, we expect exchange coupling between the Co atoms, between the Ni atoms, and possibly between the Co and Ni atoms. Such an interaction would lead to coupled magnon modes; however, we do not know a priori whether coupled or uncoupled modes exist in our system. Therefore, we can perform RT-TDDFT for various initial states and observe what modes are excited.   The four initial states prepared for the spin system of Co-Ni-Co-Ni (and their decomposition into the observed modes) are 1) Pure Co mode

4) Pure Co and pure Ni modes together
The atomic moment dynamics for each of these initial states is shown in Figure 7. The first observation is that we see both coupled and uncoupled modes for the q vectors investigated. In particular, we found decoupled modes for q ¼ AE1=2ΓX, where only the Co atoms oscillate together, and another where only the Ni atoms oscillate together. This is most clearly shown in Figure 7d, which is the same initial configuration as the þQ case in iron. If coupled modes existed for this wavevector, one would see the same frequency oscillations in both Ni and Co moments; however, very different behavior of the two elements is observed. To understand why we observe decoupled modes, we consider the Heisenberg Hamiltonian. At wavevectors q ¼ AE1=2ΓX, the effective exchange fields acting on each atom from their nearest neighbors of the other species cancel. This allows decoupled element-specific modes to form.
We can excite these element-specific modes individually, as was done for Co, as shown in Figure 7a, and for Ni, as shown in Figure 7b. In both cases, the Co/Ni atoms show transverse oscillations which are 180 out of phase with each other, whereas the x and y moments on each atom are 90 out of phase, as would be expected for magnon excitations. The interaction with the Stoner continuum is stronger for these than previously in Fe, as shown in the large deviations from simple sinusoidal behavior. The frequency of the pure Ni mode obtained by Fourier transform is 621 meV. This frequency is %200 meV higher in energy when compared with the same mode in a 4-atom nickel supercell, implying that the presence of Co has increased the beyond-nearest neighbor exchange interactions. Similar behavior is observed for the Co atoms as well. www.advancedsciencenews.com www.pss-b.com All four modes can also be excited at once by choosing a particular initial state such that it overlaps with all modes. The dynamics of this complex superposition of magnon modes is shown in Figure 7c. Here, we may also see the Goldstone and the optical modes. In multisublattice systems, the Goldstone mode must preserve the ratio of the ground-state magnetic moments. For Co 50 Ni 50 , it means that the magnitude of the transverse moments on Co must be a factor of 2.6 greater than that of Ni. The optical mode is shown in Figure 7c as a higherfrequency mode on top of pure cobalt mode oscillations. For this optical mode, the Co and Ni moments oscillate 180 out of phase with each other at a frequency of 680 meV. As this frequency is close to the pure Ni mode, it is harder to distinguish in the Ni oscillations in Figure 7c, but it can be resolved clearly if the moments on both Ni atoms are summed together, i.e., cancelling the pure Ni mode.
At thermal equilibrium, the occupations of these modes will be given by the Bose-Einstein occupation weights. In this case, as the Co modes are lower in energy, they will be more strongly occupied than the pure Ni modes, leading to different dynamics in the two magnetic sublattices. A similar situation was found [26] in Fe 50 Ni 50 which would explain the experimentally observed difference in demagnetization times between the two sublattices. [39]

Mn 3 Ga
Finally, we investigate the case of Mn 3 Ga, a ferrimagnet. In the ground state, two Mn atoms, Mn II and Mn III , are equivalent and have a magnetic moment of 2.01 μ B , whereas the other Mn atom, Mn I , is antiferromagnetically coupled with others and has a moment of À2.46 μ B .
For this system, we will investigate the magnon branches at the Γ point; thus, only a single primitive cell is required. Generally, the number of magnon branches is equal to the number of magnetic atoms in the primitive unit cell. For example, this was observed in Co 2 MnSi by Buczek et al. [31] using the linear-response TDDFT simulations, where they observed three different modes. Hence, we expect three modes due to the three Mn atoms, although very small oscillations on the Ga atom were observed for one mode.
The first mode to discuss is again the Goldstone mode. As we saw for Co 50 Ni 50 , for multisublattice systems, the Goldstone mode must preserve the ratio of the ground-state moments. For ferrimagnets, this ratio can be negative. In this case, the transverse Mn I moment must point in the opposite direction to those on the Mn II and Mn III atoms, where the sum must equal the net induced moment. To illustrate this point, say all Mn atoms initially have an induced moment of 0.005 μ B . The Goldstone mode will distribute this moment such that the Mn II and Mn III atoms have 0.0192 μ B each, whereas Mn I will have À0.0234 μ B . This distribution preserves the ratio and sums to the correct net moment. Thus, the moment on Mn I in fact reverses. However, this initial configuration must also excite another mode which oscillates around these Goldstone mode shifts.
This second mode is shown in Figure 8a,c for the dynamics of the Mn I atomic moments and the Mn II and Mn III moments, respectively. Similar to the previous optical modes, these two  The frequency of this mode is %70 meV. The Goldstone shift from an initially negative moment to the optical mode oscillating around a positive value is also shown for the m y moment in Figure 8a. The final mode is less pronounced. It is shown in Figure 8b,d, although the Mn I dynamics, or lack thereof, is shown to emphasize that this mode only involves the Mn II and Mn III atoms. This mode is similar to the Fe optical mode studied earlier where the two ferromagnetically coupled atoms oscillate 180 out of phase to each other. The striking feature of Figure 8d is how disrupted this mode is, indicating strong Stoner damping. Despite this, the signature of magnon oscillation, i.e., the phase differences between the atomic x and y moments, can still be observed. From this calculation, it is hard to assign a frequency to this mode; however, it is in the region of 200 meV. This mode likely lies in the Stoner continuum and hence would be difficult to distinguish in a LR-TDDFT calculation.

Conclusions
We began by studying a simple ferromagnet (FCC Fe) and observed the precession of the magnon modes that may be excited in a 4-atom supercell. We also saw that these magnons have q wavevectors commensurate with the supercell. The frequency of the spin-wave oscillations was compared against LR-TDDFT to validate the RT-TDDFT supercell approach. We then studied the more complex magnon modes of a multisublattice magnetic material (Co 50 Ni 50 ) where we observed that in certain cases, the Ni and Co atoms decoupled from each other, leading to element-specific modes. We also demonstrated the case where a superposition of several magnon modes can be studied altogether. Finally, we explored the magnon dynamics of a ferrimagnetic system (Mn 3 Ga), where the real-time approach revealed how the amplitude of transverse oscillations varies from atom to atom for a particular magnon mode. We also observed a strongly damped ferromagnetic Mn mode, which would be difficult to resolve in a linear-response calculation. In all cases, visualizing the magnon modes in real-time and real-space gave important insights into their behavior.
The goal of this work is to lay the foundations of RT-TDDFT as an ab-initio method for understanding and predicting magnon dynamics in many-electron systems. Such an approach is required to capture important quantum mechanical and nonlinear effects. We plan to apply this method to study how magnons behave in various dynamical situations, such as excitation of magnons via the inverse Faraday effect (IFE) where the laser frequency could be tuned for element-specific excitations, propagation of magnon wavepackets in materials, altering fundamental magnonic properties using ultrafast laser pulses as was done in the study by Singh et al., [26] or excitation of spin waves by spin current flowing across an interface. Furthermore, for finite systems such as magnetic clusters, the real-time magnon approach can be used to visualize spin-flip excitations.