The 2023 Nobel Prize in Physics awarded to Pierre Agostini, Ferenc Krausz and Anne L’Huillier ‘for experimental methods that generate attosecond pulses of light for the study of electron dynamics in matter’^[1^] and the 2018 National Academies of Science report^[2^] that spurred the ‘Brightest Light Initiative’^[3^] highlight the impactful science enabled by producing the shortest possible pulses of light^[4^–6^]. Recent advances in high-energy few-cycle pulse generation techniques^[7^,8^] allow us to probe new physical effects, including those from the carrier-envelope phase^[9^–11^], and to generate high laser intensities with moderate laser energies for applications such as relativistic plasma-based attosecond pulse generation^[12^]. These developments motivate the design of optical components with higher laser-induced damage thresholds (LIDTs) for few- to single-cycle pulses. The scaling of LIDTs with laser fluence is generally well understood for the tens of picosecond to nanosecond regimes^[13^], but is more complex for shorter pulses^[14^], especially few-cycle pulses^[15^,16^].

Previous experimental work has explored the damage and ablation of ${\mathrm{SiO}}_2$ using few-cycle approximately 800 nm wavelength pulses with durations of 5 fs^[15^,17^,18^] and 7 fs^[16^,19^–21^]. Recently Kafka et al.^[18^] determined the damage threshold for different optical components with a 5 fs pulse and Talisa et al.^[22^] explored the damage and ablation of ${\mathrm{SiO}}_2$ / ${\mathrm{HfO}}_2$ multi-layer coatings with 8 fs pulses. As highlighted in the review by Nagy et al.^[8^], post-compression techniques allow researchers around the world (e.g., Refs. [12,23–25]) to generate few-cycle pulses, but traditional optics are easily damaged by these compressed pulses.

Understanding few-cycle pulse material interactions requires advances in computational and theoretical models. There has been significant work using one-dimensional finite-difference-time-domain (FDTD) simulations to model these interactions^[26^], with recent efforts in two dimensions^[27^,28^]. FDTD simulations provide insight into the laser material interactions and can predict LIDTs, but do not capture the nonthermal nature of excited electrons, nor do they capture their ballistic motion. To expand our understanding of the interaction dynamics, we use particle-in-cell (PIC) simulations^[29^,30^]. As with FDTD simulations, PIC codes solve Maxwell’s equations on a computational grid. In addition, PIC simulations statistically represent the neutral and excited particles in the simulation with a finite number of ‘macroparticles’^[29^–31^]. The charged particles are then advanced using the Lorentz force and additional physical effects, including ionization and collisions, are often added.

PIC simulations are commonly used to study laser–plasma interactions (e.g., by Ziegler et al.^[32^]) and are increasingly being modified to simulate laser damage and related regimes. For example, Mitchell et al.^[33^] modeled crater formation in metals due to femtosecond laser ablation, Cochran et al.^[34^] modeled liquid-crystal plasma mirrors^[35^], Déziel et al.^[36^] studied laser-induced periodic surface structures and Ding et al.^[37^] and Do et al.^[38^] modeled plasmons using PIC simulations. Interactions of lasers with nano/micro-scale structures in silicon and ${\mathrm{SiO}}_2$ have been modeled with versions of the framework introduced in this work^[39^,40^] and recently Charpin et al.^[41^] developed a similar framework to explore ionization in dielectrics. Our work builds on these efforts specifically for laser damage for few-cycle pulses by including the corrected Keldysh ionization model^[28^,42^,43^] to account for ionization across a range of fluences. We are able to explore the nonthermal distribution of excited particles in the laser damage regime using PIC simulation, which improves our fundamental understanding of damage mechanisms.

In this work, we begin in Section 2 by introducing the modifications and extensions we have made to a PIC code to model laser–dielectric interactions. Then in Section 3 we discuss the material properties and simulation setup for both bulk and multi-layer targets. Next, in Section 4 we compare the predictions of this framework to existing experimental results and discuss expected damage threshold metrics. Then we apply the framework to the modeling of multi-layer mirrors in Section 5 and conclude in Section 6.

Our work extends the three-dimensional (3D) implementation of version 4.17.10 of the EPOCH^[31^] PIC code, which is designed for the study of high-energy-density physics. EPOCH is a popular open-source PIC code that can scale to run on thousands of central processing unit (CPU) cores, although we note that there are a variety of other open-source and proprietary PIC codes available with different features and implementations^[44^], such as graphics processing unit (GPU) operation (e.g., PIConGPU^[45^] and WarpX^[46^]). PIC simulations of laser–matter interactions are typically used to study only tens-to-hundreds of femtosecond timescales due to numerical instability issues and computational cost. As such, we focus on the initial laser–matter interaction and the resulting excited electron dynamics. For long-term dynamics and equilibration, one could use a final simulation state from a PIC model as the initial conditions for another model. For example, one could use the electron and ion temperatures in a two-temperature model, or consider tabulated equation of state values^[16^].

EPOCH includes multiple physics modules, but was not designed with laser–dielectric interactions in mind. Towards this goal, we have added a new photoionization model relevant for dielectric materials and a model for optical material properties, as discussed in the next sections.

Our work extends the existing ionization framework already available in EPOCH^[31^] to include the photoionization model developed by 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^[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^[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. [48,49] 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:

We may then write the ionization rate $W$ as follows:where $\kappa$ and $\epsilon$ are complete elliptic integrals of the first and second kind, $\Phi$ is the Dawson integral and

We note that this is the corrected version of the formulation, where the original contains a misprint as noted by Gruzdev^[43^]. Using the uncorrected version can result in significantly different calculations^[28^,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.^[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^[50^].

We use the Pérez/Nanbu^[51^,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: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^[31^]. EPOCH extends the model to low temperatures following an approach by Pérez et al.^[51^] and Lee and More^[53^]. More details can be found in the work of Arber et al.^[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^[31^,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^[55^], although some of the input parameters to the models are not well reported and may require fitting of output results to experimental data^[41^,56^]. Impact ionization is reduced for shorter few-cycle pulses. Models by 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^[58^] predicts this threshold to be $10\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ ^[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.

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^[46^] and the modification to EPOCH by Charpin et al.^[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^[36^,59^], although those effects are not considered here.

We begin by applying this simulation framework to a slab of fused silica ( ${\mathrm{SiO}}_2$ ) corresponding to the experiment discussed by Chimier et al.^[16^] and then apply the framework to a multi-layer mirror, as shown in Figure 1. We test a range of fluences near the reported damage threshold, where the fluence is given by $F=2{E}_\mathrm{las}/{\pi \omega}_0^2$ , with ${E}_\mathrm{las}$ being the energy of the laser pulse. The laser is introduced as a boundary condition and enters into vacuum before interacting with a slab of fused silica at normal incidence. For both the bulk and multi-layer targets, a $\lambda =800$ nm, ${\tau}_\mathrm{FWHM}=7 fs$ pulse duration, sine-squared pulse with a spot radius of ${\omega}_0=4.65\ \unicode{x3bc}$ m is modeled.

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.^[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^[60^,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^[31^].

We benchmark our framework against the experiment by Chimier et al.^[16^,21^], which finds damage with a fluence of $1.18\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ and ablation at $1.3\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ for a 7 fs FWHM pulse at normal incidence. There is some uncertainty in these thresholds. Other experiments of LIDTs for bulk silica with different experimental conditions including a shorter 5 fs pulse^[15^,18^] report thresholds from 1.5 to $1.8\ \mathrm{J}\;{\mathrm{cm}}^{-2}$ . We use a 4.65 μm spot radius in our simulations to match the experimental conditions in the work of Chimier et al.^[16^]. We note that experiments by Lenzner et al.^[15^] and Kafka et al.^[18^] used larger spot sizes, although petawatt-class laser systems^[67^] can often achieve spot sizes of the order of a few wavelengths (e.g., see Ref. [68]).

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^[69^] for free electrons or some fraction of total ionization^[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^[16^,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^[29^–31^].

Alternatively, the instability density suggested by 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.

The excited electron energy density criterion has also been suggested to predict damage^[71^,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^[71^,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. [73,75,76]). For comparison, the threshold for high-energy-density physics^[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.^[78^,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^[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.^[81^] reported $(1.177\pm 0.095)\times {10}^7\;\mathrm{J}\;{\mathrm{kg}}^{-1}$ and Khmyrov et al.^[82^] used $0.96\times {10}^7\;\mathrm{J}\;{\mathrm{kg}}^{-1}$ (from Refs. [83,84]). 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.

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.

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^[85^] library: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.

For the multi-layer dielectric mirror, ${\mathrm{HfO}}_2$ is expected to have a lower damage threshold than ${\mathrm{SiO}}_2$ ^[86^,87^] due to the lower bandgap. Therefore, the damage may be initiated in the first ${\mathrm{HfO}}_2$ layer. For example, Talisa et al.^[22^] found the damage threshold for a four-layer ${\mathrm{SiO}}_2$ / ${\mathrm{HfO}}_2$ mirror to be half that of a bulk ${\mathrm{SiO}}_2$ target. Due to the high computational cost of 3D simulations and uncertainty in material properties for ${\mathrm{HfO}}_2$ , we explore a simple mirror with a relatively small spot size to gain a better qualitative understanding of the interaction. Future validations with experiments should be coupled with more accurate material property measurements and LIDT measurements of bulk ${\mathrm{HfO}}_2$ .

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}$ .

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}$ ^[88^–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^[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.

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$ ^[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.

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.^[28^].

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.^[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.

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.

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^[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^[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^[40^].

In gases, electrons with kinetic energies exceeding $10{U}_\mathrm{p}$ are typically observed^[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.^[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.

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.

Understanding few-cycle pulse interactions with dielectric optical components and their LIDT is essential for advancing next-generation laser systems. The use of kinetic simulations in relevant regimes is becoming increasingly popular^[33^,36^,37^,41^]. Kinetic simulations allow us to capture the nonthermal nature of the initial interaction, which is important to accurately model absorption and ionization. We show that both excited electron density and energy density provide insight into the LIDT. Our framework shows good agreement with experimental LIDTs for bulk silica targets. Multi-layer mirror simulations indicate that plasma screening effects can alter the laser interaction and electron energy distribution for high fluences.

In the future, this framework can be applied to a variety of mirror and grating designs^[28^] for both near- and mid-infrared wavelengths^[39^,99^]. In addition, simulations can be inserted into an optimization algorithm to optimize the LIDT or other properties of interest for optical components^[100^]. This framework already provides a deeper qualitative understanding of the dynamics of laser damage, and we show promising quantitative agreement with experiments for few-cycle pulses. Further development of this framework, including impact ionization^[41^,54^], coupled with more precise measurements of material properties for relevant optical coating designs, can allow further validation of the framework across large ranges of laser fluence.