Open Access
Add to BibSonomy
Abstract
Modeling the laser-plasma interaction within solids is crucial in controlling ultrafast laser processing of dielectrics, where the pulse propagation and plasma formation dynamics are highly intricate. This is especially important when dealing with nano-scale plasmas where specific phenomena of plasma physics, such as resonance absorption, can significantly impact the energy deposition process. In this article, we report on adapting of a Particle-In-Cell code, EPOCH, to model the laser-plasma interaction within solids. This is performed by implementing a background permittivity and by developing and validating adapted field ionization and impact ionization modules. They are based on the Keldysh ionization theory and enable the modeling of ionization processes within solids. The implementation of these modules was validated through comparisons with a hydrodynamic code and existing literature. We investigate the necessary number of super-particles per cell to model realistic ionization dynamics. Finally, we apply the code to explore the dynamics of plasma formation in the regime of of quantized structuring of transparent films. Our study elucidates how a stack of nano-plasma layers can be formed by the interference of a pulse with its reflection on the exit surface of a high refractive index material.
© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
At high intensities, ultrashort laser pulses have the ability to generate a high-temperature and dense plasma inside solids. The control of femtosecond laser-induced plasma formation and subsequent energy deposition mechanisms in solid-state matter holds significant importance for various applications, including X-ray or extreme ultraviolet (EUV) secondary sources, dense plasma generation for studying warm dense matter, and micro-/nano-machining of transparent materials.
During the laser pulse, the ionization of the solid material occurs primarily, leading to the generation of a plasma of free electrons and holes within the material. The so-called ionization process corresponds to the promotion of electrons from the valence band to the conduction band. The plasma undergoes dynamic interactions with the laser pulse, resulting in its restructuring through phenomena such as plasma wave phenomena, resonance, surface plasmons, etc. Even after the completion of the laser pulse, the ionization dynamics can persist for several picoseconds due to avalanche and recombination processes [1].
Particle-In-Cell (PIC) codes are valuable tools for simulating the dynamics of laser-plasma interaction. These codes solve Maxwell’s equations and the equation of motion for billions of particles, enabling kinetic modeling of plasma and capturing specific plasma effects like Landau damping, which are challenging to simulate using other types of codes such as hydrodynamic models. Many PIC codes have been developed for simulating plasma physics in the gas state or the interaction of relativistic pulses with plasma surfaces [2–5]. Here, our focus is on modeling ionization and laser-plasma interaction, specifically within the bulk of transparent dielectrics, such as sapphire or fused silica, using ultrashort laser pulses of moderately high intensity ($10^{13}-10^{15}$ W.cm$^{-2}$). Unlike the conditions at the surface of dielectrics, the ionization dynamics, laser absorption, and plasma-induced defocusing significantly impact the propagation of the laser pulse and the resulting energy density distribution. Particularly interesting is the study of nano-plasmas created by ultrashort pulses, where the laser-plasma interaction plays a crucial role [6]. This physics is very rich because high energy density can be reached over nanometric transverse scales. The field-plasma interaction shows electron acceleration in the keV range, and is featured by second harmonic and THz field generation [6–10].
In our previous works using EPOCH [6,8,10], the ionization dynamics was not taken into account and the particles were placed in vacuum. To simulate the plasma formation dynamics, it is necessary to model the laser pulse propagation inside the solid, which includes accounting for dielectric permittivity and the ionization scheme specific to solids. However, EPOCH PIC code uses ionization models well suited for gases, such as Argon, Helium, but has not been developed nor tested for high-density plasmas inside solids. For field ionization, three models are implemented to provide ionization rates as a function of laser intensity: the multiphoton ionization regime, i.e. for Keldysh parameter $\mathrm {\gamma >> 1}$ is covered by the model derived by Delone et al. [11], then the tunneling ADK model derived by Ammosov et al. [12] models field ionization for moderate intensities corresponding to $\mathrm {\gamma << 1}$, and finally, for strong laser intensities exceeding ADK’s rate turning point, the barrier suppression model derived by Posthumus et al. is used [13]. Concerning impact ionization, the native model implemented in EPOCH is based on the MBELL atomic cross-section derived by Haque et al. [14].
Here, we present in section 2. the implementation of a background refractive index in the EPOCH code. PIC codes capture very well the trajectories of free electrons, but, when modeling the pulse interaction inside a solid, it is also necessary to reproduce the effect of bound electrons on electromagnetic fields, which is the dielectric permittivity. It has a significant impact on the wavelength scaling.
In Section 3, we proceed to replace the ionization models with those based on the Keldysh ionization theory of solids [15–18]. This theory provides an approximate description of the complex phenomena of solid-state ionization, encompassing both field-induced and electron impact-induced ionization processes through two distinct ionization rates. Specifically, the ionization processes involve the transfer of electrons from the valence band to the conduction band within the solid. Here, in the framework of the interaction of an ultrafast laser pulse with a wide bandgap dielectric, we consider the conduction band electrons as free particles under the parabolic-band approximation. Conversely, the holes are treated as immobile heavy ions. This is why hereafter, we will consider only plasmas of electrons and ions. In our implementation, the solid is modeled as an ensemble of neutrals immersed in a medium of relative permittivity $\varepsilon _r$ with a density that corresponds to the density of states multiplied by the crystal density. We remark that, in a more general framework, electrons and holes possess effective masses that are highly material-dependent. This can be easily implemented in the framework of particle-in-cell code, by solving the trajectories of electrons and of the holes simultaneously. We also neglect the quantum effects of Umklapp processes and the modification of the band structure. We remark that the latter holds as long as the ultrafast phase change occurs on timescales larger than the pulse duration (This time is on the order of one to several 100 fs [19]).
The validation of our implementation of the ionization model is discussed in Section 4. We investigate the influence of the number of super-particles per cell, demonstrating the necessity of a relatively high number of super-particles. The field ionization and impact ionization models are validated through comparisons with a hydro-code and existing literature. Furthermore, we calibrate the Keldysh impact ionization rate for fused silica using experimental results from the literature, where the damage threshold fluence is measured as a function of pulse duration [20].
Finally, in Section 5, we present an application of our model to the investigation of ionization dynamics in an interference configuration. This configuration has been experimentally demonstrated by Kumar et al. [21], where laser-generated nano-plasma layers induce blistering in silicon nitride.
2. Background permittivity
When modeling laser-induced plasmas in a solid, it is essential to integrate in the equations the medium permittivity $\varepsilon$, due to bound charges, since it affects the wavelength and the associated scalings. Here, the permittivity is chosen to be uniform within the simulation volume. A PIC code solves Maxwell’s equations for the electromagnetic fields together with the trajectories of particles. The laser pulse is injected via an antenna or radiative boundary condition. The fields have two origins: the input laser pulse and the currents and net charge separation due to the particle trajectories. The electric and magnetic fields are advanced in time using the Yee algorithm [4]. The boundaries necessitate a special treatment to avoid unwanted numerical reflectivity when simulating a pulse propagating in open conditions, which is Convolutional Perfectly Matched Layers (CPMLs). Therefore, inserting the information of the material permittivity within a Particle-In-Cell code requires the modification of three main parts of the code, namely the Maxwell solver, the calculation of the field at the injector and the Convolutional Perfectly Matched Layers.
2.1 Maxwell solver
In most PIC codes, only rotational equations are used in the Maxwell solver [4,5]. The permittivity $\varepsilon = \varepsilon _0 \varepsilon _r$ appears in Maxwell-Ampère equation, which reads:
2.2 Injector
In EPOCH, the fields are injected using $B$ fields. We have adapted the radiative boundary condition developed by Ramsay [22] to account for the permittivity in the injector space. For instance, for the injection in a plane $(y,z)$, the equation to update $(B_z)^n_{jkl}$ at the $x=x_{min}$ wall reads:
The permittivity appears in several places of the injector expression.
2.3 Insertion of the background permittivity in convolutional perfectly matched layers
The CPMLs were implemented in EPOCH using the method proposed by Taflove and Hagness [23]. Although this approach is valid for a medium of any index of refraction, it was implemented in EPOCH in the case of vacuum. We provide here the equations for updating the fields adapted to a homogeneous dielectric medium, with relative permeability $\mu _r = 1$ and zero electrical conductivity $\sigma = 0$, as in Ref. [24]. In the computation volume occupied by the CPMLs, the convolutional terms $\mathrm {\Psi _\textbf {E}}$ and $\mathrm {\Psi _\textbf {B}}$ are updated and injected into the fields. In the $\mathrm {y=y_{min}}$ and $\mathrm {y=y_{max}}$ plane of the simulation box, the expressions for $\mathrm {\Psi _{E_{xy}}}$ and $\mathrm {\Psi _{E_{zy}}}$ read:
The $\textbf {E}_x$ field is updated in the CPML volume, using:
Thus, the relative permittivity of the material $\varepsilon _r$ appears only in two coefficients, $\sigma _w'$ and $C_b$. It is therefore sufficient to divide them respectively by $\sqrt {\varepsilon _r}$ and $\varepsilon _r$ in EPOCH to account for the presence of the dielectric medium in the CPMLs boundaries.
2.4 Validation
The implementation validation has been performed via several straightforward simulations whose results are not shown here for the sake of conciseness: verifications of Snell’s law, Fresnel’s refraction and reflection, wave velocity and wavelength.
Here, we make a quantitative verification by simulating the electric field of an s-polarized wave obliquely impinging on a plasma ramp generated inside a medium of permittivity $\varepsilon = \varepsilon _0 \varepsilon _r$. Following the approach of Kruer, the field has the following analytical expression in the plasma ramp [25]:
In Figure 1, we show our numerical results as dots (green for $\varepsilon _r = 1$ and orange for $\varepsilon _r = 4$) and, as solid lines, the analytical expression of Eq. (10) where we highlight that permittivity appears both in the amplitude of the field and in the phase of the Airy function. The numerical values reproduce very well the analytical expectations without any rescaling. Via the interference pattern, this test allows for validating that the wavelength scaling, the field amplitude, the plasma density, and permittivity are together well reproduced. The simulation has been repeated for 5 different angles, ranging from 10 to 30$^\circ$ and different permittivities, ranging from 1 to 4, with each time a result closely following the analytical expectations. Moreover, it shows that the higher the background permittivity, the higher the maximum value reached by the field at resonance.
Fig. 1. Analytical and numerical field amplitude of an s-polarized wave (wavelength 0.8 $\mathrm {\mu}\textrm{m}$) propagating obliquely with an angle $\mathrm {\theta }=23.58^\circ$ with respect to the z-axis in a plasma ramp with density $\mathrm {n_e(z) = n_cz/L}$, where the plasma scale length is $L=10\,\mathrm {\mu}\textrm{m}$ for two different background permittivities: $\mathrm {\varepsilon _r = 1.0}$ in green and $\mathrm {\varepsilon _r = 4.0}$ in orange. The critical density is defined as $\mathrm {n_c = \varepsilon _r\omega _0^2\varepsilon _0m_0/q_0^2}$. It corresponds to a plasma density where the permittivity is zero.
3. Ionization model
3.1 Field ionization
In EPOCH, the calculation of the field ionization rate is based on the amplitude of the electric field. This rate is divided into three cases, each corresponding to a distinct regime: multiphoton absorption, tunneling, and barrier suppression ionization [4]. However, the currently implemented rates are applicable to atomic gases and not suitable for describing field ionization in solids. Accurate modeling of field ionization in solids requires a quantum description of the electron gas, such as time-dependent density functional theory [26,27]. Unfortunately, the computational cost associated with such quantum approaches is prohibitively high, especially for simulating pulses in the range of 100 fs propagating over a few tens of micrometers [6]. This computational limitation hinders the investigation of laser-plasma couplings and spatio-temporal pulse propagation in realistic scenarios.
That is why we have adopted Keldysh’s field ionization rate [15–17]. This rate is applicable to the range of electric field amplitudes typically encountered in the targeted modeling of laser-solid interaction at moderate intensities (between $10^9$ and $10^{11}$ $\mathrm {V.m^{-1}}$). Notably, Keldysh’s field ionization rate provides an estimation of the ionization rate with a relatively low computational cost [16]. Its expression is shown in Eq. (11) [28].
3.2 Collisional ionization
After being promoted to the conduction band, electrons are accelerated by the electric field. Through collisions, they can transfer their energy to valence electrons, promoting them to the conduction band. This process is known as collisional ionization [1].
The model originally implemented in EPOCH for collisional ionization was an empirical model based on an algorithm developed by Sentoku and Kemp [29], which relies on the effective cross-section of atoms. In our work, we have replaced this collisional ionization module with a description based on the Keldysh collisional ionization rate $\alpha ^{Single}$. It describes the ionization rate for a single conduction electron and is defined as follows [15,18]:
In Eq. (14), $\varepsilon _r$ denotes the relative material permittivity, $A_c$ and $A_v$ are the overlap integrals between the initial and final states of the conduction electron and the valence electron promoted to the conduction band [18].
4. Validation of ionization
In this section, we describe our validation procedure. Firstly, we assessed the impact of the number of super-particles per cell on the convergence of the electron density. We then calibrated the impact ionization rate by adjusting the values of the overlap integrals ${A_c}$ and ${A_v}$ using experimental data. The validation process was conducted separately for field and collisional ionization and subsequently for both phenomena together.
4.1 Influence of the number of super-particles per cell
In PIC codes, ionization of a super-particle results in generating a super-electron and a super-ion with the same number of real particles as the ionized super-neutral [30]. The choice of the number of super-particles per cell ($\mathrm {PPC}$) is a crucial parameter that determines the number of ionized particles during each ionization event. It directly affects the smoothness of the density evolution and has a significant impact on the simulation results. Increasing the number of super-particles per cell leads to a more discretized density evolution and higher accuracy.
Several identical 3D simulations were performed to investigate the impact of the number of super-particles per cell ($\mathrm {PPC}$) on the convergence of the solution. In these simulations, a pixel with a high density of neutral atoms ($2.2\times 10^{28}$ $\mathrm {m^{-3}}$), placed in the center of the simulation box, is excited by a plane wave. The simulation parameters are given in Table 1. This setup ensured that the ionization dynamics is decoupled from the spatial pulse propagation dynamics, enabling a comparison with the analytical formula. Figure 2 displays the results of these simulations.
Fig. 2. Evolution of the number of free-electrons generated by ionization as a function of time for different initial numbers of super-particles per cell. The black curve represents the theoretical integration of Keldysh’s ionization rate (11). The simulation parameters are provided in Table 1, and $\mathrm {t=0}$ corresponds to the time when the maximum pulse intensity reaches the pixel of neutral atoms.
Table 1. Numerical parameters used in the 3D simulations for the study of the number of super-particles per cell.
The figure clearly demonstrates that the choice of $\mathrm {PPC}$ significantly impacts the convergence of the solution. Convergence close to the theoretical number of electrons predicted by the Keldysh ionization theory (black curve) is achieved for a number of super-particles per cell equal to 256 or above and tends to reach the theoretical value as $\mathrm {PPC}$ increases. We remark that this value is also problem-dependent, depending on the steepness of the fields and density distributions. Smoother plasma distributions require less particle-per-cell.
4.2 Field ionization
To validate our implementation of the Keldysh field ionization rate, we compared our results with a Hydro code developed within our team [31,32]. We considered a 1D simulation performed by Petrov and Davis [33], where a laser pulse interacts with a dielectric solid after propagating a distance of $\mathrm {L_{vac} = c\tau _0}$ in vacuum. Figure 3 illustrates the configuration of the simulation box for this 1D setup. The parameters used in our simulations can be found in Table 2.
Fig. 3. Illustration of the simulation box at $\mathrm {t=0}$. The front of the laser pulse is positioned at the vacuum-dielectric interface, denoted as $\mathrm {x=L_{vac}}$, while the rear of the pulse is located at $\mathrm {x=0}$. The laser pulse propagates from left to right along the x-axis.
Table 2. Field ionization and collisional ionization parameters used in 1D simulations corresponding to sapphire ($\mathrm {Al_2O_3}$).
In this work, we have chosen to fix the time $\mathrm {t=0}$ when the front of the laser pulse reaches the vacuum-dielectric interface. The profile of the laser pulse is defined as follows (Eq. (16)):
The comparison between the Hydro code and the PIC code for field ionization dynamics is shown in Fig. 4. The results obtained from both codes are in good agreement and consistent with the orders of magnitude found in the literature [33–36]. With this validation, we can proceed to verify the implementation of impact ionization.
Fig. 4. Evolution of the electron density per unit area $\mathrm {\bar {n}_e(t) = \int n_e(x,t)dx}$ as a function of time for both EPOCH and Hydro codes for peak laser intensity $\mathrm {I_0 = 3\times 10^{13} \ W.cm^{-2}}$. Time $\mathrm {t=0}$ corresponds to the time when the peak of the pulse reaches $\mathrm {x=L_{vac}}$.
4.3 Impact ionization
4.3.1 Impact ionization calibration
The Keldysh impact ionization rate, given by Eq. (13), depends on various material properties such as electron/ion masses, bandgap, and the overlap integrals $\mathrm {A_c}$ and $\mathrm {A_v}$. However, calculating these overlap integrals can be challenging, as different methods can yield significantly different values (ranging from 0.01 to 0.38 for GaAs, for example) [18]. Therefore, in order to accurately simulate impact ionization, calibration of these overlap integrals is necessary.
To perform the calibration, we followed the protocol described by Déziel et al. [20]. This protocol involves determining the minimum fluence required to reach material damage, i.e. assumed to be reached at the critical density $\mathrm {n_c}$ as a function of the laser pulse duration $\tau$. According to numerical and experimental results for silicon dioxide [28,37–42], the threshold fluence follows a power law behavior, $\mathrm {F_{th}}\propto \tau ^{0.3}$, as shown in Fig. 5.
Fig. 5. Experimental and numerical evolution of the fluence threshold as a function of the pulse duration $\tau$ for fused silica (SiO2). The experimental data sets are from Refs. [28,37–42]
Table 3. Numerical parameters used in the impact ionization calibration simulations. $\mathrm {I_0}$ has been chosen for each simulation in order to fit the theoretical fluence $\mathrm {F_{th}}\propto \tau ^{0.3}$ and $\mathrm {t_0}$ is different for each value of $\tau$.
A large set of simulations (parameters provided in Table 3) were performed for different combinations of $\mathrm {A_c}$ and $\mathrm {A_v}$ values, with varying pulse durations $\tau$. For each parameter set, we varied the pulse fluence to determine the threshold at which the laser-generated plasma reaches the critical density. We plot in Fig. 5 as a black dashed line with star markers, the evolution of the threshold fluence as a function of pulse duration for $\mathrm {A_c = A_v = 0.104}$. We observe that it is in reasonable agreement with the numerical result of Déziel et al. (solid blue line) and the different experimental results of the literature. This approach could be repeated for other materials.
The numerical results obtained in EPOCH exhibit good agreement with experimental data over a significant range of pulse durations. Our implemented model in EPOCH is valid over a relatively large area of interest for the modeling of ultrashort laser pulse-solid interaction. With the successful calibration of impact ionization, the next step is to perform validation to ensure the correct implementation of our model.
4.3.2 Validation of the impact ionization
To validate the implementation of collisional ionization, we performed 1D simulations in a dielectric, similar to the field ionization validation using our modified version of EPOCH and the hydro-code mentioned previously. The simulations consider the presence of an electron-ion pre-plasma. To focus solely on the impact ionization dynamics, we used heavy ions ($\mathrm {m_i = 102\times 1836m_0}$) that transport is negligible on the simulation timescale. The initial electron and ion density distributions in the pre-plasma were defined using a centered hypergaussian function : $\mathrm {n_{e,i}} = \mathrm {n_c e^{-\left (x/w_r\right )^4}}$ and the initial neutral density profile is $\mathrm {n_0 - n_{e,i}}$ with $\mathrm {n_0=12\,n_c}$. Here, $\mathrm {n_c}$ represents the critical density, and $\mathrm {w}_\textrm{r} = 0.35\,\mathrm{\mu}\textrm{m}$ is the characteristic width. The simulation box is $\mathrm {4}\,\mathrm{\mu}\textrm{m}$ long with a discretization of 8 nm, the number of particles per cell is 512, the background permittivity is $\mathrm {\varepsilon _r=1}$ (vacuum) and the overlap integrals are $\mathrm {A_c=A_v=0.164}$.
The comparison between the results obtained from the two codes is presented in Fig. 6. The figure shows the evolution of the electron density as a function of the initial electron/ions temperature. This evolution is represented by the ratio $\mathrm {n_f/n_i}$, which indicates the ratio of the final electron number relative to the initial one. This ratio quantifies the level of ionization that occurred during the simulation. The results demonstrate that for a wide range of initial temperatures, the final electron density can exceed several times the critical density.
Fig. 6. Ratio $\mathrm {n_f/n_i}$ of the final electron density to the initial electron density as a function of the initial electron/ion temperature for both EPOCH and Hydro codes. The ratio reflects the extent of collisional ionization that occurred during the simulation.
Similar to the field ionization case, the results obtained from both codes show a close agreement, indicating that our implementation of the Keldysh impact ionization rate has been successfully carried out.
4.4 Field and collisional ionization
Continuing from the previous results, we will now investigate the evolution of the electron density by conducting 1D simulations that incorporate both field ionization and collisional ionization. These simulations are similar to those used to validate field ionization (parameters can be found in Table 2) and now include the effects of both ionization mechanisms.
As depicted in Fig. 7, the results obtained from both EPOCH and the Hydro code show a close agreement. Additionally, we have included the results obtained by Petrov and Davis [33] in the graph, which exhibit a good agreement with our simulations.
Fig. 7. Evolution of the electron density per unit area $\mathrm {\bar {n}_e(t) = \int n_e(x,t)dx}$ as a function of time for EPOCH, the Hydro code, and Petrov and Davis [33] simulations for a peak laser intensity of $\mathrm {I_0 = 3\times 10^{13} \ W.cm^{-2}}$. The simulations include both field and collisional ionization.
5. Application to quantized structuring of transparent films
In this section, we aim to numerically reproduce the experimental findings of Kumar et al. [21] on femtosecond laser-induced blistering by applying the model described above. Their experiment consisted in irradiating a thin silicon nitride ($\mathrm {Si_3N_4}$) film with a femtosecond laser pulse. The film has a thickness of $\mathrm {d=500}$ nm and a refractive index of $\mathrm {n_{film} = 1.98}$. The interaction of the laser pulse with the film exhibited Fabry-Perot interferences that led to the formation of localized plasma layers spaced by $\lambda /2$ inside the film. The subsequent plasma relaxation resulted in the creation of layered cavities and removal of nanometric layers of material [21]. Our goal is to use our simulation tool to reproduce the experimental observations.
Our simulation specifically aims to investigate plasma formation within the interference fringes. It can be expected that once the plasma layer in the fringe that is the closest to the laser source reaches the critical density, it should no longer be transparent to the field and would prevent the further increase of plasma density in the other fringes. This expectation raises the question of the dynamics of the plasma formation in such a situation, and how uniform the plasma layers can be formed.
The laser pulse used in the simulation has an intensity profile of 200 fs duration FWHM, and a central wavelength of $\lambda _0 =$ 522 nm in vacuum ($\lambda = \lambda _0/n_{\mathrm {film}}$ inside $\mathrm {Si_3N_4}$). The linearly polarized pulse has a peak intensity of $I_0 = 9\times 10^{12}$ W/cm$^2$. To mimic the Fresnel reflection of the pulse at the sample exit side, we used a reflective boundary condition at the end of the simulation box (we approximated the reflection coefficient to 1). In our simulation, we used a spatial resolution $\mathrm {d}x=0.39$ nm, i.e., 3 Debye lengths, 1024 super-particles per cell, and the bandgap of the material was set to 5.3 eV. (Note that in a 3D geometry, numerical noise is lower due to the higher number of degrees of freedom of the particle [43]. Hence, a typical resolution in 3D could be $\mathrm {d}x= 10$ nm, and $\sim 64$ super-particles per cell). The overlap integrals were set to $A_c = A_v = 0.14$ (chosen in absence of precise values in the literature).
Figure 8 illustrates the temporal evolution of the field amplitude (Fig. 8(a)) and the electron density (Fig. 8(b)) within the simulation box. In our setup, the laser pulse propagates from the $\mathrm {x=-2}\,\mathrm{\mu}\textrm{m}$ wall and is reflected back by the boundary at $\mathrm {x}=500$ nm. Neutrals are present only between $\mathrm {x=0}$ and $\mathrm {x=500}$ nm to simulate the thin silicon nitride film. From $t\sim -300$ fs to $t\sim -25$ fs, the electric field exhibits a pronounced increase in amplitude within four lobes, which are separated by a distance of $\lambda /2$. These lobes correspond to the constructive interference regions of the Fabry-Perot cavity [21]. Once the field amplitude exceeds a strength of approximately $1\times 10^{10}$ V/m, efficient ionization of the dielectric material occurs. Consequently, the plasma density rapidly increases and reaches approximately half of the critical density, defined as $\mathrm {n_c = n_{film}^2\omega _0^2\varepsilon _0m_0/q_0^2}$ for which the permittivity is zero (in Kumar et al. [21], the critical density is calculated without the $\mathrm {n_{film}^2}$ factor, which is why the density is thought over-critical).
Fig. 8. Evolution of a) field amplitude and b) electron density as a function of time along the depth of the silicon nitride layer. The critical density is $\mathrm {n_c = 1.6\times 10^{22} cm^{-3}}$ and $\mathrm {t=0}$ corresponds to the time when the maximum pulse intensity enters the simulation box. The injector is positioned at $\mathrm {x=-2}\,\mathrm{\mu}\textrm{m}$, no neutrals are present between the injector and the plane $x=0$.
In Fig. 8(a), we observe a strong absorption of the electric field around $t\simeq 0$, which is associated to the rapid rise in plasma density. Between $t=0$ and $t\simeq 50$ fs, the electric field pattern undergoes significant modifications: the field interference lobes are shifted by values between 20 and 80 nm along the propagation direction, due to the reduced refractive index caused by the presence of plasma. As time progresses further ($t\simeq 100$ fs), the electric field in the lobe situated at approximately $x=100$ nm regains a high amplitude of $1\times 10^{10}$ V/m. The electric field tunnels through the plasma layers, which are approximately 50 nm thick. We see however a significant damping of the interference lobes in the electric field pattern (a).
Due to the shift of the first interference lobe of the electric field to a position in between the first two plasma layers, the second plasma layer expands towards the left. In contrast, the two last plasma layers on the right, at $x\sim 300$ and $x\sim 450$ nm maintain a nearly constant width because of the strong attenuation of the electric field. The slight density increase in these two plasma layers is due to collisional ionization. Overall, the plasma density does not exceed the critical density during the pulse illumination itself, but might exceed it over picosecond timescales, after collisional ionization has taken place (note this does not change the overall energy density deposited by the laser pulse).
From the plasma distribution at the end of our simulation ($t=350$ fs) we expect the material to finally be ablated over the first $\sim$200 nm, where the plasma density is continuous. Depending on mechanical strength, another material layer can be ablated, corresponding to the third nanoplasma layer. Interestingly, we have numerically observed that modifying the pulse duration strongly impacts the density distribution. With a shorter pulse duration, the nano-plasma layers can be formed in a more periodic pattern, and reach a similar density.
Some numerical noise is present in our simulation because the laser field intensity is not very high ($\mathrm {10^{10}\,V.m^{-1}}$). However, we have maintained the magnitude of the numerical noise at an order of magnitude lower than the laser field by carefully tuning the number of particles per cell and grid resolution. Specifically, we employed a grid resolution three times the Debye length to mitigate the magnitude of the collisional electrostatic field. Numerical self-heating is very low at $\frac {dT}{dt}\simeq 0.02$ eV.ps$^{-1}$. The impact on the plasma density is negligible since the numerically generated field is on the order of $10^8$ to $\mathrm {10^9 \ V.m^{-1}}$, corresponding to a field ionization rate lower by at least more than 6 orders of magnitude. We have also verified the conservation of energy within less than 1%.
Overall, our simulation code effectively captures the interplay between the dynamics of nano-plasma layer formation and the reshaping of the laser field, offering valuable insights to complement the findings of Kumar et al. [21].
6. Conclusion
In conclusion, we have modified the EPOCH code by adding a background permittivity and by implementing the Keldysh ionization rates for both field-induced and collision-induced ionization in our PIC code, enabling us to model the dynamics of solid-density plasmas. The validation of our implementations involved comparisons with a previously developed hydrodynamic code, as well as with results from Petrov and Davis [33]. Additionally, a numerical study examining the influence of the number of super-particles per cell confirmed the need for a minimum of 256 super-particles per cell for convergence in the physical situations of relatively dense plasma generation by laser pulses.
It is still important to highlight that describing the ionization processes with Keldysh model is still a relatively strong approximation of the physics at play, which would require a quantum description of the evolution of the electron population within the band structure. However, the simultaneous combination of quantum and particle descriptions is computationally too heavy for fair comparison with experiments (pulse durations exceeding 100 fs, plasma lengths up to several tens of micrometers, etc) [27].
Using our modified PIC code, we investigated the evolution of the electron density in the experimental setup conducted by Kumar et al. [21]. Our simulations revealed a dynamic evolution of the free-electron plasma density coupled to the reshaping of the laser field pattern, providing valuable insights into the experimental observations. Further work will include the modeling of holes, instead of fixed ions, using the same approach. This will enable simulations of ultrafast laser pulses interaction with materials with lower bandgap and holes small effective masses [32]. Future directions could also encompass the inclusion of more refined electron dynamics in the band structure with normal and Umklapp scattering processes [44].
Funding
HORIZON EUROPE European Research Council (682032-PULSAR); Agence Nationale de la Recherche (ANR-17-EURE-0002).
Acknowledgments
Fruitful discussions with L. Froehly are gratefully acknowledged. We acknowledge the financial support of: European Research Council (ERC) 682032-PULSAR and EIPHI Graduate School ANR-17-EURE-0002. This work was granted access to the PRACE HPC resources Joliot-Curie Rome at TGCC, CEA, France under the Project "PULSARPIC" (RA5614), HPC resource Joliot-Curie Rome/SKL/KNL at TGCC, CEA, France under the projects A0070511001, A0090511001, A0130511001 and Mésocentre de Calcul de Franche-Comté.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. B. Rethfeld, D. S. Ivanov, M. E. Garcia, et al., “Modelling ultrafast laser ablation,” J. Phys. D: Appl. Phys. 50(19), 193001 (2017). [CrossRef]
2. R. Fonseca, L. Silva, F. Tsung, et al., “OSIRIS: A three-dimensional, fully relativistic particle in cell code for modeling plasma based accelerators,” (2002), pp., 342–351.
3. S. Markidis, G. Lapenta, and Rizwan-uddin, “Multi-scale simulations of plasma with iPIC3D,” Math. Comput. Simul. 80(7), 1509–1519 (2010). [CrossRef]
4. T. D. Arber, K. Bennett, C. S. Brady, et al., “Contemporary particle-in-cell approach to laser-plasma modelling,” Plasma Phys. Controlled Fusion 57(11), 113001 (2015). [CrossRef]
5. J. Derouillat, A. Beck, F. Pérez, et al., “Smilei : A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation,” Comput. Phys. Commun. 222, 351–373 (2018). [CrossRef]
6. K. Ardaneh, R. Meyer, M. Hassan, et al., “Femtosecond laser-induced sub-wavelength plasma inside dielectrics: I. Field enhancement,” Phys. Plasmas 29(7), 072715 (2022). [CrossRef]
7. J. del Hoyo, R. Meyer, L. Furfaro, et al., “Nanoscale confinement of energy deposition in glass by double ultrafast Bessel pulses,” Nanophotonics 10(3), 1089–1097 (2021). [CrossRef]
8. K. Ardaneh, M. Hassan, B. Morel, et al., “Femtosecond laser-induced sub-wavelength plasma inside dielectrics. II. Second-harmonic generation,” Phys. Plasmas 29(7), 072716 (2022). [CrossRef]
9. K. Ardaneh, R. Giust, P.-J. Charpin, et al., “Electron heating and radiation in high aspect ratio sub-micron plasma generated by an ultrafast Bessel pulse within a solid dielectric,” Eur. Phys. J. Spec. Top. 232(13), 2247–2252 (2023). [CrossRef]
10. K. Ardaneh, K.-I. Nishikawa, R. Giust, et al., “Femtosecond laser-induced sub-wavelength plasma inside dielectrics. III. Terahertz radiation emission,” Phys. Plasmas 30(1), 013301 (2023). [CrossRef]
11. N. B. Delone and V. P. Krainov, Direct (Nonresonant) Multiphoton Ionization of Atoms (Springer Berlin Heidelberg, Berlin, Heidelberg, 1994), pp. 81–117.
12. M. V. Ammosov, N. B. Delone, and V. P. Krainov, “Tunnel ionization of complex atoms and atomic ions in electromagnetic field,” in High Intensity Laser Processes, vol. 0664 J. A. Alcock, ed., International Society for Optics and Photonics (SPIE, 1986), pp. 138–141.
13. J. Posthumus, M. Thompson, L. Frasinski, et al., “Molecular dissociative ionisation using a classical over-the-barrier approach,” in Multiphoton Processes 1996, P. LambropoulosH. Walther, eds. (IOP Publishing Ltd, 1997), pp. 298–307.
14. A. K. F. Haque, M. A. Uddin, A. K. Basak, et al., “Electron impact ionization of M-shell atoms,” Phys. Scr. 74(3), 377–383 (2006). [CrossRef]
15. L. V. Keldysh, “Ionization in the field of a strong electromagnetic wave,” Zh. Eksperim. i Teor. Fiz. 47, 1 (1964).
16. T. J.-Y. Derrien, N. Tancogne-Dejean, V. P. Zhukov, et al., “Photoionization and transient Wannier-Stark ladder in silicon: First-principles simulations versus Keldysh theory,” Phys. Rev. B 104(24), L241201 (2021). [CrossRef]
17. L. Sudrie, A. Couairon, M. Franco, et al., “Femtosecond laser-induced damage and filamentary propagation in fused silica,” Phys. Rev. Lett. 89(18), 186601 (2002). [CrossRef]
18. B. K. Ridley, Quantum Processes in Semiconductors (Oxford University, 2013).
19. K. H. Bennemann, “Photoinduced phase transitions,” J. Phys.: Condens. Matter 23, 073202 (2011).
20. J.-L. Déziel, L. J. Dubé, and C. Varin, “Dynamical rate equation model for femtosecond laser-induced breakdown in dielectrics,” Phys. Rev. B 104, 1 (2021). [CrossRef]
21. K. Kumar, K. Lee, J. Li, et al., “Quantized structuring of transparent films with femtosecond laser interference,” Light: Sci. Appl. 3(3), e157 (2014). [CrossRef]
22. M. Ramsay, “Short-pulse laser interactions with high density plasma,” Ph.D. thesis, University of Warwick (2015).
23. A. Taflove and S. C. Hagness, Computational electrodynamics: the finite-difference time-domain method (Artech House, Norwood, 2005), 3rd ed.
24. K. P. Prokopidis and D. C. Zografopoulos, “Modeling plasmonic structures using LOD-FDTD methods with accurate dispersion models of metals at optical wavelengths,” J. Lightwave Technol. 35(2), 193–200 (2017). [CrossRef]
25. W. Kruer, The Physics of Laser Plasma Interactions (1st ed.), (CRC University, 2003).
26. G. Duchateau, B. Chimier, S. Coudert, et al., “Evidence of noncollisional femtosecond laser energy deposition in dielectric materials,” Phys. Rev. B 102(2), 024305 (2020). [CrossRef]
27. K. Yabana, T. Sugiyama, Y. Shinohara, et al., “Time-dependent density functional theory for strong electromagnetic fields in crystalline solids,” Phys. Rev. B 85(4), 045134 (2012). [CrossRef]
28. B. Chimier, O. Utéza, N. Sanner, et al., “Damage and ablation thresholds of fused-silica in femtosecond regime,” Phys. Rev. B 84(9), 094104 (2011). [CrossRef]
29. Y. Sentoku and A. J. Kemp, “Numerical methods for particle simulations at extreme densities and temperatures: Weighted particles, relativistic collisions and reduced currents,” J. Comput. Phys. 227(14), 6846–6861 (2008). [CrossRef]
30. C. Birdsall, “Particle-in-cell charged-particle simulations, plus monte carlo collisions with neutral atoms, pic-mcc,” IEEE Trans. Plasma Sci. 19(2), 65–85 (1991). [CrossRef]
31. B. Morel, R. Giust, K. Ardaneh, et al., “A simple solver for the two-fluid plasma model based on pseudospectral time-domain algorithm,” Commun. Comput. Phys. 29(3), 955–978 (2021). [CrossRef]
32. B. Morel, R. Giust, K. Ardaneh, et al., “Two-fluid plasma model for ultrashort laser-induced electron-hole nanoplasmas,” Phys. Rev. B 106(3), 035207 (2022). [CrossRef]
33. G. M. Petrov and J. Davis, “Interaction of intense ultra-short laser pulses with dielectrics,” J. Phys. B: At., Mol. Opt. Phys. 41(2), 025601 (2008). [CrossRef]
34. B. Rethfeld, O. Brenk, N. Medvedev, et al., “Interaction of dielectrics with femtosecond laser pulses: application of kinetic approach and multiple rate equation,” Appl. Phys. A 101(1), 19–25 (2010). [CrossRef]
35. X. Jing, Y. Tian, J. Zhang, et al., “Modeling validity of femtosecond laser breakdown in wide bandgap dielectrics,” Appl. Surf. Sci. 258(10), 4741–4749 (2012). [CrossRef]
36. B. C. Stuart, M. D. Feit, S. Herman, et al., “Optical ablation by high-power short-pulse lasers,” J. Opt. Soc. Am. B 13(2), 459–468 (1996). [CrossRef]
37. M. Lebugle, N. Sanner, N. Varkentina, et al., “Dynamics of femtosecond laser absorption of fused silica in the ablation regime,” J. Appl. Phys. 116(6), 063105 (2014). [CrossRef]
38. M. Mero, J. Liu, W. Rudolph, et al., “Scaling laws of femtosecond laser pulse induced breakdown in oxide films,” Phys. Rev. B 71(11), 115109 (2005). [CrossRef]
39. T. Q. Jia, Z. Z. Xu, X. X. Li, et al., “Microscopic mechanisms of ablation and micromachining of dielectrics by using femtosecond lasers,” Appl. Phys. Lett. 82(24), 4382–4384 (2003). [CrossRef]
40. A.-C. Tien, S. Backus, H. Kapteyn, et al., “Short-pulse laser damage in transparent materials as a function of pulse duration,” Phys. Rev. Lett. 82(19), 3883–3886 (1999). [CrossRef]
41. M. Lenzner, J. Krüger, S. Sartania, et al., “Femtosecond optical breakdown in dielectrics,” Phys. Rev. Lett. 80(18), 4076–4079 (1998). [CrossRef]
42. H. Varel, D. Ashkenasi, A. Rosenfeld, et al., “Laser-induced damage in SiO2 and CaF2 with picosecond and femtosecond laser pulses,” Appl. Phys. A 62(3), 293–294 (1996). [CrossRef]
43. J. M. Dawson, “Particle simulation of plasmas,” Rev. Mod. Phys. 55(2), 403–447 (1983). [CrossRef]
44. S. Coudert, S. Dilhaire, P. Lalanne, et al., “Unified theoretical description of out-of-equilibrium electron intraband dynamics in gold induced by femtosecond laser pulses,” Phys. Rev. B 106(13), 134308 (2022). [CrossRef]
