2.1 Keldysh photoionization

Our work extends the existing ionization framework already available in EPOCH^[ Reference Arber, Bennett, Brady, Lawrence-Douglas, Ramsay, Sircombe, Gillies, Evans, Schmitz, Bell and Ridgers ^{31 ]} to include the photoionization model developed by Keldysh^[ Reference Keldysh ^{42 ]}. The Keldysh ionization model allows us to calculate probabilities for electron transitions from the valence to the conduction band in solids due to the electric field of a laser pulse^[ Reference Gruzdev ^{47 ]}.

For a laser pulse with electric field amplitude $E$ and frequency $\omega$ interacting with a material having band gap $\Delta$ and reduced electron-hole mass ${m}^{\ast }$ , the Keldysh parameter is $\gamma =\omega \sqrt{m^{\ast}\Delta}/ eE$ , where $\gamma >>1$ is for the multiphoton ionization regime and $\gamma <<1$ is for the tunneling regime^[ Reference Keldysh ^{42 ]}. The Keldysh formulation is especially useful as it spans both regimes, which are often present when considering a laser-induced damage experiment. For example, the peak electric fields in our simulations give values of $\gamma$ ranging from about 0.3 to 0.9, which means the photoionization regime is neither multiphoton nor tunneling, suggesting that the full Keldysh formula should be utilized.

The Keldysh model is used in our simulations to calculate the probability that a given neutral macroparticle will ionize. When ionization occurs, a new electron macroparticle and corresponding ion macroparticle are created in place of the original macroparticle. Now we introduce the ionization rate equation used in our work. For brevity in the following expressions, we follow Refs. [Reference Tien, Backus, Kapteyn, Murnane and Mourou48,Reference Balling, Schou and Rep49] by introducing the variables ${\gamma}_1={\gamma}^2/\left(1+{\gamma}^2\right)$ and ${\gamma}_2=1/\left(1+{\gamma}^2\right)$ . The effective band gap is then given by the following:

(1) $$\begin{align}x=\frac{2}{\pi}\frac{\Delta}{\mathrm{\hslash}\omega}\frac{\epsilon \left({\gamma}_2\right)}{\sqrt{\gamma_1}}.\end{align}$$

We may then write the ionization rate $W$ as follows:

(2) $$\begin{align}W\left[{\mathrm{m}}^{-3}{\mathrm{s}}^{-1}\right]& =2\frac{2\omega }{9\pi }{\left(\frac{m^{\ast}\omega }{\mathrm{\hslash}\sqrt{\gamma_1}}\right)}^{3/2} \nonumber\\ &\quad \times Q\left(\gamma, x\right)\exp \left(-\pi \lfloor x+1 \rfloor \frac{\kappa \left({\gamma}_1\right)-\epsilon\left({\gamma}_1\right)}{\epsilon\left({\gamma}_2\right)}\right),\end{align}$$

where $\kappa$ and $\epsilon$ are complete elliptic integrals of the first and second kind, $\Phi$ is the Dawson integral and

(3) $$\begin{align}Q\left(\gamma, x\right)&=\sqrt{\frac{\pi }{2\kappa \left({\gamma}_2\right)}}\sum \limits_{n=0}^{\infty}\exp \left(- n\pi \frac{\kappa \left({\gamma}_1\right)-\epsilon\left({\gamma}_1\right)}{\epsilon\left({\gamma}_2\right)}\right)\nonumber\\&\quad \times \Phi \left(\sqrt{\frac{\pi^2\left(\lfloor x+1 \rfloor -x+n\right)}{2\kappa \left({\gamma}_2\right)\epsilon\left({\gamma}_2\right)}}\right).\end{align}$$

We note that this is the corrected version of the formulation, where the original contains a misprint as noted by Gruzdev^[ Reference Gruzdev ^{43 ]}. Using the uncorrected version can result in significantly different calculations^[ Reference Zhang, Menoni, Gruzdev and Chowdhury ^{28 ,} Reference Gruzdev ^{43 ]}.

For our simulation framework, the Keldysh ionization rate is evaluated in situ (with the exception of the elliptic integrals, which are read from tabulated data files with 1000 points) rather than interpolated from a tabulated form, which facilitates accuracy over a wide range of fluences. We use the first 500 terms in this infinite sum in Equation (3), which is sufficient for the regimes of interest in this paper, as illustrated in the appendix, while limiting the computational cost. An adaptive approach to a prescribed tolerance could be employed in the future. Currently each macroparticle can only be ionized once.

While the Keldysh parameter depends on the laser amplitude $E$ , the amplitudes of our pulses vary within the laser envelope. Moreover, each computational cell in a PIC simulation only considers the instantaneous electric field. To address this, we store the electric field magnitude at each time step using an array large enough to store an entire laser cycle on a rolling window, as suggested by Zhang et al. ^[ Reference Zhang, Menoni, Gruzdev and Chowdhury ^{28 ]}. Then the maximum field for the previous cycle is used to approximate the amplitude to calculate the ionization rate. This approach would encounter challenges for single-cycle pulses. In the future a more complex envelope model could be explored^[ Reference Terzani and Londrillo ^{50 ]}.

2.2 Collisional effects

We use the Pérez/Nanbu^[ Reference Pérez, Gremillet, Decoster, Drouin and Lefebvre ^{51 ,} Reference Nanbu and Yonemura ^{52 ]} binary collision module already included in the EPOCH code to account for collisions. The collision frequency in EPOCH is calculated for a charged particle $\alpha$ with charge ${q}_{\alpha }$ scattering off a charged particle $\beta$ as follows:

(4) $$\begin{align}{\nu}_{\alpha \beta}=\frac{{\left({q}_{\alpha }{q}_{\beta}\right)}^2{n}_{\beta}\ln \left(\Lambda \right)}{4\pi {\left({\varepsilon}_0\mu \right)}^2}\frac{1}{v_\mathrm{r}^3},\end{align}$$

where ${n}_{\beta }$ is the density (for particles of species $\beta$ ), $\ln \left(\Lambda \right)$ is the Coulomb logarithm, $\mu ={m}_{\alpha }{m}_{\beta }/({m}_{\alpha }+{m}_{\beta })$ is the reduced mass and ${v}_\mathrm{r}$ is a relative velocity^[ Reference Arber, Bennett, Brady, Lawrence-Douglas, Ramsay, Sircombe, Gillies, Evans, Schmitz, Bell and Ridgers ^{31 ]}. EPOCH extends the model to low temperatures following an approach by Pérez et al. ^[ Reference Pérez, Gremillet, Decoster, Drouin and Lefebvre ^{51 ]} and Lee and More^[ Reference Lee and More ^{53 ]}. More details can be found in the work of Arber et al. ^[ Reference Arber, Bennett, Brady, Lawrence-Douglas, Ramsay, Sircombe, Gillies, Evans, Schmitz, Bell and Ridgers ^{31 ]} and in the source files and documentation provided with the open-source EPOCH code (https://epochpic.github.io/). Collisions are included between all charged particles in the simulation. The Coulomb logarithm is calculated automatically with a fixed lower bound of 1.

EPOCH includes a collisional ionization routine designed for atoms^[ Reference Arber, Bennett, Brady, Lawrence-Douglas, Ramsay, Sircombe, Gillies, Evans, Schmitz, Bell and Ridgers ^{31 ,} Reference Morris, Goffrey, Bennett and Arber ^{54 ]}, but this is not well suited for impact ionization in dielectrics. A more appropriate ionization rate for our regime can be calculated with the approach by Keldysh^[ Reference Keldysh ^{55 ]}, although some of the input parameters to the models are not well reported and may require fitting of output results to experimental data^[ Reference Charpin, Ardaneh, Morel, Giust and Courvoisier ^{41 ,} Reference Déziel, Dubé and Varin ^{56 ]}. Impact ionization is reduced for shorter few-cycle pulses. Models by Petrov and Davis^[ Reference Petrov and Davis ^{57 ]} suggest collisional ionization dominates photoionization for fluences exceeding $0.4\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ , whereas the multiple-rate-equation (MRE) model by Rethfeld^[ Reference Rethfeld ^{58 ]} predicts this threshold to be $10\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ ^[ Reference Balling, Schou and Rep ^{49 ]}. For this work, we do not include impact ionization, since we model few-cycle pulses that are only 7 fs (fewer than three laser cycles) full width at half maximum (FWHM) in duration. We find good agreement with previous experiments without the need for impact ionization, but expect this to be an important consideration for longer pulses in future work.

2.3 Refraction

The optical properties of dielectrics are modeled using a spatially varying permittivity $\varepsilon$ throughout the simulation box. At the beginning of the simulation the optical properties are stored in a matrix with the same size as the simulation grid. We then modified the field solver to include a spatially dependent permittivity when advancing Maxwell’s equations, following a similar approach to that in the WarpX code^[ Reference Fedeli, Huebl, Boillod-Cerneux, Clark, Gott, Hillairet, Jaure, Leblanc, Lehe, Myers, Piechurski, Sato, Zaim, Zhang, Vay and Vincenti ^{46 ]} and the modification to EPOCH by Charpin et al. ^[ Reference Charpin, Ardaneh, Morel, Giust and Courvoisier ^{41 ]}. To easily account for arbitrary target structures, we define the shape of the optical region at the same place where particle species are initialized. There has been work including nonlinear optical effects in PIC and FDTD simulations^[ Reference Déziel, Dubé, Messaddeq, Messaddeq and Varin ^{36 ,} Reference Varin, Emms, Bart, Fennel and Brabec ^{59 ]}, although those effects are not considered here.

3.1 Grid and particle initialization

The simulation setup for bulk ${\mathrm{SiO}}_2$ is shown in Figure 1(a). The simulation box is 1.0 μm in the longitudinal ( $y$ ) direction and 9.6 μm in the transverse ( $x$ / $z$ ) directions. The target is 0.9 μm thick, leaving a 0.1 μm vacuum region before the laser interacts with the target. Simple outflow boundaries are used to allow transmission of the laser out of the simulation box with minimal reflection. The bulk silica simulations have a resolution of 5 nm in the longitudinal ( $y$ ) direction and 40 nm in the transverse directions. A higher resolution was used in the longitudinal direction to better resolve the ionization dynamics at the interface between the target surface and the vacuum region.

Then we apply this framework to multi-layer interference coatings composed of alternating fused ${\mathrm{SiO}}_2$ (yellow) and ${\mathrm{HfO}}_2$ (blue) layers, as shown in Figures 1(b) and 1(d). The simulation space is 9.6 μm by 1.533 μm by 9.6 μm long in the x, y and z directions, respectively, with a resolution of 13.8 nm in the longitudinal ( $y$ ) direction and 40 nm in the transverse directions. The thickness of the surface protective ${\mathrm{SiO}}_2$ layer is $\lambda$ /(2n), while the rest of the layer thickness is $\lambda$ /(4n), a typical Bragg quarter-wavelength mirror, where $n$ is the index of refraction for each material. The top fused ${\mathrm{SiO}}_2$ layer is placed at 0 μm and the pulse enters normally from the y-axis at –0.1 μm.

Both series of simulations are initialized with 1000 neutral ${\mathrm{SiO}}_2$ (or ${\mathrm{HfO}}_2$ ) macroparticles per cell with a temperature of 300 K. Similar to Charpin et al. ^[ Reference Charpin, Ardaneh, Morel, Giust and Courvoisier ^{41 ]}, we found that a large number of particles per cell were required for accuracy with the Keldysh ionization model in this regime. The simulations are run to a simulation time of 24 fs, which allows the pulse to make multiple passes across the simulation domain. For each simulation, we use the default time step in EPOCH of 0.95 times the Courant–Friedrichs–Lewy (CFL) limit^[ Reference Courant, Friedrichs and Lewy ^{60 ,} Reference Courant, Friedrichs and Lewy ^{61 ]}, or $0.95/\left(c\sqrt{1/\Delta {x}^2+1/\Delta {y}^2+1/\Delta {z}^2}\right)$ , where $\Delta x,\Delta y \text{ and } \Delta z$ represent the grid size in each simulation dimension^[ Reference Arber, Bennett, Brady, Lawrence-Douglas, Ramsay, Sircombe, Gillies, Evans, Schmitz, Bell and Ridgers ^{31 ]}.

3.2 Material properties

The simulations require a number of material properties as inputs to model the interaction and interpret the predictions. These include the linear refractive index, density, and band gap, for which we use standard values listed in Table 1. There is less agreement in reported values for the effective electron and hole masses ( ${m}_\mathrm{e}^{\ast }$ and ${m}_\mathrm{h}^{\ast }$ , respectively) and subsequently this results in different reduced effective masses ${m}^{\ast}=1/\left(1/{m}_\mathrm{e}^{\ast }+1/{m}_\mathrm{h}^{\ast}\right)$ . This variation leads to significant differences in predictions for the excited electron density using the Keldysh ionization model^[ Reference Zhang, Davenport, Talisa, Smith, Menoni, Gruzdev and Chowdhury ^{27 ,} Reference Gruzdev ^{65 ]}. To calculate ${m}_\mathrm{e}^{\ast }$ for fused ${\mathrm{HfO}}_2$ , we assume it is the spherically averaged effective mass around the $\Gamma$ and B points of the monoclinic ${\mathrm{HfO}}_2$ ^[ Reference Perevalov, Gritsenko, Erenburg, Badalyan, Wong and Kim ^{64 ,} Reference Zhang, Tripepi, AlShafey, Talisa, Nguyen, Reagan, Sistrunk, Gibson, Alessi and Chowdhury ^{66 ]}. We generally use material properties for m- ${\mathrm{HfO}}_2$ as those for amorphous ${\mathrm{HfO}}_2$ are less readily available. These effective electron masses are used for ionized electron particles in the simulations.

Table 1 Material properties used in simulations.

4.1 Electron density

For simulation and theoretical work, the predicted excited electron density is often used as a criterion to predict damage. Many studies use the critical electron density^[ Reference Chen ^{69 ]} for free electrons or some fraction of total ionization^[ Reference Werner, Gruzdev, Talisa, Kafka, Austin, Liebig and Chowdhury ^{70 ]}. This qualitatively makes sense, as exceeding such thresholds can result in high absorption and subsequent damage. This choice has shown good agreement with longer pulses, although recent work has suggested that this description is insufficient for modeling shorter few-cycle pulses^[ Reference Chimier, Utéza, Sanner, Sentis, Itina, Lassonde, Légaré, Vidal and Kieffer ^{16 ,} Reference Zhokhov and Zheltikov ^{71 ]}. As shown in Figure 2, the damage threshold is predicted at about $0.8\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ by the critical density criterion, which is about a 30% underestimation of the experimental LIDT. Number density is a standard output variable for PIC simulations. It is calculated by mapping the position of the macroparticles to the spatial grid based on the shape function selected for the simulation^[ Reference Birdsall and Langdon ^{29 –} Reference Arber, Bennett, Brady, Lawrence-Douglas, Ramsay, Sircombe, Gillies, Evans, Schmitz, Bell and Ridgers ^{31 ]}.

Figure 2 The electron density at the center of the x–z plane (averaged over six cells in x and y) along y at 20 fs for a series of PIC simulations at various laser fluences with a 7 fs duration pulse interacting with bulk ${\mathrm{SiO}}_2$ . The critical plasma density ${n}_\mathrm{crit}$ and electron instability density from Ref. [Reference Stampfli and Bennemann72] are labeled with dashed lines. The latter predicts damage around $1.2\;\mathrm{J}\;{\mathrm{cm}}^{-2}$ .

Alternatively, the instability density suggested by Stampfli and Bennemann^[ Reference Stampfli and Bennemann ^{72 ]} states that when the conduction band electron density reaches about 9% of the valence band electron density, the elastic shear constant will become negative and the lattice becomes unstable, which then leads directly to a very rapid melting of the crystal structure. This criterion was developed for crystals, but the fundamental physical principles – namely, the relationship between conduction band electron density and the stability of the atomic structure – are similar. By applying this criterion to the bulk fused ${\mathrm{SiO}}_2$ simulations, the damage is achieved at about $1.2\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ , which agrees with the experimental results.

4.2 Energy density

The excited electron energy density criterion has also been suggested to predict damage^[ Reference Zhokhov and Zheltikov ^{71 ,} Reference Grehn, Seuthe, Höfner, Griga, Theiss, Mermillod-Blondin, Eberstein, Eichler and Bonse ^{73 ]}. The predicted energy density is typically compared to material properties, such as the dissociation energy, or an energy barrier associated with melting or boiling^[ Reference Zhokhov and Zheltikov ^{71 ,} Reference Wang, Qi, Zhao, Wang and Shao ^{74 ]}. We compare the energy densities in our simulations to these thresholds.

Previous computational approaches typically assume some simple electron energy distribution, whereas in PIC simulation, the energy density can be calculated directly with standard outputs. For our simulations, we multiply the number density by average particle energy for a species in each cell. Alternately, this could be re-sampled to a finer or a coarser grid if the individual macroparticle positions and energies are extracted from the simulation.

Due to variations of reported material properties in the literature and uncertainty of previous simulations, the exact energy density for damage is not agreed upon. For reference, the dissociation energy of ${\mathrm{SiO}}_2$ has reported values from $54\;\mathrm{to}\;68\ \mathrm{kJ}\;{\mathrm{cm}}^{-3}$ (Refs. [Reference Grehn, Seuthe, Höfner, Griga, Theiss, Mermillod-Blondin, Eberstein, Eichler and Bonse73,Reference Jia, Xu, Li, Feng, Li, Cheng, Sun, Xu and Wang75,Reference Inaba, Fujino and Morinaga76]). For comparison, the threshold for high-energy-density physics^[ Reference Drake ^{77 ]} is approximately $100\ \mathrm{kJ}\;{\mathrm{cm}}^{-3}$ .

The energy densities related to melting or boiling are lower, where we can use the temperature-dependent heat capacity and latent heat of vaporization used by Zhao et al. ^[ Reference Zhao, Sullivan, Zayac and Bennett ^{78 ,} Reference Zhao, Cheng, Chen, Yuan, Liao, Wang, Liu and Yang ^{79 ]} to calculate an energy density of $5.7\;\mathrm{kJ}\;{\mathrm{cm}}^{-3}$ for melting and $34.7\;\mathrm{kJ}\;{\mathrm{cm}}^{-3}$ for boiling. We do have uncertainty in these values. The latent heat of vaporization has the largest contribution to the boiling criteria and there is a great deal of uncertainty in reported values. For example, the calculation above uses values from Bäuerle^[ Reference Bäuerle ^{80 ]}, who calculated the latent heat of vaporization of c- ${\mathrm{SiO}}_2$ to be $1.23\times {10}^7\;\mathrm{J}\;{\mathrm{kg}}^{-1}$ , while Kraus et al. ^[ Reference Kraus, Stewart, Swift, Bolme, Smith, Hamel, Hammel, Spaulding, Hicks, Eggert and Collins ^{81 ]} reported $(1.177\pm 0.095)\times {10}^7\;\mathrm{J}\;{\mathrm{kg}}^{-1}$ and Khmyrov et al. ^[ Reference Khmyrov, Grigoriev, Okunkova and Gusarov ^{82 ]} used $0.96\times {10}^7\;\mathrm{J}\;{\mathrm{kg}}^{-1}$ (from Refs. [Reference Grigoriev and Meilikhov83,Reference Samsonov84]). This gives a range from $28.7$ to $35.6\;\mathrm{kJ}\;{\mathrm{cm}}^{-3}$ for boiling. We expect some uncertainty in melting energy density as well.

Figure 3 shows the maximum energy density in simulations with and without collisions for a range of laser fluences. The simulations just including photoionization, (without collisions) have lower energy density at the end of the simulation and do not exceed the boiling criterion until much higher fluences are expected for the LIDT in these conditions.

Figure 3 The peak energy density for the series of PIC simulations at various fluences with a 7 fs duration pulse interacting with bulk fused silica. The experiment^[ Reference Chimier, Utéza, Sanner, Sentis, Itina, Lassonde, Légaré, Vidal and Kieffer ^{16 ]} being modeled observed damage around $1.18\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ and ablation^[ Reference Utéza, Sanner, Chimier, Brocas, Varkentina, Sentis, Lassonde, Légaré and Kieffer ^{21 ]} at $1.3\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ . The shaded horizontal line/bands indicate approximate energy density thresholds for melting, boiling and dissociation. The shaded vertical bands represent uncertainty in experimental damage and ablation thresholds. The simulations including collisions have a higher predicted energy density. Including both our Keldysh photoionization model and collisional effects shows agreement between the expected damage fluence and the dissociation energy.

In Figure 3, we see that at $1.3\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ , where ablation is observed in experiments, the simulation energy density overlaps with the reported dissociation energy values. For $1.2\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ , near the LIDT, the energy density for the simulation is around $49\ \mathrm{kJ}\;{\mathrm{cm}}^{-3}$ , exceeding the boiling threshold and near the dissociation energy.

4.3 Kinetic particle motion

The kinetic nature of PIC simulations allows us to explore the energy and motion of the excited electrons. While most previous approaches assume a thermal distribution of the excited electrons, Figure 4 shows that this is not the case during the interaction. The spectra at different times of the simulation are shown and the energy is fitted by the ${\chi}^2$ distribution using the SciPy^[ Reference Virtanen, Gommers, Oliphant, Haberland, Reddy, Cournapeau, Burovski, Peterson, Weckesser, Bright, van der Walt, Brett, Wilson, Millman, Mayorov, Nelson, Jones, Kern, Larson, Carey, Polat, Feng, Moore, VanderPlas, Laxalde, Perktold, Cimrman, Henriksen, Quintero, Harris, Archibald, Ribeiro, Pedregosa and van Mulbregt ^{85 ]} library:

(5) $$\begin{align}{\chi}^2:f\left({E}_\mathrm{e};k,\sigma \right)=\frac{E_\mathrm{e}^{k/2-1}{e}^{-\frac{E_\mathrm{e}}{2\sigma }}}{{\left(2\sigma \right)}^{k/2}\Gamma \left(\frac{k}{2}\right)},\\[-24pt]\nonumber \end{align}$$

where ${E}_\mathrm{e}$ is the energy of the electrons and fitting parameters $k$ , $\sigma$ are the degree and scale parameters, respectively. The ${\chi}^2$ distribution becomes the standard Maxwell–Boltzmann distribution when $k$ = 3 and $\sigma$ = ${k}_\mathrm{b}T$ /2. In our analysis, we primarily focus on the variation of the degrees of freedom $k$ as it is crucial in describing the main characteristics of the Maxwell–Boltzmann distribution. For both cases, with or without collisions, the excited electrons are highly nonthermal at the early stage of 4 fs, where the degree $k$ is about 0.7. Then at 10 fs, after the peak intensity passed through the target, the energy spectrum is still nonthermal with k = 1.07 when collisions are not included, as shown in Figure 4(a), while in Figure 4(b) the degree $k$ reaches 1.86. At the stable stage of 23 fs, the electrons still remain nonthermal, as shown in Figure 4(a), while in Figure 4(b), the $k$ value is about 2.56, which is approaching the Maxwell–Boltzmann distribution as indicated by the black curve giving the average energy at 23 fs. Therefore, our simulations not only show the dynamic evolution of the excited electron energy spectrum during the interaction but also show the importance of including collisions to capture particle dynamics. When collisions are included in the simulations, the energy absorbed by the electrons increases, leading to higher average energy (Figure 4) and higher maximum energy density (Figure 3) at the end of the simulation.

Figure 4 Excited electron energy distributions in a $200\ \mathrm{nm} \times 200\ \mathrm{nm} \times 900\ \mathrm{nm}$ region of the target at the center of the interaction for a $1.2\;\mathrm{J}\;{\mathrm{cm}}^{-2}$ pulse. A simulation with just field ionization is shown in (a) and a simulation with field ionization and collisions is shown in (b). A fit with the given temperature is shown. We observe the nonthermal nature, especially for early times and for simulations without collisions.

5.1 Electron density

As mentioned in Section 4, the plasma critical density may underestimate the damage threshold; it predicts the LIDT slightly above $0.2\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ in the first ${\mathrm{HfO}}_2$ layer, as shown in Figure 5. If we apply the instability density criterion to the mirror target, the damage may be achieved at about $0.5\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ in the same layer. Above $0.7\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ , the hot spots in the top two layers are almost fully ionized as the peak electron density reaches about $1.5\times {10}^{22}\ {\mathrm{cm}}^{-3}$ .

Figure 5 The electron density at the center of the x–z plane along y at 20 fs for a series of PIC simulations at various fluences with a 7 fs duration pulse interacting with a multi-layer mirror. The yellow area is ${\mathrm{SiO}}_2$ and blue area is ${\mathrm{HfO}}_2$ . Different critical electron densities are labeled in the figure with dashed lines.

5.2 Energy density

The peak electron energy density along y at the center of the x–z plane in each layer with different fluences is shown in Figure 6. The dissociation energy of fused ${\mathrm{HfO}}_2$ has not been extensively studied, particularly for amorphous samples, which are expected in these multi-layer coatings. There are reported values for the formation energy density of m- ${\mathrm{HfO}}_2$ from 48.44 to $52.64\ \mathrm{kJ}\;{\mathrm{cm}}^{-3}$ ^[ Reference Wang, Zinkevich and Aldinger ^{88 –} Reference Panish and Reif ^{94 ]}, which are similar to the dissociation energy density of ${\mathrm{SiO}}_2$ . We assume it has the same value as fused ${\mathrm{SiO}}_2$ since both of their molecules have four valence band electrons and are amorphous^[ Reference Zhang, Menoni, Gruzdev and Chowdhury ^{28 ]}, which is about $54\;\mathrm{kJ}\;{\mathrm{cm}}^{-3}$ , indicated by the red dashed line. In addition, the effects of structure on boiling and melting points are not considered, as each layer is assumed to retain the melting and boiling points characteristic of its bulk state. This criterion suggests the damage should occur at the surface at fluence of about $1.4\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ , which is even higher than the LIDT of bulk fused ${\mathrm{SiO}}_2$ discussed in Section 4.

Figure 6 The peak electron energy density at the center of the x–z plane along y at 20 fs for a series of PIC simulations at various fluences with a 7 fs duration pulse interacting with a multi-layer mirror. The yellow area is fused ${\mathrm{SiO}}_2$ and blue area is fused ${\mathrm{HfO}}_2$ . The ${\mathrm{SiO}}_2$ boiling and melting energy densities are at about $34.7$ and $5.7\;\mathrm{kJ}\;{\mathrm{cm}}^{-3}$ , respectively, indicated as yellow dashed lines, and the ${\mathrm{HfO}}_2$ boiling and melting energy densities are at about $49.9$ and $10.7\;\mathrm{kJ}\;{\mathrm{cm}}^{-3}$ , respectively, indicated as blue dashed lines.

Instead let us consider energy density thresholds for melting and boiling, as these may relate to damage within the layers of a coating. To calculate the melting and boiling energy densities, we make the following assumptions: (i) the vaporization latent heat for ${\mathrm{HfO}}_2$ is the same as that of ${\mathrm{SiO}}_2$ and (ii) the heat capacity for fused ${\mathrm{HfO}}_2$ is the same as that of m- ${\mathrm{HfO}}_2$ ^[ Reference Low, Paulson, D’Mello and Stan ^{95 ]}. These approximations give the boiling energy density of ${\mathrm{HfO}}_2$ to be about $49.9\;\mathrm{kJ}\;{\mathrm{cm}}^{-3}$ , which is much higher than that of ${\mathrm{SiO}}_2$ , $34.7\;\mathrm{kJ}\;{\mathrm{cm}}^{-3}$ , and may overestimate the actual damage threshold. Similarly, the melting energy density of ${\mathrm{SiO}}_2$ is about $5.7\;\mathrm{kJ}\;{\mathrm{cm}}^{-3}$ , and we calculated about $10.7\;\mathrm{kJ}\;{\mathrm{cm}}^{-3}$ for ${\mathrm{HfO}}_2$ .

The boiling energy density criterion predicts the LIDT at a fluence between $1.1$ and $1.2\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ on the surface of the mirror, which is close to the LIDT of bulk ${\mathrm{SiO}}_2$ . Alternately, applying the melting energy density criterion gives the LIDT to be less than $0.5\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ , and the damage site is initiated in the first ${\mathrm{HfO}}_2$ layer, as expected.

5.3 Plasma screening effects

As shown in Figure 6, the global maximum energy density at low fluences is in the first ${\mathrm{HfO}}_2$ layer, as expected. When the fluence exceeds about $1.1\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ it shifts to the top ${\mathrm{SiO}}_2$ layer and the ${\mathrm{HfO}}_2$ energy density increases slowly after the fluence reaches about $0.7\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ . This is because the excited electron density in the first ${\mathrm{SiO}}_2$ layer begins to exceed the critical plasma density (Figure 5). In the top ${\mathrm{SiO}}_2$ layer, the steady state of the dynamic simulation leads to a local maximum enhancement of the electric field at the center.

As the fluence increases beyond $0.5\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ , we observe a new peak appearing and shifting from the center to the surface in both the electron and energy density profiles. The ionization rate at the center of the top ${\mathrm{SiO}}_2$ layer is enhanced due to the strong intensity and thus the electron density reaches a maximum. The absorption and reflection will continue to increase so that the source can hardly penetrate the target by the end of the laser pulse. Therefore, the resonant pattern of the electric field is altered, and a new intensity peak appears.

To illustrate this process clearly, a two-dimensional y–z cross-section of the layers at the center of the x-axis is shown in Figure 7 at fluence of $0.7\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ . The maximum electron and energy densities are both observed in the first ${\mathrm{HfO}}_2$ layer, as shown in Figures 7(a) and 7(b), while the maximum accumulated intensity is in the first ${\mathrm{SiO}}_2$ layer in Figure 7(c). The intensity is recorded at intervals of 0.5 fs and compared across all time points to generate this figure. These results are similar to those from the FDTD simulations by Zhang et al. ^[ Reference Zhang, Menoni, Gruzdev and Chowdhury ^{28 ]}.

Figure 7 The (a) energy density and (b) electron density at the center of x on the y–z plane at 20 fs and $0.7\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ with a 7 fs duration pulse interacting with a multi-layer mirror. The mirror surface starts at y = 0, and the dashed red lines are the interfaces between layers. The white areas have no excited electrons. (c) The maximum accumulated normalized intensity over the entire simulation.

Noticeably, the electron density, energy density and intensity profiles in the protective ${\mathrm{SiO}}_2$ layer all present a bump extending from the center to the surface, which is different from in previous FDTD simulation works by Zhang et al. ^[ Reference Zhang, Menoni, Gruzdev and Chowdhury ^{28 ]}. We observed that the electron density peak first appeared at the center and then expanded at about 9 fs to the surface, and became stable at about 15 fs. These are the times when the pulse peak intensity reached the surface and when the entire pulse left the surface.

Furthermore, in Figure 6, the first ${\mathrm{HfO}}_2$ layer is more strongly affected by the lower fluence pulse. At higher fluence, there is significant growth of the energy density at the surface, while the lower layers do not show a significant increase. In addition, the energy density at $0.9\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ is the highest among all the fluences, which adds additional evidence that the first ${\mathrm{SiO}}_2$ layer has reflected more injected energy due to the plasma screening effects.

The plasma generation could also imply the breakdown threshold, which is defined as a permanent change to the optical properties. We can see from both the energy and electron density profiles that the local peak in the top ${\mathrm{SiO}}_2$ layer starts to shift at $0.5\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ and the profile is no longer symmetric at the center. This analysis is also consistent with the conclusion predicted by the instability criterion.

5.4 Particle energy

The kinetic nature of the excited electrons is shown in Figure 8. The spectrum of the electrons in the ${\mathrm{SiO}}_2$ layer is wider than that in the ${\mathrm{HfO}}_2$ layer, although both electron density and energy density reach maximum in the first ${\mathrm{HfO}}_2$ layer. For both layers, three snapshots of the energy distribution are shown. At 6 fs, in the early stage of the interaction, there are some excited electrons generated from the ionization. The electrons are highly nonthermalized since the degree $k$ is less than one. At about 10 fs, the peak intensity has passed through the target, leading to average electron energies of 24.56 and 8.61 eV for the ${\mathrm{SiO}}_2$ and ${\mathrm{HfO}}_2$ layers, respectively. The electrons are still nonthermalized as $k$ is about 1.4. At 24 fs, the pulse front has left the target for about 10 fs. The degree $k$ is above 2 and approaching 3, indicating the thermalization process is finishing towards a Maxwell–Boltzmann distribution, indicated by the black curves in Figure 8.

Figure 8 The energy histogram for the electrons in the first ${\mathrm{SiO}}_2$ layer (a) and ${\mathrm{HfO}}_2$ layer (b) at 6, 10, 24 fs for the $0.7\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ simulation with ${\chi}^2$ distribution fitted. The black curve is the Maxwell–Boltzmann distribution given the average energy at the stable stage around 24 fs.

5.5 Particle dynamics

The peak intensity of the $0.7\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ pulse is about ${10}^{14}\ \mathrm{W}\;{\mathrm{cm}}^{-2}$ , giving a theoretical ponderomotive energy^[ Reference Geoffroy, Duchateau, Fedorov, Martin and Guizard ^{96 ]} ${U}_\mathrm{p}$ of approximately 10 and 5.5 eV for ${\mathrm{SiO}}_2$ and ${\mathrm{HfO}}_2$ electrons, respectively (using the ${m}_\mathrm{e}^{\ast }$ values from Table 1). Some studies have reported keV-scale electron spectra under similar intensities but with longer pulse durations^[ Reference Geoffroy, Duchateau, Fedorov, Martin and Guizard ^{96 ]}. This enhancement is often attributed to the surface modification of the target due to ablation, which enhances the local intensity and facilitates stronger electron acceleration. In our previous work, we also observed keV-scale spectra using a 7 fs pulse by introducing pit-like defects on the target surface under comparable conditions^[ Reference Smith, Zhang, Gruzdev and Chowdhury ^{40 ]}.

In gases, electrons with kinetic energies exceeding $10{U}_\mathrm{p}$ are typically observed^[ Reference Paulus, Becker, Nicklich and Walther ^{97 ]}, and in our simulation, the highest electron energy from the first ${\mathrm{SiO}}_2$ layer reaches 150 eV. In addition, a few electrons in the first ${\mathrm{HfO}}_2$ layer achieve energies around 100 eV, as shown in Figure 9 (left). These results align with the kinetic energy predictions for excited electrons based on the Drude model. If we assume the collisional frequency is $1\ {\mathrm{fs}}^{-1}$ , then the energy would be about 118 and 65 eV for ${\mathrm{SiO}}_2$ and ${\mathrm{HfO}}_2$ , respectively, as seen in the studies by Duchateau et al. ^[ Reference Duchateau, Geoffroy, Dyan, Piombini and Guizard ^{98 ]}. Furthermore, the excited ${\mathrm{SiO}}_2$ electrons at the ${\mathrm{HfO}}_2$ – ${\mathrm{SiO}}_2$ interfaces rarely enter the ${\mathrm{HfO}}_2$ layers. There are a limited number of ${\mathrm{HfO}}_2$ electrons penetrating into the adjacent ${\mathrm{SiO}}_2$ layers for a few tens of nanometers.

Figure 9 The ${\mathrm{SiO}}_2$ and ${\mathrm{HfO}}_2$ electron energy spatial distribution for a fluence of $0.7\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ is shown on the left-hand side at 24 fs. The highest energy particles are found in the center of the first layer and at the interfaces between layers. Particle trajectories for random 2% of electrons with low (below average) final energy demonstrate the ionization dynamics in the first ${\mathrm{SiO}}_2$ layer up to 16 fs (after the pulse has left). Ionization near the surface generally occurs at later times, as indicated by the color of the tracks.

The tracks of select electrons in the first ${\mathrm{SiO}}_2$ layer are shown from 6 to 16 fs in Figure 9 (right). Most electrons near the vacuum interface are born after 10 fs, indicated by their yellowish tails. In contrast, electrons near the center of the first ${\mathrm{SiO}}_2$ layer are born earlier. This suggests that the electron density expansion in Figures 5 and 7(b) is due to direct excitation near the surface rather than displacement from the center.