Delayed-rate equations model for femtosecond laser-induced breakdown in dielectrics

TL;DR

This paper introduces a delayed-rate equations model that accurately predicts femtosecond laser-induced breakdown in dielectrics by tracking plasma density and carrier energy, offering improved physical insight and computational efficiency.

Contribution

The paper develops a simplified yet accurate delayed-rate equations model that matches complex multiple rate equations with fewer parameters, enhancing understanding and simulation of laser breakdown.

Findings

Model reproduces predictions of multiple rate equations

Limited free parameters determined by experimental data

Enables large-scale 3D simulations of laser-induced breakdown

Abstract

Experimental and theoretical studies of laser-induced breakdown in dielectrics provide conflicting conclusions about the possibility to trigger ionization avalanche on the sub-picosecond time scale and the relative importance of carrier-impact ionization over field ionization. On the one hand, current models based on single ionization-rate equations do not account for the gradual heating of the charge carriers which, for short laser pulses, might not be sufficient to start an avalanche. On the other hand, models based on multiple rate equations that track the carriers kinetics rely on several free parameters, which limits the physical insight that we can gain from them. In this paper, we develop a model that overcomes these issues by tracking both the plasma density and carriers' mean kinetic energy as a function of time, forming a set of delayed rate equations that we use to match the laser-induced damage threshold of several dielectric materials. In particular, we show that this simplified model reproduces the predictions from the multiple rate equations, with a limited number of free parameters determined unambiguously by fitting experimental data. A side benefit of the delayed rate equations model is its computational efficiency, opening the possibility for large-scale, three-dimensional modelling of laser-induced breakdown of transparent media.

Tables

Table 1. Table 1: Dielectric material parameters associated with the DRE fits given in Fig.  4 . Typical values for the linear refractive index n 0 subscript 𝑛 0 n_{0} , the bandgap ℰ g subscript ℰ 𝑔 \mathcal{E}_{g} , the recombination rate γ r subscript 𝛾 𝑟 \gamma_{r} and the molecular density ρ mol subscript 𝜌 mol \rho_{\mathrm{mol}} are gathered from various references. To estimate the molecular cross-section σ mol subscript 𝜎 mol \sigma_{\mathrm{mol}} , we have summed the individual cross-sections of the constitutive atoms, calculated as the area of a circle with a radius equal to the covalent radius. The plasma damping rate γ 𝛾 \gamma and the effective mass of the electrons m e subscript 𝑚 𝑒 m_{e} (in units of the free electron mass m 0 subscript 𝑚 0 m_{0} ) are set by fitting experimental data for the pulse-length dependence of the fluence threshold (see Figs.  4 and  5 ).

	SiO2	Al2O3	HfO2	Ta2O5	TiO2
$n_0$	1.45	1.76	2.09	2.1	2.52
${ℰ}_g$ [eV]	9.0	6.5	5.1	3.8	3.3
${\rho }_{\mathrm{mol}}\left[{10}^{28}/m^3\right]$	2.20	2.35	2.77	1.12	3.19
${\sigma }_{\mathrm{mol}}\left[{10}^{-19}/m^2\right]$	0.661	1.33	1.24	2.50	1.08
${\gamma }_r$ [ps-1]	4.0	0.0	0.0	0.0	0.0
$\gamma$ [fs-1]	2.0	1.0	0.5	0.4	0.5
$m_e$	1.0	0.8	0.4	0.5	0.3
$m_h$	1.0	1.0	1.0	1.0	1.0

Equations

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}=\underset{\mathrm{FI}}{\underbrace{{\nu }_{\mathrm{fi}}{\rho }_n}}+\underset{\mathrm{II}}{\underbrace{\alpha \rho I{\rho }_n}}-\underset{\mathrm{RE}}{\underbrace{{\gamma }_r\rho }},$$

$${\mathcal{E}}_c=\left(1+\frac{m_r}{m_h}\right)({\mathcal{E}}_g+{\mathcal{E}}_p),$$

$$\frac{\mathrm{d}{\rho }_0}{\mathrm{d}t}$$

$$\text{leavevmode}\text{leavevmode}⋮$$

$$\frac{\mathrm{d}{\rho }_j}{\mathrm{d}t}$$

$$\text{leavevmode}\text{leavevmode}⋮$$

$$\frac{\mathrm{d}{\rho }_k}{\mathrm{d}t}$$

$$\rho$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$t_{\mathrm{MRE}}$$

$${\xi }^s$$

$$=\mathrm{erfc}(r_s)+\frac{2r_s}{\sqrt{\pi }}\exp (-r_s^2),$$

$$\frac{\mathrm{d}\rho }{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\mathcal{E}}_{\mathrm{k}}^s}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\rho }^e}{\mathrm{d}t}$$

$${\rho }^e(t)={\nu }_{\mathrm{fi}}{\rho }_{\mathrm{mol}}\frac{{\mathcal{E}}_c}{{\gamma }_{\mathrm{ib}}^e\hslash \omega }\left[\exp \left(\frac{{\gamma }_{\mathrm{ib}}^e\hslash \omega }{{\mathcal{E}}_c}t\right)-1\right].$$

$$t_{\mathrm{DRE}}=\frac{1}{{\gamma }_{\mathrm{ib}}^e}\frac{{\mathcal{E}}_c}{\hslash \omega }.$$

$${\gamma }_{\mathrm{ib}}^s=\frac{\gamma }{\hslash \omega }\frac{q^2{E}^2}{2m_s({\gamma }^2+{\omega }^2)}.$$

$${\gamma }_n^s={\sigma }_{\mathrm{mol}}{\rho }_n\sqrt{\frac{2{\mathcal{E}}_{\mathrm{k}}^s}{m_s}},$$

$$E_{\mathrm{vac}}(t)=E_0\exp \left[-2\ln (2){\left(\frac{t}{\tau }\right)}^2\right],$$

$$I_{\mathrm{vac}}(t)=c{ϵ}_0|E_{\mathrm{vac}}(t){|}^2=c{ϵ}_0{E}_0^2\exp \left[-4\ln (2){\left(\frac{t}{\tau }\right)}^2\right],$$

$$F={\int }_{-\infty }^{\infty }{I}_{\mathrm{vac}}(t)\mathrm{d}t=\frac{c{ϵ}_0{E}_0^2\tau }{2}\sqrt{\frac{\pi }{\ln (2)}},$$

$$E_{\mathrm{bulk}}^2(t)=E_{\mathrm{vac}}^2(t)\frac{1-R}{n_0},$$

$$R={\left|\frac{n-1}{n+1}\right|}^2$$

$$n^2=n_0^2-\frac{{\omega }_p^2}{{\omega }^2+i\omega \gamma }.$$

$${\mathcal{E}}_{\mathrm{k}}^s<-\frac{3}{2}{\mathcal{E}}_c{\left[\ln \left(\frac{{\gamma }_{\mathrm{ib}}^s\hslash \omega }{2{\mathcal{E}}_c{\sigma }_{\mathrm{mol}}{\rho }_{\mathrm{mol}}}\sqrt{\frac{m_s\pi }{3{\mathcal{E}}_c}}\right)\right]}^{-1}.$$

$${\gamma }_e^e=\frac{4\pi {ϵ}_0}{q^2}\sqrt{\frac{6}{m_e}}{\left(\frac{2{\mathcal{E}}_{\mathrm{k}}^e}{3}\right)}^{3/2}.$$

$${\rho }_c=\left(\frac{{ϵ}_0{m}_r}{q^2}\right)n_0^2\left({\omega }^2+{\gamma }^2\right).$$

$$\frac{\mathrm{d}\tilde{i}(t)}{\mathrm{d}t}=-\gamma \tilde{i}(t)+\frac{q^2}{m}\tilde{E}(t),$$

$$\tilde{i}(t)=\frac{q^{2}E}{m({\gamma }^2+{\omega }^2)}\left[\omega \sin (\omega t)+\gamma \cos (\omega t)\right].$$

$$\tilde{P}(t)$$

$$=\frac{q^2{E}^2}{m({\gamma }^2+{\omega }^2)}\left[\omega \sin (\omega t)\cos (\omega t)+\gamma {\cos }^2(\omega t)\right].$$

$${\tilde{\mathcal{E}}}_p(t)$$

Full text

I Introduction

Computer modelling of strong-field optical phenomena in dielectrics driven by intense laser radiation is essential to understand the fundamental processes in play, e.g., during laser micro-machining, laser surgery, and high-harmonic generation in solids, to name a few. Mechanisms for laser-induced breakdown were identified and studied in various contexts Gallais et al. (2015); Jing et al. (2012); Chimier et al. (2011); Christensen and Balling (2009); Jupé et al. (2009); Jia et al. (2006); Rethfeld (2004); Kaiser et al. (2000). In the accepted picture, plasma formation in laser-driven dielectrics proceeds as follows. (1) Charge carriers are first created by field ionization (FI). (2) The charge carriers absorb energy from the laser field via inverse bremsstrahlung heating (IBH). (3) The hot charge carriers create new, cold ones through carrier-impact ionization (II). (4) The carriers created by II, in turn, gain energy from the laser field and create new carriers by II, and so on. This multiplication of charge carriers via II leads to an exponential growth of the plasma density, often referred to as an ionization avalanche. This picture applies well when the FI-II interplay has enough time to unfold, e.g., when the pulse duration is in the picosecond range or above. However, current experimental and theoretical studies of laser-induced breakdown in dielectrics provide conflicting conclusions about the relative importance of II over FI and the possibility to trigger ionization avalanche on the sub-picosecond time scale Balling and Schou (2013).

For example, a pumb-probe experiment in fused silica Lebugle et al. (2014) has shown that a significant amount of ionization can take place after the pump pulse, which cannot be described by FI alone and suggests a delayed II avalanche triggered by slowly-decaying hot plasmon excitations. In contrast, in another experiment in sapphire Guizard et al. (2010), there was no evidence of ionization avalanche. On the theory side, calculations based upon a Fokker-Planck equation in Stuart et al. (1995) lead to a strong dominance of II over FI while the simulations in Shcheblanov et al. (2012) predict kinetic energies of the charge carriers that are too low for II to be significant. It was also suggested that the condition to trigger an avalanche should be given by the laser fluence instead of the pulse duration, but the predicted threshold values differ by more than an order of magnitude (see, e.g., Rethfeld (2006); Petrov and Davis (2008)).

Actually, experiments involve different materials and laser parameters, which makes a direct comparison between them difficult. Other challenges lie in the theoretical models that are currently used to interpret the experimental observations. On the one hand, current models based on single ionization-rate equations (SRE) do not account for the gradual heating of the charge carriers which, for short laser pulses, might not be sufficient to start an avalanche. On the other hand, models based on multiple rate equations (MRE) that track the carrier kinetics on discrete energy levels rely on several free parameters, which limits the physical insight that we can gain from them. While calculation of the FI rates with the Keldysh theory Keldysh (1965) is well established, models for IBH and II can vary significantly. For example in refs. Rethfeld (2004); Christensen and Balling (2009); Bourgeade et al. (2010); Gallais et al. (2015), plasma formation was modelled within similar theoretical frameworks, but assumed different IBH rates, thus influencing directly the efficiency of II and leading to conflicting conclusions about the relative importance of II over FI and the occurrence of ionization avalanche in short pulses.

In this paper, we describe a model that overcomes the issues associated with the SRE and MRE models by tracking both the plasma density and carriers’ mean kinetic energy as a function of time, forming a set of delayed rate equations (DRE) that we use to match the laser-induced damage threshold of several dielectric materials. In particular, we show that this simplified model reproduces the predictions from the multiple rate equations, with a limited number of free parameters determined unambiguously by fitting experimental data. A side benefit of the DRE model is its computational efficiency, opening the possibilities for large-scale, three-dimensional modelling of laser-induced breakdown of transparent media.

The paper is organized as follows. First in Sec. II we present an overview of the single- and multiple-rate models. Next in Sec. III, we describe the proposed delayed-rate equations model in details. The three models (SRE, MRE, and DRE) are compared in Sec. IV. In Sec. V, we show how the DRE model can fit experimental data for the damage threshold in several dielectric materials. In Sec. VI, we discuss some of the limitations of the DRE model and, ultimately, we conclude in Sec. VII. Three Appendices gather some of the technical aspects of the model and of its implementation. All calculations are performed using a Python package that we made available online Déziel et al. .

II Overview of current rate equation models for laser-induced plasma formation in dielectrics

The modelling of the laser-induced polarization and breakdown dynamics in solid-state dielectrics is typically composed of three complementary pieces. (1) A model for the polarization density from bound electrons. (2) A model for the evolution of the conduction band population due to field ionization, impact ionization, and electron-hole recombination. (3) A model for the free-current density associated with the charge carriers (electrons and holes). This approach provides great modelling flexibility and a fair description of the underlying physics on a cycle-averaged statistical level (see Sec. VI for discussion). Below, we provide an overview of two established population-dynamics models. For reviews of bound and free currents models see, e.g., refs. Couairon et al. (2011); Balling and Schou (2013); Kolesik and Moloney (2014); Varin et al. (2018).

II.1 Single rate equation

The simplest way to describe plasma formation in dielectrics while accounting for both FI and II is the single rate equation (SRE) Stuart et al. (1995):

[TABLE]

where $\rho$ represents the carrier density. The first two terms on the right hand side —associated with the field ionization (FI) rate $\nu_{\mathrm{fi}}$ and the impact ionization (II) rate $\alpha\rho I$, respectively —weighted by the density of neutral molecules or atoms $\rho_{n}=(\rho_{\mathrm{mol}}-\rho)$ (if we account for single ionization at most). Here, $\rho_{\mathrm{mol}}$ is the molecular density, $\alpha$ is the impact rate coefficient, and $I=c\epsilon_{0}n_{0}E^{2}/2$ is the cycle-averaged laser intensity, with $n_{0}$ being the linear refractive index. The last term accounts for the recombination (RE) of electrons and holes at a rate $\gamma_{r}$.

The SRE model was first developed upon empirical observations. The linear relation between the II rate and intensity can be justified by the linearity between the heating rate of the charge carriers and the laser intensity [see Eq. (31) below]. As such, the rate at which the electrons and holes gain energy via IBH and the rate at which they give it back via II both scale linearly with intensity. However, the linear scaling between II and $\rho$ in Eq. (1) implies that all charge carriers can contribute to II, regardless of their energy. This causes an overestimation of II, especially at low fluence. In fact, an electron or a hole needs to acquire a minimum energy $\mathcal{E}_{c}$ to allow a collision where a new valence electron crosses the band gap and reaches the conduction band. To respect both energy and momentum conservation, the critical kinetic energy required for II to be possible is Kaiser et al. (2000)

[TABLE]

where $m_{r}$ is the reduced mass, $m_{h}$ is the effective hole-mass, $\mathcal{E}_{g}$ is the bandgap energy, and $\mathcal{E}_{p}$ is the ponderomotive energy. See Appendix A for the definitions associated with the mass symbols used, and Appendix B for the definition of $\mathcal{E}_{p}$.

The main difference between SRE and more advanced approaches lies in the relation between the II rate and the plasma density $\rho$. In particular, the relation should account for the gradual heating of the carriers by the laser and respect the necessity for them to reach the critical energy $\mathcal{E}_{c}$ for II to occur. Both the MRE [see Sec. II.2] and the DRE [see Sec. III] models address this issue, although with somewhat different ingredients.

II.2 Multiple rate equations

To gain insight into the dielectric breakdown process as a whole, Kaiser et al. Kaiser et al. (2000) have developed a first-principle model that accounts for the various interactions between light, phonons, and the charge carriers. It describes how FI stacks electrons in a single energy level at the bottom of the conduction band (CB), creating a sharp spike in the energy distribution at $\hbar\omega$. When subsequent photon absorption takes place, a new spike appears at $2\hbar\omega$, then another one at $3\hbar\omega$, and so on. After a few femtoseconds, these spikes broaden and disappear due to collisions (thermalization). By tracking dynamically the energy distribution, the number of charge carriers having a minimum kinetic energy of $\mathcal{E}_{c}$ [see Eq. (2)] to contribute to II is then known. II rates can then be scaled with respect to this reduced population (carriers with $\mathcal{E}>\mathcal{E}_{c}$) instead of the entire distribution as done in the SRE model.

A drawback of Kaiser et al.’s  Kaiser et al. (2000) approach is the large number of coupled differential equations that need to be solved (a few hundreds in the case of fused silica). To find a middle ground between simplicity (SRE) and completeness (Kaiser et al. Kaiser et al. (2000)), Rethfeld has developed a multiple rate equations model (MRE) by neglecting thermalization Rethfeld (2004). By doing so, the spikiness of the energy distribution is fully preserved. The energy distribution can then easily be discretized in $k=\lceil\mathcal{E}_{c}/\hbar\omega\rceil$ energy levels (plus one for the zeroth level), each separated by increments of $\hbar\omega$ and associated with an individual population. A rate equation for each level is then solved to track the entire energy distribution. For the electrons in the CB, these rate equations are

[TABLE]

The population at the zeroth level $\rho_{0}$ (the bottom of the CB) is seeded by FI via the first term of the right hand side of Eq. (3a). The next term represents the electrons that are removed from the zeroth level as they absorb photons at a rate $\gamma_{\mathrm{ib}}^{e}$. Each time an electron is removed from any of the $(j-1)$th energy level because of IBH absorption, it is added to the $j$th level, as described by Eq. (3b). After $k$ subsequent photon absorptions, electrons reach the upper $k$th level, with a kinetic energy of at least $\mathcal{E}_{c}$. At this level, IBH is artificially stopped to limit the number of rate equations. From then on, electrons can collide with neutral molecules and cause II events at a rate $\gamma_{n}^{e}$. These electrons then lose their kinetic energy and fall back to the zeroth level while bringing a second electron from the valence band (VB) to the CB (the zeroth level) [see the third term of Eq. (3a)]. Plasma relaxation at the rate $\gamma_{r}$ is also included across all energy levels.

We have extended the original model of Rethfeld  Rethfeld (2004) to account for II events caused by holes. To do so, the term $\gamma_{n}^{h}\rho_{k}^{h}$ was added to Eq. (3a) with the hole-neutral molecule collision rate $\gamma_{n}^{h}$ and the population of holes in the $k$th level $\rho_{k}^{h}$. The latter is calculated with a second set of $k+1$ rate equations that tracks the holes energy distribution. This second set is similar to Eqs. (3), but with $\rho^{h}_{j}$, $\gamma_{n}^{h}$ and $\gamma_{\mathrm{ib}}^{h}$ instead of $\rho_{j}$, $\gamma_{n}^{e}$ and $\gamma_{\mathrm{ib}}^{e}$ respectively. In the special case where electrons and holes have the same mass $m_{e}=m_{h}$, both sets of $k+1$ equations are equivalent and only one has to be solved, with $\gamma_{n}^{h}=\gamma_{n}^{e}$, $\gamma_{\mathrm{ib}}^{h}=\gamma_{\mathrm{ib}}^{e}$ and $\rho_{j}^{h}=\rho_{j}$.

Summing Eqs. (3) leads to the global plasma formation rate as:

[TABLE]

where $s=\{e,h\}$ stands for electrons and holes, respectively. Eq. (4) (MRE) and Eq. (1) (SRE) are identical, except for the second term, associated with II. For the MRE model, the II rate scales with the upper-level populations $\rho_{k}^{e}$ and $\rho_{k}^{h}$. Since the $k$th energy level is populated only after $k$ subsequent photon absorption, II is effectively delayed with respect to FI. The delay for II to unfold is approximately

[TABLE]

when accounting only for II events caused by CB electrons. To account also for VB holes, a distinct delay is set by replacing $\gamma_{\mathrm{ib}}^{e}$ in Eq. (5) by $\gamma_{\mathrm{ib}}^{h}$.

Rethfeld  Rethfeld (2004) has concluded that for pulse duration shorter than $t_{\mathrm{MRE}}$, II is negligible compared with FI because charge carriers are not heated enough to reach $\mathcal{E}_{c}$ before the end of the laser pulse. However IBH is proportional to the laser intensity, i.e., $\gamma_{\mathrm{ib}}\propto I$ (see Appendix B), which suggests that fast, sub-ps carrier heating is possible if the laser intensity is sufficiently high. Thus, a more general condition to trigger an avalanche of ionization through II is $F>It_{\mathrm{MRE}}$ Rethfeld (2006), where $F$ is the laser fluence [see Eq. (18) for definition].

III The delayed-rate equation model

We describe next a delayed-rate equation (DRE) model which addresses the lack of carrier dynamics of the SRE model, while being simpler and less computationally demanding than MRE. Numerical comparison between DRE, SRE, MRE, and experimental data will follow in Secs. IV and V.

We recall that the early energy distribution of the electrons calculated by the full kinetic approach (see ref. Kaiser et al. (2000)) exhibits sharp spikes. This has motivated the development of the MRE model that tracks the electron heating dynamics over discrete momentum levels ($\hbar\omega,2\hbar\omega,\ldots$). However, these spikes quickly broaden and disappear after only a few femtoseconds, due to collisions that drive the energy distribution towards thermal equilibrium. Following a different strategy than for MRE, that assumes that no thermalization takes place, we rely next on the approximation that on a few-laser-cycle timescale, the thermalization process can be considered as almost instantaneous (See Sec. IV for a comparison between both approaches).

Assuming a Maxwellian thermal-equilibrium energy distribution, the fraction of electrons ($s=e$) or holes ($s=h$) that have an energy higher than the critical energy $\mathcal{E}_{c}$ can be calculated analytically as

[TABLE]

where $r_{s}=\sqrt{3\mathcal{E}_{c}/2\mathcal{E}_{\mathrm{k}}^{s}}$. The ratio $\xi^{s}$ and its two contributing terms are shown in Fig. 1 as a function of the dimensionless parameter $r_{s}$. Notice that the individual contributions are almost equal at $r_{s}=0.5$, whereupon the second term rapidly becomes dominant for $r_{s}>1$. With $\xi^{s}$, the equation for the charge-carrier density can be written as

[TABLE]

which is similar to the MRE equation (4), with $\rho_{k}^{s}$ replaced by $\xi^{s}\rho$. The other terms, associated with field ionization ($\nu_{\mathrm{fi}}\rho_{n}$) and electron-hole recombination ($\gamma_{r}\rho$) are identical to both the SRE and MRE models [see Eqs. (1) and (4), respectively].

The simplicity of the DRE model comes from the possibility to track the mean kinetic energy $\mathcal{E}_{\mathrm{k}}^{s}$ of the electrons and holes, instead of the multiple level populations of MRE [see Eqs. (3)]. This is done with the single ordinary differential equation that follows:

[TABLE]

The first term on the right hand side is associated with photon absorption through IBH. The second term represents the kinetic energy lost in an II event. The final terms (in square brackets) ensure energy conservation for each ionization event and redistribute the kinetic energy among the new charge carriers generated by FI and II.

Some insight into the general behaviour of the DRE model would be useful before undertaking the numerical comparison with SRE and MRE in Sec. IV. Initially, i.e., at time $t=0$, there are no charge carriers and the ratio $\xi^{e}$ is identically zero. Field ionization will then bring electrons to the conduction band, and these electrons will be gradually heated up by the laser field, thus increasing the average electron kinetic energy and the value of $\xi^{e}$. At some point, we can expect that laser heating will be balanced by the energy losses from II, i.e., $\gamma_{\mathrm{ib}}^{e}\hbar\omega\simeq\gamma_{n}^{e}\xi^{e}\mathcal{E}_{c}$. If furthermore, we neglect recombination ($\gamma_{r}=0$) and the depletion of the valence electron population ($\rho_{n}\simeq\rho_{\mathrm{mol}}$), the electron population in the conduction band is roughly given by:

[TABLE]

whose solution is

[TABLE]

Note that $\rho^{e}(0)=0$. Eq. (12) shows an exponential increase of the free-electron density, a characteristics of an ionization avalanche. The argument in the exponential function gives the following characteristic time

[TABLE]

In the next section, we will see that Eq. (13) predicts an avalanche delay that is comparable to that obtained with the MRE model [see Eq. (5), with $k\sim\mathcal{E}_{c}/\hbar\omega$], which suggests that the plasma thermalization dynamics (as described by Kaiser et al. Kaiser et al. (2000)) has a limited impact on the avalanche process as a whole.

IV Numerical analysis of the rate models

So far, we have described three rate-equation models (SRE, MRE, and DRE) that track the temporal evolution of the charge-carrier density on a field-cycle-average, statistical level during laser-induced breakdown. We have seen that these models differ only in the way they account for impact ionization and, in particular, for the delay associated with the laser-heating process [compare Eqs. (1), (4), and (9)].

We examine the behaviour of the three rate models with respect to impact ionization by computing the ratio of the plasma density generated by impact ionization $\rho_{\mathrm{ii}}$ over the total plasma density $\rho$ when an harmonic electric field $\tilde{E}(t)=E\cos(\omega t)$ with a constant amplitude $E$ is applied. For each model, the laser intensity $I=cn_{0}\epsilon_{0}E^{2}$ is set to obtain $\rho_{\mathrm{ii}}/\rho\simeq 0.5$ after $t=100\leavevmode\nobreak\ \mathrm{fs}$ ($n_{0}$ is the refractive index of the dielectric without ionization). This intensity marks, for each model, the turning point where impact ionization becomes dominant ($\rho_{\mathrm{ii}}/\rho>0.5$). For the tests that follow, we thus define the fluence threshold for impact ionization avalanche as $F_{\mathrm{av}}=I\cdot 100\,\mathrm{fs}$.

To describe field ionization (FI), we used the Keldysh theory Keldysh (1965) that accounts for both multiphoton and tunnel ionization in a unified framework. To calculate the rate $\nu_{\mathrm{fi}}$ for solid state materials, we rely on the formalism presented in ref. (Balling and Schou, 2013). For convenience, we reproduce these equations in Appendix C, with a slightly different notation. For simplicity, we first neglect recombination ($\gamma_{r}=0$) (this contribution will be taken into account later when we compare DRE with experimental data).

In presenting the model equations in Secs. II.1,  II.2, and III, an explicit description of the laser heating rate $\gamma_{\mathrm{ib}}^{s}$ and the free-carrier-to-neutral impact rate $\gamma^{s}_{n}$ was not given. These two quantities depend on the dynamic properties of the electron-hole plasma. Assuming an harmonic laser electric field $\tilde{E}(t)=E\cos(\omega t)$, the classical Drude model leads to the following expression for the laser-heating rate [see Appendix B, in particular, Eq. (31)]

[TABLE]

In Eq. (14), the plasma damping parameter $\gamma$ accounts effectively for collisions between free carriers (e.g., $\gamma_{e}^{e}$, $\gamma_{e}^{h}$, $\gamma_{h}^{h}$, $\ldots$) and phonons. For direct collisions between charge carriers and neutral molecules $\gamma^{s}_{n}$, we used the model of ref. Balling and Schou (2013), i.e.,

[TABLE]

where $\sigma_{\mathrm{mol}}$ is the molecular impact cross-section. However, the results obtained with DRE are nearly unaffected whether we use Eq. (15) or a constant value for $\gamma_{n}^{s}$. This observation is supported by the work reported in ref. Rethfeld (2004), where it is shown that the value of $\gamma_{n}^{s}$ (or that given by the underlying model) has a small influence, as long as $\gamma_{n}^{s}\gg\gamma_{\mathrm{ib}}^{s}$.

Numerical results for a fictitious material whose properties are comparable to SiO2 are presented in Fig. 2. As expected, we see in Fig. 2(a) that the density of charge carriers generated via II predicted by MRE and DRE are delayed with respect to SRE. For both MRE and DRE, the delays follow the predicted values from Eqs. (5) and (13), $t_{\mathrm{MRE}}=83.2$ fs and $t_{\mathrm{DRE}}=76.5$ fs, respectively. As shown in Fig. 2(b), for MRE and DRE the contribution from II to the total plasma density drops sharply for fluence below the avalanche threshold ($F<F_{\mathrm{av}}$). Above threshold ($F>F_{\mathrm{av}}$), all three models show a similar trend.

We recall that MRE has been developed in the limit of an infinite thermalization time, whereas DRE was developed in the limit of an infinitesimal thermalization time. Our numerous tests reveal that the plasma formation rates are quite similar in both limits, providing compelling evidence that the thermalization time has a small impact upon the plasma formation process as a whole.

To get more insight into the DRE model, we have considered a more realistic scenario where a strong laser pulse is incident on a fictitious material similar to SiO2. The electric field envelope of the laser pulse in vacuum is modelled by a Gaussian function:

[TABLE]

where $\tau$ is the full-width at half-maximum (FWHM) duration of the pulse. The laser intensity and fluence in vacuum are then

[TABLE]

and

[TABLE]

respectively. To account for the intrinsic refractive index $n_{0}$ of the material, as well as for the laser-induced metalization, we computed the electric field in the bulk with:

[TABLE]

where

[TABLE]

and

[TABLE]

This last relation is obtained from the Drude model with a plasma frequency $\omega_{p}^{2}=q^{2}\rho/\epsilon_{0}m_{r}$, updated dynamically as the carrier density $\rho$ grows.

By solving DRE with the gaussian laser source $E_{\mathrm{bulk}}(t)$, we obtain the results displayed in Fig. 3. In the leading edge of the pulse, most of the plasma comes from FI [Figs. 3(a) and 3(d)] as is typically expected. However, as charge carriers get heated up and reach the critical energy $\mathcal{E}_{c}$, plasma growth switches to II.

We have also compared the average kinetic energy of the electrons $\mathcal{E}_{\mathrm{k}}^{e}$ to the Fermi energy $\mathcal{E}_{F}=\hbar^{2}(3\pi^{2}\rho)^{2/3}/2m_{e}$ [see Figs. 3(b) and 3(e)]. Over the entire simulations $\mathcal{E}_{\mathrm{k}}^{e}>\mathcal{E}_{F}$, which suggests that using a Fermi-Dirac distribution to get the ratios $\xi^{s}$, instead of a Maxwellian distribution, should not be necessary. This condition is respected in all the calculations performed here.

A rough estimate of an upper limit for the average kinetic energy of the charge carriers $\mathcal{E}_{\mathrm{k}}^{s}$ is obtained in the regime where $\gamma_{\mathrm{ib}}^{s}\hbar\omega\simeq\gamma_{n}^{s}\xi^{s}\mathcal{E}_{c}$ [see also the paragraph before Eq. (11)]. For moderate laser intensity, only a small fraction of carriers effectively reach the critical energy such that $\mathcal{E}_{\mathrm{k}}^{s}\ll\mathcal{E}_{c}$ at all times. In this regime, $r_{s}\gg 1$ where $\xi^{s}$ is well approximated by $\frac{2r_{s}}{\sqrt{\pi}}\exp(-r_{s}^{2})$, and it is then possible, in combination with (15), to obtain an explicit upper bound

[TABLE]

This approximation is in good agreement with the numerical results shown in in Figs. 3(b) and 3(e) (see dashed lines). Effectively, Eq. (22) predicts that the maximum average kinetic energy should not exceed the value given by the right-hand side of the inequality. For example in Fig. 3(b), the prediction from Eq. (22) is 7.214 eV and the maximum obtained from the numerical integration of DRE is 7.035 eV (a 2.54% overestimation). For the longer pulse duration case in Fig. 3(e), the predicted upper bound is 3.957 eV and the maximal value obtained in the simulation is 3.742 eV (a 5.75% overestimation).

The laser heating rate $\gamma_{\mathrm{ib}}^{e}$ and the electron-neutral collision rate $\gamma_{n}^{e}$ are shown in Figs. 3(c) and 3(f). For comparison, we display as well the electron-electron collision rate given by the following formula (see ref. Christensen and Balling (2009))

[TABLE]

It is then observed that $\gamma_{e}^{e}$ increases rapidly at the leading edge of the pulse, as the plasma gets initially build up by FI. But when II takes over FI, its value levels off to approximately 1 to 10$\leavevmode\nobreak\ \mathrm{fs}^{-1}$, which supports the hypothesis of a fast thermal relaxation and the neglect of the internal thermalization dynamics.

V Calibration of the delayed-rate equation model to experiments

The DRE model presented in Sec. III depends on a closed set of parameters. Some of them can be directly linked to material properties obtained from experimental measurements or ab initio calculations (e.g., $n_{0}$, the electron-impact cross sections, …). Below we show how effective values for the remaining parameters can be obtained by fitting the DRE model to damage-threshold data.

The laser-induced damage threshold is a common reference to benchmark laser-induced dielectric breakdown models. It is often referred to as the minimum laser fluence $F_{\mathrm{th}}$ needed to cause permanent structural modifications to the material. On the plasma formation timescale, the laser-induced damage threshold is associated with the minimum laser fluence needed to create a plasma density $\rho\gtrsim\rho_{c}$ for which the medium becomes opaque to radiation with photon energy $\hbar\omega$. Based on the complex refractive index given at Eq. (21), equating the real and imaginary parts gives the critical density that follows:

[TABLE]

To benchmark the DRE model, we have compared the results obtained by numerical integration of the underlying equations (see Sec. III) with the experimental data found in Mero et al. (2005). Computations were done as in Sec. IV while scanning both the pulse duration $\tau$ and laser fluence $F$. When the maximum carrier density reached the critical density $\rho_{c}$ [see Eq. (24)], the fluence is identified as the fluence threshold. Results are shown in Fig. 4. Fit parameters are given in Table 1.

In practice, the DRE computations shown in Fig. 4 rely only on two “free” parameters ($\gamma$ and $m_{e}$). To optimize the search for the best combination, we proceed as follows. First, we set the plasma damping rate $\gamma$ to adjust the overall scaling trend of the curve to obtain a reasonable agreement with a power-law fit of the experimental data (see below for details). Then, the effective mass parameter $m_{e}$ is chosen to fit the height of the corresponding data set. The parameters are not completely independent however [see, e.g., Eq. (24)] and it is sometimes necessary to iterate the procedure for the final set of parameters. Nevertheless, the computed curves given in Fig. 4 show that DRE succeeds at reproducing the global trend of the experimental measurements over several orders of magnitude of both pulse duration and fluence threshold. Even better fits are obtained if more free parameters are used (e.g., $\gamma$, $m_{e}$, $m_{h}$, and $\mathcal{E}_{g}$).

We emphasize that effective bandgap and mass values are typically obtained by nondestructive measurement methods, where the sample integrity is only slightly perturbed. By definition, assessing the fluence threshold implies driving the material away from the ground state and potentially inducing significant changes to its band structure. The fit values should thus be interpreted with care. Note also that effective masses are usually tensors, to account for the anisotropy of the band structure. Simulations with DRE show that the mass parameters have a significant impact on the damage threshold, which in turn suggests that the orientation of the sample with respect to the laser polarization may play an important role. This effect is likely to be more pronounced in anisotropic crystalline structures. In particular, ab initio calculations of the electronic band structure of $\mathrm{HfO}_{2}$ show that the effective masses along the different crystal planes can vary by more than an order of magnitude Garcia et al. (2004). The effective mass parameters given in Table 1 are consistent with these calculations if they are considered as effective mass values averaged over the different crystal directions.

In Fig. 5, we compare DRE with seven experimental data sets for fused silica. The typical trend across the experiments is that the fluence threshold follow a power-law dependence $F_{\mathrm{th}}\propto\tau^{\kappa}$, with $\kappa\simeq 0.3$ for $\tau<10\leavevmode\nobreak\ \mathrm{ps}$. We could reproduce that trend using DRE and the parameters for SiO2 in Table 1. Experimental data is lacking to rigorously test the model for pulse duration $<10$ fs. However, it is likely that DRE could be improved for such cases to include transient, field-cycle time scale process contributions (see Sec. VI for details). On the other hand, when neglecting laser heating (labelled as FI only), which disables impact ionization completely, the scaling agreement is lost ($\kappa\simeq 0.73$). This supports the fact that near damage threshold impact ionization plays an important role in the dielectric breakdown process, even for few-femtosecond pulse duration.

VI Discussion

We have presented the DRE model as a potential replacement of MRE to study the plasma formation dynamics during laser-induced breakdown in dielectrics. Both models improve upon the SRE model by dealing with the time delay it takes for charge carriers to gain sufficient kinetic energy from the laser field to allow the creation of new charge carriers through impact ionization and trigger an ionization avalanche. DRE and MRE predict similar delays for the first impact ionization events to occur and for a potential ionization avalanche to unfold, with characteristic values for avalanche in the 80 fs range, in agreement with trusted experiments Stuart et al. (1995); Tien et al. (1999). Extended comparison of DRE predictions with experimental data for fused silica shows that the observed damage threshold scaling ($F_{\mathrm{th}}\propto\tau^{0.3}$) can only be explained if laser-heating of the charge carriers and subsequent carrier-impact ionization is taken into account. We have shown that DRE depends on a limited number of parameters that can be unambiguously associated with effective material properties.

There are a number of technical advantages for using DRE instead of MRE. In particular, DRE requires solving less equations, offering interesting possibilities for large scale, three-dimensional calculations where computational efficiency is important. Moreover, in the three-dimensional simulations of laser induced breakdown, e.g., using the finite-difference time-domain (FDTD) or the Particle-in-cell (PIC) frameworks, it is common to see high-contrast structures in the plasma density that strongly enhance or suppress the local electromagnetic field (see, e.g., Déziel et al. (2018)). This causes significant variations in the local ponderomotive energy of the charge carriers and, in turn, of the critical energy for impact ionization [see Eq. (2)]. For MRE, this implies that numerical convergence is dictated by the number of rate equations used. This number must be chosen beforehand to account for the peak values of $\mathcal{E}_{c}$ over the entire simulation and throughout the material domain. This is an important drawback for MRE that should not be overlooked. For DRE, defined by a closed set of equations, this is not an issue.

Finally, it is important to recall that rate equation models in general describe laser-induced breakdown at the field-cycle-averaged level. Future improvements should include proper treatment of photon-assisted avalanche, often referred to as cold ionization avalanche (see, e.g., Rajeev et al. (2009)), as well as potential sub-cycle process contributions (see, e.g., McDonald et al. (2017); Zhokhov and Zheltikov (2014)).

VII Conclusion

We have provided a theoretical framework to study plasma formation during femtosecond laser-induced breakdown in dielectrics on a field-cycle average, statistical level. The model improves upon the current approaches by providing an explicit, closed-formed treatment of the charge-carrier laser-heating process that precedes the onset of carrier-impact ionization and a potential collisional ionization avalanche. In particular, we have shown that the model we propose can reproduce damage-threshold data over several orders of magnitude in both the laser pulse duration and laser fluence, while relying on a limited number of parameters related to effective material properties. A side benefit of the model is its computational efficiency that opens possibilities for large-scale, three-dimensional modelling of laser-induced breakdown and structural pattern formation in transparent media.

Appendix A Definitions for the mass symbols used in this paper

In this paper, the effective mass of the electrons in the conduction band (CB) is denoted by $m_{e}$ and the effective mass of the holes in the valence band (VB) is $m_{h}$. In some cases, the reduced mass $m_{r}^{-1}=m_{e}^{-1}+m_{h}^{-1}$ is used. For example, the IBH rate $\gamma_{\mathrm{ib}}^{e}$ for electrons is calculated with $m_{e}$ and the IBH rate $\gamma_{\mathrm{ib}}^{h}$ for holes is calculated with $m_{h}$ to give a total IBH rate $\gamma_{\mathrm{ib}}=\gamma_{\mathrm{ib}}^{e}+\gamma_{\mathrm{ib}}^{h}$, which can be calculated with $m_{r}$. The total ponderomotive energy of electrons and holes (see Appendix B) is also calculated with $m_{r}$. Finally, we refer to the free electron mass with the symbol $m_{0}$.

Appendix B Drude description of the laser-plasma dynamics

The instantaneous current $\tilde{i}(t)$ associated with the motion of a charge carrier (electron or hole) is conveniently described at a statistical-continuum level by the Drude-like single-carrier model that follows:

[TABLE]

where $\tilde{E}(t)$ is the electric field of the laser. Parameters $q$ and $m$ are the charge and mass of the charge carrier, respectively. Collisions are included phenomenologically via the damping rate $\gamma$. Given the carrier density $\rho(t)$, a current density is then defined as $\tilde{J}(t)=\rho(t)\tilde{i}(t)$.

For $\tilde{E}(t)=E\cos(\omega t)$, the steady-state solution for the single-carrier current is:

[TABLE]

Then, the power transferred instantaneously from the laser field to the charge carrier is given by

[TABLE]

The two terms in the square brackets are associated with the ponderomotive energy and inverse bremsstrahlung heating, described below.

B.1 Ponderomotive energy

The first term in the square brackets of Eq. (27) represents a carrier that gains a certain amount of energy during half of an optical cycle, before losing it during the other half, resulting in no net energy gain or loss. This is often referred to as the ponderomotive energy, whose instantaneous expression is given by the integral of the ponderomotive power, i.e., of the first term in Eq. (27), such that

[TABLE]

In general, the ponderomotive energy is expressed instead in terms of its cycle-averaged expression

[TABLE]

that reduces to the usual, free-particle expression $\mathcal{E}_{p}=q^{2}E^{2}/4m\omega^{2}$ in the limit where $\gamma=0$.

B.2 Inverse bremsstrahlung heating

The last term of Eq. (27) is associated with the absorption by the charge carrier of electrical power from the laser field resulting in a net energy gain after each optical cycle. The rate at which a quantum of light is absorbed is obtained by dividing the last term of Eq. (27) by the energy of a photon $\hbar\omega$, thus defining an instantaneous laser-heating rate as

[TABLE]

When averaged over a field cycle:

[TABLE]

Physical insight into the continuum model for the ponderomotive energy and laser-heating rate in the presence of collisions is provided in Fig. 6 [where we used a constant value for $\rho(t)$]. In the free-particle limit ($\gamma=0$), no photon is absorbed, which results in a purely ponderomotive regime ($\gamma_{\mathrm{ib}}=0$). But as $\gamma$ is increased, the amplitude of the current density decreases and the phase difference with respect to the field oscillations $\Delta\phi=\arctan(-\omega/\gamma)$ gradually shifts from $-\pi/2$ to [math], with no energy transfer to the charge carriers in the limit $\gamma\rightarrow\infty$. Optimal heating occurs when $\gamma=\omega$.

Appendix C Keldysh model for field ionization in solid-state dielectrics

The production rate of electron-hole pairs (in $\mathrm{m}^{-3}\mathrm{s}^{-1}$) induced by a strong laser field $\tilde{E}(t)=E\cos(\omega t)$ in a solid-state dielectric with bandgap energy $\mathcal{E}_{g}$ is given by the Keldysh relation (for details, see ref. (Balling and Schou, 2013), Sec. 2.3.1 and Couairon and Mysyrowicz (2007), Sec. 2.3)

[TABLE]

where

[TABLE]

with $\mathcal{K}()$ and $\mathcal{E}()$ being the complete elliptic integrals of the first and second kind, respectively, $\Phi()$ being Dawson’s integral, and $\lfloor\ldots\rfloor$ denoting the integral part of the argument. The free-particle ponderomotive energy $\mathcal{E}_{p}=q^{2}E^{2}/2m_{r}\omega^{2}$ defines the Keldysh parameter as $\Gamma=(\omega/qE)\sqrt{m_{r}\mathcal{E}_{g}}$ (see ref. (Balling and Schou, 2013)).

To get an FI rate $\nu_{\mathrm{fi}}(\gamma)$ compatible with rate-equation models [e.g., Eqs. (1), (4), and (9)], the Keldysh rate $W$ (in $\mathrm{m}^{-3}\mathrm{s}^{-1}$) is divided by the molecular density $\rho_{\mathrm{mol}}$ (in $\mathrm{m}^{-3}$) of the material. The resulting, single-molecule ionization rate $\nu_{\mathrm{fi}}=W/\rho_{\mathrm{mol}}$ is plotted in Fig. 7 for different values of the band gap energy $\mathcal{E}_{g}$.

Acknowledments

C.V. acknowledges financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) through the College and Community Innovation Program - Innovation Enhancement Grants (CCIPE 517932-17) and the Fonds de recherche du Québec - Nature et technologies (FRQNT) through the Programme de recherche pour les chercheurs et les chercheuses de collège (2019-CO-254385).