Ultrafast dynamics of exchange stiffness in Co/Pt multilayer

The exchange stiffness coefficient, Aex, represents the strength of direct exchange interactions among neighboring spins. Aex is linked to most of the magnetic properties such as skyrmion formation, magnetic vortex, magnetic domain wall width, and exchange length. Hence, the quantification of Aex is essential to understanding fundamental magnetic properties, but little is known for the dynamics of Aex on a sub-picosecond timescale. We report the ultrafast dynamcis of Aex in an ordered magnetic state in Co/Pt ferromagnetic multilayer. Time-resolved magneto-optical Kerr effect and reflectivity measurements were analyzed for various pump fluences. We reveal that the significant dynamical reduction of Aex is responsible for the dramatic increase of remagnetization time for high fluences. The analysis shows that Aex dynamically varies, strongly affecting overall ultrafast demagnetization/remagnetization process. The investigation demonstrates the possibility of Aex engineering in femtosecond timescale and thereby provides a way to design ultrafast spintronic devices. The exchange stiffness coefficient describes the interaction between neighbouring spins of a ferromagnetic system and analysing these interactions at ultrafast timescales can reveal information on their demagnetisation and remagnetisation mechanisms. Here, the authors investigate the ultrafast dynamics of the exchange stiffness coefficient of a ferromagnetic multilayer system and model the data using a generalised three temperature model.

T he exchange interaction is a fundamental aspect of ferromagnetism [1][2][3] ; this interaction, underpinning the existence of ordered magnetic states, allows the alignment of neighboring spins in a system. The origin of the direct exchange interaction is ascribed to the overlapping of electron wave functions among neighboring atoms with no classical analogy, and its strength depends sensitively on the atomic and lattice structure 4,5 , while there also exist other exchange mechanisms such as Dzyaloshinskii-Moriya interaction 6,7 , and superexchange interaction 8 . The strength of the direct exchange interaction is represented by the exchange-stiffness coefficient A ex , which is a material parameter mostly depending on the atomic and lattice structure and detailed characteristics of which are determined by the electronic band structure [9][10][11] . Fundamental mechanisms of novel spin phenomena, such as magnetic vortices 12,13 and skyrmions 14,15 , might be understood based on a quantitative analysis of A ex . It has been known that A ex is temperature-dependent 16,17 . As the temperature increases, thermal agitation reduces the degree of the ordering of neighboring spins, effectively lowering the value of A ex .
Although the temperature-dependence of A ex is relatively well recognized in static cases 17 , very little is known regarding how A ex varies on an in the ultrafast timescale. In case of ultrafast photoinduced demagnetization, as the fluence F P of a pump laser increases, the disorder in a spin system should increase so that the effective spin temperature should also increase. It is then expected that the remagnetization time of the system from an excited disordered state to a stable equilibrium state should increase as F P increases. Recent reports of the dynamics of the exchange interaction on ultrafast timescale have shown that a fundamental exchange interaction varies on scales of several tens of femtoseconds in ferromagnetic NiFe alloy 18 and antiferromagnetic KNiF 3 19 . A possibility of ultrafast control of exchange interaction by using a femtosecond pump laser has been proposed theoretically [20][21][22] . All of these results indicate that the exchange interaction or exchange splitting dramatically changes by a femtosecond pump laser. However, to our best knowledge, no systematic study has been conducted to quantify the dynamics of A ex on ultrafast timescale.
Here, we report, the dynamics of A ex on a femtosecond timescale in Co/Pt multilayers for a range of F P , demonstrating that A ex varies rapidly, affecting spin dynamics and its variation can be controlled by the pump fluence. While the electron-spin interaction strength is kept constant all the time in the conventional three-temperature model (3TM) in the study of ultrafast magnetism, the dynamic change of A ex is considered, adopting the generalized three-temperature model (G-3TM) developed by A. Manchon et al 23 .

Results
Fluence dependent remagnetization time. We performed timeresolved magneto-optical Kerr effect (TR-MOKE) measurements for [Co (6.2 Å)/Pt (7.7 Å)] 5 multilayer film, of which the magnetic properties such as perpendicular magnetic anisotropy and saturation magnetization are well known 7,24-27 . A detailed experimental configuration was reported elsewhere 28 . TR-MOKE signals were measured for 1.7 ≤ F P ≤ 28.5 mJ cm −2 for time delays of up to 30 ps. To exclude the dichroic bleaching effect, the experiments are carried out by pump beam at both 400 and 800 nm wavelength (λ pump ). TR-MOKE vs. time at different fluences is plotted in Fig. 1a for the case of λ pump = 800 nm, while no significant difference is observed in the overall trend for the case of λ pump = 400 nm. The signals were normalized by their maximum changes to compare the dynamical behaviors in remagnetization for different conditions. An external magnetic field of 1.7 kOe was applied normally to the film surface. All the measured curves exhibit clearly the dynamics of photo-induced demagnetization and subsequent remagnetization. The maximal change of demagnetization is observed around t = 0.3 ps for all the fluence cases. In case of λ pump = 800 nm, as F P was increased from 1.7 to 28.5 mJ cm −2 , the remagnetization was slowed down for F P > 9.9 mJ cm −2 . The remagnetization behavior at 1.7 ≤ F P ≤ 16.5 mJ cm −2 well-fitted with a single exponential curve (Fig. 1a, dotted lines), yielding the characteristic time τ R of remagnetization (open square in Fig. 1b). We note that τ R increased drastically as F P increased by only a factor of a few. Fitting with a single exponential curve was not valid for F P > 16.5 mJ cm −2 . A similar trend is observed for the case of λ pump = 400 nm, where τ R is fitted with a single exponential curve as well.
G-3TM analysis. For further understanding, we conducted G-3TM analysis for the TR-MOKE data 23,28,29 . G-3TM is composed of three coupled equations (see Supplementary Note 1 for more details): where T e , T l , and T s are the electron, lattice, and spin temperatures, respectively. C e , C l, and C s are the specific heats of the electron, lattice, and spin, respectively. G el , G es , and G ls are the electron-lattice, electron-spin, and lattice-spin interaction channels, respectively. P(t) is a laser source term with a Gaussian temporal profile. The term that contains K l represents lattice thermal diffusion, which is modeled to be proportional to the third power of the temperature increase of lattice system 30 . A typical relation between the magnetization and the spin temperature: ∝ (1 -(T s / T C )) 0.5 , where T C is Curie temperature of 1131 K 17 , to match the normalized TR-MOKE signal (|Δθ Kerr |/|Δθ peak |). In conventional 3TM, G el , G es , and G ls have been set to be constant over time. However, in our study, for correct analysis, the electron-spin interaction channel(G es ) was allowed to change over time. Adopting the G-3TM, G es can be written as where dt; and ð3Þ where a is a lattice constant, A ex0 an exchange-stiffness coefficient at 0 K, V unit cell volume, T F Fermi temperature, S = 3/2 spin quantum number, M magnetization and G 2 a function based on the second-order Debye function 23 (Eq. (3)). T F is Fermi temperature, chosen to be that of fcc Co (T F = 16.87 Ry/k B ) 31 . (4)) is a temperature-independent electron-spin interaction channel. In static case, it is well known that A ex / (J ex a −1 ) <S 2 > , where J ex is an exchange interaction constant, a is a lattice constant. The proportionality depends on material parameters such as a periodic lattice configuration. Since A ex can be easily measured rather than J ex , we focus on qauntifying A ex on an ultrafast timescale. When the relation between A ex and J ex is extended, we have put For simplicity, we have used the approximation, J ex~2 aA ex0 16 . In the G-3TM, G el is still assumed to be constant because the relaxation rate between the electron and the lattice is expected to be simply proportional to the temperature difference. G ls was also set to be constant throughout the simulations. The G-3TM is composed of several free parameters, so the fitting should be processed with care. First, time-resolved reflectivity R(t) data were utilized to estimate values for C e , C l , and G el , considering only the electron and lattice, based on the 2temperature model 32,33 . In the full analysis using the G-3TM, the reflectivity and MOKE data were fitted. As a constraint in the analysis, the measured values for the degree of demagnetization (D demag ) were used (Fig. 2a, d). Hysteresis loop measurement is the best way to estimate the D demag . The hysteresis loops were measured at t = −2 ps and 0.3 ps (the maximal demagnetization) using the same TR-MOKE setup with probe-beam modulation for all the F P (Methods section). An example of measurements for F P = 13.2 mJ cm −2 (Fig. 2a, inset) shows that a D demag is 70%. The excellent match has been established in all the cases (See Supplementary Note 2, where the utilization of R(t) measurement and the D demag for fitting is described).
The G-3TM fitting determines temporal evolutions of spin, electron, and lattice temperature at wavelength of pump pulse λ pump = 800 nm (Fig. 2b, c) and λ pump = 400 nm (Fig. 2e, f). The cases for very high F P corresponding to T s being very close to Curie temperature (1131 K) are not considered, where the G-3TM may not be valid. The fitted value of C e was 1.8~2.1 × 10 3 J (m 3 K 2 ) −1 and, C l was 1.8~5.0 × 10 6 J (m 3 K) −1 34,35 . (all the fitting parameters are summarized in Supplementary Note 2). C s should depend on the spin temperature T s . In the original Manchon's paper 23 which the G-3TM on, C s is determined from the numerical derivative of the spin energy. In fitting our data, we have found that the fitting becomes quite good if C s is smaller than~10 4 J (m 3 K) −1 in all cases. Thus, we used a small value of C s = 100 J (m 3 K) −1 for all cases. The upper limit of fitted C s value (~10 4 J (m 3 K) −1 ) in the present work seems to be a little bit smaller than the reported values determined from 3TM. For instance, in Ref. 35 , C s of Ni and FeCuPt are 0.2 × 10 6 J (m 3 K) −1 and 0.17 × 10 6 J (m 3 K) −1 .
The maximum values of electron temperature T e max and spin temperature T s max at t = 0.3 ps increased as F P increased; e.g., T s max at t = 0.3 ps changes from 564 to 1040 K as F P increases from 1.7 to 9.9 mJ cm −2 (λ pump = 800 nm). High F P increases the amount of energy transferred to the subsystems, so the increase of T e max and T s max is expected. The equilibrium temperature at which T e = T l = T s also increased consistently as F P increased, but it is very interesting to note that the difference between T e max and T s max got larger substantially as F P increased (Fig. 2b, c, Fig. 2e, f). In the context of the G-3TM, this observation indicates that the interaction channel G es between the electron and spin subsystem is reduced, resulting in the increase of the thermal separation of the spin system from the electron subsystem as well as thereby the increase of the time required to reach thermal equilibrium. This phenomenon may be a reason for the increase of τ R as F P increases as observed in Fig. 1b. Figure 3a is the plot of G es (Eq. (4)) as a function of T e and T s . As T s increases, G es increases then decreases for a given electron temperature. The values T e , T s , and G es determined from fitting to our experimental data at λ pump = 800 nm are shown in a gray curved line in Fig. 3a, and again in Fig. 3b for various F P . The case of T s = T e is also presented as a dotted curve for guidance in Fig. 3b. The non-monotonic nature of G es with respect to T s is a direct consequence of Eq. (3). G es increased monotonically with an increase in T s at low F P = 1.7 and 3.3 mJ cm −2 , but the increasing-then-decreasing behavior is observed at high F P > 6.6 mJ cm −2 . We suspect that diverse experimental results of ultrafast demagnetization dynamics might be originated from this different trend of G es at high F P 36-38 .
The dynamical variation of G es on a femtosecond timescale is plotted for various F P in Fig. 3c. At low F P , the simple increaseand-decrease behavior of G es is observed, with a maximum at t = 0.3 ps. The time of the maximum G es coincides with the time at which the D demag is the greatest. At high F P , G es quickly reaches the first peak right after the arrival of a pump pulse, decreased to a minimum at around the time of maximal D demag (t = 0.3 ps), and then increased again. A comparison between the behaviors of G es , G ls , and G el at t = 0.3 ps under various F P (Fig. 3d) reveals that G el is the strongest channel, G es increased at low F P , but decreases at high F P ; this trend may be the result from the feature of G es (Fig. 3a, b). On the other hand, G ls is the weakest channel (as often neglected) but becomes comparable to G es as F P increases. G ls is involved with spin-orbit coupling 23 , which might get stronger as T e or D demag 39 .
The above discussion indicates that the dynamics of the photoinduced demagnetization and remagnetization in Co/Pt spin system is mostly governed by G es and G ls . Figure 3e shows a ratio of G es to G ls at t = 0.3 ps for various fluences. G es /G ls is larger than 10 for F P < 6.6 mJ cm −2 for λ pump = 800 and 400 nm. This imbalance implies that the spin-electron interaction is dominant in this F P regime. For F P ≥ 6.6 mJ cm −2 , G es /G ls approaches unity asymptotically, indicating that spin-lattice interaction becomes increasingly important. The G-3TM fitting yields the values for G es . Equations (2) and (4) allow us to calculate A ex0 , temperatureindependent exchange-stiffness coefficient. The estimated value of A ex0 turns out to be 10.01 pJ m −1 at all the F P ; this value agrees well with a reported value for a Co/Pt multilayer 40,41 . Other analysis methods 42,43 could also reproduce the slow rate of magnetization at high fluences. It should be commented that the G-3TM based on the Hamiltonian for laser-induced demagnetization 23 allows us to separately monitor time-dependent G es and G ls as well as their ratio G es /G ls . In this work, we note that the ratio particularly seems to play an important role in determining the energy-excessive spin dynamics on a sub-ps timescale.

Discussion
The previous studies 17, [44][45][46] in static cases have shown that the temperature dependence of A ex is expressed as power of M with a scaling exponent ranging from 1.79 to 1.82 in case of Co. We set the exponent to be 1.8 and write the temporal variance of A ex as Based on Eq. (5), time-dependent A ex is plotted in Fig. 4. The increase in F P results in the reduction of A ex , as generally expected in static cases. However, the recovery of reduced A ex depends sensitively on F P . A ex decreases asymptotically as F P increases, and saturates at F P ≈ 9.9 mJ cm −2 (λ pump = 800 nm) or 12.1 mJ cm −2 (λ pump = 400 nm) without further decrease with respect to the fluence higher than this value, which is expected from the saturated behavior of the D demag at high F P . At F P > 9.9 mJ cm −2 , A ex was~1 pJ m −1 . The maximal decrease of A ex occurred at t = 0.3 ps when the maximum D demag occurs in the TR-MOKE measurement. TR-MOKE data (Fig. 1a) show a similar trend to the trend in A ex (Fig. 4). The magnetization M and A ex recovered quickly at low F P whereas the recovery becomes significantly slow for high F P .
It should be noted that Eq. (5) is valid for the steady-state case and we use the very rough assumption that considering that even in the out-of-equilibrium case, there could be a rough relation between A ex and temperature-dependent M 9 . Indeed, although we use 3TM 23,28,29 , 3 temperatures are not fundamentally well defined in the out-of-equilibrium state and only phenomenologically defined once 3TM is used. On the other hand, we consider that the M[T] might not be totally different compared to the steady-state case, since the estimated spin temperature is still below T C . The pump pulse excites the electrons around the Fermi energy so that the excited electrons occupy the allowed energy levels above the Fermi energy, while remaining electrons still follow the Fermi-Dirac distribution. Moreover, the thermal equilibrium among 3 temperatures is achieved around~10 ps and thus, the M[T] will be soon replaced back to the equilibrium case after this timescale. Therefore, we think that there might be a Possible mechanisms of A ex reduction might be involved with Stoner exchange splitting reduction 9 , where it has been reported that dynamic exchange splitting is determined by time-dependent magnetization M(t). On the other hand, magnon generation should be also an important factor 47 , where it has been reported that the magnon contribution to demagnetization is dominant only on a very short (700 fs) timescale. Thus, in our case, we consider that the magnon contribution could exist on a sub-ps timescale, while the Stoner exchange splitting reduction is lasting longer up to few tens of ps since there is still a substantial amount of demagnetized M(t) in the present work. It should be also noted that G-3TM, which our whole analysis is based on, includes the magnon generation by hot electron as a key mechanism in the model. In G-3TM, electron-spin interaction Hamiltonian intrinsically deals with the effect of magnon generation, which might be reflected in the A ex dynamics, particularly on the sub-ps timescale. The effective Stoner exchange splitting reduction is understood based on the reduced M(t) over the whole process of demagnetization and remagnetization. It should be mentioned that our film is prepared on a Si substrate with no doping, which can be approximated to be an insulator so that the spin diffusion effect could be negligible in the process of demagnetization. Fig. 3 Temporal change of interaction channels [G es (electron-spin interactio channel), G el (electron-lattice interactio channel), and G ls (lattice-spin interactio channel)] and the exchange stiffness. a 3D map of G es vs T s (spin temperature) and T e (electron temperature). Dark gray line: trajectory of G es at 9.9 mJ cm −2 for -1 to 0.3 ps. b G es vs T s for -1 to 0.3 ps at 1.7 ≤ F P (pump fluence) ≤ 9.9 mJ cm −2 (wavelength of pump pulse λ pump = 800 nm). Black dashed line: T e = T s . c G es vs. the delay time (t) between -1 to 3 ps at 1.7 ≤ F P ≤ 9.9 mJ cm −2 (λ pump = 800 nm). Vertical dashed line: t = 0.3 ps. d Interaction channels (G es , G el , and G ls ) vs F P at t = 0.3 ps for the case of 800 and 400-nm λ pump . Error bars represent that lowest standard error region (<5%) during each fitting parameter by G-3TM. G es is not fitting paramert in G-3TM. e G es /G ls vs F P at t = 0.3 ps for the case of 800 and 400-nm λ pump . Error bars represent that lowest standard error region (<5%) during each fitting parameter by G-3TM. The above analysis reveals that the significant increase of the remagnetization time (τ R ) (Fig. 1b) for high fluence is directly related to the reduction of A ex . In order to confirm how much A ex or demagnetization state affects the remagnetization process, we have carried out another series of independent micromagnetic simulations. The micromagnetic simulation was performed using the Object-Oriented Micromagnetic Framework 48 based on the Landau-Lifshitz-Gilbert (LLG) equation: where the gyromagnetic ratio γ = 2.210 × 10 5 m (A s) −1 and H eff is the effective magnetic field. Since the micromagnetic simulation does not consider the temperature variation, it is not suitable to dynamics study but still provides valuable information about the material properties for remagnetization at a fixed temperature. We set the initial degree of magnetization according to the measurement and simulated how the remagnetization proceeds for different magnetic parameters such as A ex and magnetic anisotropy. In the simulations, an external magnetic field of 1.7 kOe was applied with an angle of 0°to the surface normal of the film as in the experiments. The saturation magnetization M s of the film was set as 10 3 kA m −1 . Magnetic anisotropy constant K was set as 6 × 10 5 J m −3 . Gilbert damping constant α was set as 0.05. The cell size was 0.5 × 0.5 × 0.5 nm 3 and the sample size was 50 × 50 × 7.5 nm 3 . The initial demagnetization state is set by the experimentally-measured D demag for various F P (Fig. 1b is replotted with respect to D demag corresponding to various F P (black open squares) as shown in Fig. 5a). As the D demag increases, τ R increases, drastically at larger D demag than 60 %. To determine parameters that are most responsible for the abrupt increase of τ R with the increase of F P (or high D demag ), we performed micromagnetic simulations for A ex = 1, 8, and 15 pJ m −1 . In each simulation, A ex was fixed throughout the simulation. Simulations with A ex = 8 (Fig. 5a, green triangle) and 15 pJ m −1 (Fig. 5a, red circle) agree well with experiments at low F P (or low D demag ). The literature value 40,41 of A ex of Co/Pt multilayer for a static case is1 0 pJ m −1 , which is consistent with the range of our simulation parameter. In the simulation with A ex = 1 pJ m −1 , τ R was substantially higher than the experimental observations (low D demag ).
Magnetic anisotropy is another important parameter that might affect τ R . The magnetic anisotropy is determined mostly by the crystal structure, sample shape, and multilayer interfaces. This anisotropy produces a perpendicular magnetic anisotropy in the Co/Pt multilayer used in the present study. The micromagnetic simulation of τ R for various magnetic anisotropy constants K for various D demag from 47 to 80 % showed no significant change of τ R at 0.01 ≤ K ≤ 1.2 MJ m −3 (Fig. 5b). The measured value of K for the Co/Pt multilayer in the present study is K = 0.6 MJ m −3 , which is within our simulation range. Thus, we infer that the variation in K is not responsible for the increase of τ R . This independence is expected because the variation of K will mostly affect the total effective field without directly modifying spin-spin or spin-electron interactions.
We also systematically changed A ex in micromagnetic simulations with the D demag being fixed at 80% (Fig. 5c). As A ex was varied from 4.0 to 0.1 pJ m −1 , τ R increased from 10.8 to 30.1 ps. In particular, at A ex < 1 pJ m −1 , τ R increases rapidly in a similar manner to the experimental observation (near high F P in Fig. 1b, The corresponding values of the pump fluence are written next to the experimental data points. Yellow region: F P > 9.9 mJ cm −2 where τ R is abruptly increased. b τ R vs. K (magnetic anisotropy) for D demag of 47 (circle) and 80 % (square). A ex is changed from 1 (black) to 15 (red) pJ m −1 . c Simulated result of τ R vs. A ex . The D demag is fixed to be 80 %. Yellow region corresponds to the yellow region in b. d Simulated result of τ R vs. α (Gilbert damping constant). All error bars represent standard error of the single exponential fitting.
or near high D demag in Fig. 5a). The increase in F P can be expected to reduce A ex significantly, resulting in a large increase of τ R . The decrease in A ex is generally expected to cause the increase in τ R , because the reduced spin-spin interaction weakens the ordering among neighboring spins. This micromagnetic simulation also confirms the analysis by G-3TM.
We have carried out a simulation with the variation of damping parameter (α), as seen in Fig. 5d. A ex is set to be 11 pJ m −1 , anisotropy constant K is 0.6 MJ m −3 , and D demag is set to be 80%. τ R is found to be insensitive if α is larger than 0.005, as seen in the figure. For α smaller than 0.005, τ R drastically increases due to the significant contribution of precessional oscillation. α of Co/ Pt multilayer is reported in the range of 0.02-0.1 49,50 , which is much larger than 0.005.
In summary, we have investigated the dynamical variation of A ex on an ultrafast timescale, by TR-MOKE and reflectivity measurements in a Co/Pt multilayer for various F P . Our phenomenological analysis suggests that the ultrafast remagnetization mechanisms may be governed by the dynamically changing A ex , which is also closely related to G es , and G ls also becomes nonnegligible in case of high F P . Our comprehensive micromagnetic simulations implies that significantly reduced A ex is responsible for the large remagnetization time. These results demonstrate the possibility of engineering magnetic properties on an ultrafast timescale by modifying G es , G ls , and G el .

Methods
MOKE measurement. TR-MOKE measurements with a pump-probe stroboscope were performed on a Co/Pt multilayer. We used the femtosecond laser pulses generated by a Ti:sapphire multipass amplifier operating at a repetition rate of 3 kHz with a center wavelength of 800 nm and a pulse duration of 25 fs. We employed two pump wavelengths of 800 nm and 400 nm obtained from BBO (BaB 2 O 4 ) crystal. As probe pulses, one wavelength of 800 nm was used. F P was varied from 1.7 to 28.5 mJ cm −2 and probe fluence was 0.3 mJ cm −2 . For TR-MOKE measurements, the pump beam was modulated using a mechanical chopper at 500 Hz. An external magnetic field of 1.7 kOe was applied throughout the measurements, with an angle of 0°to the surface normal of the film to keep the initial sample condition saturated before a subsequent pump pulse.
To estimate the D demag , we conducted a series of hysteresis measurements at times of t = −2 ps and 0.3 ps under the same TR-MOKE setup with only probebeam modulation at 500 Hz, while sweeping a magnetic field from −1.7 kOe to 1.7 kOe. The hysteresis measurement at t = −2 ps gives the magnetization of an intact sample.
Sample. [Co(6.2 Å)/Pt(7.7 Å)] 5 multilayer films were deposited on Si substrates by dc magnetron sputtering, then capped by a 22-Å Pt layer to prevent the oxidation of the surface. The structure of the Co/Pt multilayers with well-defined interfaces was confirmed by a low angle X-ray diffraction and extended X-ray absorption fine structure analysis. The film had a perpendicular magnetic anisotropy (K = 0.63 MJ m −3 ) and saturation magnetization (M s = 1.04 × 10 3 kA m −1 ), which are similar to literature values [24][25][26][27] .

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