An excitonic model for the electron–hole plasma relaxation in proton-irradiated insulators

The relaxation of free electron–hole pairs generated after proton irradiation is modelled by means of a simplified set of hydrodynamic equations. The model describes the coupled evolution of the electron–hole pair and self-trapped exciton (STE) densities, along with the electronic and lattice temperatures. The equilibration of the electronic and lattice excitations is based on the two-temperature model, while two mechanisms for the relaxation of free electron–hole pairs are considered: STE formation and Auger recombination. Coulomb screening and band gap renormalisation are also taken into account. Our numerical results show an ultrafast (≪1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\ll }\,{\mathrm {1}}$$\end{document} ps) free electron–hole pair relaxation time in amorphous SiO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathrm {SiO}}_{\mathrm {2}}}$$\end{document} for initial carrier densities either below or above the exciton Mott transition. Coulomb screening alone is not found to yield the long relaxation time (≫10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathrm {\gg }}{\mathrm {10}}$$\end{document} ps) experimentally observed in amorphous SiO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathrm {SiO}}_{\mathrm {2}}}$$\end{document} and borosilicate crown glass BK7 irradiated with high-intensity laser pulses or BK7 irradiated by short proton pulses. Another mechanism, e.g. thermal detrapping of STEs, is required to correctly model the long free electron–hole pair relaxation time observed experimentally.


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 (SiO 2 ) [11].
The formation of STEs has been experimentally observed in both quartz (α-SiO 2 ) and amorphous silica (a-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 a e-mail: l.stella@qub.ac.uk (corresponding author) 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-SiO 2 as the charge density exceeds 10 20 cm −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 cm −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 pumpprobe 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 (10 18 -10 19 W/cm 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 (10 19 cm −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 carrierlattice 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 electronhole 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.

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 where q is the unit charge, 0 r gives the static dielectric permittivity of the material, φ is the electrostatic potential, n, p and n ste represent the electron, hole and STE densities (unit: cm −3 ), T n , T p and T L are the electron, hole and lattice temperatures, J n , J p and J L are the charge current densities associated with the electrons, holes and lattice, S n , S p and S L are the energy current densities associated with the electrons, holes and lattice , ρ L gives the lattice density (unit: g/cm 3 ), c L is the specific heat capacity of the lattice (unit: erg/g · K), τ n and τ 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 rad , R srh , R (n,p) aug and R ste are rates (unit: cm −3 · s −1 ) for radiative recombination, defect-assisted (Shockley-Read-Hall) recombination, Auger recombination (for electrons or holes ) and STE formation. Finally, τ ste,r and τ ste,nr give the radiative and nonradiative STE decay times, respectively (see Sect. 3).
The constitutive equations for the currents are where k B is the Boltzmann constant, M n and M p are the phenomenological electron and hole mobilities, while κ n , κ p and κ L are the electron, hole and lattice heat conductivities. It is also assumed that where the coefficients r n and r p come from the powerlaw dependence of the most effective relaxation time as a function of the kinetic energy of the electrons and holes. The power is + 3 2 for ionised impurity scattering, − 1 2 for acoustic phonon scattering and + 1 2 for optical phonon scattering [31]. If radiative recombination is neglected, Eq. (1) conserves the integral of the total energy density where the lattice energy density is defined so that ∂uL ∂t = ρ L c L ∂TL ∂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 φ = 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 aug and Q bgr (n) 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.

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 where γ 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 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/τ ste,r and 1/τ 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 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 2 it occurs in about 150 fs [12][13][14]. In this case, we can use a quasi-steadystate approximation (QSSA) for the free excitons and assume that This assumption implies that free excitons are shortlived species-or reactive intermediates-and that their density is always negligible. As a consequence of Eq. (8), we can substitute into Eq. (7) to obtain Within the domain of applicability of the QSSA, we can now distinguish between two extreme regimes: 1. When k 2 k −1 the electron-hole pairs and the free excitons reach their chemical equilibrium, given by n x = (k 1 /k −1 ) 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 Note that by using the above mentioned equilibrium condition, we have that k eff 1 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 ≈ 150 fs [12][13][14], since n x ≈ 0 -see discussion after Eq. (17). Note that this 'excitonic bottleneck' is effective only if the rate of Auger recombination is smaller than k eff 1 n 2 . Above the EMT, the electron-hole pairs are still unstable even if k eff 1 ≈ 0 and will eventually decay through Auger recombination [20].
2. When k −1 k 2 , the free excitons decay into selftrapped excitons before reaching an equilibrium with the electron-hole pair and This regime can be relevant far from the EMT and will be briefly considered at the end of Sect. 5.

Exciton dissociation equilibrium
The 'chemical' equilibrium condition for the reaction e + h ←→ x is reached when μ n + μ p = μ x [38], where μ n , μ p and μ x are the chemical potentials of electrons, holes and excitons, respectively. In the following, we will consider the general case in which μ n + μ p = 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 μ n + μ 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 and where Λ n = 2π 2 /m n k B T n and Λ p = 2π 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 2 i , where n i = 2e −Eg/2kB Tn / Λ 3 n Λ 3 p , 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 where we have used the local neutrality condition, p = n. This is the condition we assume to apply immediately (∼ 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 where Λ x = 2π 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 minimumis negligible . The analogue of the Saha's equation for semiconductors reads [39,40] where m x = (1/m n + 1/m p ) −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 − b , where b is the binding energy of the exciton.
If b > 0, the equilibrium is shifted towards bound electron-hole pairs at low electronic temperatures (T n b /k B ) and towards free electron-hole pairs at high electronic temperatures (T n 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.

Band gap renormalisation and exciton binding energy
It is experimentally observed that both the band gap and the exciton binding energy depend on the electronhole 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 pro-posed by Vashishta and Kalia (VK) [43] ε xc (r s ) = 0 b a + br s c + dr s + r 2 s , (18) where the parameters in Eq. (18)  , where E 0 = 27.2 eV is the Hartree energy, m e is the bare electron mass and r = 3.96 is the relative static dielectric permittivity of a-SiO 2 [44]. The dimensionless excitonic Wigner-Seitz radius is the ratio between the electronic Wigner-Seitz radius and the excitonic radius defined as a x = a 0 r /m x , 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 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, ΔE BGR (r s ) = μ eh (r s ) − μ eh (r s = 0), where [46] μ eh (r s ) = ε eh (r s ) + n dε eh dn = ε eh (r s ) − r s 3 dε eh dr s .
In fact, one can easily verify that only the exchange and correlation part contributes to the BGR and μ xc (r s ) = ε xc (r s ) − (r s /3) dε eh /dr s [23,47]. The renormalised excitonic binding energy also depends on r s and is defined as b (r s ) = 0 b + ΔE BGR (r s ). The value of the dimensionless excitonic Wigner-Seitz radius at which b (r s ) = 0 defines the electron-hole pair density at the EMT. For larger values of the electron-hole pair density, b gets negative and the free excitons become unstable [42]. Because of the BGR, an extra term appears in Eq. (5): Q bgr (n) = − (n∂ΔE BGR /∂n) ∂n/∂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 socalled columnar model of recombination [48] which uses the Langevin's estimate k 1 ≈ q (M n + M p ) / 0 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-SiO 2 we also have that M p 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 For electron-hole pair densities above the EMT, i.e. when b < 0, the right-hand side of Eq. (22) becomes very large in both the limits of T n | b | /k B -because of the large exponential factor-and of T n | b | /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 | b | k B T L .

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-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 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 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 twotemperature relaxation time is set to a typical value of τ = 0.1 ps [5,6], while STE formation time is set to τ ste 2 = 1/k 2 = 150 fs [12][13][14]. The lattice density and specific heat capacity are set to ρ L = 2.2 g/cm 3 and c L = 0.67 · 10 7 erg/g · 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 (0) = 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 (0) = 298.15 K or room temperature.
We consider several values of the initial electronhole pair density, n (0), below and across the EMT-the Mott density is n Mott = 1.94 · 10 19 cm −3 -but always smaller than the critical density, n c =4.25 · 10 20 cm −3 . 1 1 The critical density, nc, is defined so that ω l = ωp (nc) = ncq 2 /m x 0, where ω l is the probe frequency and ωp is the Fig. 1 Free electron-hole pair density, n (t), normalised to the initial density, n (0), for a selection of initial densities, n (0). The initial densities n (0) = 10 16 , 10 17 , 10 18 cm −3 are below, 10 19 cm −3 slightly below and 10 20 cm −3 above the EMT, respectively The initial STE density is set to zero. Auger recombination is phenomenologically modelled as R aug = (C n + C p ) n 3 , with the Auger coefficients notionally set to C p = C n = 10 29 cm 6 /s [20]. Figure 1 shows the density of free electron-hole pairs, n (t), normalised to the initial density, n (0). In all cases, the free electron-hole density decays abruptly close to 1 ps. We can rationalise the somehow counterintuitive finding that an initially larger excitation yields a slightly shorter relaxation by looking at the effective STE relaxation times shown in Fig. 2. bulk plasma frequency of the electron-hole plasma. To compute nc, we set ω l = 1.18 eV corresponding to a 1053 nm probe. If a simple Drude model is used to approximate the dielectric permittivity of the electron-hole plasma created upon laser irradiation, the plasma is expected to become strongly absorbing at ω l above the critical density. Fig. 3 Free electron-hole pair density, n (t), normalised by the initial density, n (0), as a function of the electronic temperature, Tn, for a selection of the initial density, n (0). The initial densities n (0) = 10 16 , 10 17 , 10 18 cm −3 are below, 10 19 cm −3 slightly below and 10 20 cm −3 above the EMT, respectively This time is defined as 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 (0). Note the qualitative difference between the case n (0) = 10 20 cm −3 and the remaining cases: if n (0) > n Mott , τ ste 1 (T n ) is bounded from below, i.e. there is a minimum value of the STE relaxation time. This is expected from the definition of τ ste 1 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 → 0 and T n → ∞ if b < 0, i.e. above the EMT. On the other hand, if n (0) < n Mott , τ ste 1 (T n ) 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 ≈ n (0), 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 ∝ T 3/2 n and, from Eq. (23), that τ ste 1 ∝ T 3/2 n /n (0) , 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, 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 (0) = 10 20 cm −3 and a negligible contribution in the other cases.
To further support the correlation between electronhole trapping and thermal relaxation in a-SiO 2 , in Fig. 3 we plot the free electron-hole pair density, n, Fig. 4 Electronic temperature, Tn (t), for a selection of the initial density, n (0). The initial densities n (0) = 10 16 , 10 17 , 10 18 cm −3 are below, 10 19 cm −3 slightly below and 10 20 cm −3 above the EMT, respectively as a function of the electronic temperature, T n . For n (0) < n 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 τ ste 1 (T n ) as shown in Fig. 2 which in turn depends on the Arrhenius form of the Saha equilibrium constant defined in Eq. (17). For n (0) < n Mott , τ ste 1 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 (0) = 10 20 cm −3 the minimum of τ ste 1 (T n ) in Fig. 2 occurs at T n ≈ 10 3 and that the minimal value is slightly larger than both the two-temperature equilibration time, τ , and the STE formation time, τ ste 2 . 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, τ , 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 (0).

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-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, τ ste 2 , 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-SiO 2 because of the large activation energy of about 1.5 eV [52]. This behaviour must be contrasted against the case of sapphire (Al 2 O 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-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-SiO 2 because boron is a trivalent cation (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: where k −2 (T L ) = k 0 −2 e −Ea/kB TL is the Arrhenius form of the detrapping rate, where we indicate with k 0 −2 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 k 2 case, we then obtain that The steady-state solution of Eq. (25), n ste /n 2 ≈ k eff 1 /k −2 = k eff 1 /k 0 −2 e Ea/kB TL , suggests a saturation of STEs if T L 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 τ ste 1see 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.

Conclusions
Free electron-hole pairs excited in wide gap insulators relax into STEs, if their formation is possible. This is the case for a-SiO 2 and borosilicate crown glass BK7. The STE formation occurs in about 150 fs in a-SiO 2 if the initial density of free electron-hole pairs is smaller than 10 20 cm −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-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 ( 10 ps) free electronhole pair relaxation time in a-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 ≈ 2E g /3) to suppress the formation of free excitons for any initial den-sity. 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.