Experiments

A 30 nm gold film on top of a glass substrate is considered and pump-probe transient reflectivity measurements are carried out. A 130 fs pump pulse of monochromatic light at λ 0 = 1200 nm (photon energy ħω 0 = 1.0 eV) is followed by a broadband probe pulse in the visible range. In this paper we consider two cases: the pump light incident normally from the air side and a Kretschmann configuration with the pump light incident at the appropriate Kretschmann angle (42.3°) from the glass side. The probe beam is slightly detuned from the Kretschmann angle so that the probe pulse does not excite plasmons and measures only reflectivity. The differential reflectivity, ΔR/R, is then determined as a function of time t and probe frequencies ω = 2πc/λ. The experimental sensitivity is better than a ΔR/R of 10−4, resulting in high quality spectra that are needed for precise results through the double inversion procedure. The skin depth of the Au in this spectral range is ~5 nm, so that all of the optical absorption response of the film can be measured by the ΔR/R response. Figure 1 shows a schematic of the energy band structure in gold and the intraband and interband transitions corresponding to the excitation of the pump and probe pulses, respectively. For additional experimental details, please see the Methods section.

It is important to verify that the effects discussed herein are not due to other, nonlinear optical processes as can sometimes be the case in plasmonic systems21. The analysis of the power dependence of the pump-probe signal verifies that the contribution of multi-photon absorption processes is negligible (see Supplementary Note 9 for detailed analysis of power dependence).

The double inversion procedure

For a given pump-probe set-up, the reflectivity R from the glass-gold structure corresponding to any time t and probe frequency ω can also be viewed to be a function of the metal’s instantaneous complex dielectric function \(\epsilon\) = \(\epsilon \prime\) + i\(\epsilon \prime\prime\). If we denote the unperturbed dielectric function of the metal prior to the pump as \(\epsilon_e = \epsilon \prime_{\!\!e} + i\epsilon \prime\prime_{\!\!\!\!e}\), one can write6

$$\left. {\frac{{{\mathrm{\Delta }}R}}{R} = \frac{{\partial {\kern 1pt} {\mathrm{ln}}{\kern 1pt} R}}{{\partial \epsilon \prime }}} \right|_{\mathrm{e}}{\mathrm{\Delta }}\epsilon \prime + \left. {\frac{{\partial {\kern 1pt} {\mathrm{ln}}{\kern 1pt} R}}{{\partial \epsilon \prime\prime }}} \right|_{\mathrm{e}}{\mathrm{\Delta }}\epsilon \prime\prime ,$$ (1)

where the derivatives are taken at the unperturbed dielectric constant values and \(\epsilon\) = \(\epsilon _e\) + Δ\(\epsilon\) has been assumed, which is valid for the relatively small changes that actually occur in such experiments. These derivatives (which depend on ω) are obtained numerically (see Methods). We can now use the Kramers–Kronig relation to obtain an integral equation just for \(\Delta \epsilon \prime\prime\),

$$\frac{{{\mathrm{\Delta }}R}}{R} = \left. {\frac{{\partial {\kern 1pt} {\mathrm{ln}}{\kern 1pt} R}}{{\partial \epsilon \prime\prime }}} \right|_{\mathrm{e}}{\mathrm{\Delta }}\epsilon \prime\prime + \frac{2}{\pi }\left. {\frac{{\partial {\kern 1pt} {\mathrm{ln}}{\kern 1pt} R}}{{\partial \epsilon \prime }}} \right|_{\mathrm{e}}{\mathrm{PV}}\mathop {\int}\limits_0^\infty {{\mathrm{d}}\omega {\prime}\frac{{\omega {\prime}{\mathrm{\Delta }}\epsilon \prime\prime (\omega {\prime})}}{{(\omega {\prime})^2 - \omega ^2}}}$$ (2)

where PV stands for Principal Value.

The connection between changes in the dielectric function \(\epsilon\) and changes in electron energy distribution is through the energy distribution of the joint density of states (EDJDOS) D(ħω, E)22,23 which is a 2-dimensional distribution function counting transitions due to photons of energy ħω, which excite electrons to final energy of E. Notice that the pump is in the infrared and so its photon energies of 1.03 eV are insufficient to cause interband transitions, however the probe energies assess the gold interband transitions. Thus one has a picture of possibly an initial, broad (±1.03 eV relative to the Fermi energy) initial disturbance in the carrier energy distribution that is subsequently being probed by photons involving primarily interband transitions and so the EDJDOS function needs to involve the interband transitions. This function can be evaluated explicitly for interband transition between the d-band and the s-p conduction band assuming a parabolic structure approximation of the energy bands24,25. In gold the two relevant transitions for our probe photon energies are located around the L and X points in the Brillouin zone and we include contributions from both points in our expression for the EDJDOS (see Supplementary Notes 2 and 3). Specifically,

$$\epsilon \prime\prime (\omega ) = \epsilon \prime\prime _{{\!\!\!\!\mathrm{intra}}} (\omega ) + \frac{A}{{(\hbar \omega )^2}}\mathop {\int}\limits_{E_{{\mathrm{min}}}}^{E_{{\mathrm{max}}}} {\kern 1pt} D(\hbar \omega ,E)(1 - f(E)){\mathrm{d}}E,$$ (3)

where 1 − f(E) is the probability that an upper band state of energy E is empty.

As detailed in the Supplementary Notes 1–6, we determine the parameter A by requiring that Eq. (3) approximate the empirical Johnson and Christy26 data when f(E) is taken to be f e (E), the Fermi-Dirac distribution at room temperature. \(\epsilon \prime\prime _{{\!\!\!\!\mathrm{intra}}}\) is the intraband contribution to the imaginary part of the dielectric constant which we take to be a Drude form fit to the lower frequency empirical data; for optical frequencies \(\epsilon \prime\prime _{{\!\!\!\!\mathrm{intra}}}\) is small relative to the interband contributions. The change Δf(E) = f(E) − f e (E) in the electron energy distribution after the pump is

$${\mathrm{\Delta }}\epsilon \prime\prime (\omega ) = - \frac{A}{{(\hbar \omega )^2}}\mathop {\int}\limits_{E_{{\mathrm{min}}}}^{E_{{\mathrm{max}}}} {\kern 1pt} D(\hbar \omega ,E){\mathrm{\Delta }}f(E){\mathrm{d}}E.$$ (4)

D(ħω, E) and the integral limits are obtained from explicit expressions24,25 using various energies and masses inferred from ref.27 (see also Supplementary Notes 2 and 3). Thus, given \(\Delta \epsilon \prime\prime\)(ω) inferred from solution of integral Eq. (2), integral Eq. (4) can be solved to determine Δf(E).

We should note that our approach is most readily applicable when reasonable analytical approximations for the relevant energy bands are available. Our detailed analysis (see Supplementary Notes 1–6) to obtain the EDJDOS entering in Eqs. (3) and (4) assumed parabolic approximations for the d and conduction (sp) bands and this is appropriate for noble metals Au and Ag, corresponding to Drude-like plasmonic dynamics in the conduction band and a relatively flat d band. In principle, of course, one could extend our approach beyond these metals with appropriate characterization of the relevant material’s band structure. One could also extend our approach beyond flat films to, for example, metal nanoparticles, but one would have to characterize the reflectance entering into Eqs. (1) and (2) with computational electrodynamics calculations.

Validation with model data

To validate the double inversion procedure we specified a hypothetical electron energy distribution change, Δf(E), to be a double-step-like change, Fig. 2a. We used it with Eq. (3) to generate the \(\Delta \epsilon \prime\prime\)(ω) change, Fig. 2b, and used the Kramers-Kronig relation to generate the corresponding Δ\(\epsilon \prime\)(ω). Finally, we used Fresnel calculations to find ΔR(ω)/R(ω) = [R(ω; \(\epsilon\) = \(\epsilon _e\) + Δ\(\epsilon\)) − R(ω; \(\epsilon _e\))]/R(ω; \(\epsilon _e\)), Fig. 2c. These calculations serve to provide ΔR(ω)/R(ω) data, Fig. 2c, for which the electron energy distribution change, Fig. 2a is known. We now apply our first inversion procedure to this reflectivity data, which involves numerically inverting Eq. (2) to obtain \(\Delta \epsilon \prime\prime\), Fig. 2d, and then inserting this result into Eq. (3) and numerically inverting it to obtain, Δf(E), Fig. 2e. It is clear that to a good approximation the correct energy distribution change, Fig. 2a, is in fact obtained. (The blue points in Fig. 2e are the initial distribution points of Fig. 2a). We now consider cases with experimental data for ΔR/R where the underlying Δf(E) is unknown.

Fig. 2 Validation of the double inversion method. a A hypothetical change in the electron energy distribution, Δf(E). b The corresponding \(\Delta \epsilon \prime\prime\)(ħω) due to Δf(E). c The corresponding differential reflectivity, ΔR/R(ħω) inferred with the perturbed metallic dielectric function. In our procedure we take ΔR/R(ħω) and invert the relevant integral equations to obtain \(\Delta \epsilon \prime\prime\)(ħω), d, and Δf(E), e. The latter agrees well with the original form, a. Note that in the upper and lower panels the x axis is E(eV), where zero refers to the Fermi level, while in the three middle panels the x axis is ħω, where zero refers to the gap between the d-band and the conduction band of the L-point, which is roughly 2.45 eV Full size image

Inversion of normal incidence data

Figure 3a displays experimental transient (differential) reflectivity measurements of ΔR/R for a 30 nm gold film in the normal incidence configuration where the pump and probe are incident from the air side and a glass layer is under the gold film. The increasing (red) and decreasing (blue) changes in reflectivity are clearly seen around a probe photon energy ħω ≈ 2.4 eV, on a scale of a few picoseconds. Figure 3b presents the map of our inferred \(\Delta \epsilon \prime\prime\)(ω) from (2), where the structure of positive and negative peaks is opposite to the structure of the peaks in the reflectivity. Figure 3c depicts the map of the population changes Δf(E) inferred from (4), where now the positive and negative peak structures follow the same order as those in the differential reflectivity.

Fig. 3 Inversion of experimental pump-probe data. Maps show the changes in reflectivity, dielectric function and electron occupancy changes for the normal 30 nm configuration: a ΔR/R(ħω); b \(\Delta \epsilon \prime\prime\)(ħω) solved using the integral Eq. (2); c Δf(E) solved using the second integral Eq. (4) Full size image

Figure 4 presents several spectrum cuts (vertical cuts) of the maps of Fig. 3 at t = 500, 1000, 2000, 5000 fs for the normal 30 nm configuration, where Fig. 4a shows plots of the ΔR/R data, Fig. 4b shows the inferred \(\Delta \epsilon \prime\prime\)(ħω) solved using the integral Eq. (2), and Fig. 4c presents the inferred Δf(E) solved using the second integral Eq. (4). It is seen that there is only a minor difference between the curves corresponding to t = 500 fs and t = 1000 fs, while the changes start to decay after that time. The observed distribution at these relatively late times corresponds to a Fermi–Dirac distribution at an elevated temperature relative to the asymptotic (room) temperature. We find that the maximum electronic temperature in this case is ~1000 K which is within a factor of two of simple estimates based on an extended two-temperature model10,16 (see Supplementary Note 8). The origins of this discrepancy should be analyzed further and more sophisticated models may need to be investigated, including more complete Boltzmann equation approaches28.

Fig. 4 Spectrum cuts of the data and the inverted quantities. Spectrum cuts at times 500, 1000, 2000, 5000 fs for the normal 30 nm configuration: a ΔR/R(ħω); b \(\Delta \epsilon \prime\prime\)(ħω) solved using the integral Eq. (2); c Δf(E) solved using the second integral Eq. (4) Full size image

Figure 5a, b shows electron energy distribution changes for earlier times. The results presented, from t = 70 to 500 fs show an early broad electron energy distribution change at 70 fs that is somewhat different from the double-step-like forms often assumed and then the subsequent growth of a much larger magnitude Fermi–Dirac like function at later times. Unfortunately, owing to weak/noisy signals, we could not obtain reliable results for times much <70 fs. At these early times the electron distribution has several interesting characteristics. First it is clearly not a Fermi-Dirac distribution indicating the presence of non-thermalized electron distribution in the system. The fact that it deviates from a double-step-function shape and does not extend all the way to ±1 eV (pump photon energy) relative to the Fermi energy can be attributed to the fact that thermalization process is very fast and therefore the reflectivity measurement at these early times shows some time-averaging. In other words highly excited hot carriers already start decaying within the rise time of our pump pulse. Second, there is a slight asymmetry between the electron and hole distribution which is related to the difference in electron and hole density of states, as will be explained below. From Fig. 5 one sees the hole population (negative energies) extends down to −0.65 eV whereas the excess electron population (positive energies) extends up to 0.75 eV.

Fig. 5 Changes in carrier population for the normal incidence. a Spectrum cuts of Δf(E), inferred from the second integral Eq. (4), at t = 70, 150, 250, 500 fs. The relatively broad initial electron energy distribution change at the earlier time (70 fs) evolves to a narrower and larger magnitude Fermi–Dirac distribution at later time (500 fs). b Same as a but with multiplication factors in order to better compare the various time functional forms Full size image

There are several physical origins of the permittivity and electron occupancy changes discussed here. The initial 1.0 eV photons on the ultrafast time scale appear to create an equilibrium, Fermi-smeared distribution corresponding to the holes below the Fermi energy and excess electrons above it on a 500–1000 fs time scale followed by a subsequent relaxation to phonons, on a ps time scale similar to many studies, for example, ref. 6. There may be some band shifting, possibly due to thermal expansion, that could be responsible for the slight shifting of the apparent Fermi energy or chemical potential, that is, the fact the Δf(E) in our figures is not zero at exactly E = 0.

Inversion of Kretschmann configuration data

Figure 6 displays results of applying the double inversion scheme to reflectivity data measured in the Kretschmann configuration. The Kretschmann energy distribution changes are somewhat wider than those of the normal incidence configuration in Fig. 5, although still it should be noted that they are narrower than the naive expectation, as in Fig. 2a. Figure 6 shows that the early time Kretschmann case hole population extends down to −0.75 eV and the excess electron population extends further up to 1 eV, which can be contrasted with the normal incidence values of −0.65 and 0.75 eV, indicating the broader extent of the early time electron disturbance and greater asymmetry present in the Kretschmann case.

Fig. 6 Changes in carrier population for the Kretschmann configuration. a Spectrum cuts of Δf(E), inferred from the second integral Eq. (4), at t = 136, 163, 400 fs. The relatively broad initial electron energy distribution change at the earlier time (136 fs) evolves to a narrower and larger magnitude Fermi–Dirac distribution at later time (400 fs). b Same as a but with multiplication factors in order to better compare the various time functional forms Full size image

We should note that the slight negative dips in some of the Δf cases of Figs. 5 and 6 for positive energies are likely artifacts due to imperfections in the experimental data and our double inversion procedure. In particular, the effects of data truncation and extrapolation outside the available frequency range, important considerations for the Kramers-Kronig relations, are discussed in Supplementary Note 7.

The somewhat wider distributions in the Kretschmann case (Fig. 6) are likely due to Landau damping of the plasmon29 which cannot occur in the case of normal incidence because no plasmons are excited. The overall narrower and somewhat non-symmetric shape of the early time non-Fermi energy distribution change functions relative to naive expectations is consistent with quantum mechanical calculations for the intraband hot electron generation3,30 and represents a key result of our paper. Narrower distributions are expected because of two reasons. First, high-energy electrons are created in the processes with non-conservation of linear momentum and, therefore, the probability to create electrons with high excitation energies is smaller than that for low-energy electrons; thus, the distribution tends to be closer to the Fermi level. Second, e–e relaxation times for high energy electrons are shorter and, therefore, the number of such energetic electrons should be smaller. The asymmetry of excited electrons and holes is also well expected3,30. Since the density of states depends on energy, the calculated electrons tend to have a wider energy distribution (that is, more distributed between E F and E F + ħω 0 ) as compared to the calculated holes30; such character of the distributions was found in quantum calculations. It should be also be noted that Landau damping and the resulting direct intraband transitions (diagonal arrow in Fig. 1) typically need much smaller real space confinement of fields and respective wavevectors to become resonant. However, the propagating SPPs that decay evanescently normal to the gold surface will inevitably have some wavevector components that will fall in this range. This electron decay channel can be thought of as one of the most fundamental loss mechanisms for propagating plasmons in metals. To the best of our knowledge this is the first direct observation of ultrafast Landau damping in metal plasmonic systems.

It is important to note that some of the features of the early time energy distribution changes we have inferred from the experimental results do also emerge from certain other theoretical analyses and simulations. A master equation approach31 to non-thermal electron relaxation after optical excitation in small gold and silver nanoparticles also showed significantly greater and broader responses on plasmon resonance as we have found. The Boltzmann equation simulations of ref. 28 show an evolution of the electron distribution for a silver film that is consistent with a somewhat narrower distribution resulting after a short (20 fs) pulse. Interestingly, an extended two-temperature model9 of non-thermal electron dynamics in gold films, which correlated well with experiment, also showed early time electron energy distributions with structures closer to the Fermi energy.