Establishing nonthermal regimes in pump-probe electron relaxation dynamics

Time- and angle-resolved photoemission spectroscopy (TR-ARPES) accesses the electronic structure of solids under optical excitation, and is a powerful technique for studying the coupling between electrons and collective modes. One approach to infer electron-boson coupling is through the relaxation dynamics of optically excited electrons, and the characteristic timescales of energy redistribution. A common description of electron relaxation dynamics is through the effective electronic temperature. Such a description requires that thermodynamic quantities are well-defined, an assumption that is generally violated at early delays. Additionally, precise estimation of the nonthermal window—within which effective temperature models may not be applied—is challenging. We perform TR-ARPES on graphite and show that Boltzmann rate equations can be used to calculate the time-dependent electronic occupation function $f(\epsilon ,t)$, and reproduce experimental features given by nonthermal electron occupation. Using this model, we define a quantitative measure of nonthermal electron occupation and use it to define distinct phases of electron relaxation in the fluence-delay phase space. More generally, this approach can be used to inform the nonthermal-to-thermal crossover in pump-probe experiments.

Figure 1

TR-ARPES measurements on graphite. (a) ARPES spectra from graphite measured along the $\overline{\mathrm{\Gamma }}-\overline{\mathrm{K}}$ direction before pump arrival and at zero delay. The equilibrium (before pump-arrival) temperature is 50 K. The pump pulse has a photon energy of 1.2 eV and a duration of 120 fs. (b) Momentum-integrated energy distribution curves (${\int }_k\mathrm{EDC}$) (linear scale) shows a strong deviation away from the FD function due to the compounded contribution of the dispersion and matrix element effects (see text). (c) Momentum-integrated differential energy distribution curves ($\mathrm{\Delta }{\int }_k\mathrm{EDC}$) computed from the curves in panel (b), using ${\int }_k{\mathrm{EDC}}_3$ as a reference. The blue (orange) curve is the $\mathrm{\Delta }{\int }_k\mathrm{EDC}$ at $-400$ fs (180 fs), and is labeled $\mathrm{\Delta }{\int }_k{\mathrm{EDC}}_{1,3}$ ($\mathrm{\Delta }{\int }_k{\mathrm{EDC}}_{2,3}$). The shaded region indicates the difference between two FD distributions, which bears similarity to the blue curve, but cannot describe the orange one. Accumulation of carriers within the phonon window is indicated by black arrows; dashed lines indicate the phonon window associated with the ${\mathrm{A}}_1^{\prime }$ optical phonon, with momentum $K$ and energy $\hslash {\mathrm{\Omega }}_{{\mathrm{A}}_1^{\prime }}=0.16$ eV.

Figure 2

Boltzmann simulations of electron relaxation in graphite. (a) The occupation function $f(\epsilon ,t)$ calculated from the Boltzmann rate equations (logarithmic scale). The zero delay distribution is given by the red curve in energy; two peaks associated with optical transitions and a peak created through scattering with the strongly coupled optical phonon (SCOP) are highlighted in blue and red, respectively. The dynamics of each peak are also highlighted. For visualization, we add a background of ${10}^{-5}$ to $f(\epsilon ,t)$ to represent the dynamic range of the experiment. (b) $f(\epsilon ,t)$ at the same delays in Fig. 1. (c) Differential occupation using $f(\epsilon ,t_3)$ as a reference. $\mathrm{\Delta }{f}_{1,3}$ and $\mathrm{\Delta }{f}_{2,3}$ are shown in blue and orange, respectively, reproducing the changes in sign observed in the data [Fig. 1]. The phonon window associated with the ${\mathrm{A}}_1^{\prime }$ mode is indicated by black dashed lines. (d) Electron and phonon effective temperatures extracted by fitting the distribution with a FD and inverting the Bose-Einstein distribution, respectively. Yellow shaded region indicates the 90% confidence interval of the fit. Optical (high energy) phonons are shaded in red; acoustic phonons are shaded in black. (e) The occupation in a 0.1 eV window centered at 0.6 eV is given by the red markers. The mean-squared error $[\mathrm{MSE}=1/n{\sum }_n{(y_{\mathrm{fit}}-y)}^2]$ of the fit is given by black markers. MSE is a quantitative measure of the nonthermal contributions to the electronic distribution, and it is nonzero long after the pump pulse has passed.

Figure 3

Nonthermal phases in the fluence-delay phase space. (a) Effective electronic temperature ($T_e$) extracted from a FD fit of the occupation function. The green dashed line is a contour at 1000 K; the different fluence regimes are qualitatively indicated by labels. (b) MSE obtained from the fit of the effective electronic temperature. Three phases are outlined by red, green, and magenta dashed lines, respectively; colored markers indicate the delay and fluence at which the occupation function and residuals are shown in panels (c) and (d). (${\mathrm{c}}_{1-3}$) Simulated distribution $f(\epsilon ,t)$ in phases 1, 2, and 3, respectively, are given by dotted lines matching the colored markers in panel (b). FD fits of the occupation function are shown in dashed lines in corresponding colors. The blue and red arrows in phases 1 and 2 highlight the peaks corresponding to optical excitation and phonon emission, respectively. (${\mathrm{d}}_{1-3}$) Residuals of the FD fit at delays corresponding to the markers in panel (b); despite being taken at different delays, residuals in the same phase have the same qualitative features, supporting the classification into these separate phases. (e) Cartoons indicating the primary mechanisms in the three identified phases: (phase 1) injection, cascade, and accumulation within the phonon window; (phase 2) redistribution of accumulated electrons into an energetic quasi-FD distribution; and (phase 3) cooling of the distribution inside the phonon window via phonon emission.

Figure 4

Simulation of the ${\int }_k\mathrm{EDC}$. (a) The spectral function calculated with ${\mathrm{\Sigma }}^{\prime }=0, {\mathrm{\Sigma }}^{\prime \prime }=30$ meV, along the $\mathrm{\Gamma }\text{−}\mathrm{K}$ direction for $k_z=0.2{\AA }^{-1}$. The dispersion ${\epsilon }_{\mathbf{k}}$ is shown in dashed red lines. (b) The spectral function modified by the Fermi-Dirac and convolved with a Gaussian in energy and momentum to account for system resolution effects. (c) Intensity profiles obtained by integrating the spectral function over $d{k}_x$, at $k_y=0$ (black) and $d{k}_{x}d{k}_y$ (blue). Corresponding density of states in 1D and 2D are given in red and green dashed lines, respectively. (d) Intensity profiles obtained by integrating the modified spectral function over $d{k}_x$, at $k_y=0$ (black) and $d{k}_{x}d{k}_y$ (blue).

Figure 5

Photoexcitation and e-e scattering. [(a) and (c)] Occupation function simulated with only the injection term, and the injection plus e-e scattering, respectively. Delays at [$-500$, 0, 600] fs are highlighted in blue, red, and yellow, respectively. [(b) and (d)] The excess energy in the different subsystems as a function of delay. The solid yellow line gives the electron excess energy. The solid blue line gives the energy of the coupled phonon modes. The whole phonon subsystem is given by the dashed orange line. The dashed purple line gives the total energy of the system. There is no electron-phonon coupling; thus, excess energy in the phonon system is zero.

Figure 6

Electron-phonon coupling and anharmonic phonon decay. [(a) and (c)] Simulation of $f(\epsilon ,t)$ using the injection and electron-phonon (e-ph) terms and simulation using the injection, electron-phonon, and phonon-phonon (ph-ph) terms respectively. Delays at [$-500$, 0, 600] fs are highlighted in blue, red, and yellow, respectively. [(b) and (d)] The excess energy in the different subsystems as a function of delay. In (b), e-ph coupling allows the transfer of energy from the electron subsystem to the phonon subsystem, but the absence of ph-ph scattering means the excess energy of the strongly coupled phonon (SCP) is equivalent to that of the entire phonon subsystem. In (d), turning on ph-ph scattering allows the SCP to decay into other phonon modes; thus excess energy of the SCP decreases, while the excess energy of the phonon bath remains equivalent to that in (b).

Figure 7

Simulation checks. (a) The energy of the electron and phonon subsystems in the simulation as a function of pump-probe delay are shown in blue and red solid lines. The system's total energy (magenta solid line) converges to that injected by the pump pulse with no scattering processes. (b) The left-hand side (LHS) of the Boltzmann rate equations [Eq. (3)]. Carriers are injected/depleted near 0.6 eV. The dynamics are well contained within the energy range $[-1,1]$ eV. (c) The occupation function at $-400$, 180, and 630 fs, calculated with a mesh of 5 meV (markers) and 2 meV (line). The two overlap, showing a good convergence.

Figure 8

Electron scattering rate calculated from data and experiment (a) $\partial f/\partial t$ calculated from Eq. (3) in the main text. (b) $\partial f/\partial t$ is calculated from the simulation. The output of the simulation $f(\epsilon ,t)$ is integrated into bins 50 meV wide, the tangent for each delay is extracted from a moving fit of the intensity in a 50 fs window with a first-order polynomial. (c) $dI/dt$ calculated from the data. The intensity is integrated into bins 50 meV wide, and the tangent for each delay is extracted from a moving fit of the intensity in a 50 fs window with a 1st order polynomial. The dynamics are contained within the energy range $[-1,1]$ eV. Arrows indicate nonthermal features, including the direct transition at 0.6 and 0.5 eV (blue) and the phonon-induced replica at 0.44 eV (red). Green dashed arrows mark the energy domain $\pm \hslash \mathrm{\Omega }$. In all panels, the region between 0.4 and 0.7 eV is enhanced by a factor of 20 for visualization purposes.

Figure 9

Differential maps (a, b) $\mathrm{\Delta }{\int }_k\mathrm{EDC}$ maps computed using the ${\mathrm{Ref}}_1$ (${\mathrm{Ref}}_3$), at delays -400 fs (630 fs) respectively. Blue (red) represent negative (positive) values. (c, d) $\mathrm{\Delta }f$ maps using the output of the simulation in Fig. 2 in the main text. Panel (c) and (d) are computed using $f(\epsilon ,-400\mathrm{fs})$ and $f(\epsilon ,630\mathrm{fs})$, as reference.

Figure 10

Temperature fits of the data and simulation (a) EDCs at delays [$-400$, 180, 630] fs (markers) and the phenomenological fit model (lines). (b) Temperature extracted from the fit of the data in (a) (blue markers) and from Fermi-Dirac fits of the simulation (red line). Delays shown in (a) are indicated by dashed lines in the corresponding color. (c) The mean-squared error (MSE) from the effective temperature fit of data (blue markers) and the simulation (red line).