1. Introduction

Femtosecond optical pump-probe measurements became a standard technique for accessing the ultrafast dynamics of quasiparticle excitation and relaxation in the solid state [1–3], correlated systems [4,5], and warm dense matter [6–10]. The way in which a measured transient optical response reflects the relaxation dynamics of excited matter is often indirect. The state-of-the-art approach for acquiring this knowledge is by performing direct quantum mechanical simulations of a pump-probe experiment using time-dependent density functional theory (TD-DFT) [11,12]. However, due to the significant complexity of this method, researchers often employ the two-temperature model (TTM) [13,14] or its extensions [15–18] for studying temperature dynamics, and model the optical response as a function of transient temperatures.

At low excitation energies, a linear relation between the measured signal and electron and lattice temperatures is often assumed [2]. However, in a general case of arbitrarily strong excitation, such an assumption does not necessarily hold, and a detailed understanding of how optical properties evolve with temperature is required. Analytical models allow for the treatment of simple [19,20] and noble metals [21,22], but most materials require DFT simulations at finite temperatures. The temperature-dependent optical response of systems with a bandgap is modeled via the Bethe-Salpeter equation approach, allowing for electron-hole interaction [23]. Density functional molecular dynamics (DFT-MD) combined with the Kubo-Greenwood theory has proven to be highly successful in describing systems where single-particle excitations dominate, such as metals [19,20,24] and dense plasmas [25–27].

In this work, we investigate the electron-temperature-dependent optical dielectric function of gold and ruthenium through DFT-MD simulations and the independent-particle approximation for the dielectric function. Gold has been extensively studied in this regard [7,24,28–31], which makes it an ideal material for testing our methodology at conditions typical for experiments on ultrafast laser excitation. On the example of gold we demonstrate that, although the employed approach is equivalent to the Kubo-Greenwood theory, it offers a significant advantage: an accurate description of the real part of the dielectric function without the necessity to simulate large systems of atoms, thanks to the explicit consideration of the Drude contribution to the dielectric function. It also offers an explicit treatment of electron-electron scattering, which becomes important at high temperatures. The trade-off, however, is the requirement for a priori knowledge of the damping parameter entering the Drude part.

In contrast to gold, ruthenium received less attention in previous studies. Ruthenium is a transition metal with a complicated band structure and strong electron-phonon coupling (an order-two higher than in gold [32]), which makes studying its ultrafast dynamics challenging. Extremely fast energy transfer between hot electrons and yet cold lattice makes it difficult to unravel the interplay between the two subsystems via standard optical pump-probe experiments. Ruthenium transient response was the subject of ultrafast laser experiments, e.g. [33–35], but the theoretical analysis was limited due to the absence of an accurate model for optical properties in the two-temperature regime. To the best of our knowledge, no attempts were made in the past to build such a model for ruthenium either analytically or via first-principles simulations. The resulting dielectric function of ruthenium, in combination with the aforementioned TTM, enables us to explore the transient optical response in ruthenium and compare it to recent measurements of ultrafast-laser-heated ruthenium thin films [34].

2. Model

2.1 Simulation technique

We performed first-principles simulations of the optical properties of gold and ruthenium using the Quantum Espresso package [36]. We employed the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [37] and scalar-relativistic norm-conserving pseudopotentials from the PSEUDODOJO database [38] explicitly treating 19 electrons in Au and 16 in Ru. The simulation workflow was following. First, we performed a geometry optimization procedure to get relaxed cell parameters. Next, we set up a molecular dynamics (MD) simulation of an orthogonal $2\times 2\times 2$ supercell containing 32 atoms at an ion temperature of 300 K on a coarse $2\times 2\times 2$ k-point grid. We performed MD simulation with a timestep $dt = 0.5\,\mathrm{fs}$ and a total of 2000 timesteps to ensure the supercell does not experience any significant jumps in total energy or temperature and remains structurally stable. Starting from the 1000th step, we extracted several independent ionic configurations from the MD trajectory and performed accurate simulations of their optical properties using the SIMPLE code [39], which is based on the Shirley interpolation scheme and specifically designed for simulating optical properties.

The simulation of optical properties was performed in four stages: (i) self-consistent simulation of the electron density on a $16\times 16\times 16$ k-point grid for Au, and a $24\times 16\times 16$ k-point grid for orthogonalized Ru supercells, (ii) non-self-consistent simulation of the electron wavefunctions at the $\Gamma$-point and its seven periodic images at the corners of the unit cube (see more details in [40]), (iii) construction of the optimal basis set and calculation of velocity matrix elements on the k-point grid used in step (i), (iv) interpolation of matrix elements onto a fine k-point grid twice larger in every dimension than the grid used in steps (i), (iii), and calculation of the optical dielectric function. We used 80 Ry plane wave cutoff in steps (i) and (ii), 720 bands for Au and 500 bands for Ru in step (ii). To take into account the finite electron temperature, Fermi-Dirac smearing was applied.

Additionally, we performed density functional perturbation theory (DFPT) simulations of the electron-phonon (Fan-Migdal) self-energy within the ABINIT package [41,42] using the same parameters as in our previous work [32], and used it as an intraband damping parameter in step (iv).

2.2 Dielectric function

The target quantity of our first-principles simulations is the electron-temperature-dependent complex dielectric function $\varepsilon (\omega )$, which defines optical properties such as complex refraction index $n+ik$, absorption coefficient $\alpha$, and reflectivity $R$ of a given material [43]:

For a metal, the dielectric function can be conveniently written in terms of interband and intraband (Drude) contributions [43]:

Within the Kohn-Sham quasiparticle picture and the independent-particle approximation, they are represented as follows [39]:

The dielectric function $\varepsilon (\omega )$ given by Eqs. (2)-(4) includes the dependence on electron temperature $T_e$ in several ways: energies $\epsilon _{n\mathbf {k}}$ experience a slight shift with increasing $T_e$, Fermi occupations and their derivative are functions of $T_e$, and $\gamma$ may vary with $T_e$. We consider electron-phonon (eph) and electron-electron (ee) scattering channels as the main contribution to the finite broadening of conduction electrons, $\gamma = \gamma ^{\mathrm {eph}} + \gamma ^{\mathrm {ee}}$.

The broadening arising from electron-phonon scattering can be extracted from the imaginary part of the Fan-Migdal self-energy $\Sigma ^{\mathrm {FM}}$ [44]:

Eq. (6) is written assuming thermal equilibrium between electrons and lattice. Formally, in the two-temperature case of interest $T_e \gg T_l$, one must extend the Fan-Migdal self-energy to nonequilibrium regime [45] and evaluate it for several electron temperatures, which is a computationally expensive task. In practice, the self-energy is often weakly dependent on $T_e$ [28,46], and thus can be computed once for normal conditions.

However, for electron-electron scattering, one must account for its dependence on the electron temperature. This can be done, similar to the electron-phonon contribution, by calculating the imaginary part of the GW self-energy over a wide range of electron temperatures, which again requires significant computational resources. Instead, we use temperature-dependent $\gamma ^{\mathrm {ee}}$ values obtained from the kinetic theory for Au [47], and Ru [48].

3. Optical and electronic properties at normal conditions

In this section, we present the simulated dielectric function at normal conditions $T_e = T_l = 300\,\mathrm {K}$ and compare it to available experimental data. We also discuss the electron transport properties such as electron relaxation time, DC conductivity, and resistivity, since they can be easily obtained from our simulations. Additionally, using the example of gold, we discuss the advantages of the method we used compared to the Kubo-Greenwood (KG) formalism (see, e.g., [24,49,50]) – another popular method for the first-principles simulation of optical properties derived from the real part of dynamical conductivity. While it considers the electron-ion scattering contribution without the manual inclusion of a damping parameter $\gamma$, there is a continuous discussion in the community whether the electron-electron scattering is taken into account or not [26,51–53]. The recent conclusion is that, although Kubo-Greenwood includes effects of electron-electron correlations through an exchange-correlational functional, it does not directly account for electron-electron scattering [27].

Figure 1 shows the real and imaginary parts of the optical dielectric function in gold. In the region of photon energies up to 1.5 eV, where the intraband part makes a significant contribution, our simulations are in very good agreement with experimental data. The diverging intraband tail in the $\operatorname{Im} \varepsilon$ follows by a pronounced valley arising from a gap in the positions of d-band and the Fermi level in gold. Because of this gap, dipole transitions from the d-band are suppressed under normal conditions. At higher energies, when d-s interband transitions contribute, we see a larger deviation between simulations and experiment. We associate this with an underestimation of the gap between the d-band and the Fermi level at the PBE level, resulting in the appearance of an interband hill in the $\operatorname{Im} \varepsilon$ at lower energies than in the experiment (see Fig. 1(b)). The more narrow and intense shape of this hill is also arguably due to inaccuracies of the d-band form imposed by the use of pure PBE. The usage of hybrid functionals or GW corrections for obtaining accurate band structure is beyond the scope of this work, particularly since at high electron temperatures errors associated with a choice of exchange-correlational functional become smaller.

 

Fig. 1. Real (a) and imaginary (b) parts of the dielectric function, real (c) and imaginary (d) parts of complex refractive index in gold. Red line corresponds to the results of the present simulations, blue line with circles is the experimental data of Johnson and Christy [54], purple line with squares is the measurements of Olmon et al. for a single-crystal gold sample [55].

Download Full Size | PDF

Surprisingly, the refractive index $n$ and the extinction coefficient $k$, appear to be in better agreement with the experimental values (see Fig. 1(c)-(d)). Indeed, the interband transition threshold in $n$ remains unchanged, but the absolute values of $n$ and $k$ are in much better agreement with the experiment. We believe this is due to the mutual error cancellation in the real and imaginary parts of the dielectric function.

The real and imaginary parts of the dielectric function in ruthenium are shown in Fig. 2(a)-(b). Our simulations demonstrate qualitatively the same behavior of $\operatorname{Im} \varepsilon(\omega)$ in Ru as in Au: there is a region of dominating intraband contribution below 0.5 eV followed by a strong interband transition peak around 2 eV. However, the valley in $\operatorname{Im} \varepsilon(\omega)$ separating these two regions is less pronounced in Ru than in Au, because the Fermi level crosses the overlapping d- and s-bands in Ru. Under normal conditions, d-s interband transitions in this region are not suppressed, but their relative intensity is weak due to the pseudogap-like feature in the Ru electron density of states (DOS) around the Fermi level. Simulations suggest several minor interband peaks below 1 eV and above 3 eV, which do not appear in the experiment. Unfortunately, we cannot unambiguously answer if these minor peaks are the simulation artifact or not, because the available experimental optical spectrum of Ru is quite outdated and there is no other data to compare with to the best of our knowledge. Similarly to the gold case, in Fig. 2(c) we see that the shape of refractive index $n$ is determined by the shape of $\operatorname{Im} \varepsilon(\omega)$ and hence also has minor interband peaks not resolved by the experiment. The amplitudes of $n$ and $k$ agree well with the experimental ones.

 

Fig. 2. Real (a) and imaginary (b) parts of the dielectric function, real (c) and imaginary (d) part of complex refractive index in ruthenium. Red line corresponds to the results of the present simulations, blue line with circles is the data taken from Palik’s handbook [56].

Download Full Size | PDF

Having obtained $\omega _p$ as an output of the SIMPLE simulations, we could easily calculate the DC conductivity $\sigma _0 = \tau \omega _p^2 / 4\pi$ and resistivity $\rho = 1/\sigma _0$, where the electron relaxation time within SERTA – self-energy relaxation time approximation – is simply $\tau = 1/\gamma ^{\mathrm {eph}}$. We consider that at normal conditions, the electron-electron contribution to the relaxation time is negligible. The results are shown in Table 1. The room-temperature resistivity is slightly overestimated compared to the experiment. Assuming that relaxation times are correct, we deduce that $\omega _p$ values provided by our simulations might be underestimated.

Table 1. Electron-phonon $\gamma ^{\mathrm {eph}}$ broadening, carrier relaxation time $\tau$, DC conductivity $\sigma _0$ and resistivity $\rho$ at normal conditions $T_e = T_l = 300\,\mathrm {K}$ for Au and Ru. $\tau$ values in brackets are obtained from first-principles simulations [57] (for Ru, $\tau = (2\tau _{\perp }+\tau _{\parallel })/3$), $\rho$ values in brackets are from the free-electron gas plasma frequency. Experimental resistivity $\rho _{\mathrm {exp}}$ is taken from [58].

View Table

On the other hand, the plasma frequency of a metal can be simply estimated according to the free-electron gas treatment of metals,

We benchmarked our simulations to KG-based simulations performed by Silaeva et al. [30]. The target quantity of the simulations in the KG formalism is the real part of dynamical conductivity $\sigma (\omega )$, which is directly proportional to the imaginary part of the dielectric function [59]:

The real part of the dielectric function is then restored from Kramers-Kronig (KK) relations [59]:

Figure 3 shows the comparison between our simulations and KG ones. The very good agreement of the imaginary parts of the dielectric function, shown in Fig. 3(b), confirms that both approaches yield identical ${\textrm{Im}} \varepsilon (\omega )$ / ${\textrm{Re}} \sigma (\omega )$ if one uses similar simulation parameters. However, that is not the case for the real part of the dielectric function shown in Fig. 3(a): the results reported by Silaeva differ quite significantly from ours. The KG-based ${\textrm{Re}} \varepsilon (\omega )$ has a shape closer to the interband contribution, not the total one. At first glance, we found this to be controversial: why, with ${\textrm{Im}} \varepsilon (\omega )$ being essentially the same, does the KG $\text{Re} \varepsilon$ capture mostly interband transitions?

 

Fig. 3. Total dielectric function, its interband and intraband contributions, and Kubo-Greenwood-based dielectric function taken from Silaeva et al. [30]. (a) is the real part, and (b) is the imaginary part of dielectric function. Inset in (b) shows the low-energy shoulder of ${\textrm{Im}} \varepsilon (\omega )$.

Download Full Size | PDF

The following reasoning can answer this question. While the KG approach formally includes both intra- and interband contributions to dynamic conductivity, the intraband part converges very slowly with an increase in supercell size and may require simulating up to several thousand atoms [60]. However, the usage of KK relations in Eq. (9) requires highly accurate values of the dynamical conductivity in the DC limit ($\omega \rightarrow 0$) [24], where the intraband contribution dominates. Consequently, although the Kubo-Greenwood formalism provides reasonable optical conductivity data without any a priori knowledge of the broadening parameter $\gamma$, in most cases, it necessitates a significant computational effort to obtain an accurate real part of the dielectric function for materials with a non-vanishing intraband contribution. It follows from all of the above that Silaeva et al. could not reach converged values of DC conductivity, and hence were not able to produce an accurate real part of the optical dielectric function. In contrast, the explicit inclusion of the Drude term in the dielectric function allowed us to achieve reasonably accurate $\varepsilon (\omega )$ without the costly simulations of supercells with thousands of atoms included. However, it comes at the expense of additional simulations of the $\gamma$ value, which in some cases might require a similar or even greater effort, depending on the level of theory introduced to account for broadening effects.

4. Optical properties at high $T_e$

4.1 Influence of electron-electron scattering

Before presenting the optical properties of gold and ruthenium at high electron temperatures, we will discuss the influence of the temperature-dependent damping parameter $\gamma ^{\mathrm {ee}}(T_e)$ on the dielectric function. To do so, we compare ${\textrm{Im}} \varepsilon (\omega )$ of gold at several electron temperatures obtained from our simulations with those extracted from the KG conductivity [30]. As we discussed above, the statement of whether the KG approach accounts for electron-electron scattering or not is highly debatable. Contributing to the discussion, we will consider both scenarios below. The results of the comparison are depicted in Fig. 4.

 

Fig. 4. Imaginary part of the dielectric function in gold, calculated for $T_e = 300\,\mathrm {K}, 10000\,\mathrm {K}, 25000\,\mathrm {K}$. Solid lines correspond to $T_e$-dependent $\gamma = \gamma ^{\mathrm {eph}} + \gamma ^{\mathrm {ee}}$, dashed lines to constant $\gamma = \gamma ^{\mathrm {eph}}$, and thick semi-transparent lines are the KG simulations of Silaeva et al. [30].

Download Full Size | PDF

When electron-electron scattering is disabled, the intraband shoulder aligns well with KG results, especially at not-too-high electron temperatures. At $T_e = 25 \,\mathrm {kK}$, the intraband shoulder obtained using constant damping becomes noticeably narrower than the KG one. The overall good agreement between our simulations and KG approach validates our approximation of constant $\gamma ^{\mathrm {eph}}$ in a wide range of $T_e$ considering that the KG approach does not include electron-electron scattering as well.

It is interesting to see how the inclusion of the $T_e$-dependent electron-electron scattering affects the shape of ${\textrm{Im}} \varepsilon (\omega )$. As expected, under normal conditions with $T_e = 300 \,\mathrm {K}$, electron-electron scattering contributes only insignificantly. However, as the temperature increases, its influence becomes considerable. The low-energy intraband shoulder is roughly twice as wide for $T_e = 10$ and $25 \,\mathrm{kK}$ compared to the predictions of KG.

If we come back to the previous arguments stating that the KG includes the electron-electron scattering, then we must admit that the $\gamma ^{\mathrm {ee}}(T_e)$ obtained from kinetic theory significantly overestimates the intraband contribution. There is an argument that the Lindhard screening used by Petrov et al.[47] might lead to an order of magnitude overestimation of electron-electron damping [61]. Additionally, there could be a double-counting error: coming back to Fig. 3(b), it is apparent that the low-energy shoulder is formed equally by both intra- and interband parts. We speculate that, similarly to KG approach, an interband contribution could solely capture the effect of electron-electron scattering, making its inclusion in the intraband part unnecessary.

For the sake of scientific objectivity, below we present optical properties obtained both with and without $\gamma ^{\mathrm {ee}}(T_e)$ and discuss the implications of this particular choice in comparison to the experimental results for femtosecond-laser-heated Ru in Section 5.

4.2 Optical properties of Au and Ru

Figure 5 (a) and (b) presents the electron-temperature-dependent $n$ and $k$ values in gold. The refractive index $n$ is notably influenced by the specific choice of broadening. With constant $\gamma = \gamma ^{\mathrm {eph}}$ broadening, the infrared intraband shoulder ($\omega \leq 0.5 \,\mathrm {eV}$) remains unchanged with an increase in $T_e$. Changes primarily manifest in the 0.5 – 3 eV region due to the smearing of the Fermi distribution. This smearing opens previously forbidden interband transitions at lower energies ($\omega \lesssim 2 \,\mathrm {eV}$) but reduces the probability of interband transitions at energies above 2 eV. In contrast, for the $T_e$-dependent $\gamma = \gamma ^{\mathrm {eph}} + \gamma ^{\mathrm {ee}}$ broadening, the intraband shoulder intensifies with an increase in $T_e$, causing a notable relative decrease in the d-band contribution in the 0.5-3 eV region. The extinction coefficient $k$ (Fig. 5(b)) exhibits weak dependence on $\gamma$. The $T_e$-dependent broadening results in a slightly narrower shape of the intraband shoulder compared to the constant one only in the infrared region below $1 \,\mathrm {eV}$.

 

Fig. 5. Refractive index $n$ and extinction coefficient $k$ in gold (a),(b), and ruthenium (c),(d) for different values of $T_e$. Lattice temperature is fixed at normal conditions $T_l = 300 \,\mathrm {K}$. Solid semi-transparent lines correspond to temperature-dependent damping $\gamma = \gamma ^{\mathrm {eph}} + \gamma ^{\mathrm {ee}}$, dashed lines – to constant damping $\gamma = \gamma ^{\mathrm {eph}}$.

Download Full Size | PDF

The qualitative picture remains the same for ruthenium, but the dominance of the intraband contribution in the $T_e$-dependent case in the region $\omega \lesssim 2 \,\mathrm {eV}$ is more pronounced for both $n$ and $k$, as depicted in Fig. 5(c) and (d). This is a direct consequence of the overlap between s- and d-band states, resulting in non-vanishing electron-electron scattering of s-electrons even at low temperatures and high $\gamma ^{\mathrm {ee}}$. Such a strong intraband contribution blurs the valley of weak d-s interband transitions within the 0.5-1.5 eV range.

5. Comparison with experiment for femtosecond-irradiated Ru

5.1 Fermi smearing mechanism in Ru

In our recent work [34], we investigated the mechanism of heating and degradation of ruthenium thin films irradiated by a $\omega = 1.55 \,\mathrm {eV}$ femtosecond laser pulse. We discovered that the optical response of heated Ru can be qualitatively described by the Fermi smearing mechanism [62,63], where the peak thermoreflectance signal is proportional to the change of electron occupations with increasing $T_e$:

Here, $T_0 = 300\,\mathrm {K}$, and $\Delta\epsilon$ is the difference between the energy level experiencing the optical transition and the Fermi level of Ru. We treated $\Delta\epsilon$ as the fitting parameter, independent of electron temperature, and thus having the meaning of an averaged energy level from which electrons are excited. By fitting experimental peak thermoreflectance values, we obtained $\Delta \epsilon = -1.2 \,\mathrm {eV}$.

Now, with the simulated $T_e$-dependent optical constants in ruthenium, we can assess the validity of the Fermi smearing mechanism. We calculated the relative reflectivity change using Eq. (1) and compared it to the change in the Fermi distribution given by Eq. (10). The results of the comparison are presented in Fig. 6. The Fermi smearing theory, with $\Delta \epsilon = -1.2 \,\mathrm {eV}$, suggests a steeper growth of reflectivity change with an increase in $T_e$ compared to our simulations, both accounting for and neglecting electron-electron damping. If a smaller $\Delta \epsilon = -1 \,\mathrm {eV}$ is applied, our ’eph+ee’ simulations agree well with the Fermi smearing curve up to $T_e = 5 \,\mathrm {kK}$. Constant ’eph’ damping, on the contrary, suggests slower growth at low temperatures but has the same high-temperature asymptotic and overall increasing-saturation behavior as the change in the Fermi distribution. This is because the constant damping does not influence the valley region in Ru optical spectra associated with transitions in pseudogap-like feature in Ru DOS, where the Fermi smearing mechanism can be applied.

 

Fig. 6. Change of optical reflectivity with an increasing $T_e$. Orange circles are the results of simulation with $T_e$-dependent ’eph+ee’ damping, blue squares – with constant ’eph’ damping. Solid and dashed lines represent Fermi smearing given by Eq. (10) for two different $\Delta \epsilon$.

Download Full Size | PDF

At higher temperatures, the reflectivity change obtained from simulations with enabled $T_e$-dependent electron-electron damping has a peak around $T_e = 7.5 \,\mathrm {kK}$ and then decays with further temperature increase. This deviation from the Fermi smearing theory stems from the extinction of the valley in the optical spectra at high $T_e$. The increasing-saturating-decreasing behavior of reflectivity change might, however, more accurately reflect the experimental observations. In our experiments, we observed a similar behavior in peak thermoreflectance values for a 125 nm thick film (cf. Figure 3 in [34]). We initially attributed a drop in reflectivity to film damage, but in the light of current simulations, it may be associated purely with the response of hot electrons before damage occurred.

5.2 Time-resolved optical response

Another application of our optical constants is the estimation of transient thermoreflectance after ultrafast electron heating using the TTM. This model assumes that instantly thermalized electrons transfer energy to the cold lattice via the electron-phonon scattering mechanism and relax towards equilibrium. The evolution of electron and lattice temperatures is described by two coupled heat equations:

The electron heat capacity $C_e(T_e) = \frac {\partial }{\partial T_e} \int g(\epsilon,T_e) f(\epsilon,T_e)\epsilon d\epsilon$ is calculated from the electron DOS $g(\epsilon,T_e)$ of Ru, as simulated in our recent work [32]. The lattice heat capacity $C_l$ is constant in the temperature regime we are interested in, according to the Dulong-Petit law. We employed a model for the electron thermal conductivity $k_e = k_e (T_e, T_l)$ proposed by Petrov et al. [48]. Two models for the electron-phonon coupling parameter $G$ were used: one obtained from first-principles simulations [32], and another taken from non-adiabatic tight-binding – molecular dynamics (TBMD) approach [64]. We considered a Gaussian laser pulse $S(t,x)$, taking into account the effect of multiple reflections at thin film interfaces.

Figure 7 illustrates an example of two-temperature evolution in a 30 nm Ru thin film irradiated by an incident laser fluence $F = 24.5\,\mathrm {mJ/cm^2}$. Our model predicts a sharp peak in the electron temperature for both models of the electron-phonon coupling parameter. According to the theoretical dependence of $R$ on $T_e$ (Fig. 6), one might expect to observe a similar sharp peak in thermoreflectance values in the experiment. However, experimental values exhibit a noticeably different behavior: an initial sharp increase due to ultrafast energy deposition is followed by a rather flat or slowly decaying trend. We found that the experimental trend closely aligns with the evolution of lattice temperature when the DFT-based $G$ is considered.

 

Fig. 7. Evolution of electron (solid lines) and lattice (dashed lines) temperatures on surface of 30 nm Ru thin film for an incident laser fluence $24.5\,\mathrm {mJ/cm^2}$ obtained from the TTM using two different models for electron-phonon coupling: DFT [32] and TBMD [64]. Black dash-dotted line represents the normalized lattice temperature for DFT-based electron-phonon coupling. Experimental transient thermoreflectance points are shown with gray dots.

Download Full Size | PDF

Such a good agreement between simulated lattice temperature and experimentally measured thermoreflectance trends, along with the absence of an electron temperature peak in the experiment, could be attributed to several scenarios or a combination of them. First, the experimental setup may not have allowed for the resolution of the electron temperature peak due to a timestep of $250 \,\mathrm {fs}$, while the pump pulse duration was about $85 \,\mathrm {fs}$. Second, the electron-phonon coupling might be even stronger than ab initio predictions, and with the pump-probe measurements, we may have accessed an already equilibrated electron-lattice system. Finally, in the range of excitation energies and temperatures reached, the optical properties could depend mostly on the lattice temperature. Below, we will evaluate the second and third scenarios and see if transient thermoreflectance can be successfully modeled by our first-principles optical properties depending on the equilibrium lattice temperature.

In the equilibrium case $T_e = T_l \equiv T$, we cannot ignore the influence of the lattice temperature on optical properties. To account for this, we performed several trial DFT-MD simulations at lattice temperatures up to $1500 \,\mathrm {K}$ and did not find any significant changes in optical constants compared to $T_l = 300 \,\mathrm {K}$. However, at equilibrium and not-too-high temperatures, the electron-phonon damping $\gamma ^{\mathrm {eph}}$ should vary linearly with $T$. We performed simulations of the Fan-Migdal self-energy for $T=100-300 \,\mathrm {K}$ and extracted a linear slope for this dependence, $a = 2.1517 \times 10^{-4} \,\mathrm {eV/K}$.

Figure 8 illustrates the temporal evolution of $R$ in two limits of low and high incident laser fluences. In Fig. 8(a), at a low level of excitation corresponding to peak $T \simeq 1000 \,\mathrm {K}$ in our TTM simulations, our temperature-dependent optical constants successfully overlapped with the experimental points when considering the temperature dependence of all damping channels. The agreement between the experimental data and the simulations becomes unsatisfactory if we neglect the $T$-dependence of $\gamma ^{\mathrm {eph}}$. Contrarily, in the case of high laser fluence (Fig. 8(b)), the scenario changes: optical constants simulated with a constant $\gamma ^{\mathrm {eph}}$ and temperature-dependent $\gamma ^{\mathrm {ee}}$ provide a better match with the experimental points. This discrepancy is attributed to the fact that at high excitation levels (and, consequently, high lattice temperatures, peak $T\simeq 2000 \,\mathrm {K}$), a linear temperature dependence of $\gamma ^{\mathrm {eph}}$ does not hold. Assuming a linear relation between $\gamma ^{\mathrm {eph}}$ and $T$ at high $T$ values results in an overestimation of changes in optical properties. At elevated temperatures, $\gamma ^{\mathrm {ee}}$ becomes dominant, and disregarding the temperature dependence of $\gamma ^{\mathrm {eph}}$ only introduces a minor discrepancy.

 

Fig. 8. Time-resolved thermoreflectance values in 30 nm Ru thin film for incident fluence of $19.89 \,\mathrm {mJ/cm^2}$ (a) and $51.23\,\mathrm {mJ/cm^2}$ (b). Red solid line corresponds to $R$ obtained with an account for equilibrium temperature dependence for both ’eph’ and ’ee’ damping, blue dashed line to constant ’eph’ and temperature-dependent ’ee’ damping, lilac dash-dotted line to constant ’eph’ damping. Experimental data are shown in gray dots.

Download Full Size | PDF

From the analysis provided, it becomes clear that the accurate first-principles description of ultrafast pump-probe experiments in transition metals with very large electron-phonon coupling strength $G \geq 10^{18} \,\mathrm {W/m^3 K}$ is a challenging task. It seems that the two-temperature state occurs only at very short timescales, often inaccessible in pump-probe experiments, or nonequilibrium electrons start to exchange energy with the lattice, completely missing the two-temperature regime. Accounting for the $T_e$-dependent electron-electron damping parameter $\gamma ^{\mathrm {ee}}$ is crucial for an accurate description of optical properties at high levels of excitation, whereas at low levels of excitation, knowledge of the temperature-dependent $\gamma ^{\mathrm {eph}}$ is needed. Overall, a complete understanding of the dynamics of optical properties requires their dependence on both electron and lattice temperatures. Such type of simulations are possible yet extremely demanding, as one needs to span the $(T_e, T_l)$ space using any of first-principles approaches accounting for both electron-electron and electron-phonon scattering.

6. Conclusions

We presented the simulated optical properties of gold and ruthenium in the two-temperature regime $T_e \gg T_l$, which is achievable in ultrafast laser irradiation experiments. Our calculations are carried out within the DFT framework and the independent-particle approximation for the optical dielectric function.

First, to validate our approach, we compared the simulated optical constants and derived resistivity values with experimental data available at normal conditions ($T_e=T_l=300\,\mathrm {K}$). We discussed the similarities and differences between the method we used and the Kubo-Greenwood theory. While the imaginary parts of the optical dielectric function are essentially identical in both methods, our approach provides a more accurate real part without the necessity to simulate huge supercells due to an explicit account for intraband (Drude) contribution. Second, at elevated electron temperatures, We demonstrated that the Kubo-Greenwood approach appears to be incapable of describing electron-electron scattering. In our approach, the accuracy of simulated optical properties at high $T_e$ in the low-energy region ($\lesssim 2\,\mathrm {eV}$) is determined by the accuracy of the chosen electron-electron damping parameter.

We compared simulated optical properties to experimental time-resolved thermoreflectance in Ru thin films [34]. We could reasonably match experimental data with our two-temperature optical constants at high levels of initial electron excitation, whereas at low levels, simulations with equilibrium $T_e = T_l$ were necessary. We attribute this to the very fast thermalization between electrons and lattice in Ru, requiring knowledge of optical properties at elevated lattice temperatures.