††thanks: These two authors have contributed equally.††thanks: These two
authors have contributed equally.

# The electronic and thermal response of low electron density Drude materials
to ultrafast optical illumination

Subhajit Sarkar<EMAIL_ADDRESS>Ieng Wai Un<EMAIL_ADDRESS>Yonatan Sivan<EMAIL_ADDRESS>School of Electrical and Computer
Engineering, Ben-Gurion University of the Negev and the Ilse Katz Center for
Nanoscale Science and Technology, Ben-Gurion University of the Negev, Beer
Sheva, Israel

###### Abstract

Many low electron density Drude (LEDD) materials such as transparent
conductive oxide or nitrides have recently attracted interest as alternative
plasmonic materials and future nonlinear optical materials. However, the
rapidly growing number of experimental studies has so far not been supported
by a systematic theory of the electronic, thermal and optical response of
these materials. Here, we use the techniques previously derived in the context
of noble metals to go beyond a simple electromagnetic modelling of low
electron density Drude materials and provide an electron dynamics model for
their electronic and thermal response. We find that the low electron density
makes momentum conservation in $e-ph$ interactions more important, more
complex and more sensitive to the temperatures compared with noble metals;
moreover, we find that $e-e$ interactions are becoming more effective due to
the weaker screening. Most importantly, we show that the low electron density
makes the electron heat capacity much smaller than in noble metals, such that
the electrons in LEDD materials tend to heat up much more and cool down faster
compared to noble metals. While here we focus on indium tin oxide (ITO), our
analytic results can be easily applied to any LEDD materials.

##  I Introduction

Recent years have seen a growing interest in the optical response of plasmonic
(transparent) conducting oxides and nitrides such as aluminum/gallium-doped
zinc oxide, cadmium oxide, indium tin oxide, titanium/zirconium nitride etc.
[1, 2, 3, 4, 5, 6] as alternatives to noble metals. These materials are
characterized by low electron densities and high frequency interband
transitions, such that they are usually described as Drude metals at optical
frequencies [1]. In comparison to noble metals, the real part of their
permittivities has a milder negative value, the electron density (hence,
permittivity) is highly tunable (e.g., via doping [1, 7] or post-treatment
etc. [8, 9, 10]), so that together with their CMOS compatibility, they enable
flexibility of design and implementation [1]. In what follows, we refer to
this class of materials as low electron density Drude (LEDD) materials.

Particular attention has been given to the unique epsilon-near-zero (ENZ)
point these materials have in the near infrared spectral range [2, 11, 5, 12].
While most earlier interest in ENZ materials was associated with their linear
response (e.g., in the context of supercoupling [13, 14] or shaping the
radiation pattern of a source [15, 16]), the realization that the nonlinear
optical response is inversely proportional to the (unperturbed vanishing)
permittivity attracted a range of experimental studies of the (ultrafast)
dynamics of the permittivity near the ENZ point (e.g., [2, 17, 18, 11, 4, 12,
19]). In parallel, the importance of the deviation of the band structure from
a simple parabolic dependence was realized [8, 4, 3]. From the theoretical
point of view, most attention was given to modelling the dependence of the
LEDD permittivity on the electron temperature using thermal models [20, 17, 4,
20, 21] while the temperature dynamics itself was described using the
Relaxation Time Approximation (RTA)-based Two Temperature Model (TTM) or its
extension [20, 17, 4, 11, 20, 18, 3]. However, sometimes the TTM parameters
were computed from equations suitable for parabolic bands and high density
metals. In that regard, the accuracy of the temperature modelling in LEDD
materials has not yet been verified using an underlying electronic model.

Our goal in the current work is to compute consistently the non-equilibrium
electron dynamics in LEDD materials. We employ the techniques previously used
in the context of noble metals to go beyond the state-of-the-art modelling of
LEDD materials and provide a model for their electronic and thermal dynamical
response. We focus on indium tin oxide (ITO) as a prototypical example. In
particular, in Section II, we derive the various terms in the simplest model
available for electron non-equilibrium, namely, the Boltzmann equation (BE)
without resorting to the Relaxation Time Approximation (RTA); we complement
the BE with a self-consistent treatment of the phonon system. We find that due
to the low electron density, hence, weaker screening, the $e-e$ collisions are
10 times faster than in the higher electron density noble metals. We also find
that for the same reason, the requirement of momentum conservation in $e-ph$
interactions significantly slows down the collision rate; however, due to the
opposing effect of the higher Debye energy, the $e-ph$ collision rate is
actually similar to that in noble metals. In Section III.1 we describe the
resulting ultrafast dynamics of the electron distribution, and correlate it
with the relative importance of the various underlying collision mechanisms.
In Section III.2 we derive analytic expressions for the TTM parameters and
find that the dependence of both heat capacity and chemical potential on the
electron temperature is much stronger than assumed so far, and their values
are much lower compared to noble metals. We then extract the temperature
dynamics from the BE and show an excellent match with the TTM. Our main
results are that the electron heating is much stronger and its decay is much
faster in ITO compared to noble metals due to the lower electron heat
capacity. In Section IV we conclude with further comparison to more advanced
theoretical approaches, and with a discussion of future steps.

##  II Model for the electron distribution

We determine the electron distribution in LEDD materials by solving the
quantum-like Boltzmann equation (BE). This model is in wide use for Drude
metals [22, 23, 24, 25, 26, 27, 28, 29, 30]. We focus on Indium Tin Oxide
(ITO) because it has all the unique features of a LEDD material, namely, the
low electron density, a non-parabolic conduction band and because it emerges
as a promising CMOS compatible nonlinear optical material. The energy-momentum
($\mathcal{E}-k$) relation of ITO can be expressed by the Kane quasi-linear
dispersion [31] to account for the nonparabolicity [20, 32], namely,

$\displaystyle\dfrac{\hbar^{2}k^{2}}{2m_{e}^{\ast}}=\mathcal{E}(1+C\mathcal{E}),$
(1)

where $m_{e}^{\ast}=0.3964m_{e}$ is the electron effective mass at the
conduction band minimum, and $C=0.4191$eV-1 [8] is the first-order
nonparabolicity factor. The electron density of states (eDOS) becomes

$\displaystyle\rho_{e}(\mathcal{E})=\dfrac{1+2C\mathcal{E}}{2\pi^{2}}\left(\dfrac{2m_{e}^{\ast}}{\hbar^{2}}\right)^{3/2}\sqrt{\mathcal{E}(1+C\mathcal{E})}.$
(2)

Compared with the case in which $C=0$ (no non-parabolicity), the density of
states when $C\neq 0$ is a superlinear function (instead of a radical
function) of the electron energy and is larger by a factor of
$\sqrt{1+C\mathcal{E}}(1+2C\mathcal{E})$, see Fig. 1(a). The $\mathcal{E}-k$
relation (1) allows us to represent the electron states in terms of energy
$\mathcal{E}$ rather than momentum. We also neglect interband ($d$ to $sp$)
transitions which occur in ITO only for photons having energies larger than
$3-3.5$eV [33]; these transitions, however, set the value of the background
permittivity, $\varepsilon_{b}$.

The Boltzmann equation representing electron dynamics is

$\displaystyle\frac{\partial f}{\partial t}$ $\displaystyle=$
$\displaystyle\left(\frac{\partial f}{\partial
t}\right)_{exc}+\left(\frac{\partial f}{\partial t}\right)_{e-ph\
collision}+\left(\frac{\partial f}{\partial t}\right)_{e-e\ collision},$ (3)

where $f(\mathcal{E})$ is the electron distribution function at an energy
$\mathcal{E}$, representing the population probability of electrons in a
system characterized by a continuum of states within the conduction band; this
description was shown to be suitable for particles of Drude metals of sizes as
small as a few nm [34].

The first term on the right-hand-side (RHS) of Eq. (3) describes excitation of
conduction electrons due to photon absorption, see Section II.1 below for its
explicit form. The second term on the RHS of Eq. (3) describes the population
relaxation due to collisions between electrons and phonons, see Section II.2
below for its explicit form. The third term on the RHS of Eq. (3) (see Section
II.3 below for its explicit form) represents implicitly the thermalization
induced by electron-electron ($e-e$) collisions, i.e., the convergence of the
non-thermal population into the thermalized Fermi-Dirac distribution, given by

$f^{T}\left(\mathcal{E};T_{e}\right)=\left(1+e^{(\mathcal{E}-\mu)/k_{B}T_{e}}\right)^{-1},$
(4)

where $\mu$ is the chemical potential, $k_{B}$ is the Boltzmann constant, and
$T_{e}$ is the electron temperature. Note that collisions of conduction
electrons with impurities may also be included, but while those contribute to
the permittivity, they have no effect on the electron distribution itself (on
the level of the BE (3)) [23].

Our model does not account for electron acceleration due to the force exerted
on them by the electric field (which involves a classical description, see
discussion in [35, 36]) nor for drift due to its gradients or due to
temperature gradients; these effects may be important only for nanostructures
with complex geometries. Similar simplifications were adopted in most previous
studies of electron non-equilibrium in Drude metals, e.g., [30, 37, 36].

Below, we set the electron density in ITO to be $n_{e}\sim 1.5\times 10^{27}$
m-3, such that it is characterized by a relatively low Fermi energy of
$\mathcal{E}_{F}=0.88$ eV (compared to noble metals) and a total size of
$\mathcal{E}_{\text{max}}=4$ eV [38, 39]. However, note that due to the non-
stochiometric nature of ITO (i.e., the dependence on its preparation
conditions [8, 40]), the electron density in ITO can vary from $10^{27}$ m-3
to $10^{28}$ m-3, see [39] and references therein. Similar and even lower
values are typical of other LEDD materials [1].

Similarly, the values of other parameters indicated below such as the
deformation potential, the sound velocity, Debye temperature etc. or even the
impurity density or grain size may also vary significantly between sample to
sample. Yet, the analysis presented below remains generic to all ITO and other
LEDD materials.

Figure 1: (Color online) (a) The normalized electron density of states
${\rho_{e}(\mathcal{E})}/{\rho_{e}(\mathcal{E}_{F})}$ as a function of
$\mathcal{E}/\mathcal{E}_{F}$. (b) The collision rate associated with the
$e-ph$ interaction ${1}/{\tau_{e-ph}}$ (12) with (blue solid line) and without
(orange dash-dotted line) accounting for momentum conservation and non-
parabolicity (i.e., $C\neq 0$). The blue dotted line represents the case
accounting for momentum conservation but not for non-parabolicity (i.e.,
$C=0$). The dark yellow dashed line represents the approximation of the $e-ph$
collision rate (13). (c) The thermalization rate associated with the $e-e$
interaction with (blue solid line, $C\neq 0$) and without non-parabolicity
(blue dotted line, $C=0$). The collision rates in (b) and (c) are calculated
at $T_{e}=500$ K and $T_{ph}=300$ K.

###  II.1 The quantum mechanical excitation term

The absorption of the incident light at frequency $\omega$ leads to the
excitation of electrons with initial energy
$\mathcal{E}_{\text{initial}}=\mathcal{E}$ to states with final energy
$\mathcal{E}_{\text{final}}=\mathcal{E}+\hbar\omega$. Here, we employ the
expression proposed in [29, 36] to link this process to the local electric
field. Specifically, the excitation term can be written as [29, 26, 36]

$\left(\dfrac{\partial f(\mathcal{E})}{\partial
t}\right)_{\text{exc}}=A(t)\Big{[}D_{J}(\mathcal{E}-\hbar\omega,\mathcal{E})\rho_{e}(\mathcal{E}-\hbar\omega)f(\mathcal{E}-\hbar\omega)(1-f(\mathcal{E}))\\\
-D_{J}(\mathcal{E},\mathcal{E}+\hbar\omega)\rho_{e}(\mathcal{E}+\hbar\omega)f(\mathcal{E})(1-f(\mathcal{E}+\hbar\omega))\Big{]}.$
(5)

Here, $D_{J}(\mathcal{E}_{\text{initial}},\mathcal{E}_{\text{final}})$ is the
squared magnitude of the transition matrix element for the electronic process
$\mathcal{E}_{\text{initial}}\rightarrow\mathcal{E}_{\text{final}}$ (be it
Landau damping, surface/phonon-assisted absorption, etc. [41, 42, 30]); $A(t)$
is a time-dependent coefficient ensuring energy conservation that is
proportional to the local energy density, or $|{\bf E}(t)|^{2}$ (${\bf E}(t)$
being the local electric field vector in the ITO sample)111For simplicity, we
neglect any inhomogeneity of the local electric field in the ITO sample. ; in
particular, the increase of the energy of the electron system is equal to the
energy of the absorbed pulse (via Poynting’s Theorem)

$P_{\text{abs}}\equiv\int\mathcal{E}\rho_{e}(\mathcal{E})\left(\dfrac{\partial
f(\mathcal{E})}{\partial
t}\right)_{\text{exc}}d\mathcal{E}=\dfrac{\omega}{2}\varepsilon_{0}\varepsilon^{\prime\prime}|{\bf
E}(t)|^{2},$ (6)

so that

$A(t)=\dfrac{\varepsilon_{0}\varepsilon^{\prime\prime}|{\bf
E}(t)|^{2}}{2\hbar\int
D_{J}(\mathcal{E}-\hbar\omega,\mathcal{E})\rho_{e}(\mathcal{E})\rho_{e}(\mathcal{E}-\hbar\omega)f(\mathcal{E}-\hbar\omega,t)(1-f(\mathcal{E},t))d\mathcal{E}},$
(7)

where $\varepsilon^{\prime\prime}$ is the imaginary part of the ITO
permittivity.

###  II.2 The $e-ph$ collision term

The collision term due to electron-phonon interaction via the deformation
potential scattering is given by [44],

$\displaystyle\left(\dfrac{df_{{\bf k}}}{dt}\right)_{e-ph}=\sum_{{\bf q}}$
$\displaystyle\dfrac{\pi D^{2}q^{2}}{\rho V\omega_{ph}^{{\bf
q}}}\Bigg{\\{}\Big{[}(n^{{\bf q}}_{ph}+1)(1-f_{{\bf k}})f_{{\bf k}+{\bf
q}}-n^{{\bf q}}_{ph}(1-f_{{\bf k}+{\bf q}})f_{{\bf
k}}\Big{]}\delta(\mathcal{E}_{{\bf k}+{\bf q}}-\mathcal{E}_{{\bf
k}}-\hbar\omega_{ph}^{{\bf q}})$ $\displaystyle-$ $\displaystyle\
\Big{[}(n^{{\bf q}}_{ph}+1)(1-f_{{\bf k}-{\bf q}})f_{{\bf k}}-n^{{\bf
q}}_{ph}(1-f_{{\bf k}})f_{{\bf k}-{\bf q}}\Big{]}\delta(\mathcal{E}_{{\bf
k}-{\bf q}}-\mathcal{E}_{{\bf k}}+\hbar\omega_{ph}^{{\bf q}})\Bigg{\\}},$ (8)

where ${\bf k}$ is the electron momentum, ${\bf q}$ is the phonon momentum,
$\hbar\omega_{ph}^{{\bf q}}$ is the phonon energy, $\rho=7120$ kg/m3 is the
material density [45], $D=17.2$ eV is the deformation potential 222This value
of the deformation potential was obtained by fitting the ITO permittivity (at
room temperature) calculated by the Lindhard formula [75] with that measured
experimentally [21]. and $n^{{\bf q}}_{ph}$ is the phonon distribution
function.

For simplicity, we consider only acoustic phonons since they have been found
to be responsible for the dominant scattering mechanism compared to optical
phonons [47]. We also assume that the acoustic phonons propagate at the sound
velocity such that they have a linear dispersion relation, i.e.,
$\mathcal{E}_{ph}^{{\bf q}}=\hbar\omega_{ph}^{{\bf q}}=\hbar v_{ph}q$ where
$v_{ph}=6400$ m/s [45]. Beyond the Debye energy,
$\mathcal{E}_{D}=k_{B}T_{D}\approx 0.06$eV [45], the density of phonon states
vanishes. Furthermore, we have assumed that the phonon system is in
equilibrium, so that the average phonon number satisfies the Bose–Einstein
statistics and can be characterized by the phonon (lattice) temperature
$T_{ph}$, i.e., $n^{{\bf q}}_{ph}(T_{ph})=n_{B}(\mathcal{E}_{ph}^{{\bf
q}},T_{ph})=\left(e^{\hbar\omega^{{\bf q}}_{ph}/k_{B}T_{ph}}-1\right)^{-1}$.
The two terms associated with $(n_{ph}^{{\bf q}}+1)$ correspond to the phonon
emission, whereas the two terms associated with $n_{ph}^{{\bf q}}$ correspond
to the phonon absorption. The delta-functions in Eq. (II.2) ensures energy
conservation in the various electron-phonon scattering processes.

To obtain the collision term in terms of the electron energy
$\left(\dfrac{\partial f(\mathcal{E})}{\partial t}\right)_{e-ph}$, we perform
the spherical average over the solid angle for $\left(\dfrac{df_{{\bf
k}}}{dt}\right)_{e-ph}$ in Eq. (II.2), namely, $\left(\dfrac{\partial
f(\mathcal{E})}{\partial
t}\right)_{e-ph}=\dfrac{1}{4\pi}\displaystyle\int\left(\dfrac{\partial f_{{\bf
k}}}{\partial t}\right)_{e-ph}d\cos\theta_{{\bf k}}d\phi_{{\bf k}}$, where
$\theta_{{\bf k}}$ and $\phi_{{\bf k}}$ are, respectively, the polar and
azimuthal angles of ${\bf k}$ with respect to ${\bf q}$. In particular, we
account for the phonon dispersion in the polar integral to ensure momentum
conservation,

$\displaystyle\int_{-1}^{1}d\cos\theta_{{\bf k}}\delta(\mathcal{E}_{{\bf
k}\pm{\bf q}}-\mathcal{E}_{{\bf k}}\mp\hbar\omega_{ph}^{{\bf
q}})=\begin{cases}\dfrac{m_{e}^{\ast}}{\hbar^{2}kq},&\text{if
}\left|\dfrac{m_{e}^{\ast}v_{ph}}{\hbar k}\mp\dfrac{q}{2k}\right|\leq 1,\\\
0,&\text{otherwise}\end{cases}.$ (9)

Since $m_{e}^{\ast}v_{ph}\ll\hbar q_{D}=\hbar\mathcal{E}_{D}/v_{ph}$, the
condition for the polar integral being non-zero in Eq. (9) can be simplified
to $q\lesssim 2k$. After some algebra, we arrive at

$\left(\dfrac{\partial f(\mathcal{E})}{\partial t}\right)_{e-ph}=\
\dfrac{D^{2}}{4\pi\rho(\hbar
v_{ph})^{4}}\sqrt{\dfrac{m_{e}^{\ast}}{2\mathcal{E}(1+C\mathcal{E})}}\int_{0}^{\mathcal{E}_{ph}^{\text{max}}(\mathcal{E})}(\mathcal{E}_{ph})^{2}d\mathcal{E}_{ph}\\\
\Big{[}(n_{B}(\mathcal{E}_{ph},T_{ph})+1)\Big{(}(1-f(\mathcal{E}))f(\mathcal{E}+\mathcal{E}_{ph})(1+2C(\mathcal{E}+\mathcal{E}_{ph}))\\\
-(1-f(\mathcal{E}-\mathcal{E}_{ph}))f(\mathcal{E})(1+2C(\mathcal{E}-\mathcal{E}_{ph}))\Big{)}\\\
+n_{B}(\mathcal{E}_{ph},T_{ph})\Big{(}(1-f(\mathcal{E}))f(\mathcal{E}-\mathcal{E}_{ph})(1+2C(\mathcal{E}-\mathcal{E}_{ph}))\\\
-(1-f(\mathcal{E}+\mathcal{E}_{ph})f(\mathcal{E})(1+2C(\mathcal{E}+\mathcal{E}_{ph}))\Big{)}\Big{]},$
(10)

where

$\displaystyle\mathcal{E}_{ph}^{\text{max}}(\mathcal{E})=\text{min}(\mathcal{E}_{D},2\hbar
v_{ph}k(\mathcal{E}))=\begin{cases}2\hbar v_{ph}k(\mathcal{E})&\text{for
}k(\mathcal{E})<q_{D}/2\\\ \mathcal{E}_{D}\equiv\hbar v_{ph}q_{D}&\text{for
}k(\mathcal{E})\geq q_{D}/2\end{cases}$ (11)

is the maximal energy of a phonon which can be absorbed/emitted by an electron
with energy $\mathcal{E}$ [48, 49], $\hbar q_{D}=\mathcal{E}_{D}/v_{ph}$ is
the Debye momentum and $k(\mathcal{E})$ is the momentum of an electron with
energy $\mathcal{E}$ obtained from the dispersion relation (1).

By taking the functional derivative of Eq. (10) with respect to
$f(\mathcal{E})$ [50], we obtain the collision rate associated with the $e-ph$
interaction, $\dfrac{1}{\tau_{e-ph}(\mathcal{E})}=\dfrac{\delta}{\delta
f(\mathcal{E})}\left(\dfrac{\partial f}{\partial t}\right)_{e-ph}$,

$\dfrac{1}{\tau_{e-ph}(\mathcal{E})}=\dfrac{D^{2}}{4\pi\rho(\hbar
v_{ph})^{4}}\sqrt{\dfrac{m_{e}^{\ast}}{2\mathcal{E}(1+C\mathcal{E})}}\int_{0}^{\mathcal{E}_{ph}^{\text{max}}(\mathcal{E})}(\mathcal{E}_{ph})^{2}d\mathcal{E}_{ph}\\\
\Big{[}(n_{B}(\mathcal{E}_{ph},T_{ph})+f(\mathcal{E}+\mathcal{E}_{ph}))(1+2C(\mathcal{E}+\mathcal{E}_{ph}))+\\\
(n_{B}(\mathcal{E}_{ph},T_{ph})+1-f(\mathcal{E}-\mathcal{E}_{ph}))(1+2C(\mathcal{E}-\mathcal{E}_{ph}))\Big{]}.$
(12)

The $e-ph$ collision rate (12) is plotted in Fig. 1(b). It grows monotonically
up to $k(\mathcal{E})=q_{D}/2$ ($\mathcal{E}\sim 2.8\ \mathcal{E}_{F}$),
beyond which point its energy-dependence becomes much weaker.

To gain more insight into the dependence of $\tau_{e-ph}^{-1}$ on material
parameters, we simplify Eq. (12) by expanding its integrand in a power series
in $\mathcal{E}_{ph}$ since the phonon energy is small relative to the
electron energy. After integration of the leading-order term of the integrand
($k_{B}T_{ph}\mathcal{E}_{ph}$) one gets

$\displaystyle\dfrac{1}{\tau_{e-ph}(\mathcal{E})}\approx\dfrac{D^{2}k_{B}T_{ph}}{4\pi\rho(\hbar
v_{ph})^{2}}\dfrac{m_{e}^{\ast}}{\hbar
k(\mathcal{E})}(1+2C\mathcal{E})\times\begin{cases}4k^{2}(\mathcal{E}),&\text{for
}k(\mathcal{E})<q_{D}/2\\\ q_{D}^{2},&\text{for }k(\mathcal{E})\geq
q_{D}/2\end{cases},$ (13)

and shows a descent agreement with the exact expression (12), see Fig. 1(b).
The mismatch between Eq. (12) and (13), including the dip-like feature near
$\mathcal{E}_{F}$, can be resolved by incorporating higher-order terms. Eq.
(13) shows that $\tau_{e-ph}^{-1}$ is proportional to the phonon temperature
as in noble metals [51], and that non-parabolicity increases
$\tau_{e-ph}^{-1}$ by a factor of
$k(\mathcal{E}_{F})(1+2C\mathcal{E}_{F})\sim\sqrt{1+C\mathcal{E}_{F}}(1+2C\mathcal{E}_{F})\sim
2$, similar to the electron density of states, as shown in Fig. 1(b).

As a comparison, we also plot the $e-ph$ collision rate as calculated without
accounting for momentum conservation. This shows that the $e-ph$ collision
rate can be overestimated if momentum conservation is neglected. Indeed, due
to the low electron density, the Fermi momentum ($k_{F}\approx 3.54$ nm-1) is
much smaller than the Debye momentum of ITO ($q_{D}\approx 14.3$ nm-1) so that
a substantial amount of the phonons are prohibited from interacting with the
electrons due to conservation of momentum. This can be further understood
using a phase-space argument, as detailed in Appendix A.1. In particular,
ignoring momentum conservation causes a $\sim 5$-fold over-estimate of the
$e-ph$ collision rate around the Fermi energy (see Fig. 1(b); this is the most
relevant energy regime, see Fig. 2(c) below); consequently, this causes a
30-fold over-estimate of the energy transfer rate between the electrons and
phonons (see Section III.2). This behaviour contrasts the situation in noble
metals, for which the Debye momentum (e.g., $q_{D}\approx 6.85$ nm-1 and
$q_{D}\approx 8.14$ nm-1 for Au and Ag, respectively) is smaller than their
Fermi momentum ($k_{F}\approx 11.57$ nm-1 and $k_{F}\approx 12$ nm-1 for Au
and Ag, respectively). For that reason, in noble metals almost all electrons
can interact with all phonons such that neglecting momentum conservation is
justified. Therefore, somewhat peculiarly, the contradicting effects of the
higher Debye energy in ITO (which make $\tau_{e-ph}^{-1}$ 4 times larger, see
Eq. (13)) on one hand, and the limitations on the allowed $e-ph$ scattering
events imposed by the momentum conservation (which make it $\sim 5$ times
smaller) on the other hand, make the overall magnitude of the $e-ph$ collision
rate in ITO comparable to that in noble metals.

###  II.3 The electron-electron ($e-e$) collision term

The $e-e$ interaction is responsible for the thermalization of the conduction
electrons. The corresponding collision term is derived from the screened
Coulomb interaction $U$ whose potential in momentum space is written by

$\displaystyle U({\bf
q})=\dfrac{e^{2}}{4\pi\varepsilon_{0}\varepsilon_{b}}\int
d^{3}r\dfrac{e^{-q_{\text{TF}}r}}{r}e^{-i{\bf q}\cdot{\bf
r}}=\dfrac{e^{2}}{4\pi\varepsilon_{0}\varepsilon_{b}}\dfrac{1}{q^{2}+q_{\text{TF}}^{2}},$
(14)

where ${\bf q}$ is the exchange of momentum between two electrons, and
$q_{\text{TF}}$ is the Thomas-Fermi wave vector which is given by [26]

$\displaystyle
q^{2}_{\text{TF}}=-\dfrac{4\pi^{3}e^{2}}{\varepsilon_{0}\varepsilon_{b}}\int
d\mathcal{E}\dfrac{\partial f}{\partial\mathcal{E}}\rho_{e}(\mathcal{E}),$
(15)

where $\varepsilon_{b}=4$ represents the contribution of interband transitions
to the permittivity [21, 19]. Again, following the Fermi’s golden rule
employed in [44] and accounting for the non-parabolicity of the conduction
band, we obtain the population distribution change rate

$\left(\dfrac{\partial f(\mathcal{E})}{\partial
t}\right)_{e-e}=-\dfrac{3}{16\hbar^{2}}\dfrac{1}{(2\pi)^{3}}\left(\dfrac{e^{2}}{\varepsilon_{0}\varepsilon_{b}}\right)^{2}\left(\dfrac{m_{e}^{\ast}}{\hbar^{2}q_{\text{TF}}}\right)^{3}\dfrac{\hbar}{\sqrt{2m_{e}^{\ast}(\mathcal{E}+C\mathcal{E}^{2})}}\int
d\mathcal{E}_{1}d\mathcal{E}_{2}d\mathcal{E}_{3}\\\
(1+2C\mathcal{E}_{1})(1+2C\mathcal{E}_{2})(1+2C\mathcal{E}_{3})\left[\dfrac{qq_{\text{TF}}}{q^{2}+q^{2}_{\text{TF}}}+\tan^{-1}\left(\dfrac{q}{q_{\text{TF}}}\right)\right]^{q^{\text{max}}}_{q^{\text{min}}}\delta(\mathcal{E}+\mathcal{E}_{1}-\mathcal{E}_{2}-\mathcal{E}_{3})\\\
\left[(1-f(\mathcal{E}_{3}))(1-f(\mathcal{E}_{2}))f(\mathcal{E}_{1})f(\mathcal{E})-(1-f(\mathcal{E}))(1-f(\mathcal{E}_{1}))f(\mathcal{E}_{2})f(\mathcal{E}_{3})\right],$
(16)

where $q^{\text{min}}=\max(|k-k_{2}|,|k_{1}-k_{3}|)$ and
$q^{\text{max}}=\min(k+k_{2},k_{1}+k_{3})$. By taking the functional
derivative of Eq. (16) with respect to $f(\mathcal{E})$, we obtain the
collision rate associated with the $e-e$ interaction,

$\dfrac{1}{\tau_{e-e}(\mathcal{E})}=\dfrac{3}{16\hbar^{2}}\dfrac{1}{(2\pi)^{3}}\left(\dfrac{e^{2}}{\varepsilon_{0}\varepsilon_{b}}\right)^{2}\left(\dfrac{m_{e}^{\ast}}{\hbar^{2}q_{\text{TF}}}\right)^{3}\dfrac{\hbar}{\sqrt{2m_{e}^{\ast}(\mathcal{E}+C\mathcal{E}^{2})}}\int
d\mathcal{E}_{1}d\mathcal{E}_{2}d\mathcal{E}_{3}\\\
(1+2C\mathcal{E}_{1})(1+2C\mathcal{E}_{2})(1+2C\mathcal{E}_{3})\left[\dfrac{qq_{\text{TF}}}{q^{2}+q^{2}_{\text{TF}}}+\tan^{-1}\left(\dfrac{q}{q_{\text{TF}}}\right)\right]^{q^{\text{max}}}_{q^{\text{min}}}\delta(\mathcal{E}+\mathcal{E}_{1}-\mathcal{E}_{2}-\mathcal{E}_{3})\\\
\left[(1-f(\mathcal{E}_{3}))(1-f(\mathcal{E}_{2}))f(\mathcal{E}_{1})+(1-f(\mathcal{E}_{1}))f(\mathcal{E}_{2})f(\mathcal{E}_{3})\right].$
(17)

To gain more insight into the non-parabolicity correction to the $e-e$
collision rate, we substitute the electron distribution in the intergand of
Eq. (17) by the Fermi-Dirac distribution function characterized by an electron
temperature $T_{e}$; after some algebra, we find

$\displaystyle\dfrac{1}{\tau_{\text{e-e}}(\mathcal{E})}\approx
K\left\\{\left[(\pi
k_{B}T_{e})^{2}+(\mathcal{E}-\mathcal{E}_{F})^{2}\right]+\dfrac{4}{3}\dfrac{C(\mathcal{E}-\mathcal{E}_{F})}{(1+2C\mathcal{E}_{F})^{3}}(\pi
k_{B}T_{e})^{2}+\cdots\right\\},$ (18)

where the prefactor $K$ is given by

$\displaystyle
K=\dfrac{3}{32\hbar}\dfrac{1}{(2\pi)^{3}}\left(\dfrac{e^{2}}{\varepsilon_{0}\varepsilon_{b}}\right)^{2}\left(\dfrac{m_{e}^{\ast}}{\hbar^{2}q_{\text{TF}}}\right)^{3}\dfrac{1}{k_{F}}\left[\dfrac{2k_{F}q_{\text{TF}}}{4k_{F}^{2}+q^{2}_{\text{TF}}}+\tan^{-1}\left(\dfrac{2k_{F}}{q_{\text{TF}}}\right)\right](1+2C\mathcal{E}_{F})^{3}.$
(19)

This expression provides the generalization of Fermi liquid theory [50] to
Drude materials with non-parabolic bands. Compared with the case of $C=0$, the
non-parabolicity overall causes an increase of the density of states such that
the $e-e$ thermalization rate increases by a factor of
$(1+2C\mathcal{E}_{F})^{3/2}(1+C\mathcal{E}_{F})^{-5/4}$ 333From Eq. (15), one
can deduce that the Thomas-Fermi wavevector is proportional to the square root
of the eDOS at the Fermi-energy [23]. Therefore, compared with the case of
$C=0$, $q_{TF}$ is larger by a factor of
$(1+C\mathcal{E}_{F})^{1/4}(1+2C\mathcal{E}_{F})^{1/2}$ and $k_{F}$ is larger
by a factor of $(1+C\mathcal{E}_{F})^{1/2}$. Together with the factor
$(1+2C\mathcal{E}_{F})^{3}$ coming from the increase of eDOS, the $e-e$
thermalization rate increases by a factor of
$\sim(1+2C\mathcal{E}_{F})^{3/2}(1+C\mathcal{E}_{F})^{-5/4}$. ($\approx 1.5$
for ITO), and adds a small correction term which is linear with
$\mathcal{E}-\mathcal{E}_{F}$ 444Note that the term
$(\mathcal{E}-\mathcal{E}_{F})$ in Eq. (18) should not be replaced by
$\hbar\omega$, since it involves all possible electron states, rather than
only those excited from the Fermi energy by photon absorption., see Fig. 1(c).
Eq. (19) also shows the effect of the low electron density (hence, smaller
Fermi energy) in ITO - not only it results in a narrower energy-dependence of
$\tau^{-1}_{\text{e-e}}(\mathcal{E})$, it also results in smaller Fermi
momentum, and more importantly, in weaker screening, and thus, a smaller
Thomas-Fermi wavevector 555The Thomas-Fermi wavevector is proportional to the
square root of the eDOS at the Fermi-energy [23], thus, it decreases with
electron density.. This means that $e-e$ interactions are stronger in ITO,
leading to a $\sim 10$ times faster $e-e$ collision rate than Au.

##  III Results

In the example below, we solve the BE (3) in order to obtain the electronic
(i.e., solid-state physics) response of bulk ITO systems subject to (modestly
high level) pulsed electric field illumination
$I_{\text{inc}}=I_{0}e^{-4\ln(2)(t/\tau_{p})^{2}}$, a wavelength of $1300$ nm,
duration of $\tau_{p}=220$ fs and peak intensity of $I_{0}=2.5$ GW/cm2.
However, as for the uncertainty on the various material parameters, the
results below remain generic also for other parameter values.

###  III.1 electron dynamics

Fig. 2 plots the electron distributions and the different terms corresponding
to Boltzmann equation (BE) (3) at three different time slices,
$t=-\tau_{p}/2=-110$ fs, $t=0$ fs, and $t=-110$ fs corresponding to full width
at the half maximum (FWHM) before the peak, at the peak and FWHM after the
peak of the pulse.

Fig. 2(a) shows the initial electron distribution $f(\mathcal{E})$ as a
function of normalized energy $\mathcal{E}/\mathcal{E}_{F}$. The deviation
from thermal equilibrium is clearly visible via the characteristic shoulders
above Fermi energy, corresponding to one, two etc. consecutive photon
absorption events; indeed, those were observed previously for noble metals,
e.g., in [35, 28, 55]. They originate from the structure of the excitation
term (Fig. 2(b)). However, these shoulders gradually smooth out due to (the
weak) $e-ph$ collisions, and (the much stronger) $e-e$ collisions, see Fig.
2(c) and (d), respectively. Indeed, at later times the distribution rapidly
approaches a thermal distribution (see Fig. 2(e) and (i)).

The corresponding evolution of the various terms is seen in the additional
subplots of Fig. 2. It is interesting to note the differences with respect to
the corresponding dynamics in noble metals. Specifically, due to the low Fermi
energy, there is only a single $\hbar\omega$-wide region of negative rate of
change of population
$\displaystyle\left(\dfrac{df(\mathcal{E})}{dt}\right)_{exc}$ due to photon
absorption below the Fermi energy but a corresponding multiple-shoulder
structure of positive rate above the Fermi energy; the various shoulders are
energy-dependent due to the relatively strong energy dependence of the eDOS
(see Fig. 1(a)). Moreover, near the band minimum, photon absorption is weaker
due the vanishing electron density of states. This leads to a sudden cut-off
of $\partial f/\partial t|_{exc}$ (5) near the band minimum.

The structure of the electron-phonon ($e-ph$) term is significantly different
compared to its structure in noble metals, see Fig. 2(c), (g) and (k)
corresponding to Eq. (10). The origin of this difference is the importance of
momentum conservation along with the number conserving nature of the $e-ph$
interaction which limit the possible scattering processes, see details in
Appendix A.2. Also notable is the increase in magnitude of the $e-ph$ term in
time; this is related to the rise of the overall electron energy and the rapid
increase in the number of low excess energy electrons; this effect was already
shown in [25] to lead to an acceleration in the $e-ph$ rate, which is not
captured by the RTA.

Figure 2: (Color online) (a), (e) and (i) Electron distribution (black lines),
(b), (f) and (j) photo-excitation (blue), (c), (g) and (k) electron-phonon
interaction (green), and (d), (h) and (l) electron-electron interaction
(orange) terms following illumination of ITO by a short pulse. The left column
shows the various terms at the FWHM before the peak of the pulse
($t=-\tau_{p}/2=-110$ fs), the middle column at the centre of the pulse ($t=0$
fs), and the right column at the FWHM after the peak of the pulse ($t=+110$
fs) as a function of $\mathcal{E}/\mathcal{E}_{F}$ and pulse peak intensity
$I_{0}=2.5$ GW/cm2. The inset in (b) is the zoom-in of the region between
$\mathcal{E}=1.5$ eV and $\mathcal{E}=3.5$ eV; it shows more clearly the
origin of the multiple step structure seen in (a). All terms have the same
order of magnitude because of the pulsed nature of the illumination. The
$y$-axis is in log-scale for (a), (e) and (i) and in symmetric log-scale for
all other subplots.

Finally, Fig. 2 (d), (h) and (l) show the rate of change of $f(\mathcal{E})$
due to electron-electron ($e-e$) interactions corresponding to Eq. (16). Its
shape is similar to that in noble metals (cf [36]), except near the band
minimum and for $\mathcal{E}\approx\mathcal{E}_{F}+\hbar\omega$. This is due
to the energy and number conserving nature of the $e-e$ interaction. Near the
band minimum, the positive rate of change of population is enhanced due to the
vanishing eDOS.

###  III.2 Coarse-grained dynamics - the Two Temperature Model (TTM)

Quite frequently, it is sufficient to consider the macroscopic dynamics of the
energies of the electron and phonon subsystems, respectively. This can be
achieved first by integrating over the product of the various terms in the BE
(3) with the electron energy and eDOS (see, e.g., [26, 36]), which provides
the resulting total energy of the electron subsystem. Then, a dynamic equation
for the phonon energy can be written down based on the total rate of energy
transfer from the electron subsystem. In order to make the resulting equations
more meaningful, it is customary to rewrite the resulting energies as the
product of the respective heat capacities ($C_{e}$ and $C_{ph}$) and electron
and phonon temperatures. While the latter is well-defined, it is well known
that the notion of electron temperature cannot be defined clearly in the
initial stages of the dynamics [25, 28, 56]. In this context, it is customary
to extract an effective electron temperature by calculating the electron
temperature $T_{e}$ for which the total energy of a Fermi Dirac (i.e.,
thermal) distribution $f^{T}$ (4) is the same as that of the true non-thermal
distribution $f(\mathcal{E})$ (3) (see, e.g., [28, 56, 30]), namely,

$\mathcal{U}_{e}\equiv\int\mathcal{E}\rho_{e}(\mathcal{E})f(\mathcal{E})d\mathcal{E}=\int\mathcal{E}\rho_{e}(\mathcal{E})f^{T}(\mathcal{E},T_{e})d\mathcal{E}.$
(20)

The resulting effective electron temperature emerges to be the well-defined
electron temperature once the distribution thermalizes [56]. Since $e-e$
interactions conserve the energy of the electron subsystem, the integrated BE
emerge to be

$\displaystyle C_{e}\dfrac{dT_{e}}{dt}$
$\displaystyle=P_{abs}-G_{e-ph}(T_{e}-T_{ph}),$ (21)

while the equation for the phonon temperature is

$\displaystyle C_{ph}\dfrac{dT_{ph}}{dt}$ $\displaystyle=\ \
G_{e-ph}(T_{e}-T_{ph}).$ (22)

Note that since we are interested in the ultrafast dynamics, we ignore heat
transfer to the environment (assumed to be at the environment temperature), as
this process occurs on a much longer time scale. Equations (21)-(22)
constitute the so-called “Two Temperature Model”; here, $C_{e}(T_{e})$ and
$C_{ph}$ are the electron and phonon heat capacities of ITO, and $G_{e-ph}$ is
the electron-phonon coupling. The advantage of such a coarse-grained model is
considerable, as it is significantly simpler to solve compared with the BE,
and serves as the basis for temperature-based permittivity models.

The TTM parameters are usually hard to measure directly. For example, the
phonon heat capacity was measured to be $C_{ph}=2.54\times 10^{6}$
J$\cdot$m${}^{-3}\cdot$K-1 [57] and similar values arise from a direct
calculation. More frequently, the TTM parameters are calculated theoretically
using thermal distributions at the effective electron temperature. For noble
metals, relatively simple expressions for the various emerging parameters were
obtained (e.g., [58, 56]); in particular, both electron heat capacity and
$e-ph$ coupling can be written as linear functions of the temperatures. These
parameters, however, are significantly more complicated in LEDD materials
compared to noble metals because of the non-parabolicity, and especially
because of the much smaller electron density (hence, smaller Fermi energy
$\mathcal{E}_{F}$). In principle, the TTM parameters can be evaluated through
integral expressions (see below). However, in what follows we also provide
approximate analytic expression for these parameters; these expressions are
suitable for any value of intrinsic parameters of ITO (or other LEDD
materials), which indeed vary due to fabrication/doping conditions (see e.g.,
[1, 7, 8]).

One effect of the smaller $\mathcal{E}_{F}$ is that the dependence of the
chemical potential on temperature is not negligible as in noble metals 666This
is reminiscent of what happens in semiconductors [23, 76, 77]. . To see this,
we employ the Sommerfeld expansion [60, 23] to express the total energy of the
electron subsystem $\mathcal{U}_{e}$ as a Taylor expansion in powers of
$k_{B}T_{e}$, assuming purely thermal electron distributions. In standard
textbooks, e.g. [23], the expansion is usually kept up to the second-order
only. However, since the Fermi energy of ITO is much lower than that of metals
and since the incoming illumination intensity is strong such that the electron
temperature might become non-negligible with respect to the Fermi temperature
(e.g., $\sim 10,000~{}\text{K}$ in ITO), one needs to keep the expansion at
least up to the fourth-order, namely,

$\mathcal{U}_{e}\approx\int_{0}^{\mu(T_{e})}\mathcal{E}\rho_{e}(\mathcal{E})d\mathcal{E}+\dfrac{\pi^{2}}{6}(k_{B}T_{e})^{2}\left.\dfrac{d\left(\mathcal{E}\rho_{e}(\mathcal{E})\right)}{d\mathcal{E}}\right|_{\mathcal{E}=\mu(T_{e})}+\dfrac{7\pi^{4}}{360}(k_{B}T_{e})^{4}\left.\dfrac{d^{3}\left(\mathcal{E}\rho_{e}(\mathcal{E})\right)}{d\mathcal{E}^{3}}\right|_{\mathcal{E}=\mu(T_{e})}.$
(23)

Then, the temperature-dependent chemical potential $\mu(T_{e})$ is determined
using number conservation,
$\displaystyle\int_{0}^{\infty}\rho_{e}(\mathcal{E})f^{T}(\mathcal{E},\mu(T_{e}),T_{e})d\mathcal{E}=\displaystyle\int_{0}^{\mathcal{E}_{F}}\rho_{e}(\mathcal{E})d\mathcal{E}$,
i.e.,

$\displaystyle\mu(T_{e})\approx\mathcal{E}_{F}-\dfrac{(\pi
k_{B}T_{e})^{2}}{6}\dfrac{1}{\rho_{e}(\mathcal{E}_{F})}\left.\dfrac{d\rho_{e}(\mathcal{E})}{d\mathcal{E}}\right|_{\mathcal{E}_{F}},$
(24)

where
$\displaystyle\int_{0}^{\mathcal{E}_{F}}\rho_{e}(\mathcal{E})d\mathcal{E}$ is
the number of electrons at zero-temperature (which is indeed the same as that
at 300K). This expression is plotted in Fig. 3(a) vs. the exact numerical
solution.

Similarly to $\mu$, the exact integral definition of the electron heat
capacity can be approximated as

$C_{e}(T_{e})\equiv\frac{d\mathcal{U}_{e}}{dT_{e}}=\frac{d}{dT_{e}}\left[\int\mathcal{E}\rho_{e}(\mathcal{E})f^{T}(\mathcal{E},T_{e})d\mathcal{E}\right]\\\
\approx\dfrac{\pi^{2}k_{B}^{2}}{3}T_{e}\left[\rho_{e}(\mu)+\mu\rho_{e}^{(1)}(\mu)\right]+\dfrac{7\pi^{4}k_{B}^{4}}{90}T_{e}^{3}\left[3\rho_{e}^{(2)}(\mu)+\mu\rho_{e}^{(3)}(\mu)\right]\\\
+\dfrac{d\mu}{dT_{e}}\left(\mu\rho_{e}(\mu)+\dfrac{(\pi
k_{B}T_{e})^{2}}{6}\left[2\rho_{e}^{(1)}(\mu)+\mu\rho_{e}^{(2)}(\mu)\right]+\dfrac{7(\pi
k_{B}T_{e})^{4}}{360}\left[2\rho_{e}^{(3)}(\mu)+\mu\rho_{e}^{(4)}(\mu)\right]\right),$
(25)

where
$\displaystyle\rho_{e}^{(n)}(\mu)=\left.\dfrac{d^{n}\rho_{e}(\mathcal{E})}{d\mathcal{E}^{n}}\right|_{\mu}$
and $\mu(T_{e})$ is given by Eq. (24). As shown in Fig. 3(b), up to $T_{e}\sim
1300$K, the electron heat capacity scales linearly with the electron
temperature (viz., $C_{e}\approx\gamma_{e}T_{e}$ with
$\gamma_{e}=\pi^{2}k_{B}^{2}\left[\rho_{e}(\mu)+\mu\rho_{e}^{(1)}(\mu)\right]/3$
777Due to non-parabolicity the dependence of $\gamma_{e}$ on electron density
becomes is not straightforward. In contrast, in its absence, $\gamma_{e}$
reduces to the familiar expression
$\gamma_{e}=\pi^{2}n_{e}k_{B}^{2}/2\mathcal{E}_{F}$.). Note that the value of
$\gamma_{e}$ for ITO ($12.7~{}\text{J}\cdot\text{m}^{-3}\cdot\text{K}^{-2}$)
is much smaller than for noble metals, e.g.,
$\gamma_{e}=67.6~{}\text{J}\cdot\text{m}^{-3}\cdot\text{K}^{-2}$ for gold (Au)
[62]. This smaller value of $\gamma_{e}$ in ITO is associated with the lower
electron density; indeed, compared to noble metals, the eDOS is evaluated at
the much lower chemical potential, giving rise to a smaller value for
$\gamma_{e}$. Furthermore, the cubic dependence of $C_{e}$ (emerging from the
4th-order term in the Sommerfeld expansion of the electron energy) provides
decent accuracy only up to $T_{e}\sim 3500$ K. In particular, in this regime
of electron temperatures, $C)e$ experiences a sublinear growth due to the
decrease of $\mu$ with temperature.

Lastly, the $e-ph$ coupling coefficient $G_{e-ph}$ can be obtained by
evaluating the energy transferred from the electron to the phonon subsystem,
i.e., $\left(\dfrac{\partial\mathcal{U}_{e}}{\partial
t}\right)_{e-ph}=\displaystyle\int\mathcal{E}\rho_{e}(\mathcal{E})\left(\dfrac{\partial
f}{\partial t}\right)_{e-ph}d\mathcal{E}$, where $\left(\dfrac{\partial
f}{\partial t}\right)_{e-ph}$ is given by Eq. (10). As for $C_{e}$, we
substitute the Fermi-Dirac distribution at the effective electron temperature
for $f$ and expand the integrand in a power series in $\mathcal{E}_{ph}$. In
similarity to noble metals, we find that the leading order term is
proportional to the difference between the electron and phonon temperatures,
i.e., $\left(\dfrac{\partial\mathcal{U}_{e}}{\partial t}\right)_{e-ph}\approx-
G_{e-ph}(T_{e}-T_{ph})$, where the $e-ph$ coupling coefficient is given by

$G_{e-ph}(T_{e})=\dfrac{D^{2}m_{e}^{\ast}k_{B}}{4(\pi\hbar)^{3}\rho(\hbar
v_{ph})^{4}}\int_{0}^{\infty}d\mathcal{E}\int_{0}^{\mathcal{E}_{D}}(\mathcal{E}_{ph})^{2}d\mathcal{E}_{ph}\\\
\left[(1+2C\mathcal{E})\sech^{2}\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)\tanh\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)-4Ck_{B}T_{e}\sech^{2}\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)\right]\\\
-\dfrac{D^{2}m_{e}^{\ast}k_{B}}{4(\pi\hbar)^{3}\rho(\hbar
v_{ph})^{4}}\Bigg{\\{}\int_{0}^{\infty}d\mathcal{E}\int_{0}^{\mathcal{E}^{\text{max}}_{ph}}(\mathcal{E}_{ph})^{2}d\mathcal{E}_{ph}\\\
\bigg{[}(1+2C\mathcal{E})\sech^{2}\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)\tanh\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)-4Ck_{B}T_{e}\sech^{2}\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)\bigg{]}\\\
-\int_{0}^{\mathcal{E}_{\text{min,D}}}d\mathcal{E}\left[\dfrac{8(m_{e}^{\ast}v_{ph}^{2})^{2}}{k_{B}T_{e}}\mathcal{E}^{2}(1+C\mathcal{E})(1+2C\mathcal{E})^{3}\sech^{2}\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)\right]\Bigg{\\}},$
(26)

where $\mathcal{E}_{ph}^{\text{max}}$ is given by Eq. (11).

Similar to the other TTM parameters, the $e-ph$ coupling coefficient
$G_{e-ph}$ attains a different and more complex form compared to noble metals.
Indeed, we plot the electron temperature dependence of $G_{e-ph}$ (26) in Fig.
3(c). First, we find that the magnitude of $G_{e-ph}$ of ITO ($\sim 3\times
10^{16}$ Jm-3K-1s-1) is similar to that of Au ($\sim 2.5\times 10^{16}$
Jm-3K-1s-1 [63, 26] 888For Au, Eq. (27) becomes $\displaystyle
G_{e-ph}^{(Au)}=\dfrac{D^{2}m_{e}k_{B}\mathcal{E}_{D}^{4}}{16(\pi\hbar)^{3}\rho(\hbar
v_{ph})^{4}}$; with the known parameters for Au [26], (viz., deformation
potential $D=19.3$ eV, Debye energy $\mathcal{E}_{D}=0.014$ eV, mass density
$\rho=19.3\times 10^{3}~{}\text{Kg}/\text{m}^{3}$ and phonon velocity
$v_{ph}=3240\text{m/s}$) one obtains the value above.) although the Debye
energy of ITO is 4 times larger than that of Au. The reason is that
conservation of momentum prohibits a large portion of the phonons from
interacting with the electrons, similar to the $e-ph$ collision rate in
Section II.2. Secondly, Fig. 3(c) shows that quite different from noble
metals, $G_{e-ph}$ of ITO increases with the electron temperature. This is
again a direct result of the momentum conservation constraint. In general, the
momentum conservation in $e-ph$ collision causes that electrons with momentum
$\hbar k$ can only absorb/emit phonons with momentum ranging from 0 to $2\hbar
k$ (see details in Appendix A.1). In ITO, the momentum of many of the
electrons is much smaller than the Debye momentum (due to $k_{F}=q_{D}/4$), so
the maximal momentum of a phonon which can be absorbed/emitted by an electron
increases with the momentum of the electron. Therefore, when the electron
temperature increases, more electrons occupy higher energy states; these
higher energy electrons can then interact with higher energy phonons, leading
to faster transfer of energy from the electrons to the phonon subsystem and,
thus, to an increase in $G_{e-ph}$ 999The same reasoning explains why an
increase in the phonon temperature does not have a significant effect, see Eq.
(26).. In contrast, in noble metals the momentum of electrons is much larger
than the Debye momentum so that the maximal energy of a phonon which can be
absorbed/emitted by electrons is limited by the Debye momentum instead of the
electron momentum, so that $G_{e-ph}$ is independent of the electron
temperature [26, 66, 27, 67, 36].

To gain more insight into $G_{e-ph}$, we analyze Eq. (26) in some simple
limits. First, we consider $e-ph$ coupling coefficient without accounting for
momentum conservation. This corresponds to the first term (the first two
lines) in Eq. (26), and it can be simplified to be (see Eq. (32))

$\displaystyle G_{e-ph}(T_{e})=\dfrac{D^{2}m_{e}^{\ast
2}k_{B}\mathcal{E}_{D}^{4}}{16(\pi\hbar)^{3}\rho(\hbar
v_{ph})^{4}}(1+2C\mu(T_{e}))^{2}\qquad\text{(without momentum conservation).}$
(27)

This approximation is widely used in modeling of noble metals [26, 66, 27, 67,
36], however, it is poor for ITO. Indeed, the ratio between Eq. (26) and Eq.
(27) shows that $G_{e-ph}$ is smaller than that of the case without accounting
for momentum conservation by a factor of $\left(\dfrac{\mathcal{E}_{D}}{2\hbar
v_{ph}k_{F}}\right)^{4}$ ($\sim 30$), see Fig. 3(c). Furthermore, due to the
non-parabolicity, the $e-ph$ coupling coefficient (27) is larger by an overall
factor of $(1+2C\mu)^{2}$ ($\sim 3$ for ITO) and shows an (incorrect) decrease
with the electron temperature (due to the decrease of chemical potential $\mu$
with the electron temperature $T_{e}$), see Fig. 3(c). For a more accurate
approximation, we consider the $e-ph$ coupling coefficient with momentum
conservation in the zero temperature limit. In this case, Eq. (26) becomes
(see Appendix A.3)

$\displaystyle G_{e-ph}(T_{e}\rightarrow 0\ \text{K})=\dfrac{D^{2}m_{e}^{\ast
2}k_{B}k_{F}^{4}}{\pi^{3}\rho\hbar^{3}}(1+2C\mathcal{E}_{F})^{2}\quad(k_{F}<q_{D}/2\text{
and momentum conserved)}.$ (28)

Figure 3: (Color online) (a) chemical potential $\mu$ and (b) electron heat
capacity $C_{e}$ as a function of electron temperature $T_{e}$. The blue solid
and dot-dashed lines correspond to the numerical values and analytical forms
(Eqs. (24) for $\mu$ and (25) for $C_{e}$), respectively. The dotted line in
(b) corresponds to the linear approximation of $C_{e}$ indicating the
deviation from the linear dependence above $T_{e}\sim 1300$ K. (c) The $e-ph$
coupling coefficient $G_{e-ph}$ as a function of the electron temperature (26)
(blue solid line) and its zero temperature limit (28) (blue dotted line). The
blue dashed line represents the case without accounting for momentum
conservation (27).

Having determined the various parameters appearing in the TTM equations, we
can solve them, plot the resulting (effective) electron as well as the phonon
temperature dynamics in Fig. 4, and compare it to the dynamics in noble
metals, e.g., gold (Au), for the same heat source (6). Overall, the dynamics
in these two systems is qualitatively similar, namely, the electron
temperature grows on a few 100 fs timescale (dictated by the pulse duration),
and then decays due to $e-ph$ interactions. However, the total electron
heating in ITO is much higher than in Au. The reason for that is the
difference in the values of the corresponding heat capacities; indeed, at low
temperatures the electron heat capacity is linear with the electron
temperature (i.e., $C_{e}\sim\gamma_{e}T_{e}$), and $\gamma_{e}$ for ITO is
about 5 times smaller than in Au (as already mentioned above,
$\gamma_{e}=67.6~{}\text{J}\cdot\text{m}^{-3}\cdot\text{K}^{-2}$ in Au [62]
and $\gamma_{e}=12.7~{}\text{J}\cdot\text{m}^{-3}\cdot\text{K}^{-2}$ obtained
for ITO). As a rough estimate, one can ignore $e-ph$ heat transfer in the
initial stages of the dynamics, so that the temperature rise can be easily
shown to scale as $\sqrt{\mathcal{U}_{\text{abs}}/\gamma_{e}}$ (see, e.g.,
[68]). The ratio of the maximal temperature rise in ITO and Au ($\sim 2100$ K
and $\sim 900$ K, respectively) is indeed given roughly by
$\sqrt{67.6/12.7}\sim 2.3$.

Another notable difference is that the decrease rate of the (effective)
electron temperature (and correspondingly, the rise time of the phonon
temperature) is faster in ITO compared to Au, see Fig. 4(a). To leading order,
these rates are determined by the ratio $G_{e-ph}/C_{e}$; since $C_{e}$ is
lower in ITO, but $G_{e-ph}$ is comparable in ITO and Au, the rates in ITO are
higher, see Fig. 4(b). Nevertheless, since the phonon heat capacity and the
heat absorption is similar in both ITO and Au, the eventual phonon temperature
reached is similar in the two systems (not shown).

Figure 4: (a) (Color online) Effective electron temperatures (orange lines)
and phonon temperatures (green lines) as a function of time for the pulse
parameters specified in the text. Solid and dotted lines correspond to
temperatures obtained from numerical calculation using the BE (using Eq. (20)
and (22)) and the solution of the two-temperature model (TTM) (Eqs.
(21)-(22)), respectively. Dot-dashed lines correspond to the temperature
dynamics of gold (Au) obtained from the TTM simulation with
$C_{e}=\gamma_{e}T_{e}$ (with
$\gamma_{e}=67.6~{}\text{J}\cdot\text{m}^{-3}\cdot\text{K}^{-2}$ [62]) and
electron-phonon coupling $G_{e-ph}=2.5\times
10^{16}~{}\text{J}\cdot\text{m}^{–3}\cdot\text{K}^{–1}\cdot\text{s}^{–1}$ [63,
26]. (b) The ratio $G_{e-ph}/C_{e}$ for ITO (solid line) and Au (dot-dashed
line).

##  IV Discussion and outlook

We have seen that the lower electron density along with the non-parabolicity
distinguishes the electron and heat dynamics in ITO (and more generally, LEDD
materials) from those in noble metals. In particular, we identified
significant differences in the $e-ph$ interactions, a faster $e-e$ collision
rate, a much stronger dependence of the TTM parameters on the electron
temperature and a different overall heating and dynamics due to a lower
electron heat capacity. The analytic expressions obtained for the TTM
parameters allow an easy investigation of other LEDD materials.

We have also shown that the TTM matches remarkably well the dynamics of the
effective electron temperature as well as of the phonon temperature.
Nevertheless, the TTM has known limitations in noble metals; in particular, it
assumes a-priori that the $e–e$ scattering is fast enough to establish a
thermal distribution of electrons before significant energy is transferred to
the phonons and cannot account for the accelerating rate of $e-ph$ collisions
[25]. While this is a problematic assumption for noble metals, it was shown in
[69] that this assumption holds well in light metals such as Na, Cs, Rb, and
K. In that respect, the condition of validity of the TTM in ITO does not
strictly hold, yet, the faster $e-e$ collision rate makes it is closer to be
satisfied in comparison to noble metals.

In order to improve upon the TTM, it is customary to add a dynamical equation
for the total non-thermal energy (e.g. [70, 71, 72, 36, 73] for noble metals
or [11] for ITO); this is usually done within the RTA, requiring a somewhat
ambiguous choice of an energy-averaged decay coefficients of the non-thermal
energy to the (thermal) electron and phonon subsystems. Whether such an
improvement is necessary or not requires a rigorous consideration of the
permittivity (or, e.g., reflectivity) dynamics; it might reveal differences
between the thermal and non-thermal dynamics. This complicated task is left
for a future paper. Nevertheless, even without such an analysis, we can
already say that the decrease rate of the electron temperature decay rate
($G_{e-ph}/C_{e}$) with the electron temperature implies that the internal
thermalization process between these subsystems becomes slower with increasing
illumination intensity. This explains the observation in [19] of a slower
thermalization of the reflectivity with increasing illumination intensities.

In this vein, the current work is a starting point for modelling the
permittivity dynamics of the ITO and other LEDD materials, as well as the
observed spectral broadening and fast switching, both for pulses which are a
few 100’s of femtoseconds long, as well as shorter ones. In this context, it
would be of great interest to unravel the physical mechanism underlying the
nonlinear optical response of ITO at increasingly high illumination
intensities.

###### Acknowledgements.

I. W. Un and S. Sarkar contributed equally to this work. The authors were
supported by Israel Science Foundation (ISF) grant (340/2020) and by Lower
Saxony - Israel cooperation grant no. 76251-99-7/20 (ZN 3637).

##  Appendix A The $e-ph$ collision term and the $e-ph$ coupling coefficient

###  A.1 The phase-space argument for $e-ph$ scattering

In the main text, we claimed that the $e-ph$ collision term, the corresponding
energy exchange rate and coupling coefficient can be dramatically
overestimated in ITO if momentum conservation is not taken into consideration.
However, since it is difficult to deduce how conservation of momentum is
manifested in Eqs. (9)-(12), we provide below a detailed phase-space argument.

Let us consider an electron which initially has an energy $\mathcal{E}_{1}$
interacting with phonons having energy $\mathcal{E}_{ph,1}$ or
$\mathcal{E}_{ph,2}$ as shown in green and red diamonds respectively in Figs.
5(a) and (c). Without loss of generality, we assume that the electron
initially has momentum
$(0,k_{1}=\sqrt{2m_{e}^{\ast}\mathcal{E}_{1}(1+C\mathcal{E}_{1})}/\hbar,0)$
represented by the black dot in momentum space, see Fig. 5(b). If the electron
absorbs the phonon with energy $\mathcal{E}_{ph,1}$, the energy of the
electron is changed to $\mathcal{E}_{1}+\mathcal{E}_{ph,1}$ due to the
conservation of energy. The possible final states can then be represented by a
sphere (a circle on the $k_{z}=0$ plane) whose radius is equal to the
magnitude of the momentum
$\sqrt{2m_{e}^{\ast}(\mathcal{E}+\mathcal{E}_{ph,1})(1+C(\mathcal{E}+\mathcal{E}_{ph,1}))}/\hbar$
represented by the green circle in Fig 5(b). Meanwhile, the momentum of the
electron is changed to $(0,k_{1},0)+{\bf q}_{1}$ due to conservation of
momentum, where ${\bf q}_{1}$ is the momentum of the absorbed phonon
satisfying the linear energy-momentum dispersion
$q_{1}=\mathcal{E}_{ph}^{1}/(\hbar v_{ph})$ as shown in Fig 5(c). The possible
final states satisfying momentum conservation can then be represented by a
sphere (a circle on the $k_{z}=0$ plane) centered at $(0,k_{1},0)$ and having
a radius of $q_{1}$ (the green dashed circle in Fig 5(b)) in momentum space.
Therefore, the true possible final electron states satisfying both energy and
momentum conservation can then be identified by the intersection of these two
spheres in the 3D momentum space (two circles on the $k_{z}=0$ plane), as
shown by the green dot in Fig. 5(b).

Figure 5: (Color online) (a) Color plot showing the possibility of interaction
between electrons with energy $\mathcal{E}$ and phonons with energy
$\mathcal{E}_{ph}$: (green) interacting; (red) non-interacting. The inset
shows that electrons with energy $\mathcal{E}_{1}$ can absorb a phonon with
energy $\mathcal{E}_{ph,1}$ but cannot absorb a phonon with energy
$\mathcal{E}_{ph,2}$. (b) The phase-space argument of the $e-ph$ interaction.
The black, green and red circles respectively represent the states having
energies $\mathcal{E}_{1}$, $\mathcal{E}_{1}+\mathcal{E}_{ph,1}$ and
$\mathcal{E}_{1}+\mathcal{E}_{ph,2}$ on the $k_{z}=0$ plane in the momentum
space. The green (red) dashed circles represents the possible states that can
be reached if an electron with energy $\mathcal{E}_{1}$ and momentum
$(0,k_{1},0)$ (black dot) absorbed a phonon with energy $\mathcal{E}_{ph,1}$
($\mathcal{E}_{ph,2}$). The intersection between the green solid and dashed
circles represents the final electron state of the scattering process. No
intersection between the red solid and dashed circles indicates that no final
state can be reached. (c) The linear energy-momentum relation of the acoustic
phonon. The phonon states with energy $\mathcal{E}_{ph,1}$
($\mathcal{E}_{ph,2}$) and momentum $q_{1}$ ($q_{2}$) are represented by the
green (red) diamonds.

Now, if the initial electron interacted with the phonon with energy
$\mathcal{E}_{ph,2}$ and momentum ${\bf q}_{2}$, the final electron would have
energy $\mathcal{E}_{1}+\mathcal{E}_{ph,2}$ (represented by the red circle in
Fig. 5(b)) and momentum $(0,k_{1},0)+{\bf q}_{2}$ (represented by the red
dashed circle in Fig. 5(b)). However, these two circles (two spheres in the 3D
momentum space) do not intersect each other, indicating that no final state
can be reached and thus the electrons with energy $\mathcal{E}_{1}$ do not
interact with the phonons with energy $\mathcal{E}_{ph,2}$.

Since the phonon energy is usually much smaller than the electron energy
($\mathcal{E}_{D}\ll\mathcal{E}_{1}$ in the above example), the initial and
final energy ($\mathcal{E}_{1}$ and $\mathcal{E}_{1}+\mathcal{E}_{ph,1}$,
respectively) of the electron are very close to each other (that is why the
black, green and red circles nearly overlap with each other in Fig. 5(b)).
However, the initial and final momentum of the electron can be very different
and the magnitude of the momentum difference can range from 0 to $\sim
2k_{1}$. Therefore, electrons with momentum $k_{1}$ (energy $\mathcal{E}_{1}$)
interact only with phonons having momentum smaller than $\sim 2k_{1}$ (energy
smaller then
$\sqrt{8m_{e}^{\ast}v_{ph}^{2}\mathcal{E}_{1}(1+C\mathcal{E}_{1})}$,
represented by the light green line in Fig. 5(a) and (c)) 101010This explains
why the polar integral in Eq. (9) is non-zero only when $q\lesssim 2k$.. This
happens for electrons having momentum smaller than $q_{D}/2$ (or energy
smaller than
$\mathcal{E}_{\text{min,D}}=\dfrac{1}{2C}\left(\sqrt{1+\dfrac{C\mathcal{E}_{D}^{2}}{2m_{e}^{\ast}v_{ph}^{2}}}-1\right)$),
see Fig. 5(a). For ITO, $q_{D}/2\approx 2k_{F}$ is much larger than the Fermi
momentum ($\mathcal{E}_{\text{min,D}}\approx 2.46$ eV is much larger than its
Fermi energy), such that a substantial number of phonons are prohibited from
interacting with electrons (see Fig. 5(a)), resulting in a much smaller $e-ph$
coupling coefficient than that without accounting for momentum conservation in
the electron-phonon collisions (see the comparison shown in Fig. 3(c)).

###  A.2 The function shape of the $e-ph$ collision term in the Boltzmann
equation

The conservation of momentum not only reduces the number of phonons available
for $e-ph$ collisions, but also causes $\left(\dfrac{\partial
f(\mathcal{E})}{\partial t}\right)_{e-ph}$ to exhibit a very different shape
from the shape characteristic of noble metals (see, e.g., [27, 56, 36]), as
shown in Fig. 2(c), (g), (k) and Fig. 6(a)-(c). To have a deeper understanding
of this, we simplify Eq. (10) by expanding its integrand in a power series in
$\mathcal{E}_{ph}$. After some algebra, we find that $\left(\dfrac{\partial
f(\mathcal{E})}{\partial t}\right)_{e-ph}$ is dominated by three terms,
namely,

$\left(\dfrac{\partial f}{\partial
t}\right)_{e-ph}\approx\dfrac{D^{2}}{4\pi\rho(\hbar
v_{ph})^{4}}\sqrt{\dfrac{m_{e}^{\ast}}{2\mathcal{E}(1+C\mathcal{E})}}(1+2C\mathcal{E})\\\
\Bigg{[}\dfrac{\left(\mathcal{E}_{ph}^{\text{max}}\right)^{4}}{4}\dfrac{\partial}{\partial\mathcal{E}}\left(f(\mathcal{E})(1-f(\mathcal{E}))\right)+\dfrac{\left(\mathcal{E}_{ph}^{\text{max}}\right)^{5}}{5}\dfrac{\partial^{2}}{\partial\mathcal{E}^{2}}f(\mathcal{E})\\\
+96(m_{e}^{\ast}v_{ph}^{2})^{2}\mathcal{E}(1+C\mathcal{E})(1+2C\mathcal{E})f(\mathcal{E})\left(1-f(\mathcal{E})\right)H(\mathcal{E})H(\mathcal{E}_{\text{min,D}}-\mathcal{E})\Bigg{]},$
(29)

where
$\mathcal{E}_{ph}^{\text{max}}(\mathcal{E})=\text{min}(\mathcal{E}_{D},\sqrt{8m_{e}^{\ast}v_{ph}^{2}\mathcal{E}(1+C\mathcal{E})})$
is the maximal energy of a phonon which can be absorbed/emitted by an electron
with energy $\mathcal{E}$ (the boundary between the green and red regimes in
Fig. 6(a)). The first term is proportional to the first derivative of
$f(\mathcal{E})(1-f(\mathcal{E}))$ with respect to $\mathcal{E}$. The second
term is proportional to the second derivative of $f(\mathcal{E})$ with respect
to $\mathcal{E}$. The third term is proportional to
$f(\mathcal{E})(1-f(\mathcal{E}))$; it is non-zero only for
$0<\mathcal{E}<\mathcal{E}_{\text{min,D}}$ (see the two Heaviside step
functions) and it ensures electron number conservation in the $e-ph$
interaction, i.e.,
$\displaystyle\int\rho_{e}(\mathcal{E})\left(\dfrac{\partial f}{\partial
t}\right)_{e-ph}d\mathcal{E}=0$. For noble metals,
$\mathcal{E}_{\text{min,D}}$ is much smaller than Fermi energy such that
$\mathcal{E}_{ph}^{\text{max}}(\mathcal{E})=\mathcal{E}_{D}$ for most of the
electrons and the contribution from the third term becomes negligible. In this
case, Eq. (29) reduces to the usual differential form of the $e-ph$ collision
[27].

Fig. 6(d)-(l) show the three terms in Eq. (29) as a function of electron
energy at three different time, before ($t=-110$ fs), at the centre of ($t=0$
fs) and after ($t=+110$ fs) the peak of the pulse. The shape of these three
terms can be explained by the (smeared) multi-stair-step structure of
$f(\mathcal{E})$ (see Fig. 2(a)). Since the third term in Eq. (29) has the
simplest form (proportional to $\sim f(\mathcal{E})(1-f(\mathcal{E}))$) and
the first term is proportional to the first derivative of
$f(\mathcal{E})(1-f(\mathcal{E}))$, we start with explaining the function
shape of the third term; then the first term; and lastly the second term.

Figure 6: (Color online) (a)-(c) the $e-ph$ collision term Eq. (10) (the same
as Fig. 2(c), (g) and (k)), (d)-(f) the first term, (g)-(i) the second term,
and (j)-(l) the third term of Eq. (29) following illumination of ITO by a
short pulse. The left column shows the various terms at the FWHM before the
peak of the pulse ($t=-110$ fs), the middle column at the centre of the pulse
($t=0$ fs), and the right column at the FWHM after the peak of the pulse
($t=+110$ fs) as a function of $\mathcal{E}/\mathcal{E}_{F}$ and pulse
intensity $I_{\text{inc}}=2.5$ GW/cm2. The $y$-axis is in symmetric log-scale.

To explain the function shape of the third term, we first look at the multi-
stair-step shape of $f(\mathcal{E})$ at $t=-110$ fs as shown in Fig. 2(a). The
step of $f(\mathcal{E})$ near $\mathcal{E}_{F}$ is due to the Fermi–Dirac
nature. Around this step, $f(\mathcal{E})$ changes from $\sim 1$ to $\sim 0$
while $1-f(\mathcal{E})$ changes from $\sim 0$ to $\sim 1$. This leads to a
peak in the third term of Eq. (29) ($\sim f(\mathcal{E})(1-f(\mathcal{E}))$)
near the Fermi energy, see Fig. 6(j). The small step of $f(\mathcal{E})$ near
$\mathcal{E}_{F}+\hbar\omega$ is created by the photon absorption (the non-
thermal shoulder). Around this step, $f(\mathcal{E})$ changes from $\sim
10^{-2}$ to $\sim 10^{-4}$, while $1-f(\mathcal{E})$ is nearly equal to 1.
This causes the third term of Eq. (29) ($\sim
f(\mathcal{E})(1-f(\mathcal{E}))$) to have a step-like shape near
$\mathcal{E}_{F}+\hbar\omega$ (instead of a peak), see Fig. 6(j). For $t=0$
and $t=110$ fs, due to the $e-e$ collision, the electron distribution is
smeared out (see Fig. 2(e) and (i)). This also smears out the peak and the
step of the third term in Eq. (29), see Fig. 6(k) and (l). The peak and the
step of $f(\mathcal{E})(1-f(\mathcal{E}))$ then respectively lead to a
Lorentzian dispersion shape near the Fermi-energy and a dip near
$\mathcal{E}_{F}+\hbar\omega$ in the first term of Eq. (29) since it is
proportional to the first derivative of $f(\mathcal{E})(1-f(\mathcal{E}))$,
see Fig. 6(d)-(f). Finally, since the second term of Eq. (29) is proportional
to the second derivative of $f(\mathcal{E})$ with respect to $\mathcal{E}$, it
has a Lorentzian dispersion shape near $\mathcal{E}_{F}$ and
$\mathcal{E}_{F}+\hbar\omega$, see Fig. 6(g)-(i).

Both the first and the second terms have a Lorentzian dispersion shape near
$\mathcal{E}_{F}$ but they are in opposite sign, the combination of these two
terms thus also has a Lorentzian dispersion shape. Combining this with the
peak from the third term results in the complicated shape of
$\left(\dfrac{\partial f}{\partial t}\right)_{e-ph}$ near the Fermi energy.
Finally, the dip of $\left(\dfrac{\partial f}{\partial t}\right)_{e-ph}$ near
$\mathcal{E}_{F}+\hbar\omega$ is mainly contributed from the first term in Eq.
(29), see Fig. 6(a)-(c).

###  A.3 The analytical expression of the $e-ph$ coupling coefficient
$G_{e-ph}$

In this Appendix, we provide the analytical expression of the
$T_{e}$-dependent $e-ph$ coupling coefficient $G_{e-ph}$ (Eq. (26)). We follow
the procedure mentioned in the main text, exchange the integral order, and
separate the right-hand side of Eq. (26) into two terms,
$G_{e-ph}(T_{e})=G^{(0)}_{e-ph}(T_{e})+G^{(1)}_{e-ph}(T_{e})$,

$G^{(0)}_{e-ph}(T_{e})\approx\dfrac{D^{2}}{4\pi^{3}\rho(\hbar
v_{ph})^{4}}\dfrac{m_{e}^{\ast
2}k_{B}}{\hbar^{3}}\int_{0}^{\mathcal{E}_{D}}d\mathcal{E}_{ph}(\mathcal{E}_{ph})^{3}\int_{0}^{\infty}d\mathcal{E}\dfrac{\mathcal{E}(1+2C\mathcal{E})}{(2k_{B}T_{e})^{2}}\\\
\Big{[}(1+2C\mathcal{E})\sech^{2}\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)\tanh\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)-4Ck_{B}T_{e}\sech^{2}\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)\Big{]}$
(30)

and

$G^{(1)}_{e-ph}(T_{e})\approx-\dfrac{D^{2}}{4\pi^{3}\rho(\hbar
v_{ph})^{4}}\dfrac{m_{e}^{\ast
2}k_{B}}{\hbar^{3}}\Bigg{\\{}\int_{0}^{\mathcal{E}_{D}}d\mathcal{E}_{ph}(\mathcal{E}_{ph})^{3}\int_{0}^{\mathcal{E}^{\text{min}}}\dfrac{\mathcal{E}(1+2C\mathcal{E})}{(2k_{B}T_{e})^{2}}d\mathcal{E}\\\
\Big{[}(1+2C\mathcal{E})\sech^{2}\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)\tanh\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)-4Ck_{B}T_{e}\sech^{2}\left(\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}\right)\Big{]}\\\
+\int_{0}^{\mathcal{E}_{D}}d\mathcal{E}_{ph}(\mathcal{E}_{ph})^{3}\dfrac{\mathcal{E}^{\text{min}}}{4k_{B}T_{e}}(1+2C\mathcal{E}^{\text{min}})^{2}\sech^{2}\left(\dfrac{\mathcal{E}^{\text{min}}-\mu}{2k_{B}T_{e}}\right)\Bigg{\\}},$
(31)

where
$\mathcal{E}^{\text{min}}(\mathcal{E}_{ph})=\dfrac{1}{2C}\left(\sqrt{1+\dfrac{C(\mathcal{E}_{ph})^{2}}{2m_{e}^{\ast}v_{ph}^{2}}}-1\right)$
is the minimum energy of an electron which can absorb/emit a phonon with
energy $\mathcal{E}_{ph}$. Next, we change the variables from $\mathcal{E}$ to
$x=\dfrac{\mathcal{E}-\mu}{2k_{B}T_{e}}$ and from $\mathcal{E}_{ph}$ to
$x^{\text{min}}=\dfrac{\mathcal{E}^{\text{min}}-\mu}{2k_{B}T_{e}}$ so that
$\displaystyle\int_{0}^{\infty}d\mathcal{E}=\displaystyle\int_{x_{0}}^{\infty}(2k_{B}T_{e})dx$,
$\displaystyle\int_{0}^{\mathcal{E}^{\text{min}}}d\mathcal{E}=\displaystyle\int_{x_{0}}^{x^{\text{min}}}(2k_{B}T_{e})dx$
and
$\displaystyle\int_{0}^{\mathcal{E}_{D}}d\mathcal{E}_{ph}\mathcal{E}_{ph}=\dfrac{1}{2}(8m_{e}^{\ast}v_{ph}^{2})(2k_{B}T_{e})\displaystyle\int_{x_{0}}^{x_{D}}dx^{\text{min}}(1+4Ck_{B}T_{e}(x^{\text{min}}-x_{0}))$,
where $x_{0}=-\dfrac{\mu}{k_{B}T_{e}}$,
$x_{D}=\dfrac{\mathcal{E}_{\text{min,D}}-\mu}{k_{B}T_{e}}$ and
$\mathcal{E}_{\text{min,D}}=\dfrac{1}{2C}\left(\sqrt{1+\dfrac{C\mathcal{E}_{D}^{2}}{2m_{e}^{\ast}v_{ph}^{2}}}-1\right)$.
Then, Eqs. (30) and (31) becomes

$G^{(0)}_{\text{e-ph}}\approx\dfrac{D^{2}}{4\pi^{3}\rho(\hbar
v_{ph})^{4}}\dfrac{m_{e}^{\ast
2}k_{B}}{\hbar^{3}}\dfrac{\mathcal{E}_{D}^{4}}{4}\int_{x_{0}}^{\infty}dx\Big{[}(x-x_{0})(1+4Ck_{B}T_{e}(x-x_{0}))^{2}\sech^{2}(x)\tanh{x}\\\
-4Ck_{B}T_{e}(x-x_{0})(1+4Ck_{B}T_{e}(x-x_{0}))\sech^{2}(x)\Big{]},$ (32)

and

$G^{(1)}_{\text{e-ph}}\approx-\dfrac{D^{2}}{4\pi^{3}\rho(\hbar
v_{ph})^{4}}\dfrac{m_{e}^{\ast
2}k_{B}}{\hbar^{3}}(8m_{e}^{\ast}v_{ph}^{2})^{2}(2k_{B}T_{e})^{2}\\\
\Bigg{\\{}\dfrac{1}{2}\int_{x_{0}}^{x_{D}}dx^{\text{min}}(x^{\text{min}}-x_{0})(1+2Ck_{B}T_{e}(x^{\text{min}}-x_{0}))(1+4Ck_{B}T_{e}(x^{\text{min}}-x_{0}))\\\
\int_{x_{0}}^{x^{\text{min}}}dx\Big{[}(x-x_{0})(1+4Ck_{B}T_{e}(x-x_{0}))^{2}\sech^{2}(x)\tanh{x}\\\
-4Ck_{B}T_{e}(x-x_{0})(1+4Ck_{B}T_{e}(x-x_{0}))\sech^{2}(x)\Big{]}\\\
+\dfrac{1}{4}\int_{x_{0}}^{x_{D}}dx^{\text{min}}(x^{\text{min}}-x_{0})^{2}(1+2Ck_{B}T_{e}(x^{\text{min}}-x_{0}))\\\
(1+4Ck_{B}T_{e}(x^{\text{min}}-x_{0}))^{3}\sech^{2}(x^{\text{min}})\Bigg{\\}}.$
(33)

To evaluate Eqs. (32) and (33) analytically, we define

$\displaystyle g_{n}(x)\equiv\dfrac{x^{n}}{n!}-\ln
2\dfrac{x^{n-1}}{(n-1)!}-\left(-\dfrac{1}{2}\right)^{n-1}\text{Li}_{n}(-e^{-2x}),\
n\in\mathbb{Z}$ (34)

One can verify that $g_{n}(x)=\dfrac{dg_{n+1}(x)}{dx}$,
$g_{-2}(x)=-2\sech^{2}(x)\tanh(x)$ and $g_{-1}(x)=\sech^{2}(x)$. The integral
in Eqs. (32) and (33) can thus be expressed using

$\int_{x_{0}}^{x^{\text{min}}}dx(x-x_{0})^{m}g_{n}(x)=\\\
\sum_{r=0}^{m}\left[\dfrac{(-1)^{r}m!}{(m-r)!}(x^{\text{min}}-x_{0})^{m-r}g_{n+r+1}(x^{\text{min}})\right]-(-1)^{m}g_{n+m+1}(x_{0})$
(35)

and

$\int_{x_{0}}^{x_{D}}dx^{\text{min}}(x^{\text{min}}-x_{0})^{k}\int_{x_{0}}^{x^{\text{min}}}dx(x-x_{0})^{m}g_{n}(x)=\\\
\sum_{r=0}^{m}\sum_{s=0}^{m+k-r}\left[\dfrac{(-1)^{r+s}m!}{(m-r)!}\dfrac{(m+k-r)!}{(m+k-r-s)!}(x_{D}-x_{0})^{m+k-r-s}g_{n+r+s+2}(x_{D})\right]\\\
-\sum_{r=0}^{m}\left[\dfrac{(-1)^{m+k}m!(m+k-r)!}{(m-r)!}g_{n+m+k+2}(x_{0})\right]-(-1)^{m}\dfrac{(x_{D}-x_{0})^{k+1}}{(k+1)!}g_{n+m+1}(x_{0}),$
(36)

where $m$ and $k$ are positive integers.

One can then obtain the full analytical expression of $G_{e-ph}(T_{e})$ by
substituting Eqs. (35)-(36) into Eq. (32) and (33). Since the full expression
is too long, we shall not present it here. At low electron temperatures, i.e.,
$k_{B}T_{e}\ll\mu$, we have $e^{-\frac{\mu}{k_{B}T_{e}}}\rightarrow 0$,
$\ln\left(1+e^{\frac{\mu}{k_{B}T_{e}}}\right)\rightarrow\dfrac{\mu}{k_{B}T_{e}}$
and
$\text{Li}_{2}\left(-e^{-\frac{\mu}{k_{B}T_{e}}}\right)\rightarrow-2\left(\dfrac{\mu}{2k_{B}T_{e}}\right)^{2}$,
then the sum of Eq. (32) and Eq. (33) reproduces Eq. (28).

## References

  * Naik  _et al._ [2011] G. V. Naik, J. Kim, and A. Boltasseva, Oxides and nitrides as alternative plasmonic materials in the optical range, Optical Materials Express 1, 1090 (2011).
  * Alam  _et al._ [2016] M. Z. Alam, I. D. Leon, and R. W. Boyd, Large optical nonlinearity of indium tin oxide in its epsilon-near-zero region, Science 116 (2016).
  * Diroll  _et al._ [2020] B. T. Diroll, S. Saha, V. M. Shalaev, A. Boltasseva, and R. D. Schaller, Broadband ultrafast dynamics of refractory metals: TiN and ZrN, Adv. Opt. Mater. 8, 2000652 (2020).
  * Yang  _et al._ [2017] Y. Yang, K. Kelley, E. Sachet, S. Campione, T. S. Luk, J.-P. Maria, M. B. Sinclair, and I. Brener, Femtosecond optical polarization switching using a cadmium oxide-based perfect absorber, Nat. Photonics 11, 390 (2017).
  * Kinsey and Khurgin [2019] N. Kinsey and J. Khurgin, Nonlinear epsilon-near-zero materials explained: opinion, Opt. Mater. Express 9, 2793 (2019).
  * Caspani  _et al._ [2016] L. Caspani, R. P. M. Kaipurath, M. Clerici, M. Ferrera, T. Roger, J. Kim, N. Kinsey, M. Pietrzyk, A. D. Falco, V. M. Shalaev, A. Boltasseva, and D. Faccio, Enhanced nonlinear refractive index in $\epsilon$-near-zero materials, Phys. Rev. Lett. 116, 233901 (2016).
  * Feigenbaum  _et al._ [2010] E. Feigenbaum, K. Diest, and H. Atwater, Unity-order index change in transparent conducting oxides at visible frequencies, Nano Lett. 10, 2111–2116 (2010).
  * Liu  _et al._ [2014] X. Liu, J. Park, J.-H. Kang, H. Yuan, Y. Cui, H. Y. Hwang, and M. L. Brongersma, Quantification and impact of nonparabolicity of the conduction band of indium tin oxide on its plasmonic properties, Appl. Phys. Lett. 105, 181117 (2014).
  * Pradhan  _et al._ [2014] A. Pradhan, R. Mundle, K. Santiago, J. Skuza, B. Xiao, K. Song, M. Bahoura, R. Cheaito, and P. E. Hopkins, Extreme tunability in aluminum doped zinc oxide plasmonic materials for near-infrared applications, Scientific reports 4, 1 (2014).
  * Wu  _et al._ [2021] J. Wu, X. Liu, H. Fu, K.-C. Chang, S. Zhang, and H. Y. F. nd Q. Li, Manipulation of epsilon-near-zero wavelength for the optimization of linear and nonlinear absorption by supercritical fluid, Scientific Reports 11, 15936 (2021).
  * Alam  _et al._ [2018] M. Z. Alam, S. A. Schulz, J. Upham, I. D. Leon, and R. W. Boyd, Large optical nonlinearity of nanoantennas coupled to an epsilon-near-zero material, Nat. Photonics 12, 79 (2018).
  * Bohn  _et al._ [2021] J. Bohn, T. S. Luk, C. Tollerton, S. Hutchins, I. Brener, S. Horsley, W. L. Barnes, and E. Hendry, All-optical switching of an epsilon-near-zero plasmon resonance in indium tin oxide, Nat. Commun. 12, 1017 (2021).
  * Silveirinha and Engheta [2006] M. Silveirinha and N. Engheta, Tunneling of electromagnetic energy through subwavelength channels and bends using permittivity-near-zero materials, Phys. Rev. Lett. 97, 157403 (2006).
  * Liu  _et al._ [2008a] R. Liu, Q. Cheng, T. Hand, J. J. Mock, T. J. Cui, S. A. Cummer, and D. R. Smith, Experimental demonstration of electromagnetic tunneling through an epsilon-near-zero metamaterial at microwave frequencies, Phys. Rev. Lett. 100, 023903 (2008a).
  * Alu  _et al._ [2007] A. Alu, M. G. Silveirinha, A. Salandrino, and N. Engheta, Epsilon-near-zero metamaterials and electromagnetic sources: Tailoring the radiation phase pattern, Phys. Rev. B 75, 155410 (2007).
  * Maas  _et al._ [2013] R. Maas, J. Parsons, N. Engheta, and A. Polman, Experimental realization of an epsilon-near-zero metamaterial at visible wavelengths, Nat. Photonics 7, 907 (2013).
  * Guo  _et al._ [2016a] P. Guo, R. D. Schaller, L. E. Ocola, B. T. Diroll, J. B. Ketterson, and R. P. H. Chang, Large optical nonlinearity of ito nanorods for sub-picosecond all-optical modulation of the full-visible spectrum, Nat. Commun. 7, 12892 (2016a).
  * Guo  _et al._ [2017] Q. Guo, Y. Cui, Y. Yao, Y. Ye, Y. Yang, X. Liu, S. Zhang, X. Liu, J. Qiu, and H. Hosono, A solution-processed ultrafast optical switch based on a nanostructured epsilon-near-zero medium, Advanced Materials 29, 1700754 (2017).
  * Tirole  _et al._ [2022] R. Tirole, E. Galiffi, J. Dranczewski, T. Attavar, B. Tilmann, Y.-T. Wang, P. A. Huidobro, A. Alú, J. B. Pendry, S. A. Maier, S. Vezzoli, and R. Sapienza, Saturable time-varying mirror based on an ENZ material, (2022).
  * Guo  _et al._ [2016b] P. Guo, R. D. Schaller, J. B. Ketterson, and R. P. H. Chang, Ultrafast switching of tunable infrared plasmons in indium tin oxide nanorod arrays with large absolute amplitude, Nat. Photonics 10, 267 (2016b).
  * Wang  _et al._ [2019] H. Wang, K. Du, C. Jiang, Z. Yang, L. Ren, W. Zhang, S. J. Chua, and T. Mei, Extended drude model for intraband-transition-induced optical nonlinearity, Phys. Rev. Applied 11, 064062 (2019).
  * Ziman [1972] J. M. Ziman,  _Principles of the theory of solids_ (Cambridge University Press, 1972).
  * Ashcroft and Mermin [1976] N. W. Ashcroft and N. D. Mermin,  _Solid state physics_ (Brooks/Cole, 1976).
  * Lundstrom [1990] M. Lundstrom,  _Fundamentals of carrier transport_ (Addison-Wesley, 1990).
  * Groeneveld  _et al._ [1995] R. H. M. Groeneveld, R. Sprik, and A. Lagendijk, Femtosecond spectroscopy of electron-electron and electron-phonon energy relaxation in Ag and Au, Phys. Rev. B 51, 11433 (1995).
  * Fatti  _et al._ [2000] N. D. Fatti, C. Voisin, M. Achermann, S. Tzortzakis, D. Christofilos, and F. Valleé, Nonequilibrium electron dynamics in noble metals, Phys. Rev. B 61, 16956 (2000).
  * Grua  _et al._ [2003] P. Grua, J. P. Morreeuw, H. Bercegol, G. Jonusauskas, and F. Valleé, Electron kinetics and emission for metal nanoparticles exposed to intense laser pulses, Phys. Rev. B 68, 035424 (2003).
  * Pietanza  _et al._ [2007] L. D. Pietanza, G. Colonna, S. Longo, and M. Capitelli, Non-equilibrium electron and phonon dynamics in metals under femtosecond laser pulses, Eur. Phys. J. D 45, 369 (2007).
  * Kornbluth  _et al._ [2013] M. Kornbluth, A. Nitzan, and T. Seidman, Light-induced electronic non-equilibrium in plasmonic particles, J. Chem. Phys. 138, 174707 (2013).
  * Saavedra  _et al._ [2016] J. R. M. Saavedra, A. Asenjo-García, and F. J. G. de Abajo, Hot-electron dynamics and thermalization in small metallic nanoparticles, ACS Photonics 3, 1637 (2016).
  * Kane [1957] E. O. Kane, Band structure of indium antimonide, J. Phys. Chem. Solids 1, 249 (1957).
  * Pisarkiewicz and Kolodziej [1990] T. Pisarkiewicz and A. Kolodziej, Nonparabolicity of the conduction band structure in degenerate tin dioxide, Phys. Status Solidi B 158, K5 (1990).
  * Franzen [2008] S. Franzen, Surface plasmon polaritons and screened plasma absorption in indium tin oxide compared to silver and gold, J. Phys. Chem. C 112, 6027 (2008).
  * Khurgin and Levy [2020] J. B. Khurgin and U. Levy, Generating hot carriers in plasmonic nanoparticles: When quantization does matter?, ACS Photonics 7, 547 (2020).
  * Rethfeld  _et al._ [2002] B. Rethfeld, A. Kaiser, M. Vicanek, and G. Simon, Ultrafast dynamics of nonequilibrium electrons in metals under femtosecond laser irradiation, Phys. Rev. B 65, 214303 (2002).
  * Dubi and Sivan [2019] Y. Dubi and Y. Sivan, “hot electrons” in metallic nanostructures - non-thermal carriers or heating?, Light: Sci. Appl. 8, 89 (2019).
  * Besteiro  _et al._ [2017] L. V. Besteiro, X.-T. Kong, Z. Wang, G. Hartland, and A. O. Govorov, Understanding hot-electron generation and plasmon relaxation in metal nanocrystals: Quantum and classical mechanisms, ACS Photonics 4, 2759 (2017).
  * Mryasov and Freeman [2001] O. N. Mryasov and A. J. Freeman, Electronic band structure of indium tin oxide and criteria for transparent conducting behavior, Phys. Rev. B 64, 233111 (2001).
  * Lin and Li [2014] J.-J. Lin and Z.-Q. Li, Electronic conduction properties of indium tin oxide: single-particle and many-body transport, J. Phys. Condens. Matter 26, 343201 (2014).
  * Liu  _et al._ [2008b] X. D. Liu, E. Y. Jiang, and D. X. Zhang, Electrical transport properties in indium tin oxide films prepared by electron-beam evaporation, J. Appl. Phys. 104, 073711 (2008b).
  * Brown  _et al._ [2016] A. M. Brown, R. Sundararaman, P. Narang, W. A. Goddard, and H. A. Atwater, Nonradiative plasmon decay and hot carrier dynamics: Effects of phonons, surfaces, and geometry, ACS Nano 10, 957 (2016).
  * Khurgin  _et al._ [2017] J. Khurgin, W.-Y. Tsai, D. P. Tsai, and G. Sun, Landau damping and limit to field confinement and enhancement in plasmonic dimers, ACS Photonics 4, 2871 (2017).
  * Note [1] For simplicity, we neglect any inhomogeneity of the local electric field in the ITO sample.
  * Snoke [2020] D. W. Snoke,  _Solid State Physics: Essential Concepts_ , 2nd ed. (Pearson/Addison-Wesley, Boston, MA, 2020).
  * Brinzari  _et al._ [2016] V. Brinzari, D. Nika, I. Damaskin, B. Cho, and G. Korotcenkov, Thermoelectric properties of nano-granular indium–tin-oxide within modified electron filtering model with chemisorption-type potential barriers, Physica E 81, 49 (2016).
  * Note [2] This value of the deformation potential was obtained by fitting the ITO permittivity (at room temperature) calculated by the Lindhard formula [75] with that measured experimentally [21].
  * Wang  _et al._ [2020] H. Wang, K. Du, R. Liu, X. Dai, W. Zhang, S. J. Chua, and T. Mei, Role of hot electron scattering in epsilon-near-zero optical nonlinearity, Nanophotonics 9, 4287 (2020).
  * Kim  _et al._ [2020] D. Kim, A. Aydin, A. Daza, K. N. Avanaki, J. Keski-Rahkonen, and E. J. Heller, Coherent electron dynamics in thermal lattice vibrations (2020).
  * Ridley [2013] B. K. Ridley,  _Quantum Processes in Semiconductors_ , 5th ed. (Oxford University Press, Oxford, 2013).
  * Coleman [2015] P. Coleman,  _Introduction to many body physics_ (Cambridge University Press, 2015).
  * Smith and Norris [2001] A. N. Smith and P. M. Norris, Influence of intraband transitions on the electron thermoreflectance response of metals, Appl. Phys. Lett 78, 1240 (2001).
  * Note [3] From Eq. (15), one can deduce that the Thomas-Fermi wavevector is proportional to the square root of the eDOS at the Fermi-energy [23]. Therefore, compared with the case of $C=0$, $q_{TF}$ is larger by a factor of $(1+C\mathcal{E}_{F})^{1/4}(1+2C\mathcal{E}_{F})^{1/2}$ and $k_{F}$ is larger by a factor of $(1+C\mathcal{E}_{F})^{1/2}$. Together with the factor $(1+2C\mathcal{E}_{F})^{3}$ coming from the increase of eDOS, the $e-e$ thermalization rate increases by a factor of $\sim(1+2C\mathcal{E}_{F})^{3/2}(1+C\mathcal{E}_{F})^{-5/4}$.
  * Note [4] Note that the term $(\mathcal{E}-\mathcal{E}_{F})$ in Eq. (18) should not be replaced by $\hbar\omega$, since it involves all possible electron states, rather than only those excited from the Fermi energy by photon absorption.
  * Note [5] The Thomas-Fermi wavevector is proportional to the square root of the eDOS at the Fermi-energy [23], thus, it decreases with electron density.
  * Sivan and Dubi [2021] Y. Sivan and Y. Dubi, Theory of “hot” photoluminescence from Drude metals, ACS Nano 15, 8724 (2021).
  * Stoll  _et al._ [2014] T. Stoll, P. Maioli, A. Crut, N. D. Fatti, and F. Vallée, Advances in femto-nano-optics: ultrafast nonlinearity of metal nanoparticles, Eur. Phys. J. B 87, 260 (2014).
  * Yagi  _et al._ [2005] T. Yagi, K. Tamano, Y. Sato, N. Taketoshi, T. Baba, and Y. Shigesato, Analysis on thermal properties of tin doped indium oxide films by picosecond thermoreflectance measurement, J. Vac. Sci. Technol. A 23, 1180 (2005).
  * Anisimov  _et al._ [1974] S. I. Anisimov, B. L. Kapeilovich, and T. I. Perelman, Electron emission from metal surfaces exposed to ultrashort laser pulses, Sov. Phys. JETP 39, 375 (1974).
  * Note [6] This is reminiscent of what happens in semiconductors [23, 76, 77].
  * Sommerfeld [1928] A. Sommerfeld, Zur elektronentheorie der metalle auf grund der fermischen statistik, Z. Phys. 47, 1 (1928).
  * Note [7] Due to non-parabolicity the dependence of $\gamma_{e}$ on electron density becomes is not straightforward. In contrast, in its absence, $\gamma_{e}$ reduces to the familiar expression $\gamma_{e}=\pi^{2}n_{e}k_{B}^{2}/2\mathcal{E}_{F}$.
  * Lin  _et al._ [2008] Z. Lin, L. V. Zhigilei, and V. Celli, Electron-phonon coupling and electron heat capacity of metals under conditions of strong electron-phonon non-equilibrium, Phys. Rev. B 77, 075133 (2008).
  * Allen [1987] P. B. Allen, Theory of thermal relaxation of electrons in metals, Phys. Rev. Lett. 59, 1460 (1987).
  * Note [8] For Au, Eq. (27) becomes $\displaystyle G_{e-ph}^{(Au)}=\genfrac{}{}{}{0}{D^{2}m_{e}k_{B}\mathcal{E}_{D}^{4}}{16(\pi\hbar)^{3}\rho(\hbar v_{ph})^{4}}$; with the known parameters for Au [26], (viz., deformation potential $D=19.3$ eV, Debye energy $\mathcal{E}_{D}=0.014$ eV, mass density $\rho=19.3\times 10^{3}~{}\text{Kg}/\text{m}^{3}$ and phonon velocity $v_{ph}=3240\text{m/s}$) one obtains the value above.
  * Note [9] The same reasoning explains why an increase in the phonon temperature does not have a significant effect, see Eq. (26).
  * Voisin  _et al._ [2004] C. Voisin, D. Christofilos, P. A. Loukakos, N. D. Fatti, F. Valleé, J. Lermé, M. Gaudry, E. Cottancin, M. Pellarin, and M. Broyer, Ultrafast electron-electron scattering and energy exchanges in noble-metal nanoparticles, Phys. Rev. B 69, 195416 (2004).
  * Sivan and Dubi [2020] Y. Sivan and Y. Dubi, Recent developments in plasmon-assisted photocatalysis - a personal perspective, Appl. Phys. Lett. 117, 130501 (2020).
  * Haug  _et al._ [2015] T. Haug, P. Klemm, S. Bange, and J. M. Lupton, Hot-electron intraband luminescence from single hot spots in noble-metal nanoparticle films, Phys. Rev. Lett. 115, 067403 (2015).
  * Wilson and Coh [2020] R. B. Wilson and S. Coh, Parametric dependence of hot electron relaxation timescales on electron-electron and electron-phonon interaction strengths, Commun. Phys 3, 035424 (2020).
  * Schoenlein  _et al._ [1987] R. W. Schoenlein, W. Z. Lin, J. G. Fujimoto, and G. L. Eesley, Femtosecond studies of nonequilibrium electronic processes in metals, Phys. Rev. Lett 58, 1680 (1987).
  * Carpene [2006] E. Carpene, Ultrafast laser irradiation of metals: Beyond the two-temperature model, Phys. Rev. B 74, 024301 (2006).
  * Valle  _et al._ [2012] G. D. Valle, M. Conforti, S. Longhi, G. Cerullo, and D. Brida, Real-time optical mapping of the dynamics of nonthermal electrons in thin gold films, Phys. Rev. B 86, 155139 (2012).
  * Sivan  _et al._ [2019] Y. Sivan, I. W. Un, and Y. Dubi, Assistance of plasmonic nanostructures to photocatalysis - just a regular heat source, Faraday Discuss. 214, 215 (2019).
  * Note [10] This explains why the polar integral in Eq. (9) is non-zero only when $q\lesssim 2k$.
  * Grosso and Parravicini [2014] G. Grosso and G. Parravicini,  _Solid State Physics_, 2nd ed. (Academic Press, London, 2014).
  * Sze [2006] S. M. Sze, Physics and properties of semiconductors—a review, in  _Physics of Semiconductor Devices_ (John Wiley & Sons, Ltd, 2006) pp. 5–75.
  * Sarkar  _et al._ [2022] S. Sarkar, I. W. Un, Y. Sivan, and Y. Dubi, Theory of non-equilibrium “hot” carriers in direct band-gap semiconductors under continuous illumination, New J. Phys. 24, 053008 (2022).