1. Introduction

In recent years, the rapid development of pump-probe femtosecond spectroscopy has enabled the observation of intriguing phenomena on sub-picosecond timescales^[1-7]. For example, ultrafast electron diffuse scattering provides direct access to the dynamics of the crystalline lattice and electron-phonon (e-ph) interactions with time and momentum resolution^[8-14], enabling a detailed understanding of energy transfer between hot electrons and phonons. During ultrafast processes, the system can be described as hot electrons transferring energy to a cold phonon bath via e-ph coupling. The well-known two-temperature model (TTM) provides a simple theoretical framework to describe this relaxation process, assuming that both the electron and phonon subsystems are always thermalized and characterized by distinct Fermi-Dirac (FD) and Bose-Einstein (BE) distributions at different temperatures^[15-17]. Within the TTM, the rate of energy transfer from electrons to phonons is governed by the electron-phonon coupling constant (G_ep), and an important implication is the connection between G_ep and the coupling constant λ used in the Eliashberg generalization of Bardeen-Cooper-Schrieffer (BCS) theory^[18]. This relationship allows λ to be extracted from pump-probe measurements, and interestingly, for several metals, the values of λ obtained in this way have shown good agreement with other independent methods, supporting the validity of the TTM^[19].

However, numerous theoretical calculations suggest that phonon lifetimes are typically on the order of picoseconds, which makes it unlikely for the phonon subsystem to reach thermal equilibrium on a sub-picosecond timescale^[20-26]. On the other hand, experimental results also indicate that the internal thermalization timescale for non-thermal electrons can be on the order of several hundred femtoseconds^[27-29]. These observations imply that the assumption of instantaneous thermalization of either subsystem is questionable, fundamentally challenging the basic premises of the TTM.

To gain deeper insights into the contradiction of the TTM, we combine density functional theory (DFT) calculations of e-ph and phonon-phonon (ph-ph) scattering rates with the time-dependent Boltzmann transport equation (TDBE) to investigate the time evolution of electrons and phonons under non-equilibrium conditions. Similar approaches have been applied to systems such as graphene^[30], MoS₂ monolayer^[31], FePt^[32], and other metals and semiconductors^[20,33]. These studies assume instantaneous electron thermalization and focus exclusively on the iterative evolution of the phonon distribution. This assumption is based on the fact that non-thermal electron features are observed experimentally only on femtosecond timescales, while electrons appear fully thermalized during picosecond-scale dynamics. Nevertheless, we argue that electrons can still retain localized non-thermal features under the influence of a cold phonon bath and play an important role in relaxation, because the separation between e-ph and electron-electron (e-e) scattering timescales is not always obvious. Therefore, in this work, we solve the coupled Boltzmann equations for both electrons and phonons, focusing on how non-thermal electrons and phonons together affect the ultrafast relaxation dynamics. We focus on the out-of-equilibrium dynamics in gold, silver, and aluminum. Our results reveal that even when electrons initially follow a thermal FD distribution, they can deviate substantially from the FD distribution during the relaxation process. Therefore, it becomes necessary to consider this non-thermal distribution explicitly for a quantitative understanding of the ultrafast e-ph relaxation dynamics.

2. Method

2.1 Theory of time-dependent Boltzmann transport equation

The distributions of electron (f_mk) and phonon (n_sq) evolve over time, and their dynamics are described by the TDBE^[34]:

(1)${\partial }_t{f}_{mk}=A_{mk}(1-f_{mk})-B_{mk}{f}_{mk}$

(2)${\partial }_t{n}_{sq}=A_{sq}(1+n_{sq})-B_{sq}{n}_{sq}$

where the right-hand side of the equations represents the collision integrals, and A/B correspond to the emission/absorption process. Here, e-ph scattering is considered for electrons, whereas both phonon-electron (ph-e) and ph-ph scattering are considered for phonons. The expressions can be written as:

(3)$\begin{array}{c}{A}_{mk}=\frac{2\pi }{ħN_k}\sum _{n{k}^{\prime }sq}{\left|g_{kq}^{smn}\right|}^2[f_{n{k}^{\prime }}(n_{sq}+1)\delta ({\epsilon }_{n{k}^{\prime }}-{\epsilon }_{mk}-ħ{\omega }_{sq})\\ +f_{n{k}^{\prime }}{n}_{sq}\delta ({\epsilon }_{n{k}^{\prime }}-{\epsilon }_{mk}+ħ{\omega }_{sq})]\end{array}$

(4)$\begin{array}{c}{B}_{mk}=\frac{2\pi }{ħN_k}\sum _{n{k}^{\prime }sq}{\left|g_{kq}^{smn}\right|}^2[(1-f_{n{k}^{\prime }})n_{sq}\delta ({\epsilon }_{n{k}^{\prime }}-{\epsilon }_{mk}-ħ{\omega }_{sq})\\ +(1-f_{n{k}^{\prime }})(n_{sq}+1)\delta ({\epsilon }_{n{k}^{\prime }}-{\epsilon }_{mk}+ħ{\omega }_{sq})]\end{array}$

(5)$\begin{array}{cc}& A_{sq}=\frac{2\pi }{ħN_q}\sum _{mkn{k}^{\prime }}{\left|g_{kq}^{smn}\right|}^2{f}_{n{k}^{\prime }}(1-f_{mk})\delta ({\epsilon }_{n{k}^{\prime }}-{\epsilon }_{mk}-ħ{\omega }_{sq})\\ & +\frac{\pi ħ}{4N_q}\sum _{s^{\prime }{q}^{\prime }{s}^{\prime \prime }{q}^{\prime \prime }}\frac{{\left|v_{sq{s}^{\prime }{q}^{\prime \prime }{s}^{\prime \prime }}\right|}^2}{{\omega }_{sq}{\omega }_{s^{\prime }{q}^{\prime }{\omega }_{s^{\prime \prime }{q}^{\prime \prime }}}}[\frac{1}{2}{n}_{s^{\prime }{q}^{\prime }}{n}_{s^{\prime \prime }{q}^{\prime \prime }}\delta ({\omega }_{sq}-{\omega }_{s^{\prime }{q}^{\prime }}-{\omega }_{s^{\prime \prime }{q}^{\prime \prime }})\end{array}$

(6)$\begin{array}{c}{B}_{sq}=\frac{2\pi }{ħN_q}\sum _{mkn{k}^{\prime }}{\left|g_{kq}^{sm\eta }\right|}^2(1-f_{n{k}^{\prime }})f_{mk}\delta ({\epsilon }_{n{k}^{\prime }}-{\epsilon }_{mk}+ħ{\omega }_{sq})\\ +\frac{\pi ħ}{4N_q}\sum _{s^{\prime }{q}^{\prime }{s}^{\prime \prime }{q}^{\prime \prime }}\frac{{\left|V_{sq{s}^{\prime }{q}^{\prime \prime }{q}^{\prime \prime }}\right|}^2}{{\omega }_{s{q}^{\prime }}{\omega }_{s^{\prime }{q}^{\prime \prime }}{\omega }_{s^{\prime \prime }{q}^{\prime \prime }}}[\frac{1}{2}(n_{s^{\prime }{q}^{\prime }}+1)(n_{s^{\prime \prime }{q}^{\prime \prime }}+1)\delta ({\omega }_{sq}-{\omega }_{s^{\prime }{q}^{\prime }}-{\omega }_{s^{\prime \prime }{q}^{\prime \prime }})\\ +n_{s^{\prime }{q}^{\prime }}(n_{s^{\prime \prime }{q}^{\prime \prime }}+1)\delta ({\omega }_{sq}+{\omega }_{s^{\prime }{q}^{\prime }}-{\omega }_{s^{\prime \prime }{q}^{\prime \prime }})]\end{array}$

where $g_{kq}^{smn}$ and $V_{sq{s}^{\prime }{q}^{\prime }{s}^{\prime \prime }{q}^{\prime \prime }}$ are the matrix elements for e-ph/ph-e and ph-ph scattering, respectively. N_k and N_q denote the numbers of electronic and phononic wavevectors in the Brillouin zone.

For e-e scattering, we employ the RTA:

(7)${\partial }_t{f}_{mk}={{\partial }_t{f}_{mk}|}_{e-ph}+\frac{f_{mk}-f_{mk}^0}{{\tau }_{e-e}}$

where the first term on the right-hand side is given by Eq. (1), τ_e-e is the relaxation time of e-e thermalization, and $f_{mk}^0$ is the thermalized FD distribution, which can be obtained by fitting the non-thermal distribution to a thermalized one, with energy and electron number conserved through the chemical potential:

(8)$\sum _{mk}({\epsilon }_{mk}-{\mu }^{t+\Delta t_e})f_{mk}^0(T_e^{t+\Delta t_e},{\mu }^{t+\Delta t_e})=\sum _{mk}({\epsilon }_{mk}-{\mu }^t)f_{mk}$

(9)$\sum _{mk}{f}_{mk}^0(T_e^{t+\Delta t_e},{\mu }^{t+\Delta t_e})=\sum _{mk}{f}_{mk}$

where ε_mk is the energy of electron states, μ is the chemical potential, Δt_e is the electronic iteration time step, and T_e is the average temperature of electrons. Eq. (7) is solved by employing the Lie-Trotter splitting scheme in combination with the Euler method.

Importantly, to compare with previous TDBE studies without non-thermal electrons, we consider resetting the electron distribution to the FD form immediately after each electronic time step (labeled as τ_e-e = 0). In this case, Eq. (7) is expressed as

(10)$f_{mk}=f_{mk}^0$

where $f_{mk}^0$ is determined by Eqs. (8) and (9). Moreover, to compare with the TTM, we further allow the phonon distribution to relax to the BE form after each iteration:

(11)$\sum _{sq}ħ{\omega }_{sq}(n_{sq}+\frac{1}{2})=\sum _{sq}ħ{\omega }_{sq}[n_{sq}^0\left(T_p^{t+\Delta t_p}\right)+\frac{1}{2}]$

(12)$n_{sq}=n_{sq}^0$

where ω_sq is phonon energy, T_p is the average temperature of phonons, and $n_{sq}^0$ is the thermalized BE distribution. Thus, this framework enables us to study the influence of the electron thermalization timescale on the relaxation dynamics, while at the same time offering an approach to assess the applicability of the TTM.

2.2 Electron-phonon coupling constant

In the TTM framework, the energy exchange between electrons and phonons is described by

(13)$C_e\frac{d{T}_e}{dt}=-G_{\mathrm{ep}}(T_e-T_p)$

where C_e is the electronic heat capacity, T_e and T_p are the electron and phonon temperatures, respectively. The left-hand term can be expressed as the rate of change of the total energy of electrons

(14)$C_e\frac{d{T}_e}{dt}={\partial }_t{E}_e=\sum _{mk}({\epsilon }_{mk}-\mu ){\partial }_t{f}_{mk}$

where ε_mk is the energy of electron states, μ is the chemical potential, and ${\partial }_t{f}_{mk}$ is given by Eq. (7). Thus, we can plot the energy transfer rate from electrons to phonons as a function of the average temperature difference (T_e - T_p), and the G_ep is determined as the slope of a linear fit.

2.3 Mode-dependent effective electron temperature

To quantitatively characterize the deviation of nonequilibrium electrons from the FD distribution, a mode-dependent effective electron temperature (T_mk) is introduced

(15)$f_{mk}={[\exp \left(\frac{{\epsilon }_{mk}-\mu }{k_B{T}_{mk}}\right)+1]}^{-1}$

where ε_mk is the energy of electron states and μ is the chemical potential. When the electronic distribution approaches the FD distribution, T_mk for all states converges to a single value. In contrast, away from thermal distribution, T_mk becomes broadened, and the extent of this broadening can serve as a measure of deviation from the thermalized distribution.

2.4 Workflow

As shown in Figure 1, a high-temperature FD distribution for electrons and a low-temperature BE distribution for phonons are used as the initial distributions to calculate the A_mk and B_mk, which are implemented by the EPW package^[35]. Next, we solve for the electron distribution function using the Euler method. Considering that the relaxation time of ph-ph scattering is several orders of magnitude longer than that of e-ph scattering, we first iterate the electron distribution function for 10-100 fs before starting the iteration of the phonon distribution. The calculation of the A_sq and B_sq employs both the EPW and ShengBTE packages^[36].

Display Full Size

Figure 1. Flow chart of solving the TDBT. TDBT: Time-Dependent Boltzmann Equation.

2.5 Computational details

All electronic structures were computed using the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional, using fully relativistic norm-conserving pseudopotentials, as implemented in the Quantum Espresso package^[37-39]. A 90 Ry cutoff energy was adopted, and the Brillouin zone was sampled with a 10 × 10 × 10 Monkhorst-Pack grid, with Gaussian smearing with 0.01 Ry. The phonon dispersion was computed using density functional perturbation theory and a 5 × 5 × 5 q-grid^[40]. The e-ph scattering matrix elements on a dense mesh of 20 × 20 × 20 k-points and 20 × 20 × 20 q-points were obtained using the Bloch-Wannier-Bloch transformation, with 16, 18, and 24 maximally localized Wannier functions for gold, silver, and aluminum, respectively^[41,42]. The fitted electronic and phononic band structures are shown in Figure S1,S2,S3, and the convergence of the relaxation dynamics with k-/q- grid in gold is shown in Figure S4. The electron distribution was iterated using the Lie-Trotter splitting scheme and the explicit Euler method with a time step (Δt_e) of 2 fs for gold and aluminum, and the implicit counterpart with Δt_e of 10 fs for silver. The third-order force constants were computed using the Vienna Ab initio Simulation Package and the Thirdorder package on a 3 × 3 × 3 supercell^[36,43]. An energy cutoff of 600 eV was used with PAW-PBE pseudopotentials^[37,44]. The ph-ph scattering matrix elements were calculated using the ShengBTE package^[36]. The phonon distribution was iterated using the explicit Euler method with a time step (Δt_p) of 0.02 ps for aluminum, and 0.1 ps for gold and silver.

3. Results and Discussion

3.1 Out-of-equilibrium dynamics in gold

To compare with the results of time domain thermoreflectance (TDTR) experiments, we compute the dynamics of the electron and phonon systems with an initial T_e of 2,600 K and T_p of 300 K^[45]. During the relaxation process, the concept of temperature is not directly applicable to non-equilibrium distributions. Therefore, the average electron and phonon temperatures are obtained by fitting the thermalized distributions, with energy conservation and electron number conservation enforced via the chemical potential (see Eqs. (8), (9), and (11)). Figure 2a shows the time evolution of T_e and T_p in gold under three conditions: τ_e-e = 40 fs, τ_e-e = 0 (instantaneous electron thermalization), and the TTM (instantaneous electron and phonon thermalization). The comparison reveals that whether the electron distribution is constrained to a FD distribution or the phonon distribution to a BE distribution, the overall relaxation behavior remains qualitatively similar. The main difference is the timescale for reaching equilibrium: in the TTM, the system equilibrates within approximately 12 ps, whereas in the other two cases, equilibrium is delayed to several tens of picoseconds due to the slower ph-ph scattering process (Figure S5).

Display Full Size

Figure 2. (a) Electron T_e (solid lines) and phonon T_p (dashed lines) temperatures during relaxation in gold; (b) Energy transfer rate from electrons to phonons as a function of the temperature difference. The G_ep in the TTM is determined from the slope of a linear fit. TTM: two-temperature model.

In Figure 2b, we plot the energy transfer rate from electrons to phonons as a function of the temperature difference between the two subsystems. The slope corresponds to G_ep as defined in the TTM. The fitted G_ep is 2.76 × 10¹⁶, 2.84 × 10¹⁶, and 2.79 × 10¹⁶ W/m³/K for τ_e-e = 40 fs, τ_e-e = 0, and TTM, respectively. These results are slightly larger than reported values (~2.2 × 10¹⁶ W/m³/K)^[46-48]. To examine this discrepancy, we compute the electronic heat capacity as a function of T_e during relaxation (Figure S6). The fitted electron specific heat constant (76.3 J/m³/K²) also slightly exceeds experimental measurements (67.6 J/m³/K²)^[49-51], indicating an overestimation of the electronic eigenvalues. Nevertheless, all deviations remain below 15%, supporting the accuracy of our calculations.

The difference between τ_e-e = 0 and TTM is relatively small, even though the phonons are confirmed to remain in a non-equilibrium state within 12 ps (Figure 3a), consistent with previous reports^[20,33]. Surprisingly, the results at τ_e-e = 40 fs also exhibit no significant deviation, which seems to indicate that the impact of non-thermal electrons on the dynamics may be negligible. To address this issue, we compare the rate of change of the electronic occupation and the distribution of T_mk at different times (Figure 3b,c). The rate of change in electron occupations for τ_e-e = 0 is nearly identical to that in the TTM case and slightly larger than that for τ_e-e = 40 fs. This indicates that when the electron distribution is constrained to a FD distribution, the influence of the non-thermal phonon distribution on the overall dynamics is minimal. Moreover, after the first Δt_p (0.1 ps), the distribution of T_mk already exhibits a noticeable broadening, indicating that electrons starting from a FD distribution quickly evolve into a non-thermal state through e-ph scattering. After 5 ps, the broadening becomes almost invisible, which is fully consistent with the assumptions of the TTM. These behaviors of electrons can be attributed to the weak e-ph coupling, as well as the lack of strong scattering channels. As a result, over 70% of the electrons have T_mk within 25 K of the average temperature (Figure 3d). Considering only electronic states within a ±2 eV window around the chemical potential, this value is already high, and it would increase further if a larger energy range were included. This indicates that the average temperature is representative of the collective behavior of most electrons.

Display Full Size

Figure 3. (a) Phonon distribution and (b) rate of change of the electron distribution at different times during relaxation in gold. The chemical potential is set to zero; (c) Effective temperature distribution of electronic states within ± 2 eV around the chemical potential at different times during τ_e-e = 40 fs relaxation in gold. The vertical axis represents the probability of electronic states falling within each 50 K temperature interval. The red dashed line indicates the current average temperature; (d) Time evolution of the probability that electronic states within ± 2 eV of the chemical potential fall within ± 25 K of the average temperature in gold, under the condition of τ_e-e = 40 fs.

Since gold is a case where the TTM remains valid, we assume instantaneous electron thermalization (τ_e-e = 0) to compute the conductivity, aiming to verify the accuracy of our calculations. The electron absorption and emission coefficients at each time step yield the instantaneous scattering rates, which are then used in the Boltzmann transport expression within the RTA to obtain the conductivity evolution. Experimentally, the conductivity evolution has been measured for initial T_e exceeding 16,000 K^[52]. Such high T_e in DFT calculations significantly increase the number of electronic states within the Fermi window, leading to a greater number of e-ph scattering processes and an increase in G_ep by an order of magnitude^[51,53]. At the same time, it can induce lattice melting, and many works have shown that in this regime the G_ep obtained from theory deviates markedly from experimental values^[54,55]. Thus, the initial T_e and T_p were set to 6,800 K and 300 K, respectively. The resulting time dependence of the temperatures and conductivity is shown in Figure 4. The T_p reaches the melting point of gold (~1,337 K) at 13.5 ps. The conductivity is in good agreement with a bi-exponential decay model, with relaxation times of 0.6 ps and 2.9 ps, respectively. The magnitude and trend of the conductivity are consistent with experimental observations, however, due to differences in the initial conditions, our results still require further validation^[52].

Display Full Size

Figure 4. (a) Electron (T_e) and phonon (T_p) temperatures and (b) electrical conductivity during relaxation, with initial electron and phonon temperatures of 6,800 K and 300 K, respectively. The table in (b) summarizes the fitting formula, parameters, and their standard errors.

3.2 Out-of-equilibrium dynamics in silver

The instantaneous electron thermalization reproduces the relaxation dynamics in gold owing to its weak e-ph coupling. However, this agreement is not generally applicable to other metals. We next examine the case of silver. To enable a direct comparison with gold, the same initial temperatures (2,600 K for electrons and 300 K for phonons) were assigned to the silver system. In Figure 5a, the relaxation processes for τ_e-e = 0 and the TTM are nearly identical, with the system approaching equilibrium after approximately 7 ps. In contrast, for τ_e-e = 40 fs, relaxation is much slower, reaching equilibrium around 12 ps. This indicates that when the electronic distribution is far from equilibrium, the energy transfer rate and the corresponding G_ep are reduced, as evidenced in Figure 5b. The slope remains nearly constant for the τ_e-e = 0 and TTM cases (corresponding to 4.47 × 10¹⁶ W/m³/K and 4.52 × 10¹⁶ W/m³/K, respectively), whereas for τ_e-e = 40 fs, a drop appears when the temperature difference is around 2,200 K. This suggests that using FD and BE distributions may overestimate the energy transfer rate, and as the electrons deviate from equilibrium, the G_ep converges to a lower constant value (2.27 × 10¹⁶ W/m³/K). The reported G_ep for silver (2.5-3.5 × 10¹⁶ W/m³/K^{[28,51,56,57]}) is between the results for τ_e-e = 0 and τ_e-e = 40 fs, suggesting that the actual electron thermalization rate is not instantaneous, but still faster than what is assumed in our τ_e-e = 40 fs case.

Display Full Size

Figure 5. (a) Electron T_e (solid lines) and phonon T_p (dashed lines) temperatures during relaxation in silver; (b) Energy transfer rate from electrons to phonons as a function of the temperature difference. The G_ep is determined as the slope of a linear fit.

To further understand these results, we present the phonon distribution functions, the change rate of the state-resolved electron occupations, and the distribution of T_mk in Figure 6. A large difference is observed between τ_e-e = 40 fs and τ_e-e = 0 (Figure 6a). This behavior was not seen in gold (Figure 3a). More notably, phonons across the whole frequency spectrum deviate from equilibrium in silver, whereas the deviations are primarily limited to high-frequency modes in gold. This suggests that silver possesses a broader e-ph scattering phase space and more complex scattering channels, rendering the energy transfer process more sensitive to the details of the electronic distribution. As for the electrons, the change rates of nearly all electrons in the τ_e-e = 40 fs case are around one order of magnitude smaller than those in the other two cases (Figure 3b). The comparison between τ_e-e = 0 and the TTM shows nearly identical electron occupation change rates, consistent with our observations in gold. This reinforces the conclusion that, if electrons thermalize instantaneously, the influence of the phonon non-thermal distribution on the relaxation dynamics is negligible. Moreover, when comparing τ_e-e = 40 fs to τ_e-e = 0, we observe that the sign of the change rate for some electron states even reverses, indicating a shift in the dominant e-ph scattering mechanism, from emission to absorption or vice versa.

Display Full Size

Figure 6. (a) Phonon distribution and (b) change rate of electron distribution at different times during relaxation in silver. The chemical potential is set to zero; (c) Effective temperature distribution of electronic states within ± 2 eV around the chemical potential at different times during τ_e-e = 40 fs relaxation in silver. The vertical axis represents the probability of electronic states falling within each 50 K temperature interval. The red dashed line indicates the current average temperature; (d) Time evolution of the probability that electronic states within ± 2 eV of the chemical potential fall within ± 25 K of the average temperature in silver, under the condition of τ_e-e = 40 fs.

In Figure 6c,d, the broadening of T_mk is much larger than that in gold, and a pronounced broadening is still observed even at 5 ps, by which time the TTM predicts the system to be nearly equilibrated, thereby leading to a significant deviation. The reason is that certain strongly e-ph scattering channels cause the corresponding electronic states to deviate significantly from equilibrium, and electronic thermalization cannot be completed prior to e-ph scattering due to the comparable timescales of e-ph and e-e scattering. Thus, although silver and gold exhibit similar G_ep values, indicating weak e-ph coupling in both systems, the contrasting T_mk broadening reveals that the fundamental e-ph scattering mechanisms are different.

3.3 Out-of-equilibrium dynamics in aluminum

Unlike gold and silver, aluminum exhibits much stronger e-ph coupling, with the experimental G_ep exceeding that of gold and silver by an order of magnitude^[57,58]. Moreover, the experimentally determined G_ep shows a noticeable deviation from TTM predictions^[59]. Herein, we examine the relaxation process for initial T_e of 1,900 K and T_p of 300 K (Figure 7). The relaxation processes for τ_e-e = 0 and the TTM are in close agreement, whereas τ_e-e = 40 fs exhibits a slower relaxation, consistent with the observations in silver. Owing to the stronger e-ph coupling in aluminum, the T_e approaches its equilibrium value within 0.5~1.5 ps, markedly faster than in gold and silver. The τ_e-e = 0 and TTM results agree with the calculations of Ritzmann et al., in which the T_e was found to peak at 1,800 K and approach equilibrium within ~0.7 ps^[20]. In contrast, TDTR measurements revealed that the T_e relaxed from about 1900 K to equilibrium over ~1.5 ps, consistent with our τ_e-e = 40 fs results but significantly longer than the TTM prediction^[45]. In Figure 7b, the fitted G_ep values are 2.7 × 10¹⁷, 5.8 × 10¹⁷, and 6.1 × 10¹⁷ W/m³/K for τ_e-e = 40 fs, τ_e-e = 0 and TTM, respectively. The G_ep for τ_e-e = 0 and TTM are significantly larger than that for τ_e-e = 40 fs, and the latter is consistent with the experimental result of 2.5 × 10¹⁷ W/m³/K^[58].

Display Full Size

Figure 7. (a) Electron T_e (solid lines) and phonon T_p (dashed lines) temperatures during relaxation in aluminum; (b) Energy transfer rate from electrons to phonons as a function of the temperature difference. The G_ep is determined as the slope of a linear fit.

Notably, the relaxation process in aluminum reveals that not only non-thermal electrons reduce the energy transfer efficiency, but non-thermal phonons also introduce a considerable effect compared to the TTM prediction, which is not observed in gold or silver. To understand this result, we present the phonon distribution functions and the change rate of the state-resolved electron occupations in Figure 8. For τ_e-e = 40 fs and τ_e-e = 0, high-frequency phonons are rapidly excited after 0.1 ps (Figure 8a). At 0.3 ps, the low-frequency phonons also exhibit significant deviations from a thermalized BE distribution, indicating that ph-ph scattering has already become the primary mechanism governing internal energy redistribution. At 0.6 ps, the phonon distribution over the entire spectrum exhibits an obvious non-thermal character.

Display Full Size

Figure 8. (a) Phonon distribution and (b) change rate of electron distribution at different times during relaxation in aluminum. The chemical potential is set to zero.

Moreover, Figure 8b shows that the electron distribution change for τ_e-e = 40 fs is about one order of magnitude smaller than that for τ_e-e = 0 case after 0.02 ps. The peak positions of the distribution change in the two cases remain nearly identical, indicating that the dominant scattering channels are robust. When electrons are assumed to thermalize instantaneously (τ_e-e = 0), phonon emission through these channels proceeds without restriction. In contrast, a finite thermalization rate leads to a bottleneck in the phonon excitation process, as the initial electron states become depleted while the final states become occupied. Consequently, phonon excitation is significantly faster in the τ_e-e = 0 case compared to the τ_e-e = 40 fs case (Figure 8a).

In addition, for τ_e-e = 40 fs, certain phonon modes exhibit occupations lower than the thermalized BE distribution, which is absent in the τ_e-e = 0 case. This indicates that the effect is not driven by ph-ph scattering but instead originates from the non-thermal electrons. The non-thermal electrons activate additional e-ph scattering channels that facilitate a reverse flow of energy from phonons back to electrons. At 0.6 ps, the emergence of extra high-frequency peaks in the phonon spectrum for τ_e-e = 40 fs arises from these additional scattering channels (Figure 8a). The resulting changes in high-frequency phonons further affect the low-frequency modes through ph-ph interactions. Therefore, non-thermal electrons exert a significant influence on the microscopic scattering dynamics, resulting in the distinctly different relaxation behaviors between τ_e-e = 40 fs and τ_e-e = 0 cases.

On the other hand, when comparing τ_e-e = 0 with the TTM case, Figure 8b shows that the electron distribution changes remain nearly identical at 0.02 ps, since the phonon subsystem is still close to equilibrium at this early stage. However, at 0.3 ps, the change rate of electron distribution for τ_e-e = 0 becomes smaller than that predicted by the TTM. At this time, the high-frequency phonons hold much higher occupations than those required for a thermal phonon distribution in the TTM. According to Eqs. (3) and (4), the electron emission and absorption coefficients (A_mk and B_mk) scale positively with the phonon occupation. Consequently, non-thermal high-frequency phonons enhance both coefficients. For the unoccupied electronic states above the Fermi level, where emission processes dominate, this enhancement increases the rate of change in the electron distribution. In contrast, for the occupied states at or below the Fermi level, the rate decreases. In other words, the presence of non-thermal high-frequency phonons tends to reheat electrons, thereby suppressing the net energy transfer from electrons to phonons. This result is similar to the recently reported hot electron dynamics modulated by non-equilibrium phonon excitations in semiconductors^[60]. As for gold, although non-thermal phonons also emerge during relaxation, the much weaker e-ph coupling strength limits the feedback to electrons, resulting in only minor deviations from the TTM behavior.

Finally, to further compare aluminum with gold and silver, Figure 9a shows the evolution of the distribution of T_mk. The broadening in aluminum is even more substantial than that in silver, and the probability remains below 0.3 for about 0.5 ps, corresponding to about half of the entire relaxation process (Figure 9b). This pronounced and long-lasting non-thermal electronic feature confirms that the assumption of a single electronic temperature, as adopted in the TTM, fails to capture the ultrafast relaxation dynamics in aluminum. Importantly, the deviation cannot be fully corrected by incorporating non-thermal phonons alone. Both electron and phonon non-equilibrium effects must be considered simultaneously. Our comparison across gold, silver, and aluminum demonstrates that the electron thermalization rate strongly affects the effective coupling constant. If electrons were assumed to thermalize instantaneously, the feedback from non-thermal phonons on G_ep would be largely suppressed. Therefore, explicitly treating the interplay between non-thermal electrons and phonons is essential for accurately describing ultrafast energy relaxation.

Display Full Size

Figure 9. (a) Effective temperature distribution of electronic states within ± 2 eV around the chemical potential at different times during τ_e-e = 40 fs relaxation in aluminum. The vertical axis represents the probability of electronic states falling within each 50 K temperature interval. The red dashed line indicates the current average temperature; (b) Time evolution of the probability that electronic states within ± 2 eV of the chemical potential fall within ± 25 K of the average temperature in aluminum, under the condition of τ_e-e = 40 fs.

4. Conclusion

In summary, we performed first-principles TDBE simulations that explicitly account for both non-equilibrium electrons and phonons, and compared the results with the TTM. We have demonstrated that the primary source of deviations from the TTM lies in the non-thermal electronic distribution and its thermalization dynamics, rather than in non-thermal phonons. In gold, the relaxation dynamics are well captured by the TTM, reflecting a rare case where the absence of strongly coupled scattering channels keeps the electronic distribution nearly thermalized. In silver, despite a coupling strength comparable to gold, pronounced deviations emerge once non-equilibrium electrons are included, revealing that specific e-ph scattering channels, not merely the overall coupling strength, govern the breakdown of the TTM. In aluminum, both non-equilibrium electrons and phonons strongly affect the relaxation process, rendering the TTM entirely invalid. These insights emphasize that explicitly accounting for electronic thermalization is essential for a reliable description of ultrafast relaxation processes. We further suggest that the electronic non-thermal distribution may be a central reason for the persistent discrepancies between theory and low-fluence pump-probe experiments, thereby pointing to a promising direction for future research.

Supplementary materials

The supplementary material for this article is available at: Supplementary materials.

Acknowledgments

The authors gratefully acknowledge Meng Lian, and Wen-Hao Mao for valuable discussions and helpful insights.

Authors contribution

Gui-Lin Zhu: Data analysis, conceptualization, writing-original draft, writing-review & editing.

Jing-Tao Lü: Supervision, methodology, resources, writing-review & editing.

Conflicts of interest

Jing-Tao Lü is an Editorial Board Member of Thermo-X. The other author declared that there is no conflict of interest.

Ethical approval

Not applicable.

Consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

The data and materials could be obtained from the corresponding author upon request.

Funding

This work is supported in part by the National Natural Science Foundation of China (Grant No. 22273029).

Copyright

© The Author(s) 2025.