Quasi-static strain governing ultrafast spin dynamics

The quasi-static strain (QSS) is the product induced by the lattice thermal expansion after ultrafast photo-excitation. Although the ultrafast spin dynamics driven by the QSS and thermal effects are barely distinguishable in time, they should be treated separately because of their different fundamental actions. By employing ultrafast Sagnac interferometry and the magneto-optical Kerr effect, we demonstrate quantitatively the existence of QSS and the decoupling of two effects counteracting each other in typical polycrystalline Co and Ni films. The Landau-Lifshitz-Gilbert and Kittel equations considering a magnetoelastic energy term showed that QSS, rather than the thermal energy, in ferromagnets plays a governing role in ultrafast spin dynamics. This demonstration provides a way to analyze ultrafast photo-induced phenomena. Ultrafast Sagnac interferometer provides direct access to the lattice dynamics in magnetic materials. The authors leveraged light-matter interaction to unveil the intertwined role of strain and thermal effects on the spin dynamics of ferromagnets.

M agnetoelasticity, which is the coupling of spin and strain, is a universal phenomenon in magnetic materials and enables active control of the spin states by modifying the material dimensions 1 . Thus far, it has been investigated that the spin dynamics after ultrafast spin angular momentum transfer by photo-excitation 2-7 is governed by time-dependent effective fields with origins such as magnetocrystalline 8 , dipole, Zeeman 9 , exchange [10][11][12] , terahertz 13,14 , magnetoelasticity 15,16 , and spin current [17][18][19] . These have been described mainly by electron and spin degrees of freedoms, while the lattice degree of freedom was highlighted only recently 20 .
The increase in lattice temperature by photo-excitation generates two types of strain, which are the propagating strain η p (z, t) with a few-ps temporal width and the quasi-static strain (QSS: η qss (t)) at a penetration length scale with a heat dissipation time. In the recent decade, the interaction of spins with η p has attracted considerable research attention from diverse points of view 15,16,[21][22][23] . On the other hand, QSS has not been a major focus despite its comparable amplitude to that of η p . As principal reasons, the QSS and thermal dependence of materials are inextricable, particularly in photo-induced experiments. Moreover, the QSS is difficult to quantify only with a differential reflectivity that is typical measurement data in a pump-probe experiment. Therefore, great care must be taken in determining η qss (t). In general, the three temperatures model is used to extract the timedependent temperature information of the sub-systems by fitting with the experimental curves of spin and reflectivity dynamics. These temperature profiles can account for the only temperature contribution to the reflectivity, not the strain contribution. This might deliver improper messages and hinder unveiling the new physics of ultrafast dynamics. The earlier work describes the contributions of the QSS and thermal effects to the spin dynamics of a galfenol film depending on the external magnetic field strength 24 . In addition, for an integrated understanding, the experimental quantification of the QSS and its competition with the thermal effect by various experimental conditions remain to be proven. As the QSS and the thermal effects have different contributions to dynamics (the former changes the plasmonic 25 or electronic bands 26 and the latter changes the electron populations), systematic and comprehensive measurement approaches to separate the two effects are required for the complete analysis of the fundamental mechanism of the ultrafast photo-induced spin dynamics.
Using ultrafast Sagnac interferometry (USI) and a magnetooptical Kerr effect (MOKE) instrument, we demonstrate in this work that the QSS in ferromagnetic films governs the ultrafast spin dynamics from the first ps to the ns timescale. In particular, USI is used to prove the existence of QSS by measuring the lattice expansion dynamics directly. Here, three important sets of evidence that could not be explained by the thermal effect are as follows: (i) the increase in spin precession frequency with the pump intensity, (ii) π-phase inversion of the precession, and (iii) the pronounced background distinguished from incoherent phonons and magnons. The model calculation of the Landau-Lifshitz-Gilbert (LLG) equation incorporating the time-dependent QSS strongly supports that all features mentioned are elucidated by one concept of QSS. The scenario is consolidated with full consistency using other ferromagnets (Co, Ni, and Ni x Fe 1-x ) to adjust the competition between strain and thermal effects.

Experimental configurations of MOKE and USI measurements.
The measurement geometries are shown in Fig. 1a. The bilayer structure of Al 2 O 3 (5 nm)/Co(25 nm)/Al 2 O 3 (substrate) and the trilayer structure of Co(200 nm)/Al 2 O 3 (15 nm)/Co(25 nm)/ Al 2 O 3 (substrate) were prepared by magnetron sputtering (nm in thickness are omitted hereafter). The Al 2 O 3 (5) layer in the bilayer structure was deposited to prevent oxidation. In the trilayer structure, the Al 2 O 3 (15) was used both to suppress the possible propagation of thermal magnons from the top Co layer and to match acoustic impedance with Co. The acoustic impedance Z is given by a product of mass density χ and sound velocity υ. The acoustic reflection coefficient r ac = (Z Co − Z Al2O3 )/(Z Co + Z Al2O3 ) between Co and Al 2 O 3 with χ Co = 8.9 × 10 3 kg m −3 , υ Co~5 .8 km s −1 27 and χ Al2O3 = 3.98 × 10 3 kg m −3 , υ Al2O3~1 1.2 km s −1 28 is~0.07 showing the good impedance matching.
For MOKE measurements, we used Ti:sapphire regenerative amplified laser pulses with a repetition rate of 10 kHz, a temporal width of 40 fs, and a center wavelength of 800 nm (The details of the MOKE setup are given in Methods). The frequency-doubled pump pulses (λ pu = 400 nm) excite the front sides of samples and probe pulses (λ pr = 800 nm) measure the differential Kerr rotation Δθ(t) and reflectivity ΔR(t) either of the front side for the bilayer or the back side for the trilayer, as shown in Fig. 1a-I, a-II, respectively. The strain pulse generated in the Co(200) surface in the trilayer propagates and triggers the spin precession in Co(25) layer through magnetoelastic coupling. For USI measurement, the front pump-front probe geometry we first tried partially degrades the quality of probe polarization. Therefore, in the current experiment, we make the pump pulses excite the back side of Co, and the probe pulses measure the front side ( Fig. 1a-III). In Fig. 1b, we present the simulated strain profiles η(z, t) at t = 7 (blue) and 20 ps (green) for the top Co layer in the trilayer structure for a better understanding of the time evolution of η p (z) propagating into the film depth and η qss (z) biding at a penetration length scale. For simplicity, we used η qss (t) = η(0, t) hereafter.
The main purpose of USI 29 is to measure unambiguously the lattice displacement dynamics u(z = 0, t) and extract η qss ðtÞ ¼ ∂uðz; tÞ=∂zj z¼0 in the end. The USI was built with a Ti:sapphire oscillator system operating at a repetition rate of 80 MHz, a temporal width of 40 fs, and a center wavelength of 800 nm. To avoid possible heat accumulation by a prepulse, the repetition rate is lowered down to 1 MHz by an electro-optic modulator (see "Methods" and Supplementary Note 1 for more details of the USI instrument). We confirmed that all data measured here have no heat accumulation from a prepulse by assuring zero background offset. Although ΔR(t) resulted from a change in a complex refractive indexñ is linked to η(z, t) through the relation Δñ ¼ ð∂ñ=∂TÞΔT þ ð∂ñ=∂ηÞΔη, it is relatively difficult to extract η(z, t) directly because of an unknown piezooptic property (∂ñ=∂η) of materials except for Ni, Cr, and Au at specified wavelengths 30,31 . Instead, as a circumvention, η qss (t) was obtained without introducing unknown coefficients by solving the one-dimensional wave equation incorporating the lattice temperature profile used as a driving source (see the strain calculation and parameter values in "Methods" for more details). The lattice temperature profile was calculated using a three temperatures model with proper parameter values of a bulk Co. In addition, the pump intensity calibration between the two different instruments (MOKE and USI) was performed by matching their signal levels of ΔR(t)/R. Figure 1c shows the real part of the Sagnac signal ρ = ΔR/2 R (Fig. 1c-I) and the converted lattice displacement u(0, t) = λ pr δϕ/4π (Fig. 1c-II) from the imaginary part (δϕ) of δr/r~ρ + iδϕ for the bilayer structure (blue circles for the measurement and the yellow line for the reproduction from the calculation). Here, r is defined as the complex amplitudereflection coefficient of the sample, ρ the relative change in reflectance, and δϕ the phase change induced by the pump pulse. In principle, η qss (t) can be obtained from the relation ηðz; tÞ ¼ ∂uðz; tÞ=∂z. However, as u(0, t) measured from USI is only one particular case of u(z, t), it is not possible to directly obtain η qss ðtÞ ¼ ∂uðz; tÞ=∂zj z¼0 . To bypass this problem, we first obtain a trial solution u(z, t) by calculating the wave equation combined with the three temperature models (see "Methods"). Then, by iteratively adjusting the pump intensity parameter and unknown heat coupling coefficients g ij , we find the final solution u(z, t) that best matches u(0, t) to the experimental curve measured by USI. Now that u(z, t) is known based on the measurement quantity, η qss (t) can be reproduced as shown in Fig. 1d. The positive value of η qss stands for a tensile strain resulted from the lattice expansion by a sign convention.
In contrast to the reproduced yellow curve in Fig. 1c-II, the experimental data u(0, t) shows a sharp positive peak at the first ps, indicating lattice contraction into the film direction. This may be explained by several possible candidates, such as the electronic stress by photo-excitation 32 and the magnetic stress driven by ultrafast demagnetization in FePt nanogranular structures [33][34][35] . However, the continuous film case in ref. 35 showed that the compressive strain was not produced due to the in-plane constraints upon the lattice expansion. This is a different result from our experimental data u(0, t) shown in Fig. 1c-II. However, identifying the origin is beyond the current scope of this paper, and the lattice contraction has not been considered in the reproduction process.
Variation of spin precession frequency as a function of pump intensity. The solid circles in Fig. 2a show Δθ(t) of the Co bilayer under the magnetic field H ext = 0.4 T and φ = 60°(φ: the magnetic field angle from out-of-plane direction) with various pump intensities I p . The experimental curves for other field angles (φ = 15, 30, and 45°) are presented in Supplementary Note 2. The precession frequency f was obtained by fitting the experimental curves with the following damped sinusoidal function: where τ 1 is the spin relaxation time, τ 2 the decay time of the background signal, δ 0 the initial phase of the precession, A and B the amplitudes of the first two terms on the right-hand side, and C defined as the background offset. The solid line in Fig. 2a denotes the fit curve for the highest I p . Figure 2b shows that the change in the precession frequency (Δf/f m = (f − f m )/f m ) increases with I p for all magnetic field angles (φ = 15°: black, 30°: red, 45°: green, 60°: blue circles. The error bars are marked on the data points). Here, f m is defined as the precession frequency for the minimum I p and corresponds to 13.8, 20.0, 24.4, and 27.6 GHz for the respective angles of φ (15, 30, 45, and 60°). This frequency variation observed for the polycrystalline Co film has an opposite trend to what would be expected from the conventional magnetic energy terms such as the magneto-crystalline, dipolar, and Zeeman energy. The magnitude of the effective magnetic field H eff (t) in such case would decrease with the increasing temperature of the sub-systems or the increasing I p .
In order to figure out the effect of η qss (t) on the precession frequency f, we solved the following LLG equation.
where M s is the saturation magnetization, H eff the effective magnetic field vector defined as H eff ¼ À∂F tot =μ 0 ∂M, γ s = 1.63 × 10 11 rad s −1 T −1 the gyromagnetic ratio 36   ultrafast Sagnac interferometry (USI) measurements. For MOKE, the front-pump and front-probe for the bilayer (a-I), the front-pump and back-probe schemes for the trilayer (a-II) were used. The back-pump and front-probe for the bilayer (a-III) were used for the Sagnac measurements. b Dynamic strain pulse profiles η(z, t) calculated for the top Co layer in the trilayer structure at t = 7 (blue curve) and 20 ps (green curve): the propagating bipolar pulse η p (z) and the quasi-static strain η qss (z). c Sagnac interferometric curves for Co bilayer: ρ = ΔR/2 R (c-I: ρ-relative change in reflectance induced by the pump pulse and ΔR-differential reflectivity) and the lattice displacement u(0, t) (c-II: experimental data (blue circles) and reproduced one by calculation (yellow line)). d η qss (t) profile extracted from the reproduced data u(0, t).
For the polycrystalline Co film, we used following parameter values: uniaxial magnetic anisotropy coefficient K u~0 , the magnetostriction value λ s = −62 μ 37 the mechanical stress σ s = 3(1 + ν)(1 − ν) −1 Bη qss (t), Poisson's ratio ν = 0.31, the bulk modulus B = 190 GPa 27 , and the demagnetizing factor N x = N y = 0, N z = 1. ϑ is defined as the angle between the magnetization M and σ s (out-of-plane) direction. From the time-dependent profiles of K u (t), M(t), η qss (t), and ϑ(t) calculated after solving both the three temperatures model and the wave equation based on the measurement curve u(t), the magnetization dynamics dM=dt is obtained (see "Methods" for calculation details about the three temperatures model and the wave equation).
After calculating the spin precession and extracting f from dM=dt, the simulation results are summarized in Fig. 2c. The values of Δf/f m considering η qss (t) (solid lines, on the left axis) are in good agreement with experimental data of Fig. 2b, showing an increase in Δf/f m up to 2% for a wide range of φ. In contrast, for the absence of η qss (t), Δf/f m decreased by 1.6% (dashed lines, on the right axis), as expected with conventional energy terms. We also checked the case when the spin precession was involved in neither thermal nor QSS effects. This requirement was met by the back-side measurement of the trilayer structure. The strain pulse η p (z, t) generated from the Co top layer propagates to the Co underlayer without carrying the thermal energy and contributes only to triggering spin precession through magnetoelastic coupling. From Δθ(t) with a function of I p under φ = 60°in Fig. 2d (the fit curve with Eq. (1) for the highest I p is denoted by the solid line), Δf/f m did not show a noticeable change within the experimental uncertainty (see the pink curve in Fig. 2b). Here I p = 2.6 mJ cm −2 corresponds to~70% of the burning threshold for the Co top layer.
Ferromagnetic resonance frequency dependence on T, φ, and η qss . To obtain the comprehensive picture of the QSS effect on f over a wide range of temperature T and field angle φ, the ferromagnetic resonance frequency f r was calculated using the classical Kittel equation with the magnetic free energy F, as used in the dynamics case.
where F pq ¼ ∂ 2 F=∂p∂q. The symbols ς and ξ denote the polar and azimuthal angles in a spherical coordinate system, respectively. ς eq is defined as the equilibrium angle of magnetization, and ξ = π/2 is set due to the azimuthal symmetry of the sample. As the classical Kittel equation generally treats the static regime, the QSS is assumed to be a constant supposing no thermal relaxation after a thermal equilibrium among the sub-systems. Hence, the QSS can be handled simply with thermal strain η th (T) = βΔT (thermal expansion coefficient of Co 27 : β = 13.7 μ). Figure 3a, b present the contour map of f r with the control parameters of T (x-axis) and φ (y-axis) under H ext = 0.4 T for the absence and presence of η th (T), respectively. Figure 3a shows that f r decreases monotonically with increasing T at a fixed value of φ (along the yellow line). In contrast, Fig. 3b indicates that f r increases gently with T up to~800 K (blue box region). This is because η th (T) increases almost linearly with T, unlike the barely altered M(T) owing to its high Curie temperature T c . Hence, the effect of magnetoelasticity plays a dominant role over those of the conventional magnetic energy terms. At a higher T (orange box region), the rapid drop of M(T) reduces the effective field strength, leading to the inflection point of f r around 800 K. This is not a special case only for Co but can be generalized to a broad class of magnetoelastic materials. We can expect that for low T c materials, a faster decrease of M(T) shows up at a lower T or I p , where the QSS effect on f r is comparatively weaker. The Ni can be a good tester to examine the competition of the thermal and strain effects. The Ni has low values of T c = 630 K and M = 525 emu cm −3 which are the key parameters for enhancing the thermal effect, but similar λ s , β, and B values to Co which are related to the strain effect. Figure 3c shows the contour map of f r including the QSS effect for Ni, and Fig. 3d presents Δf/f m for φ = 25, 45, and 60°selected from the contour map. In contrast to the Co case, f r decreases first at a low T, meaning that the thermal effect is dominant, and is switched to the increase with increasing T. To prove those results experimentally, Δθ(t) of Ni(270)/SiO 2 at φ = 25°and I p ≤ 2.2 mJ cm −2 was measured and plotted in Fig. 3e. As presented in the inset, after fitting the data with Eq. (1) (denoted as the solid line for the highest I p ), Δf/f m becomes first negative and then positive as I p increases (the error bars are marked on the data points). These curves reproduce the calculation (Fig. 3d) in excellent agreement supporting our analysis on thermal and QSS effects.
Governing role of QSS on ultrafast spin dynamics. We address that the QSS has a decisive role in determining the initial status of spin dynamics judged by further features that have been underrated. The first point is the phase of the spin precession. Figure 4a presents the calculated curves of ΔM z (t)/M z for the absence (red) and presence (blue curve) of η qss (z, t) under H ext = 0.54 T and φ = 60°. Here the spin precessions in Co and Ni (Figs. 2a and 3e) correspond to the blue curve, leading to a π-phase inversion compared to conventional free energy analysis.
The pictorial description of the inset explains the π-phase difference concisely. According to conventional energy analysis, the sudden decrease in M(t) after photo-excitation results in a rotation of H eff (t) (yellow solid arrow) to the out-of-plane direction (H′ eff (t)-red dashed arrow) due mainly to the decrease in dipolar energy (for polycrystalline or isotropic materials) and starts precessing around the new axis H′ eff (t). This aspect takes place when the thermal effect dominates the QSS effect for such conditions of λ s η qss ≥ 0 (λ s ≥ 0, β ≥ 0 or λ s ≤ 0, β ≤ 0) and even Unlike the Co bilayer, the sign of Δf/f m changes, as shown in the inset, confirming the competition between the thermal and the quasi-static strain (QSS) effects. All data points in the inset are presented with error bars obtained from the least-square fit. λ s η qss ≤ 0 (λ s ≤ 0, β ≥ 0) provided the magnitude of λ s is small. On the other hand, for λ s η qss ≪ 0, the QSS effect prevails over the thermal effect and drives the rotation of H eff (t) towards the inplane direction (H″ eff (t)-blue dashed arrow). Then, the spins start the precession around the H″ eff (t) axis. This orientation is determined at the initial stage of the precession, immediately after demagnetization. Ni x Fe 1-x (180)/Al 2 O 3 films were tested experimentally to verify this scenario, as plotted in Fig. 4b. As the alloys of x~0.36 and 0.8 have negligible β and λ s values, respectively, it is expected that they have low magnetoelastic energy (∝ λ s βΔT). These cases rotate H eff (t) towards H′ eff (t), inducing negative values of signals at the first half period (20-30 ps) as the conventional analysis does, as shown with the red graph in Fig. 4a.
The second feature is the pronounced decaying background (yellow dashed line in Fig. 4a). This in general has been considered as incoherent thermal phonons and magnons with nonzero wavenumbers 38 , which have the heat dissipation timescale. Rather, our analysis suggests that this needs to be interpreted as the contribution of coherent rotation of H eff (t) because of the QSS effect. This was clearly verified from the fact that Co and Ni with high magnetoelastic energy have much larger decaying offsets than those of NiFe alloys. These all features mentioned above, including the frequency increase in the spin precession, are reproduced through the model calculation incorporating only one concept of the QSS effect with the full consistency with the experimental data.

Discussion
By decoupling the thermal and the QSS effects using USI and the MOKE, we proved that the QSS has a governing role over the thermal effect on the overall behavior of ultrafast spin dynamics over a wide range of time scales. The ultrafast photo-excitation can strengthen the effective magnetic field due to the QSS effect. This leads to a higher frequency of spin precession and accounts for the π-phase inversion determined from the initial stage of the precession. Besides, the long-lived offset with the relaxation timescale of the QSS was interpreted as extra gain induced by the contribution of a coherent rotation of H eff (t). This study shows that the QSS is involved universally in a wide family of magnetoelastic materials and should be treated fundamentally.
The strain of 0.1-1% generated by photo-excitation is high enough to induce the modification of the electronic band structure 26 even at ultrafast time scales 39 . Therefore, it can be predicted that the QSS modifies the dielectric tensors bringing about features such as derivative-like changes in differential magneto-optics. And this will lead to inequivalence between its real and imaginary parts of magneto-optic signal, and involuntary extra gains in differential reflectivity as well even after thermal equilibrium timescale. To date, considerable effort has been made to identify genuine magnetism and the origin of demagnetization 3,20,[40][41][42][43][44][45][46] . Our demonstration that the QSS governs the spin dynamics provides an interpretation of ultrafast magneto-optics and therefore helps to better understand the physical mechanisms behind ultrafast phenomena by considering both strain and thermal effects.

Methods
Time-resolved pump-probe MOKE. We used Ti:sapphire regenerative amplified laser pulses with a repetition rate of 10 kHz, a temporal width of 40 fs, and a center wavelength of 800 nm. The pump pulses are frequency-doubled with a beta Barium Borate crystal and have a temporal width of 60 fs at a sample position. The pump and probe pulses were focused on the sample with diameter sizes of 150 and 30 μm, respectively, and their intensity ratio was set to~1000:1. The reflected probe pulses from the samples were split into orthogonal polarizations using a Wollaston prism and analyzed into the differential Kerr rotation Δθ(t) and differential reflectivity ΔR(t). An external magnetic field of H ext = 0.4 T was applied with various angles φ from the normal to the sample plane.
Ultrafast Sagnac interferometer. An ultrafast Sagnac interferometer based on a Ti:sapphire oscillator system issued with a center wavelength of 800 nm and a temporal width of 40 fs was employed to obtain a quantitative profile of the timeresolved lattice thermal expansion after the pump excitation. After passing through an electro-optic modulator operating at 1 MHz, the dispersed pulse was compressed to 45 fs by pairs of negative-chirped mirrors. The pump pulses were frequency-doubled by a 500-μm-thick BiB 3 O 6 crystal with a conversion efficiency of 25% and had a temporal width of 65 fs at a sample position. The probe and pump pulses were focused on the front and back sides of the samples in normal incidence. The numerical apertures of the objective lenses used in each beam line were 0.55 and 0.4, respectively. The temporal delay between p-pol. (arrives at a negative delay) and s-pol. probes (arrives at a positive delay) at the sample position were fixed to~1 ns and the s-pol. probes measure the pump-induced optical responses (ρ and δϕ) (Supplementary Note 1 for more details of the setup schematics and characteristics).
Strain calculation and parameter values. To obtain η(z, t), it is necessary to solve the three temperatures model and one-dimensional wave equation in the Co(25)/ Al 2 O 3 structure. The three temperatures model (Eq. (5)) is described as follows: where i, j = e, l, s stand for electrons, lattice, and spins, respectively. The C i and κ are the heat capacity per unit volume of bath i and the thermal conductivity. The g ij is the coupling coefficient between two baths i and j, and P(z, t) the laser source term. The ΔT ij is defined as T i − T j , the temperature difference between two baths i and j. Using boundary conditions-continuities of both the heat transfer at air/Co ð0 ¼ κ Co ∂T e ∂z j z¼0;Co Þ, Co/Al 2 O 3 ðκ Co ∂T e ∂z j z¼d;Co ¼ κ Al2O3 ∂T l ∂z j z¼d;Al2O3 Þ and the temperatures at Co/Al 2 O 3 ðT e ðz; tÞj z¼d;Co ¼ T Al2O3 ðz; tÞj z¼d;Al2O3 Þ-temperature profiles of the three baths were obtained (d is a film thickness). Here, since there is an arbitrariness for the selection of pump intensity and unknown value of C s , we carefully chose the value by matching with the experimental data of both u(z = 0, t) of the Sagnac measurement and the demagnetization Δθ L (t)/θ L (longitudinal Kerr geometry) under H ext = 0.5 T and φ = 90°. That is, T l (z, t) obtained from three temperature models calculation (Eq. (5)) is linked to u(z, t) through the onedimensional wave equation (Eq. (6)) as follows: where ρ is the mass density, υ the sound velocity, β the linear thermal expansion coefficient, and B the bulk modulus. Using boundary conditions, continuities of both the stress at air/Co ð3 1Àν 1þν B ∂u ∂z ¼ 3βBΔT l Þ, Co/Al 2 O 3 ð3 1Àν 1þν B ∂u ∂z À 3βBΔT l ¼ 3 1Àν Al2O3 1þν Al2O3 B Al2O3 ∂u ∂z Þ, and the displacements at Co/Al 2 O 3 ðuðz; tÞj z¼d;Co ¼ uðz; tÞj z¼d;Al2O3 Þ − u(z, t) were solved. Then, η qss ðtÞ ¼ ∂uðz; tÞ=∂zj z¼0 was determined by matching u(0, t) with the Sagnac measurement curve, as shown in Fig. 1c. The parameter values used in the simulation are as follows: C e = 6.6 × 10 2 T e , C l = 3.5 × 10 6 47 , C s = −W m dM 2 /dT s J m −3 K −1 48 where a molecular field prefactor W m was assumed to be 4.5 × 10 8 for Co, and C l = 3.1 × 10 6 for Al 2 O 3 49 . κ Co = 64 47 , κ Al2O3 = 20 W m −1 K −1 50 at 500 K, B = 190 GPa, ν = 0.31 for Co 27 . For dynamic heat coupling coefficients between two sub-systems, we used the following values: g el = 1.3 × 10 18 , g ls = 5.0 × 10 17 , and g es = 3.5 × 10 17 W m −3 K −1 to match with the Δθ L (t)/θ L and u(0, t) curves qualitatively. For the Curie-Weiss curve, M(T) = M s (1 − 1.058(T/T c ) α ) ζ was extracted by fitting the data in ref. 51 , here M s = 1360 emu cm −3 , T c = 1380 K, α = 3.15, and ζ = 0.50.

Data availability
The data that support the findings of this study are available from the authors upon reasonable request.