Mapping momentum-dependent electron-phonon coupling and non-equilibrium phonon dynamics with ultrafast electron diffuse scattering

Despite their fundamental role in determining material properties, detailed momentum-dependent information on the strength of electron-phonon and phonon-phonon coupling (EPC and PPC, respectively) across the entire Brillouin zone (BZ) has proved difficult to obtain. Here we demonstrate that ultrafast electron diffuse scattering (UEDS) directly provides such information. By exploiting symmetry-based selection rules and time-resolution, scattering from different phonon branches can be distinguished even without energy resolution. Using graphite as a model system, we show that UEDS patterns map the relative EPC and PPC strength through their profound sensitivity to photoinduced changes in phonon populations. We measure strong EPC to the $K$-point transverse optical phonon of $A_1'$ symmetry ($K-A_1'$) and along the entire longitudinal optical branch between $\Gamma-K$, not only to the $\Gamma-E_{2g}$ phonon as previously emphasized. We also determine that the subsequent phonon relaxation pathway involves three stages; decay via several identifiable channels to transverse acoustic (TA) and longitudinal acoustic (LA) phonons (1-2 ps), intraband thermalization of the non-equilibrium TA/LA phonon populations (30-40 ps) and interband relaxation of the LA/TA modes (115 ps). Combining UEDS with ultrafast angle-resolved photoelectron spectroscopy will yield a complete picture of the dynamics within and between electron and phonon subsystems, helping to unravel complex phases in which the intertwined nature of these systems have a strong influence on emergent properties.

The nature of the couplings within and between lattice and charge degrees of freedom is a central concern of condensed matter and materials physics. Electronphonon interactions play a dominant role in the electronic transport properties of metals [4]; they are the underlying cause of conventional superconductivity [5] and Peierls/Jahn-Teller instabilities [6]. Furthermore, they are central to our understanding of the properties of many quasiparticles including polarons [7] and phonon-polaritons [8]. Highly anisotropic (momentumdependent) EPC has been identified as a key feature of superconductivity in MgB 2 [9]. It is also intertwined with electron correlations in the iron-based superconductor FeSe [10] and has been shown to contribute to the selection of the electronic ordering vector in some charge density wave materials including ErTe 3 [11] and NbSe 2 [12]. On the other hand, PPC dictates the thermalization properties of carrier/quasiparticle excitation energy.
The subtle details of charge-lattice interactions can have an enormous impact on technologically-relevant material properties, but these interactions have so far resisted a comprehensive experimental investigation. Conventional ARPES [13], inelastic x-ray/neutron scattering [1] and Raman spectroscopy [14] provide indirect information on the EPC strength through the shifting and broadening of spectral features only over a limited part of the Brillouin zone (BZ) [1]. Detailed studies of phonon-phonon interactions and decay in materials has typically been the province of theory [15] or molecular dynamics simulations due to a lack of techniques capable of probing these interactions in any substantial detail. Time-domain approaches have recently opened new windows on the nature of EPC and PPC in materials. Time-resolved Raman spectroscopy has been used to directly measure the rate of energy exchange between photo-generated carriers and zone-center optical phonons and the subsequent relaxation of those specific nonequilibrium phonons [14,16,17]. Time-resolved inelastic x-ray scattering at synchrotron and x-ray free electron laser facilities has provided a view of non-equilibrium distributions of off zone-center phonons in InP and InSb [18] and the phonon band structure in Ge [12,19] through incoherent and coherent time-resolved diffuse scattering signals respectively.
Here we demonstrate that ultrafast electron diffuse scattering (UEDS) using radio-frequency compressed electron pulses provides a general, lab-scale, timeresolved analog of diffuse x-ray scattering [20]. This new method [21] is capable of directly determining both the relative momentum-dependent interaction strength between photogenerated carriers and phonons and the subsequent phonon-phonon interactions governing the relaxation and thermalization of the excitation energy across the entire BZ with ∼ 100 fs time-resolution [22]. In pump-probe geometry [23,24] these experiments map the transient changes to diffuse (inelastic) electron scattering patterns, which are themselves determined by the evolution of the non-equilibrium phonon distributions that follow electronic excitation. We show that this technique is particularly well suited as a probe of 2D materials using thin graphite as a model system. Specifically, optical excitation at 800 nm with 35 fs laser pulses drives vertical electronic transitions of a π − π character on the well-known Dirac cones [25] of single crystal graphite samples ( Fig. 1 a). This excitation impulsively photodopes the material with a non-equilibrium electron-hole plasma with carrier density controllable by excitation fluence. TR-ARPES experiments have found that the first stage of relaxation is for the non-equilibrium distribution of carriers to thermalize internally through carrier-carrier scattering, forming a Fermi-Dirac distribution with well defined electron temperature within ∼50 fs [17,26,27]. In this work we use UEDS patterns to determine how the energy stored in the hot electron system couples to the phonons and how the phonon system subsequently thermalizes, comparing the results-where possible-with earlier investigations using time-resolved Raman spectroscopy [14,17], pump-probe spectroscopy [28,29] and theory [15].
The photoinduced changes to the ultrafast electron scattering pattern, ∆I(q, τ ) = I(q,τ )−I(q,−∞) , are shown for several points in time following photo-excitation in Fig. 2. The evolution of these patterns from 0.5 to 100 ps is striking and encodes detailed information on changes in the phonon system. The diffuse scattering intensity at q is modulated according to population dynamics of phonon modes with momentum k [20]: The first sum is taken over j phonon branches and the second is taken over s atoms in the unit cell. f s is the atomic scattering factor, µ s is the atomic mass, and M s is the Debye-Waller factor. Also, ω j,k andê j,k are the momentum-dependent phonon frequency and polarization for branch j. Finally, n j,k is the population of the phonon mode with frequency ω j,k . The UEDS patterns ( Fig. 2) provide an ultrafast snapshot of the change in I(q) for all q proportional to the instantaneous occupancy of the ω j,k modes, multiplied by the norm of the one-phonon structure factor, F j (q). In general, the occupancy of a specific branch, n j,k , is not directly available from the inelastic scattering signal at q, since all modes j contribute to I(q); diffuse scattering is momentum resolved but energy-integrated. However, the complete phonon band structure of Si has been determined by combining modeling and thermal diffuse x-ray scattering data [20]. In addition, one can incorporate symmetryimposed inelastic scattering selection rules [30,31] which describe extinctions in F j (q) at particular points of the diffuse scattering pattern for a given phonon mode. The extinctions depend on q and the symmetry of the reduced wavevector k. If the symmetry group of the scattering vector is strict, several of the phonon modes can be inactive (F j (q) = 0), reducing the number of phonon branches that can contribute at that point. Thus, graphite and other high-symmetry 2D materials are excellent candidates for UEDS. These selection rules have been used previously in inelastic x-ray diffraction experiments [32] to measure the energy dependence of individual branches in the phonon dispersion of graphite. Here, we use the same selection rules to separate the population dynamics of individual phonon branches without energy resolution. An attractive feature of UEDS is that a discrete, strongly-coupled mode yields a peak in the differential scattering maps at the associated BZ momentum position of that phonon at short delay times; electronic excitation energy initially flows preferentially to modes with strong EPC and are the first to show an increase in diffuse scattering (Fig. 2 c). Here, diffuse scattering peaks (FWHM = 0.12Å −1 ) appear at the K−points along the reflection axes, near the {210} family of peaks (where scattering from the K − A 1 mode is allowed) and along star-like ridges joining Γ {200} − K (where LO branch scattering is allowed). The phonon dispersion relation of graphite (Fig. 1 b) shows strong softening of the LO branch near Γ and the TO branch near K due to Kohn anomalies [3,33] (Fig. 1). Earlier work suggested these strongly-coupled modes as the initial reservoir into which the electronic excitation energy flows [29] and our results confirm that hypothesis. Time-resolved Raman has previously been employed to follow the the occupancy of the zone center Γ − E 2g mode showing that it is indeed strongly-coupled [14,17]. Evidence for strong coupling to the off-zone-center K −A 1 mode has previously only been indirect and the peaks in Fig. 2 c) are the first direct observation. In addition, Fig. 2 c) indicates that coupling is strong for the entire LO branch between Γ − K not only for Γ − E 2g mode. The character of the differential diffuse scattering pattern changes dramatically through Fig. 2 c-f) as the non-equilibrium phonon distribution evolves, demonstrating their profound sensitivity to the details of the phonon occupancies. The complete timedependence of the diffuse intensity at selected points is shown in Fig. 3. Scattering from the K − A 1 mode is forbidden by symmetry at the K−points immediately proximate to the reflection axes (indicated by a green in the legend); only LO phonon scattering is observed at these points [32]. Thus, ∆I(q, τ ) at this point shows a qualita-tively distinct time-dependence (Fig. 3 b), green) versus K−points along the reflection axes at which scattering from the K − A 1 mode is allowed (Fig. 3 b), red). This includes a much slower initial rise; 730 fs (K−LO) compared to 280 fs (K−TO). Intensity near {200} (Fig. 3 b), cyan) reports on the occupancy of the strongly-coupled Γ − E 2g LO mode at early times, and exhibits a slower rise (430 fs) than the K − A 1 mode. For comparison the Γ − E 2g phonon population determined using TR-Raman [14] is shown in grey. The Raman curve is nearly identical to that shown for the K −A 1 phonon in terms of rise time and recovery, but differs from the data shown in cyan. The slower rise time observed here is likely due to the higher excitation conditions used (12 mJ/cm 2 compared to 0.2 mJ/cm 2 ) which is known to weaken the Kohn anomaly and EPC at the Γ point [14,28]. The different behavior of the cyan curve on picosecond timescales is due to scattering from the low-frequency LA and TA phonon modes involved in the dominant decay channel of the K − A 1 phonon (Fig. 1 b) and discussed  ) is almost an order of magnitude slower than that associated with the TO-K −A 1 phonon. The slow timescale decay evident in the Bragg peak [23] and reported in earlier ARPES measurements [13] does not emerge from the dynamics of any single mode, but is a composite of the decay in population of the strongly coupled optical modes (e.g. red, 1.7 ps) and the increase in population of all other modes. c) Diffuse intensity dynamics for points along the K − Γ line, and Γ − M line (inset). Rise time (blue) and amplitude (red) from single exponential fits to the early time dynamics are shown. Error bars represent covariance in fit parameters. Points close to Γ are not shown due to the interference of Bragg peak Debye-Waller dynamics of a).
further below. These LA/TA modes are not seen in the TR-Raman study, nor do they overlap with the signals measured at the K−point.
We can estimate the effective temperature of the K−point TO mode, T T O,K , using the measured increase in diffuse intensity and applying Bose-Einstein statistics to the mode population, n j,k = cosh ( ω j,k /2k B T j,k ). We determine that the effective temperature of the K−point TO mode is ∼ 1300 K by 1 ps and also that by 10 ps it has cooled back down to ∼ 450 K. The stronglycoupled TO and LO modes reach a pre-equilibrium with the laser-generated carriers in < 1 ps, while all other phonon modes remain at or near room temperature.
The diffuse scattering patterns in Fig. 2 d-f) reveal the decay channels for the generated population of stronglycoupled optical phonons as they relax. The time-scales separation between the EPC into K − A 1 and Γ − E 2g (200 − 400 fs) and the subsequent decay out of these modes (1 − 3 ps) means that the diffuse scattering pattern at 1.5 ps maps their momentum-dependent decay probability in a manner analogous to the way in which the 0.5 ps pattern indicates the relative EPC strength. The probability for the Γ − E 2g phonon to decay to two phonons of momentum k and −k was previously computed using density functional perturbation theory [15] for the distinct momentum and energy conserving channels indicated with arrows in Fig. 1 b). The results of these calculations compare well with ∆I(q, τ = 1.5 ps) in Fig. 2 d) (inset). The dominant decay channels are presented in Fig. 1 b). The hexagonal distribution of diffuse intensity approximately halfway between Γ and the BZ edges is associated with the LA-TA and LA-LA decay channels. Peaks in evident along the Γ − K lines are also in agreement with the expected location of maximum decay rate. Diffuse intensity around the K−points are associated with the TA-TA decay channel. The time-constants associated with each of these channels can be determined directly from the time-dependence of ∆I(q, τ ) at the associated k position in the BZ, shown for selected positions in Fig. 3 b). The variation in the rise time and amplitude of ∆I(q, τ ) in the BZ surrounding the {210} peak is given along the Γ − K and Γ − M directions in Fig. 3 c). The dominant TA-LA and TA-TA decay channels in the k = 0.5 K and k = 0.5 M regions have the fastest time-constants (1.0 − 1.5 ps) and largest amplitude.
The decay of the K − A 1 phonon population is evident in the relaxation of the K−TO intensity shown in Fig. 3 b) (red, 1.7 ps), is dominated by LO-LA and LA-TA channels. The LA-TA decay overlaps with that of the Γ − E 2g phonon and is not easily separated, however, the LO-LA and LA-TA channels yield an increase in intensity adjacent to Γ and K points that can be readily identified. Fig. 4 shows the time dependence of the diffuse intensity adjacent to the {200} and {210} peaks along the Γ − K direction. Adjacent to {210}, scattering from Γ − E 2g is forbidden by symmetry; thus, the fastest time-constant evident in the data, 1.5 ps (Fig. 4 (grey)), can only be assigned to the LO-LA and LA-TA decay channels of the K − A 1 phonon. Adjacent to {200}, scattering from Γ − E 2g is allowed and the early time dynamics appear approximately bi-exponential ( Fig. 4 inset (cyan)). The observed fast time-constant behaviour is associated with the strong EPC to the Γ − E 2g mode, as previously described. The slower dynamics is assigned to a composite time-scale resulting from a decrease in intensity due to the decay of Γ − E 2g phonons and an increase due to the LO-LA and LA-TA decay channels of the K −A 1 phonon. It is clear that time-resolution in UEDS can substitute for energy resolution when discriminating contributions from various phonon branches at a single momentum. On longer time-scales the dynamics evident in the phonon system remain rich. The early-time diffuse intensity bands oriented along Γ − K directions (Fig. 2 c) and d)) evolve to become intensity bands oriented along Γ − M directions. The two principle factors in these dynamics are that the decay of the strongly-coupled optical modes generates a profoundly non-equilibrium population of LA and TA phonons (Fig. 1 b), and that these branches are considerably softer along the Γ − M direction than the Γ − K in graphite. The non-equilibrium populations of LA and TA phonons thermalize within these bands by emission of lower frequency LA and TA phonons. This relaxation channel is evident in the increase in the diffuse halos surrounding the Bragg peaks in Fig. 2 e) and f) (inset), whose detailed time dependence is plotted in Fig. 4. This relaxation of the phonon distribution in the acoustic branches occurs on a time-scale of 34 ps and is a phonon analog of the carrier thermal-ization that has previously been observed in TR-ARPES experiments [17,26,27]. The lowest frequency branch of the in-plane phonon band structure consist of acoustic modes with a polarization normal to the graphene planes (ZA modes), to which the current experiments are insensitive due to electron beam orientation. Relaxation into the ZA branch and LA/TA interband relaxation is evident, however, as the decay of the diffuse halo intensity surrounding the {210} peak (Fig. 4) which occurs on the 115 ps time-scale. All of these time-scales are much faster than thermal transport of the laser deposited energy out of the probed volume, which occurs on the 10 µs timescale in the geometry of these experiments.
We have demonstrated that UEDS provides direct, momentum resolved measurements of the relative strength of EPC and PPC in graphite through the technique's ability to follow phonon population dynamics with femtosecond time resolution. Particularly notable is the unprecedented picture of phonon decay kinetics that the UEDS patterns provide, only a summary of which has been given here. UEDS is also profoundly complementary to ultrafast ARPES. Together these methods can provide a complete picture of the dynamics within and between electron and phonon sub-systems, and help unravel the physics of complex phases where the intertwined nature the electron-lattice systems determine material properties.