Exchange-striction driven ultrafast nonthermal lattice dynamics in NiO

We use femtosecond electron diffraction to study ultrafast lattice dynamics in the highly correlated antiferromagnetic (AF) semiconductor NiO. Using the scattering vector (Q) dependence of Bragg diffraction, we introduce a Q-resolved effective lattice temperature, and identify a nonthermal lattice state with preferential displacement of O compared to Ni ions, which occurs within ~0.3 ps and persists for 25 ps. We associate this with transient changes to the AF exchange striction-induced lattice distortion, supported by the observation of a transient Q-asymmetry of Friedel pairs. Our observation highlights the role of spin-lattice coupling in routes towards ultrafast control of spin order.

Due to strong correlations, ab-initio descriptions of this large band gap charge-transfer semiconductor are challenging [5][6][7][8][9][10][11][12]. The open Ni d-shell and strong superexchange leads to high-temperature antiferromagnetic (AF) order (TN ≈ 523 K), making NiO a promising candidate for room-temperature spintronic applications [13][14][15]. In this context, pioneering experiments have demonstrated coherent excitation of high-frequency magnon modes by THz pulses [16][17][18], fueling the promise of ultrafast AF spintronics. For such purposes, understanding energy transport due to coupling between the material's subsystems is of importance, in particular those to magnetic order. In equilibrium, NiO exhibits strong spin-lattice coupling and exchange striction [19], leading to a rhombohedral lattice distortion (RLD) along the [111] direction of its nominally cubic structure below [20]. Spin-phonon coupling in NiO has been recently studied in the context of magnon damping in devices [21]. However, energy transfer dynamics and couplings between the various subsystems upon optical excitation have been little studied so far [22]. In particular, to date there is no account of ultrafast lattice dynamics in NiO.
In the presence of a band gap, optically excited carriers can radiatively decay, and they can transfer energy to another subsystem, e.g. the lattice. This is typically described using coupled-heat-baths models [23][24][25], where the subsystems' transient temperatures are described by coupled rate equations. While often successfully employed [26,27], an implicit assumption is that the baths themselves remain thermalized, such that the carriers and the phonons always follow Fermi-Dirac and Bose-Einstein distributions. Whereas this is considered valid for many metals because of homogenous electron-phonon coupling and rapid electronic thermalization, this assumption has been recently challenged even for simple metallic systems [28][29][30][31]. In semiconductors, electron-phonon coupling is often strongly heterogenous [32], and the phonon dispersion is often more complicated than in metals.
Consequently, photoexcited semiconductors may experience a prolonged nonthermal lattice state, in which unconventional relaxation processes may occur. To date only a few reports describing nonthermal lattice dynamics of semiconductors have been published [33][34][35][36][37], and details about the underlying microscopic processes are scarce. In complex materials such as NiO, lattice dynamics may be substantially influenced by effects beyond electron-phonon coupling, such as spin-lattice coupling via exchange striction.
Here we use femtosecond electron diffraction (FED) to study photoinduced lattice dynamics in NiO.
Variants of FED have recently been used to study lattice dynamics in several other systems, including metals [38,39], semiconductors [34,40], heterostructures [41,42], and systems involving structural phase transitions [43,44]. We present an approach to characterize the transient lattice state, by converting Bragg reflection intensities into units of Kelvin, producing a series of temperatures associated with each scattering vector of the probe electrons (Q). Through analysis of these Q-dependent effective lattice temperatures we identify a strongly non-thermal lattice state after excitation, which lasts for tens of picoseconds. Employing the Q-dependent Ni 2+ and O 2scattering factors we gain a degree of elementsensitivity, and find that the lattice response to photoexcitation is primarily on oxygen ions, and less on Nickel. Based on the initial energy transfer time scale of ~0.3 ps, we propose a scenario in which photoexcitation perturbs the antiferromagnetically-induced RLD. This is supported by photoinduced changes in observed Q values that occur due to changes in the shape of the unit cell (i.e. the RLD), suggesting that the observed nonthermal lattice state originates from preferential occupation of phonon modes associated with the RLD. We use a 20 nm thick single crystal of NiO, epitaxially grown on a [001]-oriented NaCl crystal, which was subsequently dissolved. Diffraction and Raman measurements confirmed its bulk-like properties (see supplement). Experiments were conducted at room temperature (T0) using a compact setup [45] (Fig.   1a), with femtosecond laser excitation of ℎ = 2.16 eV and 5.2±1.3 mJ cm -2 incident fluence. This photon energy slightly exceeds half the charge-transfer gap of Δ ≈ 3.8 eV [22,46,47], and two-photon absorption is likely dominant over linear absorption (see supplement). 70 keV probe electrons transmit through the sample, producing patterns as in Fig. 1b. The response function is estimated at 200 fs.
The lattice response was measured using diffraction patterns from different pump-probe delays . The We describe Bragg intensities using structure factor calculations. Intensities of (ℎ 0) reflections depend only on the size of Q i : Here and are the scattering factors of the Ni 2+ and O 2ions, which are tabulated functions of Q [48]. and are the Debye-Waller (DW) factors for Ni and O ions, which are tabulated as functions of temperature for NiO [48]. DW factors can be expressed as = ⅔ 2 〈 2 〉 ( represents Ni i The rhombohedral distortion away from 90° reaches ~0.06° at room temperature, so a cubic structure is commonly assumed. Since and in Eq. (1) cannot be analytically separated, we adopt a temperature-based approach.
We insert the tabulated ( ) and ( ) into Eq. (1)   To demonstrate this, Fig. 3b presents as a function of Q at selected delays. Before excitation exhibits no Q dependence, as expected in a thermal state (all = 0 ). After excitation, exhibits an overall increase in temperature due to the excitation, with a continuous reduction upon increasing Q.
This marks a departure from equilibrium behavior. Given the higher sensitivity to oxygen vibrations at lower Q (Fig. 3a), this result indicates that the MSD of oxygen 〈 2 〉 initially grows disproportionately more than that of nickel 〈 2 〉, demonstrating a nonthermal lattice state at early times. This trend is subsequently suppressed, such that at 25 ps the Q-dependence of is nearly flat, indicating thermalization of the lattice. This disproportionate growth is not to be confused with the absolute difference between O and Ni MSDs, which differ also in equilibrium [49].
We quantify this nonthermal disproportionality between the O and Ni vibrational responses by empirically describing these Q-dependences as ∝ (solid lines, Fig. 3b). Fig. 4a presents ( ) at each delay, exhibiting an initial sub-picosecond reduction, and reaching a minimum at 3 ps, corresponding to the most pronounced disproportionality between O and Ni. After ~15 ps, ( ) recovers toward 0, i.e. to a thermal lattice state, which is also reflected in the difference between the curves in Fig. 2b.
( ) also enables a reliable description of the mean lattice response. Averaging -weighted data ( Fig. 2b) produces a mean lattice response ̅ ( ), presented in Fig. 4b. As in ( ), ̅ ( ) exhibits a twostep response, well-described by a biexponential. The best fit (solid line) yields a sub-picosecond process ( 1 = 0.31 ± 0.08 ps) and a second, slower one ( 2 = 4.1 ± 1.3 ps), see dashed lines. The total temperature rise is Δ = 9.6 ± 0.7 , of which 60% is the fast process. Using the NiO lattice heat capacity [50], we find that Δ corresponds to a maximal change of 16 ± 1 meV/unit cell in the lattice energy density. Careful evaluation of the data concluded that the fast process is intrinsic, while the slower process evolves as the measurement progresses (see supplement). Late delays were omitted from the fit because they exhibit recovery, similar to that of ( ) (Fig. 4a). We convert ̅ ( ) back into Bragg intensities by -weighing the ̅ ( ) fit and plugging it into Eq. 1. This produces the lines in Fig. 2a, in reasonably good agreement with the data. Disagreements exist because this washes out the distinctly nonthermal description in Fig. 4a. Our fit does not consider recovery, so at late delays the Qdependence of causes data with low/high Q to be above/below the lines.
To interpret the stronger response of the O ions, we inspect the phonon dispersion of NiO. This reveals several optical phonons that preferentially move the O ions at the Brillouin zone boundaries (e.g. along [111]), in particular the LO' mode [21]. Below , the optical modes' frequencies deviate significantly from the expected temperature dependence of crystal lattice anharmonicity [21], and they split in energy [51] due to the RLD [52,53]. Preferential coupling to such modes would lead to their enhanced population in a nonthermal phonon distribution, followed by thermalization of the excited phonon population to lower energy/momentum modes, on phonon-phonon scattering timescales [32,54]. A scenario leading to a preferential excitation of such modes could be a perturbation of the AF-induced RLD, for which a full crystallographic account was only recently reported [55]. Following Uchiyama [53], two effects contribute to this distortion, which acts along [111] (sketched in Fig. 4c). The first is a distortion caused directly by the nearest-neighbor superexchange [56]. The second originates from Coulomb forces induced by an asymmetric charge distribution around the ions, due to AF-induced band folding [52,53], predicted to be the dominant contribution [53].
This suggests two possibilities for perturbing the distortion. The first is that the excitation weakens the AF order, and therefore also the charge asymmetry, both of which then weaken the distortion. The lattice response time ( 1 ≈ 0.31 ps) should then reflect magnetic excitation times. Reported opticallyexcited magnon data indeed suggest similar times scales [57][58][59]. However, a magnetic diffraction experiment has disputed these results [60]. The second possibility is that the charge asymmetry is directly perturbed, without involving magnetism. An excitation above the gap could cause a displacive phonon excitation. However, time scales associated with electronic excitations in NiO are much shorter than 1 [61], and displacive phonon excitations in similar binary oxides occurs on much shorter time scales [62,63], so we conclude that this possibility is unlikely.
This brings about the following scenario. Electrons are optically excited above the gap and magnetic order rapidly weakens. This reduces the asymmetric charge distribution, triggering lattice motions that weaken the RLD. Due to this preferential electron-phonon coupling, modes that lead to a reduction of the RLD are preferentially occupied, which is observed as disproportionally higher growth of 〈 2 〉. The lattice subsequently reaches an elevated thermal state via phonon-phonon coupling within ~25 ps. This scenario should directly affect other observables, such as the phononic energy gap or the rhombohedral angle itself. To support this scenario, we consider how varying the RLD would affect the observed scattering vectors, i.e. the positions of the spots (Fig. 1b). We divide them into Friedel pairs of the form (ℎ 0) and (ℎ ̅ ̅ 0), and extract the change in peak position ( ) = ( ) − 0 for each spot individually, as well as the average of each pair. The average curves of all pairs closely reproduce each other (black symbols in Fig. 4c), and exhibit a slow decrease of Q indicating isotropic lattice expansion (i.e. a = b at all delays).
Combining with (Fig. 4b) produces an expansion coefficient of ~10 -5 K -1 , in agreement with literature [65]. However, the individual spots in every pair deviate symmetrically around this mean (shown for the (ℎ00) family in Fig. 4c). To quantify this, we introduce the asymmetry Λ = 2 (ℎ 0) − 2 (ℎ ̅ ̅ 0) (accounting for Ewald's sphere curvature, see supplement). A nonzero Λ represents deviations from an orthonormal unit cell, and scales as Λ ∝ (ℎ + ) (see supplement). The inset presents Λ( ) All intensities and peak positions in our diffraction experiments were extracted from best fits to twodimensional peak functions. Figure S1 -schematic of the sample preparation process: NiO is evaporated on NaCl; the NaCl is dissolved using water, and the NiO film is placed on a Cu TEM Grid. Top: microscope image of the sample.

Raman Scattering
Raman spectroscopy measurements were performed to assess whether the sample's phonon dispersion agrees with that of bulk NiO. Measurements were conducted in the backscattering configuration with an optical excitation wavelength of 405 nm. The incident laser light was focused by a microscope objective onto the sample surface (spot diameter ~2 µm). The backscattered light was collected by the same objective (confocal configuration), spectrally dispersed by an 80 cm spectrograph (LabRam HR Evolution, Horiba/Jobin Yvon) and detected with a liquid nitrogen cooled CCD. For Rayleigh light suppression, a band pass filter with ultra-narrow spectral bandwidth was used. Figure S2 presents spectra taken at room temperature, with vertical lines indicating positions of Raman peaks due to second-order phonon scattering observed for bulk NiO. Figure S2 -Raman spectra collected from the sample at room temperature. Two-phonon features known from bulk NiO are indicated.
Here we detail the laser excitation. The probe beam was imaged on the sample, and its size was estimated using a 2000 lines/inch grid. We assume its profile to be a flat top ellipsoid of full widths (140 ± 13 × 114 ± 13) 2 , and we label its perimeter as for convenience.
The pump beam was imaged at a point equivalent to the sample position. Its best fit to a two-dimensional Pseudo-Voigt profile has a full width at half maximum (FWHM) of (400 ± 20 × 245 ± 12) 2 .
The pump energy was = 7.6 ± 0.2 per pulse, at a repetition rate of 4 kHz.
To estimate the incident fluence, we consider the fraction of this energy that overlaps with the probe pulse in an approach similar to the one presented by Harb et al. [1]. This is done by taking the ratio of two integrals over the pump profile: once only up to the limits of the probe pulse ( ), and once over all space. The incident fleunce is then given by: Here is the probe spot area. R is the reflectivity of the sample at the pump photon energy (2.16 eV), taken as R=0.16 [2]. Using Eq. (S1) we reach the final value of = 5.2 ± 1.3 mJ cm -2 .
Next we estimate the linear and two-photon absorption. The linear absorption coefficient at the pump photon energy is 260 cm -1 [2]. Considering , this produces a negligibly small absorbed fluence of 1 ℎ = 0.003 mJ cm -2 (with ~20% error). For estimating two-photon absorption, we follow the steps described by Zeuschner et al. [3]. We make the same thin-layer approximation to reach an intensity of the form: Here z is the film thickness, 0 is the incident intensity, and is the two-photon absorption coefficient.
We also assume that the pump pulse is a Gaussian envelope with a FWHM of = 50 ± 10 fs. It then follows that the absorbed fluence from two-photon absorption is The available literature value for two-photon absorption at twice our pump energy is β = 0.12 cm/MW [4].
This produces 2 ℎ = 0.091 ± 0.047 mJ cm -2 . This is over 30 times larger than the linearly absorbed fluence, indicating that the excitation is dominated by two photon absorption. These values take into account repeated internal reflections of the pump within the sample, which add ~4% more fluence.
Importantly, we note that the value of was measured with 10 ns pulse lengths. As our pump pulse is ~50 fs, we regard 2 ℎ as a lower limit for nonlinear absorption of our pump pulse.

Use of tabulated values
In this work we used tabulated values for the scattering factors and for the Debye-Waller factors. Highenergy electron diffraction from NiO was studied in detail in Ref. [5]. The authors parametrized the scattering factors for Ni 2+ and O 2using a modified version of the Mott-Bethe formula of the form: Here the X-ray scattering form factors are replaced by a parametrized sum of Gaussians, and the momentum transfer s is related to our definition of Q by = /2. To Calculate ( ) and ( ), we used Δ = ±2, which is the most conservative choice when considering the element sensitivity provided by (i.e. this value produces the smallest variation in across the measured Q range). We note that even when the Mott-Bethe formula [6,7] is used with tabulated X-ray scattering factors, the results of our analysis remain nearly unchanged. A more complete account of this approach was given by Peng in Ref. [8].

Sample evolution during the measurement
In this section we discuss the evolution of the sample during the measurement, which could indicate degradation. The mean lattice temperature, ̅ , is presented in Fig. 4b of the main text. Its response to photoexcitation exhibits two processes: a sub-picosecond process and a second, longer process of a few picoseconds. The experimental data presented in this manuscript were accumulated over a period of ~24h. Figure S3a presents ̅ , as extracted at different stages of accumulation. Clearly the data spread improves with time. However, it is also apparent that the sub-picosecond process remains the same as time progresses (albeit better-resolved), while the slower process does not. The slower process is nearly absent at the beginning of the accumulation time, and grows in prominence as accumulation progresses.
From this we conclude that the sub-picosecond process is intrinsic to the system, while the second process is induced during the experiment, such as due to increased absorption from photodoping.
The Q-dependent behavior described in the main text is not affected by this. To demonstrate this, Figure   S3b presents the parameter ( ) (as in Fig. 4a of the main text) at the same acquisition times as in a. Figure S3 -(a) Mean lattice temperature and (b) the ( ) parameter, as functions of pump-probe delay. The data were extracted in the same way as the data in Fig. 4a and 4b of the main text. The different curves are generated from data taken after different acquisition times up to 24h (the same times are used in both panels). The data in (a) emphasize that the sub-picosecond process is always present, but the slower process evolves during the measurement. The data in (b) demonstrate that the Q-dependent behavior is present and does not qualitatively change.

Friedel pair analysis
In this section we present background regarding the asymmetry analysis and deformation of the unit cell.
In a rhombohedral system the size of a scattering vector reads: 2 (ℎ ) = (ℎ 2 + 2 + 2 ) sin 2 + 2(ℎ + + ℎ )(cos 2 − cos ) 2 (1 − 3 cos 2 + 2 cos 3 ) (S6) From inspection we find that = 90° recovers the cubic functional dependence. Asymmetry is defined in the main text as Inspecting Eq. (S6), we conclude that Λ is always zero. However, an experiment probes a cut of reciprocal space given by Ewald's sphere, which is not flat. Therefore, reflections are not cut precisely through their center, but rather near their center, providing an observed maximum. This scenario is depicted in Figure   S4 for a Friedel pair of the form (ℎ 0). In this depiction it is clear that both reflections in the pair acquire a small out of plane component of the same sign. Therefore, an assumed Λ(ℎ 0) is actually Λ(ℎ | |), which can be expressed as Λ(ℎ 0) = 4 ( + ℎ) cos 2 (cos − 1)(1 + 2 cos ) ≡ −2 (ℎ + ) ( ) (S8) When plugging in = 0°, we find Λ(ℎ 0) = 0, so even if exists (experimentally it always does), a nonzero Λ can only occur if the system is not orthonormal. This description is valid if = , as in the cubic case. In our experiment we do not identify any deviation from this. Figure S4 -sketch of the experimental observation of Friedel pair with = 0. The curvature of Ewald's sphere causes both reflections to acquire a small out-of-plane component = .
In the following we present examples of Λ. Importantly, we note that asymmetry will occur in our data also due to experimental artifacts related to the experimental geometry (e.g. sample tilts). These effects would constitute a constant baseline for Λ. To avoid this, we consider only the photoinduced changes to Λ, and not the absolute quantity, by subtracting the baseline. The top panels of Figure S5 present Λ( ) for six Friedel pairs along the two primary axes: ±(ℎ 0 0) and ±(0 0). Since the data here are based on the positions of Bragg peaks on our 2D detector, we avoid any possible artifacts associated with unit conversion by presenting in units of detector pixels instead of Å -1 . We find that dynamics observed for all six pairs along a given axis are similar, but are larger in magnitude as ℎ or grows.
To understand this, we reexamine Eq. (S8). Since = = , photoinduced dynamics of Λ occur only due to , the rhombohedral angle. Furthermore, the photoinduced change of Λ should scale in magnitude with (ℎ + ). To test this, the lower panels of Figure S5 present the same data, divided by ℎ or . We find that this causes the Λ( ) curves to fall onto each other. This demonstrates that the photoinduced changes of Λ that we observe are in good agreement with Eq. (S8), pointing to a change in the unit cell's shape.
We note that while the changes originate from the rhombohedral angle , the change in may not be continuous. Within the observed time frame it could occur that ≠ ≠ , in which case ( ) in Eq. (S8) would be ( , , ).