Time- and momentum-resolved phonon population dynamics with ultrafast electron diffuse scattering

TL;DR

This paper develops a quantitative framework to analyze ultrafast electron diffuse scattering data, enabling the extraction of phonon populations and coupling constants across the Brillouin zone, demonstrated on graphite.

Contribution

It introduces a detailed method for interpreting UEDS signals to determine phonon mode occupancies and coupling constants, advancing the understanding of lattice-charge interactions.

Findings

Extracted momentum- and mode-dependent electron-phonon and phonon-phonon coupling constants.

Determined the A1' phonon mode-projected electron-phonon coupling strength in graphite.

Validated the method by comparing results with other techniques and simulations.

Abstract

Interactions between the lattice and charge carriers can drive the formation of phases and ordering phenomena that give rise to conventional superconductivity, insulator-to-metal transitions, and charge-density waves. These couplings also play a determining role in properties that include electric and thermal conductivity. Ultrafast electron diffuse scattering (UEDS) has recently become a viable laboratory-scale tool to track energy flow into and within the lattice system across the entire Brillouin zone, and deconvolves interactions in the time domain. Here, we present a detailed quantitative framework for the interpretation of UEDS signals, ultimately extracting the phonon mode occupancies across the entire Brillouin zone. These transient populations are then used to extract momentum- and mode-dependent electron-phonon and phonon-phonon coupling constants. Results of this analysis are presented for graphite, which provides complete information on the phonon-branch occupations and a determination of the $A_1'$ phonon mode-projected electron-phonon coupling strength $\langle g_{e,A_1'}^2 \rangle = 0.035 \pm 0.001$ eV$^2$ that is in agreement with other measurement techniques and simulations.

Tables

Table 1. Table 1: Coupling strength between electronic system, the A 1 ′ superscript subscript 𝐴 1 ′ A_{1}^{\prime} phonon, and the lattice system. Uncertainty is derived from fit covariances.

	Coupling strength [$W{m}^{-3}{K}^{-1}$]
$G_{e,A_1^{\prime }}$	$\left(6.8\pm 0.3\right)\times {10}^{17}$
$G_{A_1^{\prime },l}$	$\left(8.0\pm 0.5\right)\times {10}^{17}$
$G_{e,l}$	$\left(0.0\pm 6.0\right)\times {10}^{15}$

Equations

$${\mathrm{a}}_1=a{\mathrm{e}}_1,{\mathrm{a}}_2=\frac{a}{2}\left(\sqrt{3}{\mathrm{e}}_2-{\mathrm{e}}_1\right),{\mathrm{a}}_3=c{\mathrm{e}}_3,$$

$$I(\mathrm{q},t)=I_0(\mathrm{q},t)+I_1(\mathrm{q},t)+...,$$

$$I_0(\mathrm{q},t)$$

$$I_1(\mathrm{q},t)$$

$$|F_{1j}(\mathrm{q},t){|}^2={\left|\sum _s{e}^{-W_s(\mathrm{q},t)}\frac{f_s(\mathrm{q})}{\sqrt{{\mu }_s}}(\mathrm{q}\cdot {\mathrm{e}}_{j,s,\mathrm{k}})\right|}^2$$

$$W_s(\mathrm{q},t)=\frac{1}{4{\mu }_s}\sum _{j,\mathrm{k}}|a_{j,\mathrm{k}}(t){|}^2|\mathrm{q}\cdot {\mathrm{e}}_{j,s,\mathrm{k}}{|}^2$$

$$|a_{j,\mathrm{k}}(t){|}^2=\frac{2\hslash }{N_c{\omega }_{j,\mathrm{k}}(t)}\left(n_{j,\mathrm{k}}(t)+\frac{1}{2}\right).$$

$$\frac{\Delta I(\mathrm{q},t)}{N_c{I}_e}=\sum _j\frac{\Delta n_{j,\mathrm{k}}(t)}{{\omega }_{j,\mathrm{k}}(t_0)}|F_{1j}(\mathrm{q},t_0){|}^2$$

$${\mathrm{I}}_{\mathrm{k}}(t)={\mathrm{F}}_{\mathrm{k}}{\mathrm{n}}_{\mathrm{k}}(t)$$

$${\mathrm{I}}_{\mathrm{k}}(t)$$

$${\mathrm{n}}_{\mathrm{k}}(t)$$

$${\mathrm{F}}_{\mathrm{k}}$$

$$\left[\begin{array}{ccc}|F_{11}(\mathrm{k}+{\mathrm{H}}_1,t_0){|}^2& \dots & |F_{1N}(\mathrm{k}+{\mathrm{H}}_1,t_0){|}^2\\ ⋮& \ddots & ⋮\\ |F_{11}(\mathrm{k}+{\mathrm{H}}_M,t_0){|}^2& \dots & |F_{1N}(\mathrm{k}+{\mathrm{H}}_M,t_0){|}^2\end{array}\right]$$

$$C_e(T_e)\frac{\partial T_e}{\partial t}=\sum _i{G}_{ep,i}(T_e-T_{ph,i})+f(t-t_0)$$

$$\Bigg{\{}C_{ph,j}(T_{ph,j})\frac{\partial T_{ph,j}}{\partial t}=\\ \sum_{i\neq j}G_{ep,i}(T_{e}-T_{ph,i})+G_{pp,ij}(T_{ph,j}-T_{ph,i})\Bigg{\}}_{j=1}^{N}$$

$$n_{j,\mathrm{k}}(t)\propto {\left[\exp \left(\frac{\hslash {\omega }_{j,\mathrm{k}}}{k_B{T}_{ph,j}(t)}\right)-1\right]}^{-1}.$$

$$n_{j,\mathrm{k}}(t)\propto \frac{k_B{T}_{ph,j}(t)}{\hslash {\omega }_{j,\mathrm{k}}}-\frac{1}{2}+O\left(T_{ph,j}^{-1}(t)\right).$$

$$C_e(T_e)$$

$$+1.16\times 10^{-3}{T}_e^2$$

$$+2.6\times 10^{-7}{T}_e^3.$$

$$\begin{gathered} C_{ph,j}(T_{ph,j})= \\ {k}_B{\int }_0^{{\omega }_D}d\omega D(\omega ){\left(\frac{\hslash \omega }{k_B{T}_{ph,j}}\right)}^2\frac{e^{\hslash \omega /k_B{T}_{ph,j}}}{{\left(e^{\hslash \omega /k_B{T}_{ph,j}}-1\right)}^2} \end{gathered}$$

$$C_e(T_e)\frac{\partial T_e}{\partial t}$$

$$-G_{e,A_1^{\prime }}(T_e-T_{A_1^{\prime }})$$

$$-G_{e,l}(T_e-T_l)$$

$$C_{A_1^{\prime }}(T_{A_1^{\prime }})\frac{\partial T_{A_1^{\prime }}}{\partial t}$$

$$-G_{A_1^{\prime },l}(T_{A_1^{\prime }}-T_l)$$

$$C_l(T_l)\frac{\partial T_l}{\partial t}$$

$$+G_{A_1^{\prime },l}(T_{A_1^{\prime }}-T_l)$$

$${\mathrm{P}}_{j,\mathrm{k}}={\left[\begin{array}{cccc}{\omega }_{j,\mathrm{k}}& {\mathrm{e}}_{j,s=1,\mathrm{k}}& \dots & {\mathrm{e}}_{j,s=M,\mathrm{k}}\end{array}\right]}^T$$

$$\Vert {\mathrm{P}}_{i,\mathrm{k}}-{\mathrm{P}}_{j,{\mathrm{k}}^{\prime }}\Vert =|{\omega }_{i,\mathrm{k}}-{\omega }_{j,{\mathrm{k}}^{\prime }}{|}^2+\sum _s\Vert {\mathrm{e}}_{i,s,\mathrm{k}}-{\mathrm{e}}_{j,s,{\mathrm{k}}^{\prime }}\Vert .$$

$$\frac{\partial T_{e}}{\partial t}-\sum_{j}\frac{\partial T_{ph,j}}{\partial t}=\sum_{j}\Bigg{[}\frac{G_{ep,j}}{C_{e}}~{}(T_{e}-T_{ph,j})\\ -\sum_{i}\bigg{(}\frac{G_{ep,i}}{C_{ph,j}}~{}(T_{e}-T_{ph,i})+\frac{G_{pp,ij}}{C_{ph,j}}~{}(T_{ph,i}-T_{ph,j})\bigg{)}\Bigg{]}$$

$$\frac{\partial T_e}{\partial t}-\sum _j\frac{\partial T_{ph,j}}{\partial t}$$

$$-\sum _{i,j}\frac{G_{ep,i}}{C_{ph,j}}(T_e-T_{ph,i}).$$

$$\dot{\lambda }(t)-a(t)\lambda (t)=0$$

$$a(t)=\sum _j\left(\frac{G_{ep,j}}{C_e}-\sum _i\frac{G_{ep,i}}{C_{ph,j}}\right).$$

$$\lambda (t)=\exp \int dt\left[a(t)\right].$$

$$\frac{1}{\tau }\approx \sum _j\left(\frac{G_{ep,j}}{C_e}-\sum _i\frac{G_{ep,i}}{C_{ph,j}}\right).$$

Full text

I Introduction

Elementary excitations and their mutual couplings form the fundamental basis of our understanding of diverse phenomena in condensed matter systems. The interactions between collective excitations of the lattice system (phonons) and charge carriers, specifically, are known to lead to superconductivity, charge-density waves, multiferroicity, and soft-mode phase transitions Eliashberg (1960); Hur et al. (2004). These carrier-phonon interactions are also central to our understanding of electrical transport, heat transport, and energy conversion processes in photovoltaics and thermoelectrics Zhao et al. (2016). Phonons can themselves be intimately mixed in to the very nature of more complex elementary excitations, as they are in polarons or polaritons, or intertwined with electronic, spin, or orbital degrees-of-freedom, as it now seems is the case for the emergent phases of many strongly-correlated systems that exhibit complex phase diagrams like high-$T_{c}$ superconductors Miyata et al. (2017); Kim et al. (2012); Lanzara et al. (2001).

Our inability to fully characterize the nature of elementary excitations and to quantify the strength of their momentum-dependent interactions has been one of the primary barriers to our understanding of these phenomena, particularly in complex anisotropic materials. Ultrafast pump-probe techniques provide an opportunity to study couplings between elementary excitations rather directly. Photoexcitation can prepare a non-equilibrium distribution of quasiparticles or other selected modes whose subsequent relaxation dynamics and coupling to other degrees of freedom can be followed in the time-domain. Under favourable circumstances, qualitatively distinct channels can be disentangled by their associated spectra (response functions) and time-scales. This field has evolved rapidly over the last decade, both from the perspective of the selectivity of the initial excitation and the ability to probe the subsequent dynamics over a broad range of frequencies. For example, spectroscopic pump-probe techniques in the terahertz range have been used to interrogate the link between electrons/holes and optical phonons in hybrid lead halide perovskites Lan et al. (2019) and to investigate the time-ordering of phenomena behind charge-density waves in titanium diselenide Porer et al. (2014).

The low-photon momentum associated with optical frequencies, however, prevents the most commonly applied optical photon-in, optical photon-out techniques from providing a full characterization of the wavevector-dependent interactions between elementary excitation. Time-resolved Raman and Brillouin spectroscopies, for example, are limited to the interrogation of zero-momentum (zone-center) phonons for this reason Tsen and Ferry (2009); Yan et al. (2009); Yang et al. (2017).

In recent years, non-optical ultrafast techniques have been developed to probe wavevector-dependent dynamics. The most mature of these approaches is time- and angle-resolved photoemission spectroscopy (trARPES), which has been used to assemble a complete picture of the dynamics of the electronic and spin excitations following the photoexcitation of materials Johannsen et al. (2013); Gierz et al. (2015); Stange et al. (2015); Yang et al. (2017); Rohde et al. (2018).

The ability to directly interrogate wavevector-dependent dynamics within the phonon system, on the other hand, is an extremely recent development. Ultrafast X-ray diffuse scattering Trigo et al. (2010); Zhu et al. (2015); Wall et al. (2018) is one technique that has the potential to reveal lattice excitation dynamics across the whole Brillouin zone. This approach leverages the remarkable brightness of the beams available from the current generation of X-ray free-electron laser facilities to measure the time dependence of the diffuse (phonon) scattering from materials following photoexcitation.

At the laboratory scale, there has been similar progress made in furthering ultrafast electron beam brightness which has enabled equivalent diffuse scattering experiments to be performed. Ultrafast electron diffuse scattering (UEDS) has the potential to be transformative in that it can provide both a wavevector-resolved view of the coupling between electron and lattice systems Harb et al. (2016); Chase et al. (2016); Waldecker et al. (2017); Stern et al. (2018); Konstantinova et al. (2018) and the wavevector dependence of the interactions within the phonon system itself. The large scattering cross-section of electrons, combined with the relative flatness of the Ewald sphere, potentially allows for the simultaneous measurement of both the average lattice structure (via Bragg scattering) and lattice excitation dynamics (via diffuse scattering) in specimens as thin as a single atomic layer.

In this work we provide a description of the signals contained in UEDS measurements and a comprehensive and broadly applicable computational method for UEDS data reduction based on density functional perturbation theory (DFPT). Specifically, we present a procedure to recover phonon population dynamics as a function of the phonon branch and wavevector, and a determination of wavevector-dependent (or mode-projected) electron-phonon coupling constants from those phonon population measurements. This method uses only the measured time-resolved UEDS patterns and DFPT determinations of the phonon polarization vectors as inputs. The application of this approach to the case of photodoped carriers in the Dirac cones of thin graphite is demonstrated. The electron-phonon coupling strength to the strongly-coupled $A_{1}^{\prime}$ phonon at the $\bm{\mathrm{K}}$-point of the Brillouin zone, and the nonequilibrium optical and acoustic phonon branch populations as a function of time following excitation across the whole Brillouin zone are all determined from the UEDS measurements.

II Experimental and Computational methods

The change in experimental electron scattering intensity of graphite, photoexcited with $35\text{\,}\mathrm{fs}$ pulses of $800\text{\,}\mathrm{nm}$ light at a fluence of $12\text{\,}\mathrm{mJ}\text{\,}{\mathrm{cm}}^{-2}$, are presented in Figure 2 for a few representative time-delays. This section provides details on the experimental parameters, data processing steps, and computational techniques that are used in this work.

II.1 Data acquisition

UEDS measurements are pump-probe experiments in which an ultrafast laser pulse is used to photoexcite a thin single-crystal specimen at $t=t_{0}$, followed by probing the specimen with an ultrafast electron pulse at $t=t_{0}+\tau$, resulting in the acquisition of a transmission electron scattering pattern. By scanning across time-delays $\tau$, the dynamics in the $\mathrm{fs}$ – $\mathrm{ns}$ ($10^{-15}\textrm{--}10^{-9}$ $\mathrm{s}$) range can be recorded. UEDS data can be acquired coincidentally during ultrafast electron diffraction (UED) experiments with state-of-the-art detection cameras, although UEDS intensities are empirically $10^{4}$ to $10^{6}$ times less intense than those of Bragg diffraction. UEDS measurements are inherently statistical. Pump-probe experiments sample many decay processes, incoherently in time. Hence, all possible decay channels are represented, proportionally to their statistical likelihood.

Scattering measurements presented in this work use bunches of $10^{7}$ electrons, accelerated to $90\text{\,}\mathrm{keV}$, at a repetition rate of $1\text{\,}\mathrm{kHz}$. A radio-frequency cavity is used to compress electron bunches to $\approx$150\text{,}\mathrm{fs} at the sample, as measured by a home-built photoactivated streak camera Kassier *et al.* ([2010](#bib.bib24)). More detailed descriptions of this instrument are given elsewhere Chatelain *et al.* ([2012](#bib.bib25)); Morrison *et al.* ([2013](#bib.bib26)); Otto *et al.* ([2017](#bib.bib27)). Analysis of static diffraction patterns indicate a momentum resolution of $0.06\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$, while the range of visible reflections is consistent with a real-space resolution of $<$1\text{\,}\mathrm{\SIUnitSymbolAngstrom}. $35\text{\,}\mathrm{fs}$ pump laser pulses of $800\text{\,}\mathrm{nm}$ ($1.55\text{\,}\mathrm{eV}$) light are used to photoexcite a single-crystal flake of freestanding single crystal natural graphite, provided by Naturally Graphite. The flakes were mechanically exfoliated to a thickness of $70\text{\,}\mathrm{nm}$.

The interrogated film area covers $500\text{\times}500\text{\,}\mathrm{\SIUnitSymbolMicro m}$, with a pump spot of $1\text{\times}1\text{\,}\mathrm{mm}$ full-width half-max (FWHM) ensuring nearly uniform illumination of the probed volume. The film was pumped at a fluence of $12\text{\,}\mathrm{mJ}\text{\,}{\mathrm{cm}}^{-2}$, resulting in an absorbed energy density of $8\text{\,}\mathrm{J}\text{\,}{\mathrm{m}}^{-3}$. The scattering patterns are collected with a Gatan Ultrascan 895 camera; a $2.54\text{\times}2.54\text{\,}\mathrm{cm}$ phosphor screen fiber-coupled to a $2048\text{\times}2048\text{\,}\mathrm{p}\mathrm{x}$ charge-coupled detector (CCD) placed $25\text{\,}\mathrm{cm}$ away from the sample. The experiment herein consists of time-delays in the range of $-40\text{\,}\mathrm{ps}680\text{\,}\mathrm{ps}$. Per-pixel scattering intensity fluctuations over laboratory time reveals a transient dynamic range of $1:10^{8}$, allowing the acquisition of diffraction patterns and diffuse scattering patterns simultaneously 111The intensity fluctuations of pixel values across scattering patterns acquired before photoexcitation are $10^{8}$ times smaller than the brightest Bragg reflection.

II.2 Processing and corrections

Scattering intensity patterns are inherently redundant due to the point-group symmetry of the scattering crystal. When this symmetry is not broken by photoexcitation or the dynamical phenomena itself, it is possible to use this redundancy to enhance the signal to noise ratio of a UEDS data set. No dynamical phenomena breaking point-group symmetry was observable within the raw signal to noise of the current measurements, so the measured patterns have been subject to a six-fold discrete azimuthal average based on the $D_{6h}$ point-group of graphite. This discrete rotational average effectively increases the signal-to-noise ratio of this data set by a factor of $\sqrt{6}$ and is therefore an effective data processing step given the weakness of the diffuse scattering signals. “Scattering intensities” is henceforth implied to mean six-fold averaged scattering intensities.

Scattering from a few samples with varying thicknesses ($10\text{\,}\mathrm{nm}90\text{\,}\mathrm{nm}$) was acquired. There was no quantitative difference in the Bragg scattering dynamics, indicating that scattering from these samples is kinematical to a good approximation. The expected effects of multiple scattering on the diffuse scattering intensity distribution will be discussed further below.

II.3 Computational methods

Structure relaxation was performed using the plane-wave self-consistent field program PWscf from the Quantum ESPRESSO software suite Giannozzi et al. (2017). The graphite structure was fully relaxed using a $18\text{\times}18\text{\times}10$ $\bm{\mathrm{k}}$-point mesh centered at $\bm{\mathrm{\Gamma}}$ and force (energy) threshold of $1\text{\times}{10}^{-8}\text{\,}\mathrm{R}\mathrm{y}\mathrm{/}\mathrm{B}\mathrm{o}\mathrm{h}\mathrm{r}$ ($1\text{\times}{10}^{-14}\text{\,}\mathrm{R}\mathrm{y}$). The dynamical matrices were computed on $5\text{\times}3\text{\times}3$ $\bm{\mathrm{q}}$-point grid using a self-consistency threshold of $1\text{\times}{10}^{-18}\text{\,}\mathrm{R}\mathrm{y}$. The resulting graphite structure has the following lattice vectors:

[TABLE]

where $a=$2.462\text{,}\mathrm{\SIUnitSymbolAngstrom}, $c=$6.837\text{\,}\mathrm{\SIUnitSymbolAngstrom}, and $\left\{\bm{\mathrm{e}}_{i}\right\}$ are the usual Euclidean vectors. Graphite has four atoms in the unit cell, with two groups of two atoms forming lattices rotated with respect to each other. This structure respects the symmetries of the $D_{6h}$ point group 222The space group of this structure is $P6_{3}/mmc$ (Hermann-Mauguin symbol) or $D_{6h}^{4}$ (Schoenflies symbol)..

The phonon frequencies $\left\{\omega_{j,\bm{\mathrm{k}}}\right\}$ and polarization vectors $\left\{\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\right\}$ were computed using the PHonon program in the Quantum ESPRESSO software suite, using the B86b exchange-coupled Perdew-Burke-Ernzerhof (B86bPBE) generalized gradient approximation (GGA) and the projector augmented-wave (PAW) method Becke (1986); Perdew et al. (1996); Blöchl (1994). The cutoff-energy of the wavefunction was set to $100\text{\,}\mathrm{R}\mathrm{y}$, while the cutoff energy for the charge density was set to $1200\text{\,}\mathrm{R}\mathrm{y}$, and a Fermi-Dirac smearing of $0.06\text{\,}\mathrm{R}\mathrm{y}$ was applied. To include the dispersion energy between the carbon layers, the exchange-hole dipole moment (XDM) method was used Becke and Johnson (2007).

III Theory

Similar to X-ray scattering, under the kinematical approximation the measurement of the total scattering intensity at scattering vector $\bm{\mathrm{q}}$ and time $t$, $I(\bm{\mathrm{q}},t)$, of an electron bunch interacting with a thin film of crystalline material, can be separated as follows:

[TABLE]

where the intensity $I_{n}$ represents the scattered intensity of an electron that interacted with $n$ phonons. Specifically, $I_{0}$ represents diffraction, or Bragg scattering, and $I_{1}$ represents the first-order diffuse scattering intensity. The experimentally-observed ratio $I_{0}/I_{1}$ ranges between $10^{4}$ – $10^{6}$ 333Detector counts for the brightest Bragg peak reaches as much as 20 000 counts, while the average diffuse feature shown in Figure 2 is 0.2 counts.. Higher-order terms have much smaller cross-sections, hence much lower contribution to scattering intensity, and are therefore ignored in this work. The expressions for the intensities $I_{0}$ and $I_{1}$ are given below:

[TABLE]

where $N_{c}$ is the number of diffracting cells, $I_{e}$ is the intensity of scattering from a single event, $\bm{\mathrm{q}}$ is the wavevector (or scattering vector), $\bm{\mathrm{k}}=\bm{\mathrm{q}}-\bm{\mathrm{H}}$ is the reduced wavevector associated to $\bm{\mathrm{q}}$ with respect to the nearest Bragg reflection $\bm{\mathrm{H}}$ (see Figure 1), $s$ are indices associated with atoms in the crystal unit cell, $\bm{\mathrm{R}}_{s}(t)$ is the real-space atomic position of atom $s$, $W_{s}(\bm{\mathrm{q}},t)$ is the Debye-Waller factor of atom $s$, $f_{s}(\bm{\mathrm{q}})$ are the atomic form factors, $j\in\{1,2,...,N\}$ runs over phonon modes, and $\left\{n_{j,\bm{\mathrm{k}}}(t)\right\}$ and $\left\{\omega_{j,\bm{\mathrm{k}}}(t)\right\}$ are the population and frequency associated with phonon mode $j$, respectively Wang (1995); Xu and Chiang (2005).

The diffuse scattering intensity contribution of each mode $j$ is weighted by a factor, called the one-phonon structure factor $|F_{1j}(\bm{\mathrm{q}},t)|^{2}$:

[TABLE]

where $\mu_{s}$ is the mass of atom $s$, and $\left\{\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\right\}$ are the wavevector-dependent polarization vectors associated with phonon mode $j$ for atom $s$. The one-phonon structure factors represent the contribution of phonon mode $j$ on the overall intensity at a specific scattering vector $\bm{\mathrm{q}}$ and time $t$. $|F_{1j}(\bm{\mathrm{q}},t)|^{2}$ are a measure of two things: the locations in Brillouin zone where the phonon mode polarization vectors $\left\{\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\right\}$ are aligned in such a way that they will contribute to diffuse scattering intensity on the detector, expressed via the terms $\left\{\bm{\mathrm{q}}\cdot\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\right\}$; and the strength of the contribution of a single scattering event, including the effect of the instantaneous disorder in the material, expressed via the quantities $\left\{e^{-W_{s}(\bm{\mathrm{q}},t)}f_{s}(\bm{\mathrm{q}})/\sqrt{\mu_{s}}\right\}$.

Instantaneous disorder is described by the Debye-Waller factors $\left\{W_{s}(\bm{\mathrm{q}},t)\right\}$. The general expression of the anisotropic Debye-Waller factor for atom $s$, $W_{s}(\bm{\mathrm{q}},t)$, is given below:

[TABLE]

where $a_{j,\bm{\mathrm{k}}}(t)$ is the phonon mode vibration amplitude for mode $j$ at reduced wavevector $\bm{\mathrm{k}}$ Xu and Chiang (2005):

[TABLE]

Debye-Waller factors describe the reduction of intensity at scattering vector $\bm{\mathrm{q}}$ due to the effective deformation of a single atom’s scattering potential that results from the collective lattice vibrations in all phonon modes. Wavevector-specific information is in general impossible to extract from the transient changes to the Debye-Waller factors that result from photoexcitation and the non-equlibrium phonon populations that such excitation produces.

The expression for $I_{1}(\bm{\mathrm{q}},t)$ in Equation (2) (and related quantities) apply rigorously under single electron scattering conditions. The most probable type of multiple scattering event affecting $I_{1}(\bm{\mathrm{q}},t)$ is diffuse scattering followed by secondary Bragg scattering, not multiple or consecutive diffuse scattering events Cowley and Fields (1979); Wang (1995). A secondary Bragg scattering event only changes the electron wavevector by a reciprocal lattice vector; thus, this type of multiple scattering results in a redistribution of diffuse intensity from lower-order to higher-order Brillouin zones (further from $|\bm{\mathrm{q}}|=0$). However, the wavevector dependence of experimental diffuse intensities is not strongly influenced, even under experimental conditions where such dynamical effects are important. The strength of the scattering selection rules implied by the $\left\{\bm{\mathrm{q}}\cdot\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\right\}$ terms, and described further below, are reduced as the proportion of multiple scattering increases.

IV Results and Discussion

In this section a comprehensive procedure for ultrafast electron diffuse scattering data reduction will be presented. This approach recovers the time- and wavevector-dependent phonon population dynamics in each of the phonon branches. In addition, the population dynamics of the $A_{1}^{\prime}$ phonon, a strongly-coupled optical phonon in graphite, is used to demonstrate the extraction of wavevector-dependent (mode-projected) electron-phonon coupling constants from ultrafast electron diffuse intensities.

IV.1 Calculation of phonon polarization vectors and eigenfrequencies

The quantitative connection between the observed diffuse intensity $I_{1}(\bm{\mathrm{q}},t)$ and the phonon populations, $\left\{n_{j,\bm{\mathrm{k}}}\right\}$, is provided by the one-phonon structure factors, $|F_{1j}(\bm{\mathrm{q}},t)|^{2}$. A determination of $|F_{1j}(\bm{\mathrm{q}},t)|^{2}$ requires phonon polarization vectors $\left\{\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\right\}$ and associated frequencies $\left\{\omega_{j,\bm{\mathrm{k}}}\right\}$. Density functional perturbation theory (DFPT) is a widely used, readily-available method to compute these phonon properties.

The separation of frequencies and polarization vectors into modes is key to the calculation of one-phonon structure factors $|F_{1j}(\bm{\mathrm{q}},t)|^{2}$. The phonon frequencies $\left\{\omega_{j,\bm{\mathrm{k}}}\right\}$ and polarization vectors $\left\{\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\right\}$ were computed independently at every $\bm{\mathrm{k}}$; however, diagonalization routines have no way of clustering eigenvalues and eigenvectors into physically-relevant groups (i.e. phonon branches). Association between atomic motions (given by polarization vectors) and a particular mode is only well-defined near the $\bm{\mathrm{\Gamma}}$-point Paulatto et al. (2013). Clustering of phonon properties into phonon branches is described in detail in Appendix A.

The polarization vectors and frequencies, calculated for irreducible $\bm{\mathrm{k}}$-points 444The coverage of irreducible $\bm{\mathrm{k}}$-points is important. Only computing phonon properties along high-symmetry lines is fraught with peril, given that polarization vectors can vary significantly not only along high-symmetry lines, but over the entire Brillouin zone., were extended over the entire reciprocal space based on crystal symmetries using the crystals software package René de Cotret et al. (2018).

IV.2 Debye-Waller calculation

A key component of the computation of one-phonon structure factors $|F_{1j}(\bm{\mathrm{q}},t)|^{2}$ is the calculation of the Debye-Waller factors $\left\{W_{s}(\bm{\mathrm{q}},t)\right\}$, representing the instantaneous disorder of the material. Based on Equation (4), terms of the form $\left\{W_{s}(\bm{\mathrm{q}},t)\right\}$ are not sensitive reporters on the wavevector dependence of nonequilibrium phonon distributions because their value at every $\bm{\mathrm{q}}$ involves a sum of all mode, at every reduced wavevector $\bm{\mathrm{k}}$. Phonon population dynamics can only affect the magnitude of the Debye-Waller factors. The potential time dependence of the Debye-Waller factors was investigated, via the time dependence of mode populations $\left\{n_{j,\bm{\mathrm{k}}}\right\}$. Profoundly non-equilibrium distributions of phonon modes were simulated, with all modes populated equivalently to a temperature of $300\text{\,}\mathrm{K}$ except one mode at high temperature 555A maximum of $5000\text{\,}\mathrm{K}$ for optical modes, and $1000\text{\,}\mathrm{K}$ for acoustic modes. The discrepancy between maximum temperatures represents the fact that the heat capacity of acoustic modes is much higher.. These extreme non-equilibrium distributions increased the value the terms $\sum_{s}W_{s}(\bm{\mathrm{q}},t)$ by at most $1.5\text{\,}\mathrm{\char 37\relax}$ for optical modes, and $5\text{\,}\mathrm{\char 37\relax}$ for acoustic modes. Since these fractional changes are constant across $\bm{\mathrm{q}}$, wavevector-dependent changes in UEDS signals are not impacted significantly by transient changes to the Debye-Waller factors and any time dependence of the one-phonon structure factors $|F_{1j}(\bm{\mathrm{q}},t)|^{2}$ that result from the Debye-Waller factors themselves can be ignored to a good approximation.

IV.3 One-phonon structure factors

Using the computed Debye-Waller factors from the previous section, the calculation of $|F_{1j}(\bm{\mathrm{q}},t)|^{2}$ was carried out, from Equation (3), for the eight in-plane phonon modes of graphite: the longitudinal modes LA, LO1 – LO3, and the transverse modes TA, TO1 – TO3 666The calculation of $|F_{1j}(\bm{\mathrm{q}},t)|$ is trivial for out-of-plane modes ZA, ZO1 – ZO3 because $\bm{\mathrm{q}}\cdot\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\equiv 0$ for these modes.. The resulting one-phonon structure factors of a few in-plane modes, with occupations equivalent to a temperature of $300\text{\,}\mathrm{K}$, are shown in Figure 3. One-phonon structure factors $|F_{1j}(\bm{\mathrm{q}},t)|^{2}$ display striking scattering vector dependence (selection rules) based on the nature of phonon polarization vectors $\left\{\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\right\}$. Specifically, near $\bm{\mathrm{\Gamma}}$, the one-phonon structure factor for longitudinal modes is highest in the radial direction, because the polarization of those modes is parallel to $\bm{\mathrm{q}}$. On the other hand, $|F_{1j}(\bm{\mathrm{q}},t)|^{2}$ for transverse modes is highest (near $\bm{\mathrm{\Gamma}}$) in the azimuthal direction for transverse modes, because the polarization of those modes is perpendicular to $\bm{\mathrm{q}}$.

An alternative view of one-phonon structure factors is presented via weighted dispersion curves, an example of which is shown in Figure 4. This presentation allows easy comparison of the relative weights of the one-phonon structure factors along high-symmetry lines for different phonon branches.

A cursory inspection of the weighted dispersion curves in Figure 4 suggest that there are regions in the Brillouin zone where diffuse intensity is strongly biased towards a single mode (strong scattering selection rule) based on the relative intensities of one-phonon structure factors. Careful analysis reveals that there are very few wavevectors $\bm{\mathrm{q}}$ for which a particular phonon mode’s one-phonon structure factor is strongly dominant. Figure 5 presents a comparison of the relative intensities of one-phonon structure factors weighted by phonon frequency, which is indicative of mode population as per Equation (2). Only $37\%$ of wavevectors visible in measurements shown in Figure 2 have a mode that contributes over $50\%$ of the quantity $\sum_{j}|F_{1j}(\bm{\mathrm{q}})|^{2}/\omega_{j,\bm{\mathrm{k}}}$; only $1\%$ of wavevectors have a phonon mode that contributes more than $75\%$.

The results of Figure 5 show that quantitative answers regarding phonon dynamics from UEDS measurements cannot generally be obtained by inspection; at almost any wavevector $\bm{\mathrm{q}}$, at least two phonon modes contribute significantly to the transient diffuse scattering intensity. Therefore, a more robust procedure, presented in the next section, must be employed to extract wavevector- anndmode-dependent phonon populations from UEDS intensities.

IV.4 Population dynamics across the Brillouin zone

Transient electron diffuse intensity has been used elsewhere Chase et al. (2016); Waldecker et al. (2017); Stern et al. (2018); Konstantinova et al. (2018) as an approximation to the population dynamics of particular modes. However, one can extract the transient wavevector-dependent phonon population dynamics $\left\{\Delta n_{j,\bm{\mathrm{k}}}(t)\right\}$ by combining the measurements of $\Delta I(\bm{\mathrm{q}},t)$ with the calculations of one-phonon structure factors and associated quantities presented above.

For many materials (including graphite), the temperature dependence (and hence time dependence) of the phonon mode vibrational frequencies is negligible, because such dependence is proportional to anharmonic couplings between branches Calizo et al. (2007); Judek et al. (2015); hence, $\omega_{j,\bm{\mathrm{k}}}(t)\equiv\omega_{j,\bm{\mathrm{k}}}(t_{0})$. Moreover, as is discussed in Section IV.2, the temperature dependence (and hence time dependence) of the Debye-Waller factors — and therefore the one-phonon structure factors — has a much smaller magnitude that the variations due to other terms in Equation (2). Therefore, we have $|F_{1j}(\bm{\mathrm{q}},t)|^{2}\equiv|F_{1j}(\bm{\mathrm{q}},t_{0})|^{2}$ for all times. In this case, transient scattering intensity at the detector, $\Delta I(\bm{\mathrm{q}},t)$, can be expressed as follows:

[TABLE]

for $\bm{\mathrm{k}}$ away from $\bm{\mathrm{\Gamma}}$, where there might be interference with elastic scattering signals.

At every reduced wavevector $\bm{\mathrm{k}}$ and time-delay $t$, there are $N$ different values $\left\{\Delta n_{j,\bm{\mathrm{k}}}(t)\right\}$ that must be determined — one for each phonon mode. Since the one-phonon structure factors $|F_{1j}(\bm{\mathrm{q}},t_{0})|^{2}$ vary over the total wavevector $\bm{\mathrm{q}}$, the transient diffuse intensity for at least $N$ Brillouin zones must be considered so that Equation (6) can be solved numerically. A linear system of equations must be solved at every reduced wavevector $\bm{\mathrm{k}}$.

Let $\{\bm{\mathrm{H}}_{1},...,\bm{\mathrm{H}}_{M}\mid M\geq N\}$ be the chosen reflections from which to build the system of equations. Then, the transient phonon population of mode $j$ at every $\bm{\mathrm{k}}$ and time $t$, $\Delta n_{j,\bm{\mathrm{k}}}(t)$, solves the linear system:

[TABLE]

where

[TABLE]

for $\bm{\mathrm{k}}$ vectors away from the $\bm{\mathrm{\Gamma}}$ point. These linear systems of equations can be solved numerically, provided enough experimental data ($M\geq N$).

The choice to solve for $\bm{\mathrm{n}}_{\bm{\mathrm{k}}}(t)$, rather than for the change in population, is related to the degree of confidence that should be placed in the calculation of phonon polarization vectors and frequencies. The phonon polarization vectors are mostly affected by the symmetries of the crystal. On the other hand, phonon vibrational frequencies might be influenced by non-equilibrium carrier distributions. Solving for the ratio of populations to frequencies, rather than populations, is more robust against the uncertainty in the modelling, because the one-phonon structure factors only take into account the polarization vectors.

We also note that the procedure presented above can be easily extended to (equilibrium) thermal diffuse scattering measurements, where the phonon populations are known at constant temperature, but the phonon vibrational spectrum is unknown. Therefore, using pre-photoexcitation data of a time-resolved experiment, one could infer the phonon vibrational frequencies, which are then used to determine the change in populations using the measurements after photoexcitation. This scheme only relies on the determination of phonon polarization vectors.

Wavevector-dependent phonon population dynamics in graphite

We applied this general formalism for wavevector-dependent phonon population decomposition to the transient diffuse intensity patterns of photoexcited graphite shown in Figure 2. Since the diffraction patterns have been symmetrized, one would expect that using intensity data for reflections related by symmetry would be redundant. However, better results were achieved by using the entire area of the detector. We expect this is due to minute misalignment of the diffraction patterns and uncertainty in detector position which are averaged out when using all available data. The Brillouin zones associated with all in-plane Bragg reflections $\bm{\mathrm{H}}$ such that $|\bm{\mathrm{H}}|\leq$12\text{,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1} were used, for a total of fourty-four Brillouin zones ($M=44$), many more than the minimum required for the eight in-plane phonon modes of graphite ($N=8$). The physical constraint that $\Delta n_{j,\bm{\mathrm{k}}}(t)>0~{}\forall t$ was applied 777The positivity constraint $\Delta n_{j,\bm{\mathrm{k}}}(t)>0~{}\forall t$ means that the phonon population of a branch cannot drop below its equilibrium level. While not necessary, it leads to more reliable solutions. via the use of a non-negative approximate matrix inversion approach Lawson and Hanson ([1987](#bib.bib47)) to solve Equation ([7](#S4.E7)) at every reduced wavevector $\bm{\mathrm{k}}$ and time-delay $t$. Stable solutions were found for reciprocal space points where $|\bm{\mathrm{k}}|>$0.45\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}, where there is no interference between elastic $I_{0}(\bm{\mathrm{q}},t)$ and diffuse $I_{1}(\bm{\mathrm{q}},t)$ signals. Figure 6 presents the direct decomposition of diffuse intensity into wavevector-dependent, transient phonon population changes $\left\{\Delta n_{j,\bm{\mathrm{k}}}(t)\right\}$, for a few in-plane modes that are particularly relevant to graphite. The discussion of physical processes that explain the wavevector-dependent transient phonon populations follows.

In graphite, optical excitation creates a nonthermal phonon distribution, increasing population primarily in two strongly-coupled optical phonons (SCOP): $A_{1}^{\prime}$, located near the $\bm{\mathrm{K}}$ point, and $E_{2g}$, at $\bm{\mathrm{\Gamma}}$ Kampfrath et al. (2005). Dynamics measured at the earliest time scales ($<$5\text{,}\mathrm{ps}$$) are discussed qualitatively in a previous publication Stern et al. (2018). Using the measured population dynamics of Figure 6, we can track the transfer of energy across the Brillouin zone quantitatively.

By conservation of both momentum and energy, the anharmonic decay of the two SCOP transfers population into mid-Brillouin zone acoustic modes. The early time-points presented in Figure 6 confirm that quickly after photoexcitation ($<$500\text{,}\mathrm{fs}$$), the transverse optical mode TO2 is strongly populated at $\bm{\mathrm{K}}$, indicative of the expected strong electron-phonon coupling to the $A_{1}^{\prime}$ phonon. The transfer of energy away from the TO2 mode is already well underway at $5\text{\,}\mathrm{ps}$, associated with an increase in acoustic modes along the $\bm{\mathrm{\Gamma}}$ – $\bm{\mathrm{M}}$ line. This is in accordance with the phonon band structure, where the mid-point along the $\bm{\mathrm{\Gamma}}$ – $\bm{\mathrm{M}}$ line favors occupancy of the TA mode (Figure 4). This behaviour intensifies from $1.5\text{\,}\mathrm{ps}$ to $100\text{\,}\mathrm{ps}$.

The initial increase ($<$5\text{,}\mathrm{ps}$$) of TA population along the $\bm{\mathrm{\Gamma}}$ – $\bm{\mathrm{M}}$ line is in excellent agreement with predicted anharmonic decay probabilities from the $E_{2g}$ phonon Bonini et al. (2007). The (small) increase at $500\text{\,}\mathrm{fs}$ in LA population at $\tfrac{1}{2}\bm{\mathrm{K}}$ is also in line with calculated decay probabilities from the $E_{2g}$ phonon by anharmonic coupling 888On the other hand, a monotonic increase in LA population at $\bm{\mathrm{\Gamma}}$ is expected from the decay of the other SCOP, $A_{1}^{\prime}$, which has a high probability to decay into a pair of LA-LO modes. However, we were not able to measure population changes close enough to $\bm{\mathrm{\Gamma}}$..

Over longer time scales ($>$25\text{,}\mathrm{ps}$$), the TA population has pooled significantly at $\tfrac{1}{3}\bm{\mathrm{M}}$ and $\bm{\mathrm{M}}$. There are no three phonon anharmonic decay processes that start in a purely transverse mode; the only allowed interband transitions are L $\leftrightarrow$ T + T and L $\leftrightarrow$ L + T, where L (T) represents a longitudinal (transverse) mode Lax et al. (1981); Khitun and Wang (2001); therefore, a build-up of population in the TA mode is expected. At $\tfrac{1}{2}\bm{\mathrm{M}}$, computed lifetimes predict that both LA and TA phonons will favor decay processes into out-of-plane phonons (ZA) Paulatto et al. (2013) that are not visible (have zero one-phonon structure factor) in the [001] zone-axis geometry in which these UEDS experiments were conducted. Similarily, the ZA phonons predominantly decay back into in-plane phonons, implying that the phonon thermalization in the acoustic branches occurs through a mechanism that exchanges in-plane and out-of-plane modes. Additionally, the computed LA and TA anharmonic lifetimes are predicted to significantly drop at $\tfrac{1}{2}\bm{\mathrm{K}}$ and $\tfrac{1}{2}\bm{\mathrm{M}}$, respectively, due to the activation of Umklapp scattering to the ZA phonons. Our measurements corroborate these predictions, as can be seen by TA and LA population at the mid-Brillouin zone (dashed white hexagon) being relatively lower than average. The confirmation of those predictions fundamentally relies on UEDS’ ability to probe the entire Brillouin zone at once.

The robustness of such an analysis must be emphasized. The decomposition of transient diffuse intensity change via Equation (7) admits no free parameter. Given sufficient data, a single optimal solution exists.

IV.5 Wavevector-dependent electron-phonon coupling

The flow of energy between electronic and phononic subsystems is typically crudely modelled using the two-temperature model Allen (1987). This model assumes that the electronic system and the phononic system can each be associated with temperatures $T_{e}$ and $T_{ph}$, throughout the dynamics. Effectively this approximation assumes that the internal thermalization dynamics of each system is much more rapid than any processes that couple the two system. It is evident from the earlier description of the UEDS data from graphite that this assumption is (rather generally) quite a poor one; the idea that the phononic subsystem is internally thermalized does not hold on the timescales typically associated with energy flows between the electron and phonon systems following photoexcitation. On these timescales the phonon occupations are generally very far from being thermalized. UEDS allows to move beyond the two-temperature approximation; by leveraging momentum-resolution, mode-dependent electron-phonon and phonon-phonon couplings can be extracted from the transient change in mode populations $\left\{\Delta n_{j,\bm{\mathrm{k}}}(t)\right\}$. Specifically, wavevector-dependent phonon population dynamics determined in the previous section will now be used to determine the electron-phonon and phonon-phonon coupling strength of the $A_{1}^{\prime}$ phonon.

The formalism of the two-temperature model can be extended to the non-thermal lattice model (NLM) model Waldecker et al. (2016), where every phonon branch $j$ has its own molar heat capacity $C_{ph,j}$, and temperature $T_{ph,j}$:

[TABLE]

[TABLE]

where $f(t-t_{0})$ is the laser pulse profile, and $C_{e}$ and $T_{e}$ are the electronic heat capacity and electron temperature, respectively 999The electronic system thermalizes in approximately $100\text{\,}\mathrm{fs}$ Stange et al. (2015), and hence after we can consider the electronic system to be well-described by a single temperature $T_{e}$.. This model accounts for discrepancies in coupling between the electronic system and certain phonon modes, which occurs for example in graphite — where some modes are strongly-coupled to the electron system via Kohn anomalies Piscanec et al. (2004).

Observations of transient changes in mode populations are related to mode temperatures $T_{ph,j}(t)$ via the Bose-Einstein distribution:

[TABLE]

We can decompose the above expression with a Laurent series Wunsch (2005) to show explicitly that the mode population is proportional to temperature, for appropriately high $T_{ph,j}$:

[TABLE]

Hence, in the case of measurements presented herein, $\Delta n_{j,\bm{\mathrm{k}}}(t)\propto\Delta T_{ph,j}(t)$, where the initial temperature is known to be $300\text{\,}\mathrm{K}$.

We now use the NLM to extract the couplings to the $A_{1}^{\prime}$ mode from population measurements. The differential $A_{1}^{\prime}$ phonon population $\Delta n_{A_{1}^{\prime}}(t)$ is obtained by integrating over the region of the wavevector-dependent TO2 phonon population, in a circular arc centered at $\bm{\mathrm{K}}$ ($|\bm{\mathrm{k}}-\bm{\mathrm{K}}|\leq$0.3\text{,}{\mathrm{\SIUnitSymbolAngstrom}}^{-1}$$). This location is shown in Figure 6. In order to correlate the mode population measurements with the NLM, the heat capacities of the electronic system and every phonon mode must be parametrized.

The electronic heat capacity $C_{e}$ is extracted from experimental work by Nihira and Iwata (2003):

[TABLE]

Over this range of time-delays, thermal expansion (or contraction, in the case of graphite) has not yet occurred Chatelain et al. (2014). No changes in Bragg peak positions — indicative of lattice parameter changes — is observed within the experimental range of time-delays $\tau\leq$680\text{,}\mathrm{ps}$$. We can therefore calculate the heat capacity of each graphite mode $j$ as the heat capacity at constant volume Ziman (1979):

[TABLE]

where $k_{B}$ is the Boltzmann constant, $\omega_{D}$ is the Debye frequency, and $D(\omega)$ is the phonon density of states. Momentum resolution of UEDS allows for a simplification, where a single frequency contributes to the heat capacity in the $A_{1}^{\prime}$ mode — $D(\omega)=\delta(\omega)$. Moreover, we can reduce the number of coupled equations in Equations (11) and (12). Simultaneous conservation of momentum and energy during the decay of an $A_{1}^{\prime}$ phonon can only be satisfied in a few reciprocal space locations. Using first-principles calculations, it is possible to determine the decay probabilities. One such calculation, reported by Bonini et al. (2007), allows us to define an effective heat capacity into which the $A_{1}^{\prime}$ population drains, $C_{l}$ 101010This effective heat capacity $C_{l}$ is composed of 9% chance to decay into two TA modes, 36% change to decay into a TA mode and an LA mode, and 55% chance to decay into either an LA and TA, or LO and LA.. Therefore, the energy dynamics at $\bm{\mathrm{K}}$ can be specified in terms of a system of three equations:

[TABLE]

where $G_{e,A_{1}^{\prime}}$, $G_{A_{1}^{\prime},l}$, and $G_{e,l}$ are constants. Solving this system of equations gives the temperature evolution of each of the subsystems. The evolution in the $A_{1}^{\prime}$ population can be used as a proxy for the mode temperature $T_{A_{1}^{\prime}}(t)$; minimizing the difference between observed population dynamics and modelled temperature changes yields the coupling constants $G_{e,A_{1}^{\prime}}$, $G_{A_{1}^{\prime},l}$, and $G_{e,l}$. The resulting temperature transients are presented in Figure 7. The extracted coupling constants have been listed in Table 1. This model correctly identifies the strong electron-phonon coupling of the $A_{1}^{\prime}$ mode, $G_{e,A_{1}^{\prime}}$, as compared with the rest of the relevant modes, $G_{e,l}$.

From the coupling constant $G_{e,A_{1}^{\prime}}$, mode-projected electron-phonon coupling value $\langle g_{e,A_{1}^{\prime}}^{2}\rangle$ can be determined. In the case of the coupling between the electron system and the $A_{1}^{\prime}$ phonon, the heating rate of $G_{e,A_{1}^{\prime}}=$(6.8\pm 0.3)\text{\times}{10}^{17}\text{,}\mathrm{W}\text{,}{\mathrm{m}}^{-3}\text{,}{\mathrm{K}}^{-1} (Table [1](#S4.T1)) corresponds to a mode-projected electron-phonon coupling value of $\langle g_{e,A_{1}^{\prime}}^{2}\rangle=$0.035\pm 0.001\text{\,}{\mathrm{eV}}^{2} (see appendix B for details). These values are in agreement with recent trARPES measurements and simulations Johannsen et al. (2013); Stange et al. (2015); Rohde et al. (2018); Na et al. (2019).

V Conclusion

UEDS provides direct access to wavevector-resolved, non-equilibrium phonon populations and is, in this sense, a lattice-dynamical analog of trARPES. A robust and generally-applicable UEDS data reduction method has been described that provides detailed information on transient changes in phonon populations across the entire Brillouin zone that follow photoexcitation in single-crysalline materials. This method takes only the observed UEDS patterns and computed one-phonon structure factors as inputs and can easily be extended with minimal alterations to ultrafast X-ray diffuse scattering. A procedure for computing the required phonon properties using DFPT, and their potential time dependence via the Debye-Waller factors was described in detail. This method was demonstrated for the case of photodoped carriers in the Dirac cones of thin graphite, where the phonon populations were tracked. Finally, the mode dependence of couplings between electron and phonons have been demonstrated at a specific point in the Brillouin zone, where the strongly-coupled optical phonon $A_{1}^{\prime}$ is located. Mode-projected electron-phonon coupling value for the $A_{1}^{\prime}$ phonon was extracted, using the non-thermal lattice model, and corroborated with numerous other experiments and simulations.

Direct determination of wavevector-dependent, transient phonon populations holds great promise for the study of phenomena that emerge primarily due to the coupling of electronic and lattice degrees of freedom, and specifically those involving strongly anisotropic interactions. In particular, with sufficient time-resolution, the applicability of the Kramers-Heisenberg-Dirac theory to Raman scattering measurements in graphene/graphite could be explored, via the detection of early-times (<$$50\text{\,}\mathrm{fs}) phonon populations in the strongly-coupled optical phonon $E_{2g}$ Heller et al. (2016). Another potential extension concerns influence of non-equilibrium carrier distributions on phonon vibrational frequencies. Many systems, charge-density wave materials in particular, exhibit phonon modes that harden or soften at high temperatures and selective electronic excitation which can be used to explore such phenomena in greater depth.

Acknowledgements

L. P. R. de C. thanks H. Seiler for illuminating discussions regarding the role of the Debye-Waller effect on diffuse intensity, and M. X. Na for providing insights regarding the relationship between heat rates and mode-projected electron-phonon coupling constants.

This research was enabled in part by support provided by Calcul Quebec (www.calculquebec.ca) and Compute Canada (www.computecanada.ca)

This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC), the Fonds de Recherche du Québec - Nature et Technologies (FRQNT), the Canada Foundation for Innovation (CFI), and the Canada Research Chairs (CRC) program.

Author contribution

M. S. and B. J. S. conceptualized the work. L. P. R. de C. performed the research. L. P. R. de C. and B.J.S. wrote the manuscript. J.-H. P. performed the DFPT calculations. All authors helped edit the article.

Appendix A Clustering of phonon eigenvalues and eigenvectors into branches

This section describes the clustering of phonon polarization vectors $\left\{\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\right\}$ and frequencies $\left\{\omega_{j,\bm{\mathrm{k}}}\right\}$ into physically-relevant categories, i.e. branches. The general idea behind the procedure is that phonon properties are continuous. The variation of a property should not display any discontinuity along any path in the Brillouin zone.

Let $\bm{\mathrm{P}}_{j,\bm{\mathrm{k}}}$ be an abstract vector representing the polarization vectors and frequency of branch $j$ at reduced wavevector $\bm{\mathrm{k}}$. We represent $\bm{\mathrm{P}}_{j,\bm{\mathrm{k}}}$ as the following vector:

[TABLE]

where the index $s$ runs for all $M$ atoms in the unit cell ($M=4$ in the case of graphite). We define the metric between two abstract vectors $\bm{\mathrm{P}}_{i,\bm{\mathrm{k}}}$ and $\bm{\mathrm{P}}_{j,\bm{\mathrm{k}}}$ as follows:

[TABLE]

A one-dimensional path $\gamma(\bm{\mathrm{k}})$ connecting all $\bm{\mathrm{k}}$-points was defined, starting at $\bm{\mathrm{\Gamma}}$. At $\bm{\mathrm{\Gamma}}$, polarization vectors are associated with a mode based on geometry and oscillation frequency. For example, polarization vectors $\left\{\bm{\mathrm{e}}_{j,s,\bm{\mathrm{k}}}\right\}$ parallel to their wavevector $\bm{\mathrm{k}}$ for all atoms $s$ is a longitudinal mode; if the associated frequency is $\approx$0\text{,}\mathrm{THz}$$, this mode can be labelled longitudinal acoustic. Then, polarization vectors and frequencies at any point along the path were assigned to modes that optimized continuity. That is, the assignment of phonon branches $j$ at $\gamma(\bm{\mathrm{k}}+\bm{\mathrm{\bm{\mathrm{\Delta}}}})$, $\bm{\mathrm{P}}_{j,\gamma(\bm{\mathrm{k}}+\bm{\mathrm{\bm{\mathrm{\Delta}}}})}$ based on the assignment at $\gamma(\bm{\mathrm{k}})$, $\bm{\mathrm{P}}_{i,\gamma(\bm{\mathrm{k}})}$, minimized the distance $\lVert\bm{\mathrm{P}}_{i,\gamma(\bm{\mathrm{k}})}-\bm{\mathrm{P}}_{j,\gamma(\bm{\mathrm{k}}+\bm{\mathrm{\bm{\mathrm{\Delta}}}})}\rVert$.

The procedure described above, adapted for numerical evaluation, is part of the scikit-ued software package René de Cotret et al. (2018).

Appendix B Calculation of mode-projected electron-phonon coupling from heating rates

Consider the coupled equations of the non-thermal lattice model in Equations (11) and (12). These coupled first-order ordinary differential equations will admit solutions for $T_{e}(t)$ and $\left\{T_{ph,j}(t)\right\}$. After photoexcitation ($f(t-t_{0})\to 0$), the appropriate summations of those equations yields the following single equation:

[TABLE]

where the temperature dependence of $C_{e}$ and $\left\{C_{ph,j}\right\}$ has been omitted for brevity. In the case of graphite, at early times ($<$5\text{,}\mathrm{ps}$$), phonon-phonon coupling $G_{pp,ij}$ is much weaker at the $\bm{\mathrm{K}}$-point (Table 1). Therefore, we may simplify the above equation to a more manageable system:

[TABLE]

By performing a substitution $\lambda=T_{e}-\sum_{j}T_{ph,j}$, the equation above simplifies to a familiar situation:

[TABLE]

where

[TABLE]

The time dependence comes from the time-evolution of the individual temperatures. In the case of phonon temperatures, the phonon population dynamics are directly related to temperature dynamics according to Equation (14). Equation (24) is a separable equation with solution:

[TABLE]

For a slow-varying integrand $a(t)\approx a$, then $a=1/\tau$, where $\tau$ is a compound variable representing the relaxation of the system. This leads to the following form:

[TABLE]

As a specific example, the above expression reduces nicely in the case of the two-temperature model, where all phonon modes are considered to be thermalized with each other, with isochoric heat capacity $C_{ph}$:

[TABLE]

and we see that $\tau$ physically represents the relaxation time of the electronic system into the lattice. Equation (27) can be thought of as a sum of relaxation times between the electronic subsystem and specific modes $\tau_{e,j}$:

[TABLE]

The final state in relating heating rates to their mode-projected coupling values requires knowledge about density of states. Because the measurements herein consider only in-plane interactions, we use an approximate electronic density of states for graphene close to the Dirac point Neto et al. (2009):

[TABLE]

where $A$ is the unit cell area and $v_{F}=$9.06\text{\times}{10}^{5}\text{,}\mathrm{m}\text{,}{\mathrm{s}}^{-1}$$ is the Fermi velocity 111111Note that factors of $\hbar$ are often ignored, including in the accompanying reference.. The electronic density of states is related to the mode-projected electron-phonon coupling $\langle g_{ep,j}^{2}\rangle$ as follows Na et al. (2019):

[TABLE]

where $\hbar\omega_{\nu}$ corresponds to the optical excitation energy ($1.55\text{\,}\mathrm{eV}$ or $800\text{\,}\mathrm{nm}$).