## Abstract

To gain control over magnetic order on ultrafast time scales, a fundamental understanding of the way electron spins interact with the surrounding crystal lattice is required. However, measurement and analysis even of basic collective processes such as spin-phonon equilibration have remained challenging. Here, we directly probe the flow of energy and angular momentum in the model insulating ferrimagnet yttrium iron garnet. After ultrafast resonant lattice excitation, we observe that magnetic order reduces on distinct time scales of 1 ps and 100 ns. Temperature-dependent measurements, a spin-coupling analysis, and simulations show that the two dynamics directly reflect two stages of spin-lattice equilibration. On the 1-ps scale, spins and phonons reach quasi-equilibrium in terms of energy through phonon-induced modulation of the exchange interaction. This mechanism leads to identical demagnetization of the ferrimagnet’s two spin sublattices and to a ferrimagnetic state of increased temperature yet unchanged total magnetization. Finally, on the much slower, 100-ns scale, the excess of spin angular momentum is released to the crystal lattice, resulting in full equilibrium. Our findings are relevant for all insulating ferrimagnets and indicate that spin manipulation by phonons, including the spin Seebeck effect, can be extended to antiferromagnets and into the terahertz frequency range.

## INTRODUCTION

In solids, vibrations of the crystal lattice have a significant impact on the orbital dynamics of the electrons (Fig. 1A). They strongly modify properties such as electrical conductivity and may even cause insulator-to-metal transitions (*1*). Likewise, the interplay between phonons and electron spins (Fig. 1A) is relevant for equally drastic phenomena including colossal magnetoresistance, femtosecond magnetization control (*2*–*6*), and the spin Seebeck effect (*7*, *8*). Recently, ultrafast optical techniques have provided new insights into the coherent coupling of individual phonon and spin modes (*9*–*14*).

Despite this progress, the microscopic origins of spin-phonon interactions remain an intriguing problem. Even the equilibration of crystal lattice and electron spins, arguably the conceptually simplest collective process, is far from being understood. This notion is highlighted in the model system yttrium iron garnet Y_{3}Fe_{5}O_{12} (YIG), which is one of the best-studied magnetically ordered solids (*15*, *16*) and ubiquitous in the field of magnonics (*8*). Only estimates of the time constant of spin-phonon equilibration exist, and they extend over as many as six orders of magnitude from ~1 μs (*17*) to ~250 ps (*18*, *19*) and down to ~1 ps (*20*). Here, we introduce an approach based on resonant phonon excitation that allows us to directly probe the interactions between the crystal lattice and electron spins over multiple time scales.

### Experiment

A schematic of our terahertz-pump magneto-optic–probe experiment is shown in Fig. 1A and detailed in Materials and Methods and fig. S1A. An incident, intense, ultrashort terahertz pump pulse (*21*) (photon energy of ~0.1 eV, duration of ~250 fs; see fig. S1B) selectively excites the crystal lattice by resonantly driving infrared-active transverse optical (TO) phonons. The impact on the sample’s magnetic order is monitored from femtoseconds to milliseconds by measuring the magneto-optic Faraday rotation of a time-delayed probe pulse (see Fig. 1A and fig. S1A). In addition, we measure isotropic changes in the optical transmittance of the sample on ultrafast time scales.

As samples, we choose model systems for spin-wave dynamics in magnetic insulators (*7*, *8*, *16*, *22*, *23*): pure YIG and bismuth/gallium-substituted YIG (BiGa:YIG). In these ferrimagnets, magnetic Fe^{3+} ions at a and d sites in the unit cell (Fig. 1B) comprise two nonequivalent, ferromagnetic sublattices with magnetization *M*_{a} and *M*_{d}, respectively, which couple antiferromagnetically. The 2:3 ratio of a to d sites results in a nonzero net magnetization *M*_{a} + *M*_{d} below the Curie temperature *T*_{C}. The Faraday rotation θ of the probe pulse is determined by (*24*)(1)where *a*_{a} and *a*_{d} are the sublattice magneto-optic coefficients. Compared to YIG, the Faraday effect of BiGa:YIG is enhanced by about one order of magnitude and of opposite sign (see Materials and Methods and fig. S1C) (*24*).

In both materials, we expect that the analysis of the pump-induced dynamics is simplified due to the following features. First, the sublattice magnetizations *M*_{a} and *M*_{d} almost entirely arise from the spin rather the orbital angular momentum of the electrons. At room temperature, the *g* factor of YIG differs by only 0.13% from that of the free electron (*25*). Second, the sizeable electronic band gap (2.85 eV for YIG) implies that the electronic orbital degrees of freedom remain in their ground state, thereby significantly reducing the complexity of the pump-induced processes.

Figure 2A displays the absorptance spectrum of a 15-μm-thick BiGa:YIG film from 10 to 45 THz, where absorption is known to be due solely to infrared-active phonons (*26*). Our ab initio calculations show that the pump pulse centered around 21 THz (red spectrum in Fig. 2A) predominantly excites long-wavelength TO(Γ) lattice normal modes characterized by an asymmetric Fe-O stretch vibration (see Fig. 1B, fig. S4, and Materials and Methods).

## RESULTS AND DISCUSSION

### Ultrafast spin dynamics

Figure 2B shows the relative pump-induced change Δθ/θ_{0} = θ(*t*)/θ_{0} − 1 in the Faraday rotation as a function of the delay *t* since excitation. Here, θ_{0} = θ(−2 ps) refers to the equilibrium case. When pumping off the TO(Γ) phonon resonances, at around 38 THz (blue pump spectrum in Fig. 2A), a relatively small signal is found (blue curve in Fig. 2B). In marked contrast, we witness a response more than one order of magnitude stronger for resonant phonon excitation around 20 THz (red pump spectrum in Fig. 2A): an ultrafast single-exponential drop of the Faraday signal with a time constant as short as τ_{fast} = 1.6 ps (red curve in Fig. 2B). On much longer time scales, an additional exponential reduction of θ with a time constant of τ_{slow} = 90 ns is found (Fig. 2C). Recovery back to the initial state occurs over about 1 ms (fig. S2D).

Because the pump-probe signal grows linearly with the pump fluence (inset of Fig. 2C), excitation is dominated by one-photon absorption, whereas strong-field effects such as field or impact ionization are negligible (*27*). This notion is further supported by Fig. 2A, which demonstrates that the sample absorptance and transient Faraday rotation are found to depend on the pump frequency in a very similar way.

We emphasize that almost identical dynamics are observed when a pure YIG sample, instead of the BiGa:YIG film, is used (see fig. S2, B to D). Additional control experiments confirm that, as soon as the pump pulse has left the sample, the transient Faraday signal Δθ(*t*) reliably reflects the dynamics of the sublattice magnetizations *M*_{a} and *M*_{d} with time-independent coefficients *a*_{a} and *a*_{d} in Eq. 1 (see Materials and Methods and fig. S3).

Figure 2B also shows the relative optical transmittance change of the BiGa:YIG sample for the case of resonant pumping. The signal starts with a peak-like feature whose width coincides with that of the terahertz pump pulse. Already for delays *t* > 1 ps, the signal is almost constant, changing by less than 10% in the subsequent evolution. This relaxation to a quasi-steady state is substantially faster than seen for the transient Faraday rotation, which still doubles its value at *t* = 1 ps in the following 5 ps (Fig. 2B). Furthermore, the transmittance signal is found to be independent of the sample magnetization (fig. S3C). These features suggest that the transmittance change predominantly monitors the redistribution dynamics of the pump-deposited energy in the crystal lattice.

Although the relative signal changes are small and still in the perturbative regime, our results in Fig. 2 show a proof of concept that resonant phonon excitation provides an ultrafast manipulation of magnetic order. It only involves the crystal lattice and electron spins, yet no electronic orbital degrees of freedom (Fig. 1A). The picosecond spin dynamics observed here (Fig. 2B) are unexpected because they are five orders of magnitude faster than the spin coherence lifetime (>0.1 μs) of YIG, which is known to be one of the longest among magnetically ordered materials (*8*, *15*).

### Temperature dependence

To characterize the transient states established on the fast (τ_{fast} = 1.6 ps) and slow (τ_{slow} = 90 ns) time scales, respectively, we increase the sample temperature *T*_{0} from 300 to 420 K. Simultaneously, we measure the equilibrium Faraday rotation θ_{0} as well as its pump-induced change Δθ at ultrashort (*t* = 10 ps; Fig. 2B) and long delays (*t* = 1 μs; Fig. 2C) after phonon excitation.

The resulting θ_{0} versus *T*_{0} has the typical shape (*15*) of a ferrimagnet’s static magnetization curve (Fig. 3A). Its slope ∂θ_{0}/∂*T*_{0} steepens with rising *T*_{0} until the transition into the paramagnetic phase occurs at the Curie temperature *T*_{C} = 398 K of BiGa:YIG. In contrast, the pump-induced Faraday signal at *t* = 1 μs (Fig. 3B) increases with *T*_{0} and reaches a maximum right below *T*_{C}, reminiscent of the derivative of the static curve (Fig. 3A). Indeed, we find that Δθ(1 μs) versus *T*_{0} closely follows (∂θ_{0}/∂*T*_{0})Δ*T* (Fig. 3B). Here, Δ*T* = 0.39 K is the increase in equilibrium temperature as calculated from the energy density deposited by the 1-μJ pump pulse (see Materials and Methods). The agreement of both curves shows that ~1 μs after pumping, the BiGa:YIG film is in full thermodynamic equilibrium characterized by temperature *T*_{0} + Δ*T*.

Figure 3B reveals that the changes in the Faraday signal at 10 ps are systematically yet nonuniformly smaller than at 1 μs, in agreement with Fig. 2 (B and C). While for temperatures *T*_{0} < 380 K, Δθ(10 ps) amounts to roughly 80% of Δθ(1 μs), this ratio decreases strongly when *T*_{0} is increased. For example, slightly above 380 K (where both curves reach their maxima) and slightly below 400 K, Δθ(10 ps) has reduced to only 50 and 20% of Δθ(1 μs), respectively. Consequently, Δθ(10 ps) versus *T*_{0} in Fig. 3B cannot be made to agree with Δθ(1 μs) versus *T*_{0} by a simple global rescaling of the Δθ(10 ps) values. In particular, the shape of Δθ(1 μs) versus *T*_{0} agrees much better with a suitably scaled derivative ∂θ_{0}/∂*T*_{0} of the static Faraday rotation (black curve in Fig. 3B) than Δθ(10 ps) versus *T*_{0}. Therefore, at time *t* = 10 ps after pump excitation, the spin system is in a state that is significantly different from the equilibrium state found at *t* = 1 μs.

### Analysis of spin couplings

To understand the microscopic mechanism driving the ultrafast change in magnetic order after resonant phonon excitation (Fig. 2B), we note that solids exhibit only three fundamental spin couplings. They can be understood as effective magnetic fields exerting torques on spins. In the following, we discuss all of them in terms of their capability to modify the sublattice magnetizations *M*_{a} and/or *M*_{d} on ultrafast time scales.

First, spin-orbit (SO) coupling is often considered to dominate the ultrafast demagnetization of laser-excited ferromagnetic metals in the absence of transport (*28*). In our experiment, however, we find identical ultrafast dynamics for pure YIG and BiGa:YIG (see fig. S2B), although the Bi ions are known to increase SO coupling. This notion is manifested in the magneto-optic effects, which are enhanced by more than one order of magnitude (*24*). Therefore, SO coupling plays a negligible role in the ultrafast spin dynamics seen in Fig. 2B.

The same conclusion can be drawn for the second type of coupling, spin-spin magnetic dipole (SSMD) interaction, which has a comparable strength to SO coupling in YIG (*15*). We note that SO and SSMD coupling also account for magnetic anisotropy, which can undergo strong photo-induced changes in the garnet YIG:Co (*29*). Ultrafast laser-induced modification of the anisotropy field of a YIG:Co film was even shown to drive precessional magnetization switching (*29*). The reported time constant of this process (~20 ps) is, however, still more than one order of magnitude longer than the ultrafast dynamics (τ_{fast} = 1.6 ps) seen here. In addition, we do not observe any sign of coherent precession, which is typical of anisotropy changes in the perturbative regime. Consequently, neither SO nor SSMD coupling can explain the spin dynamics observed on the τ_{fast} scale here.

We finally consider the only remaining fundamental spin coupling, isotropic exchange, which is the strongest spin interaction in most magnets (*2*). In YIG, it is responsible for the ferrimagnetic order and the high magnon frequencies extending to more than 20 THz at room temperature (*15*, *30*). Therefore, exchange interaction may well account for the picosecond spin dynamics observed here (Fig. 2B). Because it conserves the total spin angular momentum, the magnetization changes of a and d sublattices cancel each other. Despite Δ*M*_{a} + Δ*M*_{d} = 0, the resulting change in Faraday rotation is nonzero because the magneto-optic coefficients *a*_{a} and *a*_{d} of YIG differ significantly (see Eq. 1) (*24*). In summary, our analysis of the fundamental spin couplings shows that only isotropic exchange interaction is capable of explaining the picosecond spin dynamics seen in our experiment.

### Model

We suggest the following scenario of exchange-mediated spin-phonon coupling: The resonantly excited TO(Γ) phonons decay into other modes significantly faster than the time constant τ_{fast} = 1.6 ps of the measured ultrafast spin dynamics (*26*), thereby heating up the crystal lattice. The assumption of ultrafast phonon relaxation is consistent with our transient transmittance data (Fig. 2B), which indicate that the phonons reach an approximately stationary state within less than 1 ps. As the heat capacity of the crystal lattice of YIG is two orders of magnitude higher than that of the spin system (*31*), the lattice acts akin to a bath whose temperature is suddenly increased from *T*_{0} to *T*_{0} + Δ*T*.

As a consequence of the ultrafast lattice heating, the O^{2−} ions, which are the lightest, will undergo additional random deflection Δ*u*(*t*) (Fig. 4A). This perturbation modulates the superexchange (*32*) of adjacent a-Fe^{3+} and d-Fe^{3+} spins and the associated coupling constant *J*_{ad} by (*16*, *33*, *34*)(2)We now put this model to test by conducting calculations of both dynamic and equilibrium states and by their comparison to our experimental observations on the ultrafast (τ_{fast}; Fig. 2B) and the slower time scale (τ_{slow}; Fig. 2C).

### Dynamics on the τ_{fast} scale

To determine the ultrafast rate of change of *M*_{a} und *M*_{d}, we perform atomistic spin-dynamics simulations based on ~10^{6} coupled spin equations of motion (see Materials and Methods) (*30*). We include time-dependent exchange parameters through Eq. 2, where Δ*u*(*t*) is assumed to be random with a variance given by the pump-induced temperature increase Δ*T* of the ionic lattice.

Simulation results are shown in Fig. 4B. At times *t* < 0, the magnetizations of both sublattices fluctuate around their constant means *M*_{a0} and *M*_{d0}. However, when fluctuations Δ*J*_{ad}(*t*) of the exchange parameter are switched on at *t* = 0, *M*_{a} decreases linearly with time until Δ*J*_{ad}(*t*) is switched off. We find the opposite behavior for Δ*M*_{d}. The sum curve Δ*M*_{a} + Δ*M*_{d} has more than one order of magnitude smaller amplitude. It arises from the much weaker random fields of the bath due to spin couplings other than isotropic exchange (see Materials and Methods). Figure 4B demonstrates that thermal modulation of *J*_{ad} induces demagnetization of the two spin sublattices by the same amount, and thus transfers energy into the spin system.

The slope of the simulated Δ*M*_{a}(*t*)/*M*_{d0} (Fig. 4B) can be compared directly to the initial slope of the pump-probe signal Δθ(*t*)/θ_{0} (Fig. 2B), which amounts to approximately −0.1% ps^{−1}. Agreement between experiment and theory is obtained by choosing ∂*J*_{ad}/∂*u*~10 *J*_{ad} Å^{−1} (see Materials and Methods). This value and the theory of superexchange (*32*) imply that the overlap integral of the ground-state wave functions of adjacent O^{2−} and Fe^{3+} ions changes by 100% when their distance changes by ~0.4 Å. As the distance between the two ions is one order of magnitude smaller than the lattice constant *a* of YIG (*15*), we estimate ∂*u*/∂*a*~0.1 and obtain ∂*J*_{ad}/∂*a*~(∂*J*_{ad}/∂*u*)(∂*u*/∂*a*)~1 *J*_{ad} Å^{−1}. This result is in agreement with recent ab initio calculations (*34*) of the exchange coupling constants as a function of *a*, which yielded ∂*J*_{ad}/∂*a*~0.7 *J*_{ad} Å^{−1}. Therefore, phonon-modulated exchange coupling successfully and consistently explains the rate (∂Δθ/∂*t*)/θ_{0} at which magnetic order of YIG is observed to reduce immediately after resonant phonon excitation.

According to our model, ∂*M*_{a}/∂*t*, ∂*M*_{d}/∂*t*, and, thus, ∂Δθ/∂*t* (Eq. 1) are proportional to the difference between lattice and spin temperature (Eq. 10), in agreement with the linear pump fluence dependence of the transient Faraday rotation Δθ seen in the experiment (inset of Fig. 2C). Consequently, the measured exponential decay of Δθ to a constant value (Fig. 2B) indicates that the electron spins approach the lattice temperature *T*_{0} + Δ*T* with the time constant τ_{fast} = 1.6 ps.

Note that the subsequent slower dynamics evolve with a time constant of τ_{slow} = 90 ns. We conclude that for all pump-probe delays *t* ≪ τ_{slow}, the transient changes in *M*_{a} and *M*_{d} arise solely from inter-sublattice exchange coupling. Because this interaction conserves the total spin angular momentum, the dynamics are under the constraint Δ*M*_{a} + Δ*M*_{d} = 0 throughout that time window.

### Dynamics on the τ_{slow} scale

As indicated by our temperature-dependent measurements (Fig. 3), the garnet film approaches full equilibrium at temperature *T*_{0} + Δ*T* with the time constant τ_{slow}. The total magnetization *M*_{a} + *M*_{d} of this state is lower than for *t* ≪ τ_{slow}, where the change Δ*M*_{a} + Δ*M*_{d} is still zero. Thus, transfer of angular momentum from the spin system to the lattice has certainly occurred for times *t* ~ τ_{slow} and above.

The only microscopic spin couplings capable of changing the total spin and, thus, magnetization *M*_{a} + *M*_{d} are SO and SSMD coupling. If SO coupling were dominant, then one would expect a difference in the spin dynamics in pure YIG and in BiGa:YIG in which SO coupling is enhanced by partial Bi substitution. However, the signal-to-noise ratio of the Faraday signal seen on the τ_{slow} scale of pure YIG (fig. S2B) currently does not allow us to draw a conclusion about which of the two couplings is dominant.

### Constrained spin state

As discussed above, the spin system is under the constraint of constant magnetization *M*_{a} + *M*_{d} = *M*_{a0} + *M*_{d0} for times *t* ≪ τ_{slow}, whereas it is in full, unconstrained equilibrium for *t* ≫ τ_{slow} (Fig. 2C). Assuming that for *t* ≫ τ_{fast}, the spins and crystal lattice are characterized by a single temperature *T*_{0} + Δ*T*, both the constrained and unconstrained spin state can be fully described by equilibrium statistical physics (see Materials and Methods).

Using a mean-field approximation, we obtain the change (∂*M*_{l}/∂*T*_{0})Δ*T* in the sublattice magnetizations versus *T*_{0} for heating with and without the constraint of conserved *M*_{a} + *M*_{d} (Fig. 5A). The constrained ∂*M*_{a}/∂*T*_{0} of minority spin-sublattice a (green line in Fig. 5A) closely follows its unconstrained counterpart. In contrast, the constrained ∂*M*_{d}/∂*T*_{0} of majority spin-sublattice d (blue line in Fig. 5A) is systematically smaller than the unconstrained ∂*M*_{d}/∂*T*_{0} because sublattice d can only demagnetize as much as sublattice a (Δ*M*_{d} = −Δ*M*_{a}). In particular, the difference between constrained and unconstrained ∂*M*_{d}/∂*T*_{0} increases considerably when the temperature *T*_{0} approaches the Curie point (Fig. 5A). We emphasize that this behavior is in excellent qualitative agreement with the measured Faraday signals Δθ(10 ps) and Δθ(1 μs) versus *T*_{0} (Fig. 3B).

It is interesting to note that in the equilibrium formalism used here, the constant total magnetization is reinforced by a Lagrange multiplier, which can be interpreted as a virtual homogeneous magnetic field (see Eq. 13). One could also term this field “spin pressure,” in analogy to the pressure that is required to keep a gas of particles in a bottle of constant volume while the gas temperature is increased. In our experiment, the spin pressure builds up on the time scale given by τ_{fast} and is released by angular momentum transfer to the crystal lattice on the time scale τ_{slow}.

### Picture of spin-phonon equilibration

In summary, our dynamic and equilibrium calculations are fully consistent with the time scales, fluence dependence, temperature dependence, and magnitude of the transient Faraday rotation found in the experiment. This agreement leads us to the following picture of the flow of energy and angular momentum of phonon-pumped YIG (Fig. 5B): [1] the pump pulse excites zone-center TO(Γ) phonons that [2] decay within the pump duration, thereby increasing the crystal lattice temperature. The additional thermal modulation of a-d exchange induces [3a] transfer of angular momentum between a- and d-spin sublattices and implies [3b] energy transfer from the phonon to the spin system with the time constant τ_{fast} = 1.6 ps (Fig. 2B). The excess magnetization of this constrained state decays by [4] transfer of angular momentum and energy between crystal lattice and spins with the time constant τ_{slow} = 90 ns, resulting in full, unconstrained equilibrium (Fig. 2C). The angular momentum transfer is mediated by SO and/or SSMD interactions (*22*), which do not conserve *M*_{a} + *M*_{d}.

## CONCLUSION

We use a terahertz pump magneto-optic–probe experiment to investigate spin-lattice equilibration from the femtosecond to the microsecond time scale. Combined with a spin-coupling analysis and atomistic spin-dynamics simulations, our results reveal that the speed of spin-phonon relaxation in ferrimagnetic insulators critically depends on whether one refers to energy or angular momentum.

On a picosecond scale, spins and phonons reach a quasi-equilibrium through phonon-modulated exchange interaction. It induces redistribution of energy among phonons and spins, whereas angular momentum is only transferred between spins. On a much longer time scale in the nanosecond range, the additional transfer of angular momentum between spins and crystal lattice results in full equilibration. The quasi-equilibrium persists on the intermediate time scale and can be understood as a transient spin state with elevated temperature yet unchanged net magnetization. It could be considered as a realization of magnon populations with nonzero chemical potential, which were recently introduced in the theory of magnon transport by the spin Seebeck effect (*19*).

Transfer of angular momentum between spin sublattices is also a key process in magnetic switching of ferrimagnetic metallic alloys by optical femtosecond laser pulses (*35*–*38*). This process is considered to follow the preferential ultrafast demagnetization of one of the two spin sublattices. Note that the flow of angular momentum can proceed through various mechanisms such as direct exchange torque (*36*) or nonlocal processes (*38*) including the transport of spin-polarized electrons. Our results reveal a new mechanism for manipulating the exchange interaction on ultrafast time scales and highlight the important role phonons can play in the direct exchange dynamics between spin sublattices.

In terms of applications, our results suggest that resonant phonon excitation is a new pathway to the preparation of spin states with increased temperature yet unchanged total magnetization. This route may be particularly interesting for the manipulation of antiferromagnetic order (*2*), where spin angular momentum is inherently conserved. Finally, the ultrafast spin-lattice coupling of YIG implies that the magnon temperature follows the phonon temperature with a delay on the order of only 1 ps, thereby shifting the cutoff frequency of the bulk spin Seebeck effect (*7*) to the terahertz range.

## MATERIALS AND METHODS

### Samples

Our iron garnet films were grown by liquid-phase epitaxy on substrates of gadolinium gallium garnet Gd_{3}Ga_{5}O_{12} (GGG). To prevent pump absorption by the substrate, we transferred several iron-garnet films on diamond windows after mechanically removing the GGG. Test measurements confirmed that GGG pump absorption was irrelevant to the ultrafast dynamics.

We studied two types of garnet samples: pure yttrium iron garnet YIG and bismuth/gallium-substituted iron garnet Bi_{x}Y_{3−x}Ga_{y}Fe_{5−y}O_{12} (*x* = 1.53, *y* = 1.33, BiGa:YIG). They have, respectively, a Curie temperature of 545 and 398 K, and in-plane and out-of-plane magnetic anisotropy. Film thicknesses covered a range from 7 to 20 μm. The chemical composition of the garnet films was determined by x-ray fluorescence measurements. Substitution of Y by Bi in the BiGa:YIG sample enhanced the magneto-optic Faraday effect by one order of magnitude (*24*). Typical hysteresis loops obtained by measuring the static Faraday rotation of the probe beam are shown in fig. S1C.

Terahertz spectra of the sample absorptance (which equals one minus the reflected and transmitted power, both normalized to the incident power) were obtained by combined reflection and transmission measurements with a broadband terahertz time-domain spectrometer. Figure 2A shows that maximum absorptance occurs at 21 THz, which is 2 THz above the highest-frequency TO(Γ) phonon resonance (*26*). According to Hofmeister and Campbell (*26*), absorption in this frequency range solely arises from TO phonons, while crystal-field and charge-transfer gap transitions are located at much higher frequencies (*24*).

### Experimental design

A terahertz pump pulse was used to resonantly excite long-wavelength TO phonons of the iron-garnet film under study. The instantaneous magnetic state of the sample was determined by a subsequent optical probe pulse measuring the magneto-optic Faraday effect. In this way, spin dynamics were monitored over a large range of pump-probe delays ranging from femtoseconds to milliseconds. Details of our setup are shown in fig. S1.

To generate intense terahertz pump pulses, we used an amplified Ti:sapphire laser system delivering pulses (energy of 15 mJ, center wavelength of 800 nm, duration of 40 fs, and repetition rate of 1 kHz) to drive two optical parametric amplifiers that, in turn, generated intense infrared pulses (pulse energy of 1.5 and 1.2 mJ, center wavelength of 1280 and 1400 nm, and duration of 50 and 50 fs, respectively). By difference-frequency mixing of the two infrared pulses in a nonlinear-optical GaSe crystal (thickness of 1 mm) (*21*), we obtained intense terahertz pulses (tunable from 15 to 60 THz, energy of typically 3 to 20 μJ, pulse duration around 250 fs) with a stable carrier-envelope phase. For sample excitation, the terahertz pulse was focused onto the sample surface to a spot with a diameter of 180 μm. Its transient electric field *E*(*t*) was measured by electro-optic sampling (see below). The pulse spectrum was obtained by Fourier transformation of *E*(*t*) and shown in Fig. 2A.

Low-noise probe pulses (8 fs, 800 nm, 0.5 nJ) were derived from the pulse train of the Ti:sapphire seed laser (repetition rate of *f*_{rep} = 80 MHz). They traversed the sample collinearly with the pump after a delay *t*. To measure the probe’s transient polarization rotation Δθ(*t*), we used a polarimeter consisting of a Wollaston-prism polarizer followed by two fast photodiodes. Measurement of the probe ellipticity Δη(*t*) was accomplished by putting a quarter-wave plate before the Wollaston prism. The resulting photodiode current was a train of temporally isolated electrical pulses (duration of 5 ns each, pulse-to-pulse distance of 1/*f*_{rep} = 12.5 ns), which was sampled using a computer-controlled fast analog-digital converter card (*39*). In this way, we were able to measure signals at delays of *t = τ + j*/*f*_{rep} (with integer *j*) between −100 ns and ~1 ms within a single pump shot and with a delay spacing of 1/*f*_{rep} = 12.5 ns. The ultrashort offset τ was set from −2 to 200 ps by a variable mechanical delay stage.

Pump-probe signal traces Δθ_{±}(*t*) were taken for the sample magnetization saturated in opposite directions using an external magnetic field *B*_{ext} of ±15 mT, respectively. To measure samples with in-plane or out-of-plane anisotropy, we oriented the magnetic field 45° with respect to the incident pump and probe beams. To vary the sample temperature, we mounted the sample on a resistive heater. The sample temperature was measured with two thermocouples and an imaging infrared thermometer, all of them yielding consistent temperature values. To measure the terahertz electric field, we substituted the sample by a GaSe crystal and measured the transient ellipticity Δη(*t*)∝*E*(*t*) induced by the electro-optic effect (*21*).

Optionally, our setup also permitted measurement of the transient isotropic transmittance of the sample. For this purpose, we removed all polarization-sensitive components such as the Wollaston prism and wave plate and detected the probe-pulse energy with one of two fast photodiodes of the polarimeter. Considering the small attenuation length (*26*) of ~1 μm of the resonant pump pulse into the iron-garnet sample, we estimate that the temporal blurring of the pump-probe signals due to pump-probe velocity mismatch is smaller than 10 fs and, therefore, negligible.

### Signal analysis

Up to first order in magnetization, the Faraday rotation θ (and likewise the ellipticity η) of the probe’s outgoing polarization is a weighted average (*24*) of the magnetization of the two spin sublattices *l =* a, d, that is(3)Here, *a*_{l} are the local magneto-optic coupling coefficients (Verdet constants) of YIG, *M*_{l} = **e**_{pr} ⋅ **M**_{l} is the magnetization of sublattice *l* projected on the propagation direction unit vector **e**_{pr} of the probe, *d* is the sample thickness, and *b* is an offset arising from any magnetization-independent optical anisotropy. The angular brackets 〈.〉 denote averaging over the probed sample volume.

In equilibrium, the *a*_{l} are constant throughout the probed magnetic volume, 〈*a*_{l}*M*_{l}〉 = *a*_{l}〈*M*_{l}〉. The impact of the pump pulse makes all quantities in Eq. 3 time-dependent, resulting in relatively small changes θ(*t*) – θ_{0} of the signal θ_{0} = θ(−2 ps) measured 2 ps before arrival of the pump pulse. According to Eq. 3, the nonmagnetic offset *b* cancels by calculating the asymmetric part(4)The first term on the right-hand side scales with a weighted average of the sublattice magnetizations along **e**_{pr}, which is the quantity we are interested in. The second term, however, makes a contribution independent of the Δ*M*_{l} and arises from a possible pump-induced modulation Δ*a*_{l} of the Verdet constants. Figure S3 shows that the second term is negligible and that Δθ reflects the true magnetization dynamics, Δθ = ∑*a*_{l0}〈Δ*M*_{l}〉, once the pump pulse has left the sample.

### Excitation profile and heating

Because the absorption length of the terahertz pump pulse (*26*) (~1 μm) is smaller than the sample thickness *d* (~10 μm), pump excitation is inhomogeneous along the *z* axis normal to the film plane. Because the pump-induced signal Δθ(*t*) depends linearly on the absorbed pump fluence *F*_{abs} (inset of Fig. 2C), the pump-induced change in the local Faraday rotation Δθ(*z*,*t*) is proportional to the pump energy density *w*_{abs}(*z*) absorbed in the vicinity of the plane at *z*. That is, one has and Δθ(*z*,*t*) = *C*(*t*)*w*_{abs}(*z*), where *C*(*t*) captures the temporal dynamics. As a consequence, the total pump-induced signal(5)is independent of the shape of the absorption profile *w*_{abs} (*z*). Therefore, without loss of generality, we can consider the absorption profile to be homogeneous with *w*_{abs} = *F*_{abs}/*d*. Our argumentation is valid as long as transport of heat from the YIG film into the substrate is negligible. Figure S2D shows that this assumption is fulfilled for pump-probe delays up to at least 1 μs.

Assuming that the sample has reached (local) thermal equilibrium at *t* = 1 μs, we estimate the average temperature increase Δ*T* (1 μs) of the BiGa:YIG film by dividing *F*_{abs} by the pumped sample volume [*d* × (pump diameter)^{2} × π/4], by the mass density, and by the specific heat capacity (*31*) of YIG. We obtain a value of Δ*T* (1 μs) = 0.39 K for an incident pump energy of 1 μJ, which can be interpreted as the temperature increase of a homogeneously excited YIG film.

We checked the effects due to accumulative heating of the sample by reducing the repetition rate of the pump pulses from 1000 to 500 Hz by means of a mechanical chopper. The magnitude of pump-induced signal that persists after 1 ms was found to be smaller than 10% of the pump-induced signal at *t* ~ 1 μs. It corresponds to a static homogeneous temperature increase of the probing volume of only ~39 mK, which has a negligible influence on our results, including the temperature dependence seen in Fig. 3B.

### YIG ab initio calculations

Our ab initio phonon calculations were performed using the finite-displacement method implemented in the open-source package phonopy (*40*). The electronic-structure calculation and equilibrium geometry search of YIG were performed within the density functional theory framework, using the Vienna Ab initio Simulation Package (*41*) and adopting the generalized gradient approximation (*42*) of the exchange-correlation functional. Our calculations and resulting phonon dispersion relations are detailed in the Supplementary Materials.

From the element-resolved vibrational density of states (fig. S4), we found that a continuum of modes from ~7 to 20 THz contributes to the displacement of the O^{2−} ion with approximately frequency-independent weight. This result implies that the vibration of the O^{2−} ion has a mean frequency of Ω_{O}/2π ~ 14 THz and a correlation time of τ_{O} ~ 10 fs.

### Dynamic spin model

Our consideration of all possible spin-coupling mechanisms (SO, SSMD, exchange; see main text) leads us to the view that the ultrafast regime of the spin dynamics of YIG after phonon excitation arises from isotropic exchange interaction. More precisely, the thermal vibrations of the ionic lattice modulate the a-d exchange interaction (Fig. 4A), thereby transferring spin angular momentum from the a- to d-spin sublattice and quenching each sublattice magnetization at the same rate.

To test this hypothesis, we calculated both the rate of spin angular momentum transfer between sublattices (based on numerical simulations) and the temperature dependence of the sublattice magnetizations (based on an analytical treatment). The starting point is the Heisenberg spin Hamiltonian of YIG (*15*)(6)where is the spin angular momentum operator of the total electron spin of an Fe^{3+} ion (spin quantum number *S* = 5/2) located at site *j* in spin-sublattice *l*. Orbital contributions to the magnetic moments are negligibly small (*25*). In YIG, there are two different crystallographic sites *l* = a, d for the Fe^{3+} ions, which are, respectively, surrounded by six O^{2−} ions (octahedral sites) and four O^{2−} ions (tetrahedral sites). The exchange constants were determined as follows (*15*): If for a given *l*, *l*′, the *j* and *j*′ represent next neighbors, then . Otherwise, . Thus, the exchange part of the Hamiltonian is fully determined by the three constants *J*_{aa}, *J*_{dd}, and *J*_{ad} = *J*_{da}.

Following Eq. 2, the phonon-induced modulation of the exchange interaction is modeled as variation(7)of the exchange coupling constant between *l* = a and *l*′ = d spin sublattices. Here, ∂_{u} = ∂/∂*u*, and Δ*u* is the pump-heat-induced deflection of the O^{2−} ion mediating the superexchange of adjacent a- and d-type Fe^{3+} ions, as depicted in Fig. 4A. This assumption is justified because *J*_{ad} predominantly involves the a-Fe^{3+} and d-Fe^{3+} ions and the O^{2−} ion in between. Because the O^{2−} ion is, by far, the lightest ion in the YIG unit cell, its motion is expected to modulate the a-d superexchange most strongly.

### Atomistic spin-dynamics simulations

To calculate the rate of spin angular momentum transfer between sublattices, we conducted numerical simulations based on Eqs. 6 and 7. For the implementation, we assumed an ensemble of classical spins and replaced each spin operator by . Here, **s**_{i} is a vector of length one, and the compound index (*lj*) is summarized in one index *i*. Likewise, the Hamilton operator *Ĥ* turns into the Hamilton function *H*. For the exchange constants *J*_{aa}, *J*_{dd}, and *J*_{ad}, we used the values 0.33, 1.15, and −3.43 meV, respectively (*15*).

By applying Langevin theory, a stochastic Landau-Lifshitz-Gilbert–type equation is derived for each spin (*30*, *43*, *44*)(8)where γ_{i} is the gyromagnetic constant of a spin with magnetic moment μ_{i}**s**_{i}. Damping due to coupling to nonspin degrees of freedom (the so-called bath) was quantified by λ_{i} = 2 × 10^{–5}. The effective field **B**_{eff,}_{i} acting on **s**_{i} has three contributions: a deterministic part due to the exchange fields of the adjacent spins for temporally constant and two stochastic components **ξ**_{i}(*t*) and Δ**B**_{exch,i}(*t*) that arise from coupling of the spins to the bath. The first field **ξ**_{i}(*t*) transfers spin angular momentum and energy between the spin system and the bath having temperature *T*_{0}. According to the fluctuation-dissipation theorem, in equilibrium, the fluctuations **ξ**_{i}(*t*) are balanced by a friction force, the second, Gilbert-type term in Eq. 8, thereby driving the spin system into an equilibrium state with temperature *T*_{0} (*30*, *43*, *44*). The time constant of this relaxation process is set by the damping parameters λ_{i} and found to be irrelevant on the picosecond scale considered here.

The second stochastic field arises from the phonon-induced modulation of the exchange interaction (Eq. 7). It is distinctly different from **ξ**_{i}(*t*) as it causes energy transfer into the spin system but leaves the total spin angular momentum unchanged. Because the noise field at temperature *T*_{0} is already accounted for by **ξ**_{i}(*t*), the variance of Δ*J*_{ij}(*t*) scales with the temperature increase Δ*T* of the crystal lattice induced by the pump pulse. More precisely, owing to the equipartition theorem, the pump-induced deflection Δ*u* of the O^{2−} ion obeys , where Ω_{O}/2π~14 THz is the mean O^{2−}-ion vibration frequency (see above).

To determine the initial rate of phonon-driven angular momentum transfer between the spin sublattices, we assume phonons and spins to be thermalized at temperatures *T*_{0} + Δ*T* and *T*_{0}, respectively (see main text). Following our ab initio results (see above), Δ*u*(*t*) is assumed to have a correlation function with a width of τ_{O} ~ 10 fs, which is accordingly modeled as(9)The equation of motion (Eq. 8) was integrated numerically for *N* = 64^{3} = 2.6 × 10^{5} unit cells at equilibrium temperature *T*_{0} = 300 K. At each time step *t*, we calculated the two sublattice magnetizations with *l* = a, d. In our numerical implementation, the pump-induced thermal noise of the exchange constants (Eqs. 7 and 9) can be switched on or off over arbitrary time intervals and with variable strength (∂_{u}*J*_{ad})^{2}*ΔT*.

The results of our atomistic simulations are shown in Fig. 4B and fig. S5. We found the demagnetization rate of each sublattice scales according to(10)with constant *f =* 2 × 10^{16} m^{2} J^{–2} K^{–1} ps^{–1}. On the other hand, our experiment shows that the Faraday signal Δθ and thus |∂_{t}*M*_{l}/*M*_{l0}| scale with the pump fluence (see inset of Fig. 2C) and thus the initial Δ*T* as well. In other words, we have(11)where *g* = 1.9 × 10^{−3} K^{−1} ps^{−1} was determined from the initial slope of Δθ/θ_{0} ~ Δ*M*_{l}/*M*_{l0} (Fig. 2B) seen in our experimental results. Comparison of Eqs. 10 and 11 yielded the estimate Å^{−1}.

According to the theory of superexchange (*32*), *J*_{ad} is proportional to the fourth power of the overlap integral of the Fe^{3+} and O^{2−} ions. Therefore, upon changing the distance of the two ions, the overlap integral undergoes relative changes on the order of 2.5 Å^{−1}.

### Mean-field sublattice magnetization

Our previous results strongly indicate that spins and phonons reach quasi-equilibrium for times *t* ≫ τ_{fast} after phonon excitation. As exchange interaction leaves the total spin angular momentum of the electrons unchanged, this thermal state is constrained by the boundary condition of conserved total spin angular momentum. Such constraints are ubiquitous in thermostatics (such as for a gas contained in a closed bottle) and can be treated by equilibrium statistical physics (*45*).

Our goal is to determine the magnetization of the d- and a-spin sublattices of YIG in thermal equilibrium for two cases: the unconstrained (UC) situation without boundary condition and the constrained (C) situation with the boundary condition of constant total spin angular momentum(12)In both cases, we need to calculate the statistical operator of the system, from which the expectation value of an observable *Â* (such as a sublattice magnetization) can be derived. In thermal equilibrium and for both the UC and C cases, the statistical operator is given by (*45*)(13)where β = 1/*k*_{B}*T*, with *k*_{B} being the Boltzmann constant. The partition function *Z* is determined by the normalization condition . In the UC case, we set **p** = 0 in Eq. 13, thereby resulting in the standard canonical statistical operator (*45*). In the C case, the constraint of conserved spin angular momentum (Eq. 12) has to be fulfilled by proper choice of the Lagrange multiplier **p** in Eq. 13. We note that **p** can be interpreted as a virtual homogeneous magnetic field that is adjusted to reinforce the boundary condition of Eq. 12. It is analogous to the chemical potential that keeps the number of particles in a given macroscopic system constant.

The general Heisenberg spin Hamiltonian *Ĥ* of Eq. 6 is too complex for an analytical calculation of . Consequently, we apply the mean-field approximation (*32*), , omit the number , and allow for alignment exclusively along the *z*-axis unit vector **e**_{z}, , as set by the external magnetic field **B**_{ext} = *B*_{ext}**e**_{z}. In addition, we assume the spin order homogeneously on each sublattice *l*, that is, [sublattice approximation (*32*)]. By evaluating the of YIG (see above), we finally obtain the mean-field (MF) Hamiltonian(14)Here, is defined as above, and denotes the number of the next neighbors a site in sublattice *l* has in sublattice *l*′. For YIG, one has *N*_{aa} = 8, *N*_{dd} = 4, *N*_{ad} = 6, and *N*_{da} = 4 (*15*). We can now write the statistical operator of Eq. 13 in a very compact form as(15)where the effective field(16)captures the mean field (first and second terms) and the scalar virtual field *p* = **p** ⋅ **e**_{z} imposed by the constraint of Eq. 12. By taking advantage of the fact that all spin operators in Eq. 15 commutate with each other, we obtain the expectation value(17)for any Fe^{3+} spin of sublattice *l*. Evaluation of this equation leads to an implicit relationship for (18)where *S* = 5/2, and B_{S} is the Brillouin function (*32*). The magnetization *M*_{l} of sublattice *l* is proportional to , where *v*_{l} is the number of *l*-Fe^{3+} ions per unit cell with *v*_{d} = 12 and *v*_{a} = 8.

For the UC case (*p* = 0), Eq. 18 is a system of two coupled equations and solved numerically, thereby yielding the unconstrained for each temperature *T*_{0}. For the C case, the boundary condition of Eq. 12, rewritten as(19)adds a third equation to Eq. 18, which determines *p*. Note that in our experiment, we start from an UC equilibrium given by the at temperature *T*_{0}. Subsequently, the pump pulse increases *T*_{0} by a small amount Δ*T* while keeping the total spin angular momentum constant. Because we are only interested in small relative changes of the , we linearized Eqs. 18 and 19 in terms of Δ*T* and solved analytically for the small changes .

To be consistent with our atomistic spin-dynamics model (see above), we chose the ratio of the coupling parameters *J*_{ad}, *J*_{aa}, and *J*_{dd} to be identical to the values given in (*15*), but rescaled by a factor of 0.36 to fit the experimentally observed critical temperature of *T*_{c} = 398 K. For comparison with our experiment, we finally determined the Faraday rotation resulting from the calculated sublattice magnetizations *M*_{l} (see Eq. 1). The results of this procedure are shown and discussed in Fig. 5A and fig. S6.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/7/eaar5164/DC1

Supplementary Methods

1. Details of YIG ab initio calculations.

Fig. S1. Experimental details.

Fig. S2. Supporting pump-probe data.

Fig. S3. Nonmagnetic signal contributions.

Fig. S4. Ab initio phonon calculations.

Fig. S5. Supporting data of atomistic spin-dynamics simulations.

Fig. S6. Calculated equilibrium sublattice magnetizations.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank G. E. W. Bauer and S. M. Rezende for stimulating discussions.

**Funding:**T.K. acknowledges funding through the grants European Research Council H2020 CoG TERAMAG (grant no. 681917) and DFG SFB/TRR 227 “Spin Dynamics” (project A05). I.R. acknowledges funding through German Federal Ministry of Education and Research (BMBF) grant 05K16BCA “Femto-THz-X.” P.M. and P.M.O. acknowledge financial support from the Swedish Research Council (VR), the K. and A. Wallenberg Foundation (grant no. 2015.0060), and the Swedish National Infrastructure for Computing. M.G. thanks for funding through BMBF grant no. 05K10KEB. R.V.P. acknowledges financial support from the Russian Scientific Foundation (grant no. 16-12-10456). M.W. acknowledges funding by the Max Planck Society.

**Author contributions:**I.R. conceived the original experimental idea. S.F.M., I.R., and T.K. designed the experiment. S.F.M. built the setup, performed measurements, and analyzed data. I.R., A.M.K., and R.V.P. prepared the samples, which were characterized by S.F.M., I.R., A.M.K., A.P., and M.G. The theoretical model was developed by T.K. and S.F.M., with contributions from J.B., P.M., A.P., and P.M.O. Atomistic spin-dynamics simulations were conducted by J.B. Equilibrium magnetization curves were calculated by S.F.M. and T.K. Ab initio calculations of electronic and vibrational properties were performed by P.M. and P.M.O. The manuscript was written by T.K., J.B., and S.F.M., with contributions from A.P., P.M.O., P.M., M.W., and I.R. All authors contributed to the discussions of the results and commented on the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).