Anisotropic ultrafast spin dynamics in epitaxial cobalt

Ultrafast quenching of magnetic order at sub-picosecond time scales triggered by femtosecond laser pulses and its subsequent recovery were observed in a ferromagnetic nickel thin film in the pioneering experiment by Beaurepaire et al.¹ Since then, many experiments have confirmed the occurrence of this phenomenon in metallic thin-film ferromagnets.^2–19 Significant theoretical progress^20–26 has been made toward finding the fundamental microscopic mechanisms that are able to explain how angular momentum is lost and recovered at these ultrafast time scales, orders of magnitude faster than expected, e.g., by the textbook Landau–Lifshitz theory.²⁷ Despite a two-decade-long quest, a complete understanding of the phenomenon is still lacking.

The role of the lattice in ultrafast magnetism has been discussed since the early years following the pioneering experiment. The Elliott–Yafet-type spin-flip scattering was put forward as a possible mechanism through which the angular momentum can be transferred from the spin system to the lattice, although the efficiency of this mechanism has been debated.^{22,23,28–30} Surprisingly, only very few experimental studies^31–35 have investigated epitaxial systems, where the crystalline structure of the sample can be properly modeled. Recently, using femtosecond x-ray diffraction, it was observed that a femtosecond optical pulse can trigger ultrafast coherent terahertz longitudinal acoustic phonons (up to 4 THz) in an epitaxial iron thin film,³⁴ disproving the common assumption that the lattice cannot respond coherently on ultrafast time scales. Even more recently, another ultrafast x-ray experiment on a similar iron film³⁵ suggested the possibility of the ultrafast version of the Einstein–de Haas experiment, where the demagnetization of the material is compensated by a coherent mechanical rotation of the body, in this case, driven by the generation of transverse acoustic phonons at terahertz frequencies. Notwithstanding, the unambiguous detection of the involvement of the lattice structure in ultrafast magnetization dynamics is still to be achieved.

In this Letter, we investigate a different model system, an epitaxial hcp(1$1¯$00) cobalt thin film, which we probe with a femtosecond optical pump–probe setup. Using the time-resolved magneto-optical Kerr effect (TR-MOKE), we measure the demagnetization and recovery of the sample magnetization on the femto- and picosecond time scales. The magnetization is set along the easy or hard magnetization axis, which corresponds to two orthonormal lattice directions in the thin-film plane. Surprisingly, there do not yet exist systematic studies on possible anisotropic ultrafast spin dynamics in epitaxial model systems with strong magneto-crystalline anisotropy, which would allow us to pin-down the role of the orientation-dependent electron–phonon coupling. In the following, we show that our measurements can reveal distinct magnetization dynamics coupled to the anisotropic lattice structure, even without the atomic resolution given by an x-ray probe, and that is consistently explained by ab initio calculations of the anisotropic electron–phonon interaction.

A 15-nm-thick epitaxial hcp(1$1¯$00)-cobalt thin film was grown as Co[1$1¯$00] on a MgO(110) substrate and a Cr(211) seed layer. The Co layer was capped with a 3-nm-thick Pt layer. The easy axis of magnetization is along the $c−$ axis [0001] and lies in the plane of the film, as shown in Fig. 1(a). The hard axis of magnetization [11$2¯$0] is perpendicular to it and also in the sample plane. This strong in-plane magnetic anisotropy of the film enabled us to measure the ultrafast demagnetization along two different crystalline orientations by a simple in-plane rotation of the sample, as shown schematically in Fig. 1(b). Magneto-optical loops along the easy (orange) and hard (blue) magnetization axes measured with the longitudinal MOKE are shown in Fig. 1(c). The loops are qualitatively similar to the vibrating sample magnetometer data presented in the supplementary material. The probing configuration and laser setup are the same as for the time-resolved data shown below. In this way, the time-resolved data can be directly normalized with respect to the magnetization loops.

The pump–probe experiments were performed with an amplified Ti:Sapphire laser, with a pulse duration of approximately 40 fs, a repetition rate of 1 kHz, and a central wavelength of 800 nm. As a pump, we used the 400-nm optical pulses generated by frequency doubling the fundamental of the laser, using a β-barium borate crystal. As a probe, we used the fundamental of the laser at 800 nm. The pump was incident at an angle $θpump≈10$° and the probe at an angle $θ2≈55$°, as shown in Fig. 1(a). In this configuration, we are optimized to measure the longitudinal MOKE, proportional to the in-plane component of the magnetization.³⁶ In addition, the different pump and probe energies allow us to suppress coherent optical artifacts due to the formation of transient gratings in the film.³⁷ The intensity autocorrelator trace measured at the sample position returned an 80 fs Gaussian, which is the actual resolution of the measurement. This is consistent with the convolution of two 55-fs long pulses which have acquired finite dispersion along the beam path.

Figures 2(a) and 2(b) show the TR-MOKE measurements performed at selected fluences for the easy and hard magnetization axes. The delay traces are calculated as the difference of the delay traces recorded using magnetic fields of equal magnitude but opposite sign and then normalized by the amplitude of the magnetization loops shown in Fig. 1(c). The applied magnetic fields were ±400 mT and ±1000 mT for the easy and hard magnetization axes, respectively, enough to reach saturation. These values were chosen to result in the same effective field, as demonstrated by the measurement of the sample ferromagnetic resonance (FMR), see the supplementary material. However, we have checked that none of the observations reported below depend on the magnitude of the external field, given that the sample is saturated. Plotting the difference of opposite fields allows for isolating the pure transient magnetic signal and removes the contribution from the transient reflectivity signal. The transient reflectivity is shown in the supplementary material and has a maximum relative change of 0.2% for the highest fluence. This value has to be compared with the maximum relative change in Kerr rotation of about 20%. The much larger variation indicates that the Kerr signal is measuring genuine magnetization dynamics and not optical artifacts, i.e., $ΔθK(t)/θK0=ΔM(t)/M0$⁠.³⁸ Furthermore, we checked the dynamical Kerr rotation and ellipticity of all combinations of s and p pump and probe polarizations,^31,39,40 and the MOKE response stayed the same in shape and amplitude, within the experimental uncertainty.

The demagnetization curves shown in Fig. 2 illustrate the typical response observed in this experiment: a rapid quench of the magnetization on a timescale of the order of $∼100$ fs followed by a fast recovery on the timescale of $∼1$ ps, and finally a much slower recovery on tens of ps time scales. The figure also shows that the maximum demagnetization increases monotonically with the absorbed fluence. In order to accurately determine the demagnetization and recovery time constants, we fitted the ultrafast demagnetization data with the equation given in Ref. 41. In those equations, the dynamics is described by the decay time τ_m for the ultrafast demagnetization, $τR1$ for the fast recovery, and $τR2$ for the slow recovery, which we can extract by a careful fitting procedure described in the supplementary material.

For the three fluences presented in Fig. 2, there is only a negligible difference (i.e., within the error bars) in the demagnetization amplitude at each fluence for the two different magnetization orientations. The change of demagnetization time constant τ_m for these fluences is below the resolution of our measurement,⁹ and we obtained the best fit with $τm=120$ fs for all these measurements. By using this value for τ_m, we could reliably extract the fast recovery time constant $τR1$ and the slow recovery time constant $τR2$ for all the measurements. We also note that the slow recovery time $τR2$ is a coarse approximation of the dynamics and that excludes the full response of the magnetization, including the ferromagnetic resonance. Hence, we do not discuss it further in the following.

We, instead, focus on the fast recovery time $τR1$⁠. Figure 3 shows the extracted $τR1$ for both orientations as a function of the absorbed fluence. The fast recovery time increases with increasing fluence, consistent with previous reports.^9,42,43 In addition, we also observe that the fast recovery of the magnetization along the hard axis orientation is always faster than for the easy axis orientation. The data can be fitted assuming a linear dependence and forcing the fit to go through the origin. The slopes of the lines are approximately 378 fs cm²/mJ for the easy magnetization axis and 270 fs cm²/mJ for the hard magnetization axis. Therefore, this is one of the main experimental findings of this work; the ultrafast magnetization dynamics in hcp Co recovers systematically faster along the hard magnetization axis than along the easy axis. In the absorbed fluence range up to 2 mJ/cm², the ratio of the easy/hard slopes is approximately 1.40 ± 0.02, i.e., the recovery is approximately 40% slower along the easy magnetization axis. One could argue that if this is due to an intrinsic non-equilibrium spin-scattering mechanism within the material, one would expect not only the magnetization recovery to be faster but also the quenching. However, the expected change of the demagnetization time constant at these fluences is comparable with the resolution of our measurement, and such difference may not be measurable within our experimental resolution.

In order to test this hypothesis, we looked at slower demagnetization by increasing the maximum absorbed fluence by a factor of two, i.e., 4 mJ/cm², close to the sample damage threshold, observed as a permanent change in the sample reflectivity and magneto-optical signal. The data recorded at such fluence are shown in Fig. 4 for the two magnetization axes, where we also show the corresponding transient reflectivity. With the demagnetization process slowed down, we can now resolve the different time constants of the quenching, with $τm,hard≈120$ fs and $τm,easy≈165$ fs, and show that indeed the hard magnetization axis has an overall faster dynamics than the easy axis. We also note that in this measurement, $τm,easy/τm,hard≈1.37$ and $τR1,easy/τR1,hard≈1.32±0.02$⁠, suggesting that the same microscopic mechanism may be involved in governing the quenching and relaxation processes. This also shows that the dynamics along the easy magnetization axis is substantially slower than the one along the hard axis even at much larger fluences and where this trend can be observed directly from the experimental data.

In order to check the robustness of our results, we performed several measurements with the goal of falsifying this hypothesis. At first, we set the fluence to 1.5 mJ/cm² and repeated the paired easy and hard axes measurements in three different random sample spots at least a few hundred micrometers apart. This test aimed at checking the robustness of our result across the sample surface. Then, we set the fluence to 2 mJ/cm² and chose three additional sample spots. In this case, first, we measured the hard axis and then the easy axis, in order to exclude that the anisotropic relaxation is related to a slow sample degradation and hence that it is an effect of the measurement sequence. Finally, we performed a measurement at a new sample spot in order to check one intermediate fluence of 1.7 mJ/cm². We observe some variation of the absolute values of the extracted time constants, but in all cases we found that the magnetization dynamics along the easy axis recovers always slower than the hard axis. Taking the average value of all these ratios indicates that the recovery time along the easy axis is $33±16%$ longer than the hard axis. A detailed summary of the measurement protocol and of the data analysis is given in the supplementary material.

We now turn to the transient reflectivity data shown in Fig. 4. Along both the [0001] (easy) and [11$2¯$0] (hard) axes, there is a measurable delay between the maximum change in reflectivity and the maximum change in magnetization, by approximately 250 fs for the easy magnetization axis and 150 fs for the hard magnetization axis. This is similar to what was observed in Ref. 6 for bcc iron pumped with 800 nm pulses and probed at shorter wavelengths (⁠$500–540$ nm). As stated in that work, this indicates that the spin dynamics follows the onset of a non-equilibrium electronic distribution.

In order to explain the magneto-crystalline anisotropy of the de- and re-magnetization time constants, first, we consider the overall picture for energy transfer from the laser-excited conduction electrons. The conduction electrons thermalize within 100 fs and transfer their energy to the cold phonons and magnons due to electron–phonon and electron–magnon coupling. On a timescale longer than the initial electron thermalization, this energy transfer process is reasonably well described by the two-temperature⁴⁴ or three-temperature model.¹ The energy that the magnon and phonon systems can receive is proportional to their heat capacities C_s and C_ph. As the former is small, it is often approximated as zero,⁹ i.e., the dominant energy flow occurs to the lattice. The rate of increase in the lattice temperature is then given by $Cph ∂Tph/∂t=G(Te−Tph)$⁠, where T_e (T_ph) is the electron (lattice) temperature and G is the electron–phonon coupling constant.^45,46 This quantity can be computed ab initio; it is given by⁴⁵ 

where $g(εF)$ is the electronic density of states at the Fermi level ε_F and $α2F(Ω)$ is the Eliashberg function, with Ω the phonon frequency variable (for explicit expressions, see Refs. 26 and 30).

To investigate the dependence of G on the magnetization axis, we calculated the phonon spectra and electron–phonon matrix elements of hcp Co self-consistently, using the ELK full-potential code.⁴⁷ Note that the spin–orbit interaction was included, which allows us to examine the influence of the magnetization axis on G. We find a significant magneto-crystalline dependence of G: for magnetization along [0001], the calculated coupling was $G=1.8×1018$ W m⁻³ K⁻¹, while for magnetization along [11$2¯$0], it was $2.8×1018$ W m⁻³ K⁻¹. Hence, the electron–phonon coupling for M along the hard axis is about 50% larger as for M along the easy axis.

Analyzing next where the magneto-crystalline anisotropy in G comes from, we found that the differences between the phonon dynamical matrices and phonon spectra computed for the hard and easy magnetization axes are small (see the supplementary material for the computed phonon spectra). However, we found that the electron–phonon interaction is much more sensitive to the magnetization direction, due to small shifts of the energy levels near the Fermi energy induced by spin–orbit coupling. The biggest changes in the electron–phonon interaction were obtained for high-energy phonons.

The implication of the larger electron–phonon coupling G for the hard magnetization axis is a stronger transfer of energy from hot electrons to cool phonons, and thus, a faster electron cooling leads to a faster remagnetization, which is fully consistent with the magnetization dynamics measured for t > 0.5 ps. The high magneto-crystalline anisotropy of the remagnetization rate corroborates that the recovery trend is driven by a mechanism that depends strongly on spin–orbit interaction, as the Elliott–Yafet electron–phonon spin-flip scattering. Note that this does not exclude the involvement of magnons in the demagnetization signal. Transverse spin excitations have been detected in demagnetizing ferromagnetic films^6,19 and are expected to be responsible for the MOKE signal here as well, but for the magnetization recovery, the energy flow to the lattice is responsible. Furthermore, we observe that, also, the Gilbert damping α is anisotropic with respect to the magnetization axis (see the supplementary material). In particular, we estimate $αha≈0.11$ when the magnetization is along the hard axis and $αea≈0.085$ when it is aligned parallel to the easy axis. The ratio $αha/αea≈1.29$ is approximately the same as the ratio of the inverse of the magnetization relaxation time constants for the two respective axes and similar to the calculated ratio of 1.55 for electron–phonon coupling constants. This suggests that the spin–orbit coupling, responsible for the anisotropic coupling to the lattice, is also responsible for the anisotropic damping. We do not speculate further on this observation and leave it to future works, given that a much larger anisotropic Gilbert damping has been recently reported.⁴⁸ Note that a possible strong coupling between magnetic and elastic properties in Co has already been pointed out as an explanation for the unusual pressure dependence of the sound velocity in Co.⁴⁹ 

An accurate description of the demagnetization dynamics in the first few hundred fs is, however, a more complex issue. Note that it is observed only at very high fluences. For the here-used Co sample with a 3-nm Pt cap layer, it is conceivable that optically induced superdiffusive spin currents^15,16,24 travel immediately with ballistic Fermi velocities (∼1 nm/fs) into the Pt layer and will contribute to the initial demagnetization. These superdiffusive spin currents j_s arise from the quenching of the Co atomic moment, $∂M/∂t∝js$⁠, which is practically isotropic and independent of spin–orbit interaction. Considering next an electron–phonon picture, the transfer of spin angular momentum from the electrons to the phonons is given by the Elliott–Yafet electron–phonon spin-flip scattering,⁹ which is exactly described by the spin-flip Eliashberg function (see Refs. 26 and 30). This quantity has a very similar spectral dependence as the conventional $α2F(Ω)$⁠, but it is about 40 times smaller.²⁶ It has nevertheless the same magneto-crystalline anisotropy as the common $α2F$⁠. The electron–phonon spin-flip scattering for M along the hard axis is, thus, larger, which would imply a faster magnetization decay in the first few hundred fs, consistent with our measurements. It needs to be emphasized, though, that in this time interval, there will be nonthermal electron populations that depend on the used fluence and, as mentioned, non-equilibrium processes as superdiffusion will be involved, which limit the validity of the two-temperature model as well as of the here-used quasi-equilibrium electron–phonon scattering description. We can, therefore, only conclude that the right trend is given on the very short timescale.

In conclusion, we performed ultrafast magneto-optical pump–probe experiments on epitaxial hcp cobalt, in order to measure the magnetization dynamics along the easy and hard magnetization axes. We observed a systematic 33% slower quenching and relaxation dynamics along the easy magnetization axis. Our ab initio calculations reveal a large magneto-crystalline anisotropy in the electron–lattice coupling and Elliott–Yafet spin-flip scattering, which explains the observed anisotropic magnetization dynamics. The interplay between the different temporal scales could be important for the overall demagnetization or switching processes. Furthermore, our study introduces a different approach to probe, using wavelengths in the optical range, the role of the lattice anisotropy in ultrafast magnetism. We envision that future experiments that mimic our approach will be able to explore other crystalline materials with well-defined lattice structures. The investigation of model systems, as opposed to polycrystalline ones, allows, moreover, for theoretical models to be tested to a greater accuracy. We anticipate that such studies may give important hints toward completely solving the question of the dissipation of angular momentum at ultrafast time scales, which is yet not settled after more than two decades of research.

See the supplementary material for (i) fitting procedure; (ii) static characterization of the sample: crystalline, magneto-crystalline, and magneto-optical measurements; (iii) dynamic characterization: transient reflectivity dynamics and FMR; and (iv) ab initio calculation of the phonon spectra.

The authors gratefully acknowledge B. Wehinger for useful discussion. V.U. and S.B. acknowledge support from the European Research Council, starting Grant No. 715452 “MAGNETIC-SPEED-LIMIT.” R.M. and E.E.F. were supported by U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Contract No. DE-SC0018237. P.M.O. acknowledges the support from the Swedish Research Council (VR), the Knut and Alice Wallenberg Foundation (Grant No. 2015.0060), the CRC/TRR 227 “Ultrafast Spin Dynamics,” and the Swedish National Infrastructure for Computing (SNIC). K.C. acknowledges support from the Czech Science Foundation (Grant No. 18-07172S) and The Ministry of Education, Youth and Sports Large Infrastructures for Research, Experimental Development and Innovations project “e-Infrastructure CZ-LM2018140.”