1 Introduction

The passage of swift ions through semiconductors or insulators generates a large density of electron–hole pairs, which eventually recombine either radiatively or non-radiatively. Most of the energy deposited by the electron–hole pairs is exchanged with the atomic lattice and can yield permanent radiation damage, such as point or extended defects [1]. The pathway to radiation damage involves multiple time and length scales. Three microscopic mechanisms have been proposed for the early stages of this pathway: Coulomb explosion [2, 3], thermal spikes [4,5,6,7], and excitonic [1, 8,9,10]. The excitonic mechanics requires the formation of self-trapped excitons (STEs), which is known to occur in several wide gap insulators, including alkali halides and oxides such as silica (\({{\mathrm {SiO}}_{2}}\)) [11].

The formation of STEs has been experimentally observed in both quartz (\(\alpha \)-\({{\mathrm {SiO}}_{2}}\)) and amorphous silica (a-\({{\mathrm {SiO}}_{2}}\)) [12,13,14]. These groundbreaking femtosecond pump–probe experiments found that free electron–hole pairs relax into STEs in about 150 fs. Subsequent measurements indicate a relaxation time between 50 and 220 fs, suggesting a longer free electron–hole pair relaxation time as the laser intensity is increased [15, 16]. Grojo et al. [17] also observed an abrupt increase of the free electron–hole pair relaxation time in a-\({{\mathrm {SiO}}_{2}}\) as the charge density exceeds \(10^{20}\) \({\mathrm {cm}}^{\mathrm {-3}}\). Critical examination of the experimental conditions and modelling suggests that relaxation to STEs may not be the dominant mechanism for electron–hole densities as large as \(10^{22}\) \({\mathrm {cm}}^{\mathrm {-3}}\) [18]. Since for large electron–hole pair densities Coulomb screening is known to hinder the formation of free excitons—which are the precursor of STEs—Auger recombination will become dominant [19, 20]. The transition from a gas of free excitons to an electron–hole plasma driven by Coulomb screening is known as the exciton Mott transition (EMT) [21,22,23].

Along with the well-established femtosecond pump–probe laser approach, it is now possible to ‘pump’ electron–hole pairs in transparent insulators with ultrashort (about 3.3 ps [24]) proton pulses generated by target normal sheath acceleration [25] and then probe the electron–hole plasma optically [24, 26,27,28]. Both laser and proton irradiation can achieve large power densities (\({{\mathrm {10}}^{18}}\)–\({\mathrm {10}}^{\mathrm {19}}\) \({\mathrm {W/cm}}^{\mathrm {3}}\)) [5]. Protons deposit the energy in a cylinder of only a few nanometres radius along their trajectories, while a laser beam is typically a few micrometers across. It is then possible to locally excite a large density of electron–hole pairs (\({\mathrm {10}}^{\mathrm {19}}\) \({\mathrm {cm}}^{\mathrm {-3}}\)) without permanently damaging the sample on a micrometric scale [28].

Initial experiments with an ultrafast proton ‘pump’ measured the transient opacity of a borosilicate crown glass BK7 using an optical streaking technique [24]. The transient opacity lasted for an unexpectedly long time (hundreds of ps), much longer than the carrier–lattice thermal equilibration time predicted from a typical electron–phonon timescale of about 0.1 ps [5, 6]. This long relaxation time has been confirmed by subsequent experiments [27, 28] and points towards the presence of a ‘cold’—i.e. of temperature comparable with the lattice temperature—persistent electron–hole plasma which can still absorb in the near infrared (1053 nm) over hundreds of ps.

In this article, we investigate the effectiveness of Coulomb screening in determining long free electron–hole pair relaxation times. Our conclusions are based on the numerical solution of a simplified set of hydrodynamic equations for the coupled evolution of the electron–hole and STE densities, alongside with the electronic and lattice temperatures. In particular, we found that Coulomb screening alone does not yield a substantial increase of the relaxation time and in fact it causes a slight reduction of it. On the other hand, long free electron–hole pair relaxation times are possible as a consequence of a combination of Coulomb screening and an effective thermal detrapping of STEs.

2 Hydrodynamic model

Our phenomenological description is based on a simplified version of a set of hydrodynamic equations, commonly known as energy-transport model, originally developed to model ‘hot’ electrons and holes in semiconductor devices [6, 29]. The full set of equations reads

$$\begin{aligned}&\nabla \cdot \left[ \frac{\epsilon _{0}\epsilon _{r}}{q}\nabla \phi \right] =-\left( p-n\right) \;,\nonumber \\&\frac{\partial n}{\partial t}-\nabla \cdot \frac{{\varvec{J}}_{n}}{q}\nonumber \\&\quad =-\left( R_{\mathrm {rad}}+R_{\mathrm {srh}}+R_{\mathrm {aug}}^{\left( n\right) }+R_{\mathrm {aug}}^{\left( p\right) }+R_{\mathrm {ste}}\right) \;,\nonumber \\&\frac{\partial p}{\partial t}{+}\nabla \cdot \frac{{\varvec{J}}_{p}}{q} =-\left( R_{\mathrm {rad}}{+} R_{\mathrm {srh}}+R_{\mathrm {aug}}^{\left( n\right) }+R_{\mathrm {aug}}^{\left( p\right) }+R_{\mathrm {ste}}\right) \;,\nonumber \\&\frac{\partial n_{\mathrm {ste}}}{\partial t} =R_{ste}-\left( \frac{1}{\tau _{\mathrm {ste},r}}+\frac{1}{\tau _{\mathrm {ste},nr}}\right) n_{\mathrm {ste}}\;,\nonumber \\&\frac{3}{2}k_{B}\frac{\partial \left( nT_{n}\right) }{\partial t}+\nabla \cdot {\varvec{S}}_{n}\nonumber \\&\quad =-{\varvec{J}}_{n}\cdot \nabla \phi -\frac{3}{2}k_{B}n\left( \frac{T_{n}-T_{L}}{\tau _{n}}\right) \nonumber \\&\qquad -\frac{3}{2}k_{B}T_{n}R_{\mathrm {rad}}+R_{\mathrm {aug}}^{\left( n\right) }E_{g}\;,\nonumber \\&\frac{3}{2}k_{B}\frac{\partial \left( pT_{p}\right) }{\partial t}+\nabla \cdot {\varvec{S}}_{p}\nonumber \\&\quad =-{\varvec{J}}_{p}\cdot \nabla \phi -\frac{3}{2}k_{B}p\left( \frac{T_{p}-T_{L}}{\tau _{p}}\right) \nonumber \\&\qquad -\frac{3}{2}k_{B}T_{p}R_{\mathrm {rad}}+R_{\mathrm {aug}}^{\left( p\right) }E_{g}\;,\nonumber \\&\rho _{L}c_{L}\frac{\partial T_{L}}{\partial t}+\nabla \cdot {\varvec{S}}_{L} \nonumber \\&\quad =\frac{3}{2}k_{B}n\left( \frac{T_{n}-T_{L}}{\tau _{n}}\right) +\frac{3}{2}k_{B}p\left( \frac{T_{p}-T_{L}}{\tau _{p}}\right) \nonumber \\&\qquad +R_{\mathrm {srh}}E_{g}+R_{\mathrm {ste}}\left( E_{g}-E_{x}\right) \;, \end{aligned}$$
(1)

where q is the unit charge, \(\epsilon _{0}\epsilon _{r}\) gives the static dielectric permittivity of the material, \(\phi \) is the electrostatic potential, n, p and \(n_{\mathrm {ste}}\) represent the electron, hole and STE densities (unit: \({\mathrm {cm}}^{\mathrm {-3}})\), \(T_{n}\), \(T_{p}\) and \(T_{L}\) are the electron, hole and lattice temperatures, \({\varvec{J}}_{n}\), \({\varvec{J}}_{p}\) and \({\varvec{J}}_{L}\) are the charge current densities associated with the electrons, holes and lattice, \({\varvec{S}}_{n}\), \({\varvec{S}}_{p}\) and \({\varvec{S}}_{L}\) are the energy current densities associated with the electrons, holes and lattice , \(\rho _{L}\) gives the lattice density (unit: \(\mathrm {g/cm^{3}})\), \(c_{L}\) is the specific heat capacity of the lattice (unit: \(\mathrm {erg/g\cdot K})\), \(\tau _{n}\) and \(\tau _{p}\) are phenomenological relaxation times based on the two-temperature model of the carrier–lattice thermal equilibration [30], \(E_{g}\) is the energy gap, \(E_{x}\) is the STE energy and \(R_{\mathrm {rad}}\), \(R_{\mathrm {srh}}\), \(R_{\mathrm {aug}}^{\left( n,p\right) }\) and \(R_{\mathrm {ste}}\) are rates (unit: \(\mathrm {cm^{-3}\cdot s^{-1}})\) for radiative recombination, defect-assisted (Shockley–Read–Hall) recombination, Auger recombination (for electrons or holes ) and STE formation. Finally, \(\tau _{\mathrm {ste},r}\) and \(\tau _{\mathrm {ste},nr}\) give the radiative and non-radiative STE decay times, respectively (see Sect. 3).

The constitutive equations for the currents are

$$\begin{aligned} \begin{aligned}\frac{{\varvec{J}}_{n}}{q}&=-M_{n}n\nabla \phi +\frac{k_{B}M_{n}}{q}\nabla \left( nT_{n}\right) \;,\\ \frac{{\varvec{J}}_{p}}{q}&=-M_{p}p\nabla \phi -\frac{k_{B}M_{p}}{q}\nabla \left( pT_{p}\right) \;,\\ {\varvec{S}}_{n}&=-\kappa _{n}\nabla T_{n}-\frac{5}{2}k_{B}T_{n}\frac{{\varvec{J}}_{n}}{q}\;,\\ {\varvec{S}}_{p}&=-\kappa _{p}\nabla T_{p}+\frac{5}{2}k_{B}T_{p}\frac{{\varvec{J}}_{p}}{q}\;,\\ {\varvec{S}}_{L}&=-\kappa _{L}\nabla T_{L}\;, \end{aligned}\nonumber \\ \end{aligned}$$
(2)

where \(k_{B}\) is the Boltzmann constant, \(M_{n}\) and \(M_{p}\) are the phenomenological electron and hole mobilities, while \(\kappa _{n}\), \(\kappa _{p}\) and \(\kappa _{L}\) are the electron, hole and lattice heat conductivities. It is also assumed that

$$\begin{aligned} \begin{aligned} \kappa _{n}&=\left( \frac{5}{2}+r_{n}\right) \frac{k_{B}^{2}}{q}T_{n}M_{n}n\;,\\ \kappa _{p}&=\left( \frac{5}{2}+r_{p}\right) \frac{k_{B}^{2}}{q}T_{p}M_{p}p\;, \end{aligned}\nonumber \\ \end{aligned}$$
(3)

where the coefficients \(r_{n}\) and \(r_{p}\) come from the power-law dependence of the most effective relaxation time as a function of the kinetic energy of the electrons and holes. The power is \(+\frac{3}{2}\) for ionised impurity scattering, \(-\frac{1}{2}\) for acoustic phonon scattering and \(+\frac{1}{2}\) for optical phonon scattering [31]. If radiative recombination is neglected, Eq. (1) conserves the integral of the total energy density

$$\begin{aligned} u_{tot}&=\frac{1}{2}\epsilon _{0}\epsilon _{r}\left( \nabla \phi \right) ^{2}+\frac{E_{g}}{2}\left( n+p\right) +\frac{3}{2}nk_{B}T_{n}+\frac{3}{2}pk_{B}T_{p}\nonumber \\&\quad +n_{\mathrm {ste}}E_{\mathrm {ste}}+u_{L}\;, \end{aligned}$$
(4)

where the lattice energy density is defined so that \(\frac{\partial u_{L}}{\partial t}=\rho _{L}c_{L}\frac{\partial T_{L}}{\partial t}\).

The original hydrodynamic equations have been obtained by taking the moments of the Boltzmann transport equation for semiconductors [32]. The equations for the moments form an exact, yet infinite, hierarchy which in practice needs to be truncated. Here, we have followed the approach introduced in Ref. [33], which approximates the moments of third or higher order using equilibrium Boltzmann statistics, i.e. non-degenerate electrons and holes which have fully equilibrated at temperature \(T_n\) and \(T_p\), respectively. The validity of this approximation can be checked a posteriori by looking, e.g. at the degeneracy factor defined as the argument of the exponential on the right-hand side of Eq. (13) or (14). We found that Boltzmann statistics is indeed accurate away from the EMT, i.e. when the electron–hole pairs are either completely dissociated or bounded. Extensions to the degenerate case are known [34, 35], but are not considered in this article for the sake of simplicity. It is also worth noting that the hydrodynamic equations have been obtained within the relaxation time approximation. However, a different relaxation time is assumed for each of the moments [33]. To this extent, the approach is fundamentally phenomenological and its accuracy relies on the way the phenomenological relaxation times are obtained. A complete discussion of the known limitations of the hydrodynamic model of transport in semiconductors can be found in Ref. [29].

Since in this article we are concerned with the picosecond response of excited transparent insulators, we carry on by neglecting radiative and defect-assisted (Shockley–Read–Hall) recombination. These two processes typically characterise the nanosecond response of semiconductors and insulators [36]. We will also assume local charge neutrality so that we can take \(\phi =0\), \(p=n\) and \(T_{p}=T_{n}\). This is usually a good approximation after a few femtoseconds [6]. To keep the presentation simple and uncluttered, we also neglect the spatial dependency of the fields. This is usually not a good approximation, but it will not affect our discussion of the electron–hole pair density relaxation—see Sect. 5.

The set of simplified hydrodynamic equations considered in this article reads

$$\begin{aligned} \frac{\partial n}{\partial t}&=-\left( R_{\mathrm {aug}}+R_{\mathrm {ste}}\right) \;, \nonumber \\ \frac{\partial n_{\mathrm {ste}}}{\partial t}&=R_{\mathrm {ste}}\;,\nonumber \\ 3k_{B}\frac{\partial \left( nT_{n}\right) }{\partial t}&=-3k_{B}n\left( \frac{T_{n}-T_{L}}{\tau }\right) +R_{\mathrm {aug}}E_{g}+Q_{\mathrm {bgr}}\left( n\right) \;,\nonumber \\ \rho _{L}c_{L}\frac{\partial T_{L}}{\partial t}&=3k_{B}n\left( \frac{T_{n}-T_{L}}{\tau }\right) +R_{\mathrm {ste}}\left( E_{g}-E_{\mathrm {ste}}\right) \;, \end{aligned}$$
(5)

where \(\tau =\left( \tau _{n}+\tau _{p}\right) /2\), \(R_{\mathrm {aug}}=R_{\mathrm {aug}}^{\left( n\right) }+R_{\mathrm {aug}}^{\left( p\right) }\) and \(Q_{\mathrm {bgr}}\left( n\right) \) is a correcting term due to band gap renormalisation, see Sect. 3.2. In the second line of Eq. 5, we have also neglected the STE decay as it occurs on a timescale typically longer than that considered in this article.

3 Excitonic processes

In this section, we use an analogy with chemical reactions, in which the ‘chemical’ species are electrons, holes, free excitons and STEs. According to this analogy, we can write down the following set of reactions

$$\begin{aligned} \begin{aligned}e+h&\longleftrightarrow x\\ x&\longrightarrow \mathrm {ste}\\ \mathrm {ste}&\longrightarrow \gamma \quad \mathrm {(radiative)\;,}\\ \mathrm {ste}&\longrightarrow \mathrm {FP}\quad \mathrm {(non-radiative)\;,} \end{aligned} \end{aligned}$$
(6)

where \(\gamma \) stands for photon and ‘FP’ for Frenkel pair of point defects. The last two reactions give the two STE decay pathways [10]. We initially assume that the STEs cannot ‘untrap’ and postpone the discussion of the thermal detrapping of STEs to the end of Sect. 5. For the sake of simplicity, we also neglect the impact dissociation of STEs.

The corresponding system of kinetic equations is

$$\begin{aligned} \begin{aligned}\frac{dn}{dt}&=-k_{1}n^{2}+k_{-1}n_{x}\;,\\ \frac{dn_{x}}{dt}&=k_{1}n^{2}-k_{-1}n_{x}-k_{2}n_{x}\;,\\ \frac{dn_{\mathrm {ste}}}{dt}&=k_{2}n_{x}-\left( \frac{1}{\tau _{\mathrm {ste},r}}+\frac{1}{\tau _{\mathrm {ste},nr}}\right) n_{\mathrm {ste}}\;, \end{aligned} \nonumber \\ \end{aligned}$$
(7)

where \(k_{1}\) and \(k_{-1}\) are the reaction rate constants of the forward and backward reactions in the first line of Eq. (6), \(k_{2}\) is the reaction rate constant in the second line of Eq. (6) and \(1/\tau _{\mathrm {ste},r}\) and \(1/\tau _{\mathrm {ste},nr}\) are the reaction rate constants in the third and fourth lines of Eq. (6). These last two rates are also found in the fourth line of Eq. (1). Neglecting the radiative and non-radiative decay of STEs, we have that \(n+n_{x}+n_{\mathrm {ste}}\) is conserved. Rate constants can depend on both electronic and lattice temperatures, see Sec. 3.2

At very large electronic temperatures, the equilibrium of the first reaction is shifted towards free electron–hole pair and there will be very few free excitons. This condition is generally fulfilled immediately after irradiation since the electronic temperature is large [6]. Exciton trapping can be a very fast (possibly barrierless [37]) process, e.g. in a-SiO\(\mathrm {_{2}}\) it occurs in about 150 fs [12,13,14]. In this case, we can use a quasi-steady-state approximation (QSSA) for the free excitons and assume that

$$\begin{aligned} \frac{dn_{x}}{dt}=k_{1}n^{2}-k_{-1}n_{x}-k_{2}n_{x}\approx 0\;. \end{aligned}$$
(8)

This assumption implies that free excitons are short-lived species—or reactive intermediates—and that their density is always negligible.

As a consequence of Eq. (8), we can substitute

$$\begin{aligned} n_{x}\approx \frac{k_{1}n^{2}}{\left( k_{-1}+k_{2}\right) } \end{aligned}$$
(9)

into Eq. (7) to obtain

$$\begin{aligned} \begin{aligned}\frac{dn}{dt}&\approx -\frac{k_{1}k_{2}n^{2}}{\left( k_{-1}+k_{2}\right) }\;,\\ \frac{dn_{\mathrm {ste}}}{dt}&\approx \frac{k_{1}k_{2}n^{2}}{\left( k_{-1}+k_{2}\right) }-\left( \frac{1}{\tau _{\mathrm {ste},r}}+\frac{1}{\tau _{\mathrm {ste},nr}}\right) n_{\mathrm {ste}}\;. \end{aligned} \end{aligned}$$
(10)

Within the domain of applicability of the QSSA, we can now distinguish between two extreme regimes:

  1. 1.

    When \(k_{2}\ll k_{-1}\) the electron–hole pairs and the free excitons reach their chemical equilibrium, given by \(n_{x}=\left( k_{1}/k_{-1}\right) n^{2}\), before the free excitons decay into STEs. In this case, \(k_{2}\) at the denominator of the right-hand side of the first line of Eq. (10) can be neglected, leading to

    $$\begin{aligned} \frac{dn}{dt}\approx -k_{2}\left( \frac{k_{1}}{k_{-1}}\right) n^{2}=-k_{1}^{\mathrm {eff}}n^{2}\;. \end{aligned}$$
    (11)

    Note that by using the above mentioned equilibrium condition, we have that \(k_{1}^{\mathrm {eff}}n^{2}=k_{2}n_{x}\). This condition is fulfilled close to the EMT because \(k_{-1}\) gets large—see Eq. (22). Hence, close to the EMT the free electron–hole pairs can become long-lived, even if \(1/k_{2}\approx 150\) fs [12,13,14], since \(n_{x}\approx 0\) — see discussion after Eq. (17). Note that this ‘excitonic bottleneck’ is effective only if the rate of Auger recombination is smaller than \(k_{1}^{\mathrm {eff}}n^{2}\). Above the EMT, the electron–hole pairs are still unstable even if \(k_{1}^{\mathrm {eff}}\approx 0\) and will eventually decay through Auger recombination [20].

  2. 2.

    When \(k_{-1}\ll k_{2}\), the free excitons decay into self-trapped excitons before reaching an equilibrium with the electron–hole pair and

    $$\begin{aligned} \frac{dn}{dt}\approx -k_{1}n^{2}\;. \end{aligned}$$
    (12)

    This regime can be relevant far from the EMT and will be briefly considered at the end of Sect. 5.

3.1 Exciton dissociation equilibrium

The ‘chemical’ equilibrium condition for the reaction \(e+h\longleftrightarrow x\) is reached when \(\mu _{n}+\mu _{p}=\mu _{x}\) [38], where \(\mu _{n}\), \(\mu _{p}\) and \(\mu _{x}\) are the chemical potentials of electrons, holes and excitons, respectively. In the following, we will consider the general case in which \(\mu _{n}+\mu _{p}\ne 0\), while \(n=p\). This condition corresponds to a locally neutral mixture of electrons and holes which has not yet fully equilibrated, being the ‘chemical’ equilibrium between electrons and holes characterised by the equation \(\mu _{n}+\mu _{p}=0\).

If the electrons and holes are not degenerate, Boltzmann statistics applies. For a simple spin-unpolarised parabolic two-band model of the insulator, we can then write that

$$\begin{aligned} n=\frac{2}{\Lambda _{n}^{3}}\exp \left( -\frac{E_{g}-\mu _{n}}{k_{B}T_{n}}\right) \;, \end{aligned}$$
(13)

and

$$\begin{aligned} p=\frac{2}{\Lambda _{p}^{3}}\exp \left( \frac{\mu _{p}}{k_{B}T_{n}}\right) \;, \end{aligned}$$
(14)

where \(\Lambda _{n}=\sqrt{2\pi \hbar ^{2}/m_{n}k_{B}T_{n}}\) and \(\Lambda _{p}=\sqrt{2\pi \hbar ^{2}/m_{p}k_{B}T_{n}}\) are the thermal wavelengths of electrons and holes, respectively. The effective masses of the electrons and holes are indicated by \(m_{n}\) and \(m_{p}\), respectively. The top of the valence band has been taken as the reference energy.

For fully equilibrated electrons and holes, the law of mass action \(np=n_{i}^{2}\), where \(n_{i}=2e^{-E_{g}/2k_{B}T_{n}}/\sqrt{\Lambda _{n}^{3}\Lambda _{p}^{3}}\), immediately follows from the ‘chemical’ equilibrium condition. The intrinsic electron–hole pair density \(n_{i}\), is negligible for wide band gap insulators at room temperature. More in general we have that

$$\begin{aligned} np=n^{2}= & {} \frac{4}{\Lambda _{n}^{3}\Lambda _{p}^{3}}\exp \left( -\frac{E_{g}-\mu _{n}-\mu _{p}}{k_{B}T_{n}}\right) \nonumber \\= & {} \frac{4}{\Lambda _{n}^{3}\Lambda _{p}^{3}}\exp \left( -\frac{E_{g}-\mu _{x}}{k_{B}T_{n}}\right) \;, \end{aligned}$$
(15)

where we have used the local neutrality condition, \(p=n\). This is the condition we assume to apply immediately (\(\sim \) 10 fs) after irradiation, when electrons and holes have thermally equilibrated among themselves through carrier–carrier scattering, but they still have to recombine.

By assuming non-degenerate free excitons at the same temperature, \(T_n\), of electrons and holes, we can write that

$$\begin{aligned} n_{x}=\frac{4}{\Lambda _{x}^{3}}\exp \left( -\frac{E_{x}-\mu _{x}}{k_{B}T_{n}}\right) \;, \end{aligned}$$
(16)

where \(\Lambda _{x}=\sqrt{2\pi \hbar ^{2}/m_{x}k_{B}T_{n}}\) is the thermal wavelength of the excitons and \(E_{x}\) is the energy of the free excitons at rest. The bare mass of the free exciton is \(m_{x}=m_{n}+m_{p}\), and we assume that the difference between the bare and effective exciton mass—defined as the inverse curvature of an excitonic band close to its minimum—is negligible . The analogue of the Saha’s equation for semiconductors reads [39, 40]

$$\begin{aligned} \frac{n_{x}}{n^{2}}= & {} \left( \frac{\Lambda _{n}\Lambda _{p}}{\Lambda _{x}}\right) ^{3}\exp \left( \frac{E_{g}-E_{x}}{k_{B}T_{n}}\right) \nonumber \\= & {} \left( \frac{2\pi \hbar ^{2}}{m_{x}^{\star }k_{B}T_{n}}\right) ^{\frac{3}{2}}\exp \left( \frac{\epsilon _{b}}{k_{B}T_{n}}\right) \end{aligned}$$
(17)

where \(m_{x}^{\star }=\left( 1/m_{n}+1/m_{p}\right) ^{-1}\) is the reduced mass of the exciton. Equation (17) is obtained by dividing Eq. (16) by Eq. (15) and setting \(E_{x}=E_{g}-\epsilon _{b}\), where \(\epsilon _{b}\) is the binding energy of the exciton.

If \(\epsilon _{b}>0\), the equilibrium is shifted towards bound electron–hole pairs at low electronic temperatures (\(T_{n}\ll \epsilon _{b}/k_{B}\)) and towards free electron–hole pairs at high electronic temperatures (\(T_{n}\,\gg \,\epsilon _{b}/k_{B}\)). In fact, the exciton binding energy is not fixed and depends on the density of electron–hole pairs and free excitons, see next section. In particular, the exciton binding energy can become negative for large electron–hole pair densities, leading to an instability of the free excitons even at low temperature.

3.2 Band gap renormalisation and exciton binding energy

It is experimentally observed that both the band gap and the exciton binding energy depend on the electron–hole pair density [38]. This is due to the extra Coulomb screening which is provided by the excited carriers which reduces the effective attraction between electrons and holes (free or bound). As a consequence, both the band gap and the exciton binding energy tend to decrease as the electron–hole pair density increases [41, 42] . The two effects compensate, and it is found experimentally that the exciton absorption frequency does not depend on the electron–hole pair density [38]. Both the band gap renormalisation (BGR) and exciton binding energy reduction can be estimated by using the electron–hole exchange and correlation functional proposed by Vashishta and Kalia (VK) [43]

$$\begin{aligned} \varepsilon _{xc}\left( r_{s}\right) =\epsilon _{b}^{0}\left( \frac{a+br_{s}}{c+dr_{s}+r_{s}^{2}}\right) \;, \end{aligned}$$
(18)

where the parameters in Eq. (18) are \(a=-4.8316\), \(b=-5.0879\), \(c=0.0152\) and \(d=3.0426\). The bare (Wannier–Mott) excitonic binding energy is defined as \(\epsilon _{b}^{0}=E_{0}m_{x}^{\star }/2m_{e}\epsilon _{r}^{2}\), where \(E_{0}=27.2\) eV is the Hartree energy, \(m_e\) is the bare electron mass and \(\epsilon _{r}=3.96\) is the relative static dielectric permittivity of a-\({{\mathrm {SiO}}_{2}}\) [44]. The dimensionless excitonic Wigner–Seitz radius

$$\begin{aligned} r_{s}=\left( \frac{4\pi }{3}n\right) ^{-\frac{1}{3}}/a_{x} \end{aligned}$$
(19)

is the ratio between the electronic Wigner–Seitz radius and the excitonic radius defined as \(a_{x}=a_{0}\epsilon _{r}/m_{x}^{\star }\), where \(a_{0}\) is the Bohr radius [45].

The VK exchange and correlation functional has been obtained by fitting calculations for Ge and Si performed within the fully self-consistent approximation of Vashishta and Singwi [43]. It is remarkable that the fit is almost independent of the detailed band structure and just depends on \(r_s\).

In the case of a locally neutral electron–hole fluid, the contribution from classical electrostatics is zero and the full electron–hole functional reads

$$\begin{aligned} \varepsilon _{eh}\left( r_{s}\right) =\frac{3\epsilon _{b}^{0}}{5}\left( \frac{9\pi }{4}\right) ^{\frac{2}{3}}\frac{1}{r_{s}^{2}}+\varepsilon _{xc}\left( r_{s}\right) \;, \end{aligned}$$
(20)

where the first term accounts for the kinetic energy of electrons and holes [46]. The BGR is obtained by computing the variation of chemical potential, \(\Delta E_{BGR}\left( r_{s}\right) =\mu _{eh}\left( r_{s}\right) -\mu _{eh}\left( r_{s}=0\right) \), where [46]

$$\begin{aligned} \mu _{eh}\left( r_{s}\right) =\varepsilon _{eh}\left( r_{s}\right) +n\frac{d\varepsilon _{eh}}{dn}=\varepsilon _{eh}\left( r_{s}\right) -\frac{r_{s}}{3}\frac{d\varepsilon _{eh}}{dr_{s}}\;. \end{aligned}$$
(21)

In fact, one can easily verify that only the exchange and correlation part contributes to the BGR and \(\mu _{xc}\left( r_{s}\right) =\varepsilon _{xc}\left( r_{s}\right) -\left( r_{s}/3\right) d\varepsilon _{eh}/dr_{s}\) [23, 47]. The renormalised excitonic binding energy also depends on \(r_{s}\) and is defined as \(\epsilon _{b}\left( r_{s}\right) =\epsilon _{b}^{0}+\Delta E_{BGR}\left( r_{s}\right) \). The value of the dimensionless excitonic Wigner–Seitz radius at which \(\epsilon _{b}\left( r_{s}\right) =0\) defines the electron–hole pair density at the EMT. For larger values of the electron–hole pair density, \(\epsilon _{b}\) gets negative and the free excitons become unstable [42]. Because of the BGR, an extra term appears in Eq. (5): \(Q_{\mathrm {bgr}}\left( n\right) =-\left( n\partial \Delta E_{BGR}/\partial n\right) \partial n/\partial t\).

From the first two lines of Eq. (7) in the limit of \(k_{2}=0\), we obtain that \(k_{-1}/k_{1}=n^{2}/n_{x}\). Following Sekiguchi and Shimano [42], we assume that \(k_{1}\) only depends on the electronic temperature, \(T_{n}\), and the lattice temperature, \(T_{L}\). The assumption is consistent with the so-called columnar model of recombination [48] which uses the Langevin’s estimate \(k_{1}\approx q\left( M_{n}+M_{p}\right) /\epsilon _{0}\epsilon _{r}\), where \(M_{n}\) and \(M_{p}\) are the electron and hole mobilities, respectively. The mobilities can depend on both \(T_{n}\) and \(T_{L}\) [29]. If the lattice temperature is close to room temperature, in the case of a-\(\mathrm {SiO_2}\) we also have that \(M_{p}\ll M_{n}\) because the holes tend to localised and their transport is activated [49]. As a consequence of the Saha dissociation equilibrium, Eq. (17), the reaction rate constant for the free exciton dissociation can be written as

$$\begin{aligned} k_{-1}=k_{1}\left( \frac{2\pi \hbar ^{2}}{m_{x}^{\star }k_{B}T_{n}}\right) ^{-\frac{3}{2}}\exp \left( -\frac{\epsilon _{b}}{k_{B}T_{n}}\right) \;. \end{aligned}$$
(22)

For electron–hole pair densities above the EMT, i.e. when \(\epsilon _{b}<0\), the right-hand side of Eq. (22) becomes very large in both the limits of \(T_{n}\ll \left| \epsilon _{b}\right| /k_{B}\)—because of the large exponential factor— and of \(T_{n}\,\gg \,\left| \epsilon _{b}\right| /k_{B}\). —because of the large prefactor, while the exponential factor gets close to one. The applicability of the QSSA approximation is then confirmed in these two limits, which generally include the condition of very large electronic temperature immediately after irradiation and the thermal equilibrium between the lattice and the electron–hole pairs, if \(\left| \epsilon _{b}\right| \,\gg \, k_{B}T_{L}\).

4 Results

We have solved numerically the simplified hydrodynamic equations (5) to assess the consequences of Coulomb screening and BGR in the evolution of the electron–hole plasma generated upon proton irradiation of a-\({{\mathrm {SiO}}_{2}}\). The band gap energy (bare, i.e. not renormalised) is set to \(E_{g}=8.7\) eV and the STE energy to \(E_{\mathrm {ste}}=E_{g}-5.6\) eV [16]. The values of the effective mass of the electron (0.5 \(m_e\)), light (0.5 \(m_e\)) and heavy (5 \(m_e\)) holes are taken from Ref. [50]. The static relative permittivity is set to \(\epsilon _{r}=3.96\). The corresponding Wannier–Mott exciton energy is \(E_{x}=E_{g}-0.367\) eV, and the excitonic Bohr radius is \(a_{x}=4.95\) Å. The two-temperature relaxation time is set to a typical value of \(\tau =0.1\) ps [5, 6], while STE formation time is set to \(\tau _{2}^{\mathrm {ste}}=1/k_{2}=150\) fs [12,13,14]. The lattice density and specific heat capacity are set to \(\rho _{L}=2.2\) \(\mathrm {g/cm^{3}}\) and \(c_{L}=0.67\cdot 10^{7}\) \(\mathrm {erg/g\cdot K}\), respectively.

By assuming that the average energy needed to generate an electron–hole pair is approximately three times the band gap—the so-called Klein’s rule [6, 51]—the initial electronic temperature is estimated as \(T_{n}\left( 0\right) =2E_{g}/3k_{B}=67,300\) K. Note that the initial electronic temperature is independent of the initial electron–hole pair density. The initial lattice temperature is \(T_{L}\left( 0\right) =298.15\) K or room temperature.

We consider several values of the initial electron–hole pair density, \(n\left( 0\right) \), below and across the EMT—the Mott density is \(n_{\mathrm {Mott}}=1.94\cdot 10^{19}\) \({\mathrm {cm}}^{\mathrm {-3}}\)—but always smaller than the critical density, \(n_{c}{=}4.25\cdot 10^{20}\) \({\mathrm {cm}}^{\mathrm {-3}}\).Footnote 1 The initial STE density is set to zero. Auger recombination is phenomenologically modelled as \(R_{\mathrm {aug}}=\left( C_{n}+C_{p}\right) n^{3}\), with the Auger coefficients notionally set to \(C_{p}=C_{n}=10^{29}\) \(\mathrm {cm^{6}/s}\) [20].

Fig. 1
figure 1

Free electron–hole pair density, \(n\left( t\right) \), normalised to the initial density, \(n\left( 0\right) \), for a selection of initial densities, \(n\left( 0\right) \). The initial densities \(n\left( 0\right) =10^{16},10^{17},10^{18}\) \({\mathrm {cm}}^{\mathrm {-3}}\) are below, \(10^{19}\) \({\mathrm {cm}}^{\mathrm {-3}}\) slightly below and \(10^{20}\) \({\mathrm {cm}}^{\mathrm {-3}}\) above the EMT, respectively

Full size image

Figure 1 shows the density of free electron–hole pairs, \(n\left( t\right) \), normalised to the initial density, \(n\left( 0\right) \). In all cases, the free electron–hole density decays abruptly close to 1 ps. We can rationalise the somehow counter-intuitive finding that an initially larger excitation yields a slightly shorter relaxation by looking at the effective STE relaxation times shown in Fig. 2.

Fig. 2
figure 2

Effective STE formation time, \(\tau _{1}^{\mathrm {ste}}\left( t\right) \), as a function of the electronic temperature, \(T_{n}\), for a selection of initial densities, \(n\left( 0\right) \). The initial densities \(n\left( 0\right) =10^{16},10^{17},10^{18}\) \({\mathrm {cm}}^{\mathrm {-3}}\) are below, \(10^{19}\) \({\mathrm {cm}}^{\mathrm {-3}}\) slightly below and \(10^{20}\) \({\mathrm {cm}}^{\mathrm {-3}}\) above the EMT, respectively

Full size image

This time is defined as

$$\begin{aligned} \tau _{1}^{\mathrm {ste}}=1/k_{1}^{\mathrm {eff}}n=\left( k_{-1}/k_{1}\right) /k_{2}n\;, \end{aligned}$$
(23)

where we have used Eq. (11) in the last equality, and shows a strong dependence on both the electronic temperature , \(T_{n}\), and the initial electron–hole pair density, \(n\left( 0\right) \). Note the qualitative difference between the case \(n\left( 0\right) =10^{20}\) \({\mathrm {cm}}^{\mathrm {-3}}\) and the remaining cases: if \(n\left( 0\right) >n_{\mathrm {Mott}}\), \(\tau _{1}^{\mathrm {ste}}\left( T_{n}\right) \) is bounded from below, i.e. there is a minimum value of the STE relaxation time. This is expected from the definition of \(\tau _{1}^{\mathrm {ste}}\) and the behaviour of \(k_{-1}/k_{1}\) according to Eq. (22). In particular, the right-hand side of Eq. (22) becomes large in both limits of \(T_{n}\rightarrow 0\) and \(T_{n}\rightarrow \infty \) if \(\epsilon _b<{0}\), i.e. above the EMT. On the other hand, if \(n\left( 0\right) <n_{\mathrm {Mott}}\), \(\tau _{1}^{\mathrm {ste}}\left( T_{n}\right) \) is a monotonically increasing function of \(T_{n}\). In all cases, for large electronic temperature we have that the electron–hole pairs are still fully dissociates, i.e. \(n\approx n\left( 0\right) \), while the exponential on the right-hand side of Eq. (22) is close to one. As a consequence, we also have that \(k_{-1}/k_{1}\propto T_n^{3/2}\) and, from Eq. (23), that \(\tau _{1}^{\mathrm {ste}}\propto T_{n}^{3/2}/n\left( 0\right) \) , which explains why the relaxation shown in Fig. 1 is slightly faster for larger initial densities.

Coulomb screening and BGR are irrelevant during the initial part of the relaxation because \(k_{B}T_{n}\) largely exceeds the exciton binding energy, \(\epsilon _{b}\). At such high electronic temperature, all free excitons would be dissociated also in the absence of the BGR. We also found that Auger recombination gives only a minor contribution (about 6% of the electron–hole pair recombination) for \(n\left( 0\right) =10^{20}\) \({\mathrm {cm}}^{\mathrm {-3}}\) and a negligible contribution in the other cases.

Fig. 3
figure 3

Free electron–hole pair density, \(n\left( t\right) \), normalised by the initial density, \(n\left( 0\right) \), as a function of the electronic temperature, \(T_{n}\), for a selection of the initial density, \(n\left( 0\right) \). The initial densities \(n\left( 0\right) =10^{16},10^{17},10^{18}\) \({\mathrm {cm}}^{\mathrm {-3}}\) are below, \(10^{19}\) \({\mathrm {cm}}^{\mathrm {-3}}\) slightly below and \(10^{20}\) \({\mathrm {cm}}^{\mathrm {-3}}\) above the EMT, respectively

Full size image

To further support the correlation between electron–hole trapping and thermal relaxation in a-\({{\mathrm {SiO}}_{2}}\), in Fig. 3 we plot the free electron–hole pair density, n, as a function of the electronic temperature, \(T_{n}\). For \(n\left( 0\right) <n_{\mathrm {Mott}}\) a threshold effect is evident as the electron–hole pairs make an abrupt transition from free to self-trapped at a sharp value of \(T_{n}\). This is due to the dependence of \(\tau _{1}^{\mathrm {ste}}\left( T_{n}\right) \) as shown in Fig. 2 which in turn depends on the Arrhenius form of the Saha equilibrium constant defined in Eq. (17). For \(n\left( 0\right) <n_{\mathrm {Mott}}\), \(\tau _{1}^{\mathrm {ste}}\) goes very steeply to zero as the electronic temperature decreases. The transition becomes less abrupt as the Mott density is approached. Note that in the case of \(n\left( 0\right) =10^{20}\) \({\mathrm {cm}}^{\mathrm {-3}}\) the minimum of \(\tau _{1}^{\mathrm {ste}}\left( T_{n}\right) \) in Fig. 2 occurs at \(T_{n}\approx 10^{3}\) and that the minimal value is slightly larger than both the two-temperature equilibration time, \(\tau \), and the STE formation time, \(\tau _{2}^{\mathrm {ste}}\).

Fig. 4
figure 4

Electronic temperature, \(T_{n}\left( t\right) \), for a selection of the initial density, \(n\left( 0\right) \). The initial densities \(n\left( 0\right) =10^{16},10^{17},10^{18}\) \({\mathrm {cm}}^{\mathrm {-3}}\) are below, \(10^{19}\) \({\mathrm {cm}}^{\mathrm {-3}}\) slightly below and \(10^{20}\) \({\mathrm {cm}}^{\mathrm {-3}}\) above the EMT, respectively

Full size image

Figure 4 shows the electronic temperature, \(T_{n}\), as a function of time. The initial part of the relaxation shows a negligible dependence on the initial electron–hole density, and it is determined by the two-temperature equilibration time, \(\tau \), only. Eventually, the electronic and lattice temperatures become equal, i.e. the lattice and electronic degrees of freedom thermally equilibrate. However, without a diffusive process in Eq. (5) which can equilibrate the local lattice temperature to that of the environment, the difference between the final and initial lattice temperatures shows an artificial proportionality to \(n\left( 0\right) \).

5 Discussion

In this section, we discuss the validity and implications of the main assumptions of our simplified model. Neglecting the spatial dependency of the fields and the associated diffusive processes is obviously a severe approximation which was explicitly made in order to single out the effects on the time evolution due to Coulomb screening, only. From a qualitative point of view, re-introducing the diffusive terms will lead to a further decrease of the free electron–hole pair density in the core region of the track. Excitons in principle are more mobile than charge carriers, but their ultrafast trapping will strongly limit their diffusion. This is an efficient mechanism to localise the energy initially transferred to the electronic degrees of freedom and generate lattice defects (Frenkel pairs) in a-\({{\mathrm {SiO}}_{2}}\) [1, 8]. With the qualitative effect of carrier diffusion in mind, the results of the simplified hydrodynamic model for the free electron–hole pair density can be safely taken as an upper limit, i.e. the free electron–hole pair relaxation can be even faster than what shown in Fig. 1, but not slower. As a consequence, re-introducing the diffusive terms and solving the full hydrodynamic model, Eq. (1), are not expected to demonstrate an increased free electron–hole pair relaxation time driven by Coulomb screening and BGR, alone.

Assuming a constant STE formation time, \(\tau _{2}^{\mathrm {ste}}\), is also an approximation that has been recently criticised by Jürgens et al. [16]. These authors suggest that the long free electron–hole pair relaxation time observed by increasing the laser intensity can be explained by thermal detrapping of STEs. However, the effect of detrapping is minor in a-\({{\mathrm {SiO}}_{2}}\) because of the large activation energy of about 1.5 eV [52]. This behaviour must be contrasted against the case of sapphire (\({\mathrm {Al}}_{\mathrm {2}}{\mathrm {O}}_{\mathrm {3}}\)) for which a much longer (about 100 ps) free electron–hole pair relaxation time has been reported [53], along with a much smaller activation energy for detrapping, between 0.04 an 0.51 eV [54]. The case of borosilicate crown glass BK7 is also relevant because long (exceeding 10 ps) free electron–hole pair relaxation times—longer than for a-\({{\mathrm {SiO}}_{2}}\)—have been reported both using optical [15] and proton [24, 26,27,28] probes. Borosilicate glasses present more non-bridging oxygen atoms than a-\({{\mathrm {SiO}}_{2}}\) because boron is a trivalent cation (\({\mathrm {B}}^{{3+}}\)) which disrupts the silica tetrahedral network and can create localised soft modes [55]. Although STE formation has been observed in BK7, it is possible that the presence of a large number of non-bridging oxygen atoms both hinder the STE formation and facilitate their thermal detrapping [56].

The addition of a STE detrapping mechanism into Eq. (7) yields the following set of kinetic equations:

$$\begin{aligned} \begin{aligned}\frac{dn}{dt}&=-k_{1}n^{2}+k_{-1}n_{x}\;,\\ \frac{dn_{x}}{dt}&=k_{1}n^{2}-k_{-1}n_{x}-k_{2}n_{x}+k_{-2}n_{\mathrm {ste}}\;,\\ \frac{dn_{\mathrm {ste}}}{dt}&=k_{2}n_{x}-k_{-2}n_{\mathrm {ste}}-\left( \frac{1}{\tau _{\mathrm {ste},r}}+\frac{1}{\tau _{\mathrm {ste},nr}}\right) n_{\mathrm {ste}}\;, \end{aligned}\nonumber \\ \end{aligned}$$
(24)

where \(k_{-2}\left( T_{L}\right) =k_{-2}^{0}e^{-E_{a}/k_{B}T_{L}}\) is the Arrhenius form of the detrapping rate, where we indicate with \(k_{-2}^{0}\) the frequency factor—of the order of the Debye frequency of the material—and with \(E_{a}\) the activation energy.

By applying the QSSA in the \(k_{-1}\ll k_{2}\) case, we then obtain that

$$\begin{aligned} \frac{dn}{dt}\approx -k_{1}^{\mathrm {eff}}n^{2}+k_{-2}n_{\mathrm {ste}}\;. \end{aligned}$$
(25)

The steady-state solution of Eq. (25), \(n_{\mathrm {ste}}/n^{2}\approx k_{1}^{\mathrm {eff}}/k_{-2}=\left( k_{1}^{\mathrm {eff}}/k_{-2}^{0}\right) e^{E_{a}/k_{B}T_{L}}\), suggests a saturation of STEs if \(T_{L}\,\gg \, E_{a}/k_{B}\) and a persistent population of free electron–hole pairs if the initial density is larger than the Mott density. For initial densities smaller than the Mott density, there is no minimal value of \(\tau _{1}^{\mathrm {ste}}\)—see Fig. 2—and the Saha dissociation equilibrium will eventually favour the formation of free excitons as the electronic and lattice degrees of freedom approach thermal equilibrium. Hence, Coulomb screening and BGR, along with an efficient thermal detrapping of STEs, can provide an explanation for the observed long free electron–hole pair relaxation times in BK7. Note that the QSSA is not expected to hold far below the EMT as the STE formation rate, \(k_{2}\), can become much larger than the exciton dissociation rate, \(k_{-1}\)—see Eq. (22). As a consequence, thermal detrapping of STEs alone cannot yield a longer free electron–hole pair relaxation time below the EMT.

6 Conclusions

Free electron–hole pairs excited in wide gap insulators relax into STEs, if their formation is possible. This is the case for a-\({{\mathrm {SiO}}_{2}}\) and borosilicate crown glass BK7. The STE formation occurs in about 150 fs in a-\({\mathrm {SiO}}_{2}\) if the initial density of free electron–hole pairs is smaller than \(10^{20}\) \({\mathrm {cm}}^{\mathrm {-3}}\). This ultrafast mechanism has been previously observed with pump–probe experiments and recently confirmed by optically probed samples irradiated with ultrashort proton pulses generated by target normal sheath acceleration. At variance with a-\({{\mathrm {SiO}}_{2}}\), BK7 samples display a longer free electron–hole pair relaxation time of about hundreds of ps after irradiation by either protons or high-intensity laser pulses. Two microscopic mechanisms have been suggested in the literature to explain this long relaxation time: (i) suppression of the formation of free excitons—precursor of the STEs—for densities above the EMT due to Coulomb screening and (ii) thermal detrapping of STEs if the STEs are only shallowly trapped. In this article, we have tested the first mechanism by solving a simplified set of hydrodynamic equations for the coupled evolution of the electron–hole and STE densities, along with the electronic and lattice temperatures. Our numerical results do not show a long (\({\gg }\,{\mathrm {10}}\) ps) free electron–hole pair relaxation time in a-\({{\mathrm {SiO}}_{2}}\) even for initial densities exceeding the Mott density. In fact, a slightly shorter relaxation time was observed in this case. This rather counter-intuitive result is rationalised by noting that, immediately after irradiation, the electronic temperature is large enough (\(k_{B}T_{n}\approx 2E_{g}/3\)) to suppress the formation of free excitons for any initial density. Coulomb screening is then irrelevant for the exciton formation in such a ‘hot’ electron–hole plasma. We then conclude that another mechanism, e.g. thermal detrapping of STEs, is required to correctly model the long free electron–hole pair relaxation time observed in some irradiated material, including BK7. The effect of detrapping is expected to be negligible in a-SiO2 because of the large activation energy, but in the case of BK7 the presence of a large number of non-bridging oxygen atoms can indeed facilitate the thermal detrapping of STEs. Our results are based on a simplified set of hydrodynamic equations in which we have neglected the spatial dependency of the fields and the associated diffusive processes. This is a severe approximation which was deliberately introduced to single out the effects due to Coulomb screening, only, on the time evolution. Despite this severe approximation, our results can be safely taken as an upper limit for the free electron–hole pair relaxation time and we do not expect the main conclusion of this article to be affected by the re-introduction of the diffusive terms. Modelling based on the full set of hydrodynamic equations is currently in progress.