Coherent helicity-dependent spin-phonon oscillations in the ferromagnetic van der Waals crystal CrI3

The discovery of two-dimensional systems hosting intrinsic magnetic order represents a seminal addition to the rich landscape of van der Waals materials. CrI3 is an archetypal example, where the interdependence of structure and magnetism, along with strong light-matter interactions, provides a new platform to explore the optical control of magnetic and vibrational degrees of freedom at the nanoscale. However, the nature of magneto-structural coupling on its intrinsic ultrafast timescale remains a crucial open question. Here, we probe magnetic and vibrational dynamics in bulk CrI3 using ultrafast optical spectroscopy, revealing spin-flip scattering-driven demagnetization and strong transient exchange-mediated interactions between lattice vibrations and spin oscillations. The latter yields a coherent spin-coupled phonon mode that is highly sensitive to the driving pulse’s helicity in the magnetically ordered phase. Our results elucidate the nature of ultrafast spin-lattice coupling in CrI3 and highlight its potential for applications requiring high-speed control of magnetism at the nanoscale.

S ince its first demonstration 1 , the ultrafast manipulation of magnetism has been a major topic of research, with significant efforts aimed at unraveling the nature of dynamic demagnetization 2,3 , the excitation of collective magnetic modes 4,5 , and the realization of all-optical switching of magnetic phases [6][7][8] . The recent emergence of two-dimensional (2D) van der Waals (vdW) materials hosting intrinsic long-range magnetic order 9-11 presents a new, relatively unexplored platform for investigating these phenomena in systems where the interplay between structural order and exchange interactions plays a pivotal role. Moreover, the inherently strong light-matter interactions typical of vdW materials open exciting new possibilities for exploring the fundamental properties of magneto-structural coupling at the nanoscale under the influence of intense optical excitation [12][13][14] . This in turn has the potential to yield revolutionary leaps in low-dimensional magneto-optical devices for future storage and quantum information applications, especially when used in heterostructured device architectures 15 .
CrI 3 is a prototypical example of a magnetic vdW material hosting out-of-plane ferromagnetic (FM) order down to the monolayer limit 10 stabilized by uniaxial anisotropy, which opens a gap in the low energy zone-center magnon spectrum 16 . Even more intriguing is the fact that the magnetic phase is highly dependent on the number of atomic layers. This is due to the nature of the fundamental Cr-Cr superexchange interactions 17 , which stabilize an antiferromagnetic (AFM) phase with broken inversion symmetry in even-layered crystals. This gives rise to unique nonlinear optical phenomena 18 and highlights the strong interdependence between the various crystal structures (Fig. 1a), magnetic orders, and optical responses present in CrI 3 . Additionally, the presence of phonon-induced lattice distortions that can modulate the length scale of spin-spin interactions may unlock the possibility for dynamic light-driven coherent coupling between the lattice and spin degrees of freedom 19 .
Here, we uncover the nature of this coupling by measuring the transient magnetization and coherent vibrational dynamics in CrI 3 after femtosecond optical photoexcitation. We find that the demagnetization dynamics are driven by spin-flip scattering, while coherent pump helicity-dependent oscillations in the timeresolved polarization rotation (TRPR) signal highlight the strong influence of magnetic order on the c-axis A 1g optical phonon mode. This can be interpreted within the framework of a dynamic spin-lattice coupling mechanism, which opens a unique pathway for manipulating magnetic order through the vibrational degree of freedom, with significant implications for future nanoscale opto-magnetic applications.

Results
Ultrafast Demagnetization. In our experiments, 1:55 eV (800 nm), 85 fs pump and probe pulses (measured at the sample position) were focused at near-normal incidence onto as-grown bulk-like flakes of CrI 3 (Curie temperature, T c ¼ 61 K) (Fig. 1b). Figure 2 shows the resulting TRPR signal obtained in the FM phase at T ¼ 15 K under right circularly polarized (σ þ ) pumping. We observe a rapid sub-picosecond decrease in the magnetization, a further reduction over tens of picoseconds, a 16 GHz oscillation due to a coherent strain wave, and an eventual nanosecond timescale recovery. This two-step demagnetization process, referred to as type-II dynamics, has been observed in other ferromagnets 20,21 and occurs when demagnetization is not completed before electron-phonon equilibration 2 .
More insight can be obtained by considering the optical transitions in CrI 3 associated with 1:55 eV photoexcitation. The absorption spectrum of bulk CrI 3 shows a weak resonance peak at 1:5 eV 22 , attributed to transitions between the partially filled t 2g and unfilled e g levels resulting from the octahedral ligand fieldinduced splitting of the Cr 3+ d-orbital, as shown in Fig. 1c. Despite the even parity of these states, this transition is allowed due to mixing with various odd-parity states 22 , and may also be Schematics of a the crystal structure of the rhombohedral (R 3) and monoclinic (C2=m) phases of CrI 3 , b the femtosecond TRPR experiment, where QWP is a quarter wave plate, WP is a Wollaston prism, and BD is a balanced detector, and c the relevant energy levels in CrI 3 , where the red arrow denotes optical transitions at $ 1:5 eV that generate excited electrons and holes. Dashed lines indicate the dominant atomic orbital character of the ligand-field split levels. Transitions between the e g , t 1u , and t 2u levels significantly contribute to coherent phonon generation. Fig. 2 Ultrafast photoinduced demagnetization in CrI 3 . Time-resolved polarization rotation signal under σ þ pumping at T ¼ 15 K, where the blue dots are experimental data and the green trace is a fit using the M3TM (see the SI for more details). The oscillatory component in the M3TM fit curve was obtained by separately fitting the oscillatory component of the experimental data with a decaying sinusoidal function. The inset shows the magnetization dynamics on shorter timescales under σ þ (blue), σ À (red), and linearly polarized (green) pumping at T ¼ 15 K. associated with a low-lying bright charge-transfer exciton 23 . Regardless, below T c the 1:55 eV pump pulse drives transitions predominantly involving majority-spin states, due to significant spin polarization of the valence and conduction bands 23,24 . This leads to a strong preferential absorption of σ þ light in spin-up domains at 1:55 eV 23 .
The helicity-dependent absorptivity of CrI 3 suggests that the ultrafast magnetic response should be nearly negligible for left circularly polarized (σ À ) pump pulses. However, as seen in the inset of Fig. 2, we observe that the TRPR signal under σ À pumping has a sign opposite to the σ þ case. This helicitydependent response is preserved for all temperatures T < T c . We attribute this to the presence of multiple domains within the photoexcited region 15,25 ; i.e., spin-down domains preferentially absorb σ À photons, leading to their subsequent demagnetization. Nevertheless, horizontally polarized pumping yields a negative signal ( Fig. 2 inset), implying that the volume fraction of spin-up domains is higher in the probed region.
Accordingly, upon absorption of a σ þ pulse, majority-spin electrons and holes are excited in CrI 3 . Due to spin-orbit coupling (SOC), the wave functions of these carriers are a mixture of pure spin states 26 . This allows for a finite probability of Elliott-Yafet spin-flip scattering processes 27 , which mediate spin relaxation, particularly in the hole population due to the comparatively small valence band exchange splitting 23 . Furthermore, the photoexcited carrier density of $10 19 cm À3 exceeds the Mott density criterion, given the nanometer-scale exciton radius 23 , leading to a transient quasi-metallic state. This makes it possible to describe ultrafast dynamics in CrI 3 using the microscopic three-temperature model (M3TM), in which the excited electronic system supplies the energy for demagnetization, while interactions with the lattice allow for angular momentum dissipation 2 . A fit to the demagnetization dynamics with the M3TM, shown in Fig. 2, accurately reproduces both demagnetization steps with a spin-flip probability of a sf ¼ 0:175, consistent with other materials demonstrating type-II dynamics 20,21 . More detail is included in the supplementary information (SI).
Coherent dynamics. We now turn our attention to the coherent dynamics. The top panel of Fig. 3a shows the TRPR signal at T ¼ 75 K > T c under σ þ and σ À pumping after subtracting the non-oscillatory background. Above T c , the background signal is relatively small (as shown in Supplementary Fig. 4), and most likely originates from optical effects such as pump-induced birefringence, while below T c , the background is large and primarily due to the ultrafast demagnetization phenomena discussed in the previous section. We observe pronounced coherent oscillations in the time domain, and the corresponding power spectral densities (PSDs) plotted in the bottom panel reveal two distinct modes at~2:37 THz and~3:87 THz, corresponding to in-plane (A 1 1g ) and c-axis (A 2 1g ) Raman-active phonons, respectively [28][29][30][31] . These modes are excited via impulsive stimulated Raman scattering (ISRS), which is permitted by their A 1g symmetry 30,31 .
Notably, the PSD of both modes is nearly identical for both pump helicities above T c . However, as shown in Fig. 3b, this symmetry is broken below T c , where σ À pumping leads to a significantly smaller A 2 1g amplitude. In contrast, the amplitude of the A 1 1g mode remains relatively invariant to the choice of pump polarization. The temperature and pump-helicity dependent asymmetry in the phonon amplitude can be clearly seen in Fig. 3c, where we plot the σ þ =σ À ratio of the integrated spectral peaks of both oscillatory modes vs. T. Here, the A 2 1g ratio below T c consistently follows a FM order-parameter-like function, / ffiffiffiffiffiffiffiffiffiffiffiffiffiffi T c À T p . This highlights the strong sensitivity of the coherent amplitude of the A 2 1g mode to the underlying magnetic order in the system. In contrast, the A 1 1g mode ratio shows minimal variation with temperature. This difference is striking and suggests that the helicity dependence of the A 2 1g mode does not originate from a dominant optical effect (i.e., the larger absorptivity of σ þ light in CrI 3 below T c ), as this would impact the two modes in a similar manner. Additionally, we note that time-reversal symmetry breaking associated with magnetic ordering can lead to non-zero off-diagonal tensor elements for the A 1g modes 32 , which in turn can lead to helicity-dependent asymmetries in the equilibrium spectral response 31,33 . However, this is unlikely to be the origin of the pump helicity-dependent effects we observe in our experiments. First, the influence of offdiagonal tensor elements should be present for all modes of A 1g symmetry, in stark contrast to our results where only the c-axis mode is sensitive to the pump polarization; this strongly implies that mode symmetry alone cannot account for our observations. Second, recent ab initio calculations show that the magnitude of the off-diagonal elements is highly sensitive to the photon energy, becoming non-trivial only in the vicinity of the dominant c The σ þ =σ À ratio of the integrated Fourier transform peaks of the measured signal at the A 1 1g (orange) and A 2 1g (green) mode frequencies, and the simulated helicity-dependent ratio at the A 2 1g frequency (pink) vs. normalized temperature. The solid green line is a fit using a FM order-parameter-like function / ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi T c À T p and the solid orange line is a guide to the eye. Error bars for the experimental data in c were obtained from a bootstrap sampling analysis. NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-31786-3 ARTICLE ligand-metal charge transfer resonances above 2 eV 34 . Our photon energy is far below these resonances, further supporting an alternative mechanism for the helicity-dependent effects observed here. Finally, we note that fitting our time-domain data with a multi-component decaying sinusoidal model reveals the presence of an additional oscillatory mode below T c , consistent with recent reports of AFM order dependent modes at 3:73 THz, close to the A 2 1g mode frequency 31,33,35 . However, this mode does not show a pronounced helicity dependence in our measurements (more detail is provided in Supplementary Fig. 5 and the associated discussion).
To understand the clear sensitivity of the A 2 1g mode to the underlying magnetic order, we recall that this mode corresponds to a c-axis oscillation of iodine atoms (shown schematically in the lower panel of Fig. 3a) leading to an oscillatory trigonal distortion of the CrI 6 octahedra. This distortion can transiently modulate spin exchange interactions through anisotropy-or exchangemediated pathways 19,36 . We confirmed this using density functional theory (DFT) to calculate ΔE ¼ E FM À E AFM , the energy difference per spin between FM and AFM configurations, which is a measure of the overall exchange strength. Tracking its dependence on the average equilibrium lattice displacement, X, of the A 2 1g phonon yields ∂ΔE=∂X % À32 meV=Å per Cr atom, signaling the importance of spin-lattice coupling in CrI 3 . Here, X > 0 denotes trigonal compression of the CrI 6 octahedra.
To explain the helicity-dependent TRPR, we then developed a spin-phonon model for CrI 3 . This describes interacting S ¼ 3=2 Cr moments on a layered honeycomb lattice that are coupled to the local A 2 1g Einstein phonon displacements. The total Hamilto- Here X r and P r are the displacement and momentum of the Einstein phonon at lattice site r, M % 6:3 10 À25 kg is the mass of the three iodine atoms per Cr spin, and Ω % 3:87 THz is the A 2 1g phonon frequency. The spin Hamiltonian is where rr 0 h i γ denotes nearest-neighbor (NN) sites within a honeycomb layer with bond direction γ, rr 0 h i h i and rr 0 h i h i h i denote pairs of second-NN and third-NN, respectively, and rr 0 h i c denotes NN sites between two neighboring honeycomb layers. The exchange constantsJ H1;rr 0 , J H2 , J H3 , J K , and J D quantify NN Heisenberg, second-NN Heisenberg, third-NN Heisenberg, Kitaev, and Dzyaloshinskii-Moriya (DM) interactions, respectively, with the DM (unit) vectorsD rr0 as specified in ref. 37 . The inter-layer exchange constant J Hc quantifies the interaction between neighboring layers of the honeycomb lattice. The outof-plane single-ion anisotropy is given by where J A is the anisotropy energy. The A 2 1g lattice distortion modifies the Heisenberg exchange asJ H1;rr 0 ¼ J H1 1 þ ϕ X r þ X r 0 À Á =2 À Á , with spin-phonon coupling parameter ϕ. We neglect phonon coupling to the remaining exchange terms, since they are small in comparison. Guided by previous ab initio calculations [38][39][40] , we fix ðJ H1 ; J H2 ; J H3 ; J K ; J D ; J A ; J Hc Þ ¼ À2:4; À0:1; 0:1;0:9; À0:2; ð À0:2; À0:59Þ meV to capture the spin-wave dispersion 16,37 and the ordering temperature T c (see the SI for more details). Setting ϕ ¼ 1:98 Å À1 then reproduces our DFT results for ∂ΔE=∂X. For T ( T c , our model predicts an equilibrium displacement X h i $ 6 10 À4 Å, comparable to the experimentally estimated magnetostriction 41 . Table 1 summarizes the quantities relevant to the dynamical spinphonon model and static spin exchange parameters.
Using this spin-phonon model, we simulate the impact of the pump pulse as an instantaneous helicity-dependent lattice distortion, X r; where m r is the local projection of the spin onto the equilibrium magnetization axis, and the ± sign corresponds to the σ ± pump helicity. ξ 1 determines the overall strength of the distortion due to the Raman process, and ξ 2 represents its helicity dependence; the latter can arise from the spin selectivity of the helicity-dependent photoexcitation, which transiently enhances or suppresses the local Cr moment 42 . The functional form of the lattice distortion X r; ± t ð Þ can be rationalized as follows: from a mean-field perspective, the equilibrium lattice displacement scales as X r ¼ αm 2 r , for some proportionality constant α (see the SI for more details). Assuming that the pump pulse induces a pump helicity-dependent magnetization change ±δ, the phonon equilibrium position is shifted to X 0 r; ± ¼ αðm r ± δÞ 2 , which is a special case of the aforementioned functional form, with ξ 1 ¼ αδ 2 and ξ 2 ¼ 2=δ, that only depends on the single parameter δ (see the SI for details).
We solve for the subsequent dynamics of the spin-phonon system by numerically integrating the coupled equations of motion (see Methods). For T > T c , we find persistent oscillations in the average phonon displacement X ± t ð Þ ¼ 1 N ∑ r X r; ± ðtÞ at the phonon frequency Ω, while the uniform magnetization m ± t ð Þ does not exhibit coherent dynamics (see Supplementary Fig. 9d). Remarkably, in the FM phase (T < T c ), the distortion leads to coupled, coherent oscillations in X ± t ð Þ and the magnetization m ± t ð Þ at Ω, as shown in Fig. 4. The coherent m ± t ð Þ oscillations lead to oscillations in the polarization rotation, creating an additional temperature and helicity-dependent contribution to the TRPR signal below T c . Setting ξ 2 ¼ 0:18=μ B , we find that m þ Ω ð Þ=m À Ω ð Þ plotted vs. T=T c shows excellent agreement with the experimental TRPR ratio below T c (Fig. 3c). Furthermore, in our model, the choice of ξ 2 also fixes the overall distortion strength, ξ 1 ¼ 0:01 Å. We can use this to calculate that the resulting coherent oscillations in the magnetization are of magnitude Δm σ ðtÞ=m $ Oð10 À4 Þ (c.f. Fig. 4b), consistent with our experimental observations.

Discussion
Ultimately, the overarching physical mechanism underlying the helicity-dependent response of the A 2 1g mode centers upon the coherent coupling between the magnetization and the A 2 1g phonon. Above T c , the femtosecond optical pump drives transitions involving the Cr-like e g level. The partial occupation of this level leads to a strongly Jahn-Teller active ion, causing the system to undergo an ultrafast trigonal distortion that triggers the coherent phonon mode. We observe these oscillations above T c in the polarization sensitive detection configuration due to the imperfect cancellation between the p-and s-polarized components of the probe pulse, an effect that has been observed for A 1g phonons in other systems 43,44 . This gives us a baseline signal in the nonmagnetic phase that serves as a comparison to the response below T c . Below T c , the modulation of the spin exchange by this phonon mode leads to the intertwining of lattice vibrations and coherent spin oscillations. Moreover, symmetry-allowed virtual transitions from strongly iodine 5p-like states (e.g., t 1u and t 2u , see Fig. 1c), which occur simultaneously with the real t 2g to e g transitions described earlier, lead to helicity-dependent changes in the local single ion moment, transiently enhancing (suppressing) it under σ þ (σ À ) photoexcitation in a majority spin-up photoexcited volume due to optical selection rules. For example, for a spin-up domain, a right circularly polarized pulse would excite a spin-up electron from the t 1u or t 2u level to the partially occupied e g level. This transiently enhances the local single-ion moment, prior to the system returning to the ground state through phonon emission. As shown above, this contributes a helicity-dependent component to the impulsive lattice distortion, akin to a dynamic magnetostriction process. Accordingly, this gives rise to a measured oscillatory signal that is an admixture of the non-temperaturedependent vibrational contribution (analogous to the T > T c case) and a temperature-dependent polarization modulation driven by the magneto-optical Kerr effect through spin-phonon coupling. While there is no symmetry constraint preventing this from also occurring for the A 1 1g mode, our DFT calculations reveal that the exchange modulation is approximately a factor of 2 smaller than for the A 2 1g mode. Therefore, we hypothesize that the lack of appreciable helicity dependence for this mode is due to the combined effects of the smaller spin-phonon coupling and smaller lattice distortion associated with the A 1 1g mode as compared to the A 2 1g mode, which can be inferred from the smaller Raman peak intensity observed via spontaneous Raman scattering 31 and its comparatively smaller oscillatory amplitude, as observed in our experiments.
Our work demonstrates several key findings in the burgeoning field of the ultrafast control of vdW magnets. First, the inherently strong spin-phonon coupling in these materials provides a fruitful route for the control of both the vibrational and magnetic dynamics via photon helicity, allowing for a new degree of control over their transient properties. Second, while previous spontaneous Raman results required resonant excitation, our observation of coherent spin-phonons driven via ISRS with photons well below the dominant ligand-metal charge transfer resonances demonstrates that femtosecond optical spectroscopy may be a significantly more sensitive probe of vibrational and magnetic excitations in CrI 3 . Moreover, it may provide a pathway to more completely characterize the vibrational modes of CrI 3 , especially with regard to differences in the tensors governing stimulated and spontaneous Raman scattering from phonons 45,46 . Finally, our observation of an additional 3:73 THz vibrational mode in bulk CrI 3 is in stark contrast to past spontaneous Raman measurements, where it was only observed in few-layer flakes. The persistence of this mode in thick samples, which was previously linked exclusively to near-surface AFM order in the few-layer crystal 30,33 , suggests that it is far more robust than previously believed. Moreover, in our experiments, any sort of a near-surface AFM region would make up an extremely small percentage of the overall photoexcited volume due to the large optical penetration depth, belying the clear signature of this mode below T c . Instead, our simultaneous observation of modes linked to both FM and AFM order in truly bulk samples is in excellent agreement with a recent report demonstrating the coexistence of monoclinic and rhombohedral structural phases (Fig. 1a) in bulk CrI 3 over a wide temperature range 47 . Together with our results, this suggests that FM and AFM exchange couplings and phases may be present throughout the volume of the crystal 47 , rather than restricted to the near-surface region, as was previously hypothesized.
We note that the implied magnetization changes in our theoretical picture are in excess of the saturation magnetization, suggesting that the microscopic mechanisms which drive the lattice displacement are not fully captured by our phenomenological model. Their precise nature thus remains to be uncovered and would be a fruitful avenue for future research. Another crucial area of future experimental efforts will be to study the impact of sample thickness on non-equilibrium magnetization dynamics, which would be critical for future low-dimensional device applications. Finally, a fundamental challenge to overcome is the relatively low damage threshold of CrI 3 , especially under focused ultrafast illumination, which makes fluence-dependent measurements a distinct challenge. This may be circumvented, however, by driving the system at longer wavelengths, where the absorption is significantly lower. Here, nonlinear effects are easier to drive, due to the ability to use higher pump fluences, and may prove successful in more efficiently inducing transient magnetic phenomena in CrI 3 .
In conclusion, our measurements and simulations highlight the inextricable link between magnetic order and structure in CrI 3 , due to the intimate relationship between the strength of exchange interactions, trigonal lattice distortions, and the orbital character of the states involved in the photoexcitation process. This in turn allows for greater flexibility in the manipulation of magnetization and vibrational dynamics in vdW materials, especially using all-optical techniques that exploit nonlinear processes to drive coherent phenomena. Finally, our results shed new light on the immense potential for CrI 3 and similar 2D magnetic materials in the next generation of opto-magnetic technologies and provide insight into enabling ultrafast optical control of magnetism at the atomic limit.

Methods
Time-resolved experiments. Bulk-like flakes of CrI 3 , grown using chemical vapor transport, were deposited onto an Au coated sample holder and placed in a variable temperature liquid helium flow optical cryostat. Ultrafast pulses were supplied by a regeneratively amplified Ti:Sapphire laser, which generated 1:55 eV pulses (800 nm central wavelength) with an 85 fs duration (measured at the sample position) at a 100 kHz repetition rate. The beam was split into pump and probe arms using a plate beam splitter. The pump beam was sent through a mechanical delay line, and both beams were passed through achromatic wave plates and subsequently focused onto optically flat regions of the sample at near-normal incidence using a 20X near-IR optimized apochromatic objective. The nominal 1=e 2 focal spot diameter of the probe was approximately 10-15 μm, and the pump was approximately 20 μm, yielding a pump fluence of $0:64 mJ=cm 2 . The latter was just below the damage threshold of the sample (i.e., yielding no observable signal degradation or photodarkening of the sample), while also being the lowest fluence we could utilize while still maintaining a sufficient ratio with respect to the probe fluence. We performed experiments where the pump polarization was horizontally polarized (i.e., parallel to the optical table surface), right circularly polarized, and left circularly polarized. The probe beam was linearly polarized, and upon reflection, the beam was decomposed into its two orthogonal linear components using a Wollaston prism. These components were then differentially detected using a balanced Si photodiode detector and the output signal was fed into a lock-in amplifier to isolate the TRPR signal. Both the pump and probe beams were modulated using a 7/5 slotted optical chopper, allowing the lock-in amplifier to be referenced to the sum intermodulation frequency (1:2 kHz) to eliminate pump scattering effects.
Density functional theory calculations. Our ab initio simulations were carried out within density functional theory. We used both the pseudopotential plane wave method implemented in the Vienna ab initio simulation package (VASP) 48,49 and in Quantum Espresso (QE), and the full-potential all-electron linearized augmented plane wave (FP-LAPW) method implemented in the Wien2k code 50 . The local density approximation was used for the exchange-correlation functional throughout all our ab initio simulations. For the VASP calculations, we chose an energy cutoff of 500 eV for the plane-wave basis set and used a 10 × 10 × 10 (10 × 10 × 3) Γ-centered k-point mesh to sample the bulk (monolayer) Brillouin zone to perform the structure optimization, Γ-point phonon, Raman tensor (aided with the vasp_raman python script 51 ), and electronic structure calculations. The QE calculations used a plane-wave basis set with a cutoff energy of 60 Ry. A 10 × 10 × 10 k-point mesh was used to calculate the electronic structure of bulk CrI 3 and to support a calculation of the phonon dispersion at the Γ-point. In the Wien2k simulations, we focused on calculating the electronic band structure for bulk CrI 3 by taking a 15 × 15 × 15 k-point mesh with a muffin-tin radius of 2:50a 0 (Cr) and 2:35a 0 (I), where a 0 is the Bohr radius.
We also used VASP to perform total energy calculations for the FM and AFM states of a CrI 3 monolayer, from which a Heisenberg model with nearest-neighbor exchange interactions was fit to obtain the exchange interaction. The spin-lattice coupling strength was then determined by distorting the equilibrium lattice according to the displacement field for a chosen vibration mode. The electronic structure obtained in our simulations on both bulk and 2D CrI 3 is consistent with those reported in the literature 28,52,53 . We note that SOC reduces the band gap. However, since SOC does not qualitatively change the lattice dynamics 28 , we did not include it in the calculation of the Raman tensor and estimate of the spin-lattice coupling strength.
Additionally, we performed LDA+U calculations using a Hubbard parameter of U eff ¼ 2:2 eV on the Cr 3d orbitals 52 to compare with the LDA calculations. Our results are in good agreement with those reported by Liu et al. 54 . We found the lattice constants change slightly, which suggests a minor change of phonon modes. This is consistent with the results found by similar studies 39,55 . For the electronic properties, as shown in Supplementary Fig. 13 of the SI, we found the spin down Cr-3d bands are pushed to higher energy whereas the energy gap remains at about 1 eV, close to that with LDA. This observation is also consistent with that reported by Liu et al. for monolayer CrI 3 54 . This band structure change could affect the photoinduced carrier density and the initial relaxation. However, it is not likely to change the long-time recombination and relaxation, and the physics of the spinlattice coupling driven relaxation process.
Monte Carlo simulations. We performed classical Monte Carlo simulations on the coupled spin-phonon system described by the Hamiltonian H ¼ H ph þ H sp specified in Eqs. (1)-(3) to generate a thermal ensemble of states. We studied lattices of 2 × L × L × L c spins, with system sizes up to L ¼ 32 and L c ¼ 8, with periodic boundary conditions. To improve the statistical convergence of the simulations, we employed a parallel tempering scheme 56 to simulate 144 logarithmically spaced temperature points between T min % 9 K and T max % 120 K in parallel. Each simulation was equilibrated for 10 6 sweeps before taking measurements of the specific heat and the equilibrium magnetization for an additional 1:5 10 7 sweeps; a single sweep is defined as one attempted update per spin or phonon degree of freedom. As shown in the SI, the chosen coupling constants reproduce the experimental spin-wave dispersion 37 . In addition, the computed temperature-dependent specific heat and magnetization yield T c ¼ 42:8 K which shifts to 71 K when considering quantum corrections by rescaling to the effective spin length S eff ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi , with S ¼ 3=2 for Cr (see the SI for more details). This value is in reasonable agreement with the experimental result T c ¼ 61 K.
Dynamical Simulations. We modeled the impact of the pump laser as a helicitydependent local lattice distortion, given by X r; ± t ¼ 0 ð Þ¼X r þ ξ 1 1 ± ξ 2 m r À Á at time t ¼ 0. Here, X r is the equilibrium lattice displacement at lattice site r, m r is the local magnetization of the spin at site r along the global equilibrium magnetization direction, and the sign ± distinguishes σ þ and σ À polarized pump pulses. Such a distortion could arise from an ultrafast trigonal splitting of a photoexcited Jahn-Teller active Cr e g level, with the splitting being sensitive to the local Cr spin due to Hund's coupling. We applied this ultrafast distortion to individual configurations selected from our Monte Carlo ensemble with L ¼ 32 and L c ¼ 8. The post-distortion dynamics were described using the Landau-Lifshitz equations for spins given by and dP r dt In Eq. (4), r γ denotes the NN lattice site of r along the in-plane bond direction γ and r 0 (r 00 ) denotes the second (third) in-plane NN; the inter-plane NN direction is denoted by r c . These equations were numerically integrated using a ninth-order Runge-Kutta algorithm 57 and adaptive time steps with a local relative error tolerance of ε rel ¼ 10 À13 . This yielded converged results up to timescales of approximately 50 ps in the magnetically ordered phase. We averaged the resulting dynamical observables over 10000 Monte Carlo configurations with positive magnetization (when projected onto the twofold degenerate polarization axis). In particular, we extracted the time-dependent deviation of the magnetization, Δm ± t ð Þ ¼ hm ± t ð Þi À m, where m is the mean magnetization in thermal equilibrium and hm ± t ð Þi is the time-dependent magnetization after the ultrafast σ ± -helicity-dependent distortion was applied. We determined the ratio of Fourier components m þ Ω ð Þ=m À Ω ð Þ by calculating the Fourier transform m ± ω ð Þ ¼ FT Δm ± t ð Þ È É and fitting a Gaussian profile to the peak found at the A 2 1g phonon frequency Ω. The ratio m þ Ω ð Þ=m À Ω ð Þ is independent of the overall distortion ξ 1 , while a helicity-dependent splitting ξ 2 % 0:18=μ B was found to reproduce our experimental results below T c . Additional results on the time-resolved magnetization and phonon displacement are shown in the SI. Finally, our results are qualitatively robust upon tuning the Hamiltonian parameters, as long as the Kitaev exchange remains finite. In particular, our results also hold for a single-layer honeycomb model as well as for an alternate model with dominant Kitaev interactions, which has been proposed in Ref. 16 .

Data availability
The data that support the findings of this study are available from the corresponding authors upon reasonable request, in order to comply with Los Alamos National Laboratory policy on data security.

Code availability
The Vienna ab initio simulation package is available at https://www.vasp.at. The Quantum Espresso package is available at https://www.quantum-espresso.org. The Wien2k package is available at http://susi.theochem.tuwien.ac.at. The Monte Carlo and dynamical simulation codes are available from the corresponding authors upon reasonable request.