Experimentally validated modeling of the optical energy deposition in highly ionized ambient air by strong femtosecond laser pulses

The deposition of femtosecond laser optical energy in gases leads to the emission of secondary electromagnetic and acoustic radiation. These optoacoustic components have a complex nonlinear dependency on the laser beam characteristics, such as the pulse energy, duration, wavelength and the focusing conditions, as well as on the optical and elastic characteristics of the gaseous medium. The initial interaction times are governed by the high electronic excitation and ionization. These phenomena result in a self-modulation of the laser pulse, significantly affecting the optical energy deposition on the medium. Such complex nonlinear phenomena are very difficult to be studied via analytical equations. To address this, a multiphysics Particle-In-Cell model is applied for the evaluation of the optical energy deposition and plasma generation from tightly focused femtosecond pulses in ambient air. The computational domain of the model is built to describe optical energy deposition in its full spatiotemporal scale. The model is validated by experimental results of the absorbed energy. The agreement between the computational and experimental results provides the basis for the future development of an advanced microstructural Finite Element Method model, which, combined with the Particle-In-Cell model, will have the ability of delivering detailed insights for all the sub-domains and timescales varying from nano- to femto-seconds of the laser-induced breakdown phenomenon.

Energy deposition in ambient air or other gases from ultrashort laser pulses is a highly complex phenomenon that involves a multitude of linear and nonlinear physical processes. Overall, the chain of processes starts with the deposition of optical energy in a confined area of the targeted medium, where the conditions for the generaa e-mail: kkaleris@upnet.gr (corresponding author) tion of free electrons are met. In this region, the optical intensity surpasses a minimum threshold beyond which, the medium becomes opaque for the particular laser wavelength. It is known that the threshold depends on various parameters, such as the duration and wavelength of the laser pulse [16][17][18]. For femtosecond laser pulses with sufficient energy, the dominant mechanisms of free electron generation are tunneling and above threshold ionization. The free electrons interact with the heavy particles, namely the ions, atoms and molecules, in the ionized region. From these interactions, optical radiation is emitted due to radiative recombination of free electrons and de-excitations of bound electrons to less excited states. At the same time, energy is transferred from the electrons to the heavy particles, leading to a local thermalization of the medium. The increased heat in the confined volume causes a rapid expansion that generates pressure at the boundary with the surrounding air. The expansion takes place at supersonic velocities, essentially forming a shock wave and leading to the progressive de-excitation of the volume. When the pressure in the expanded region is equilibrated by the pressure of the surrounding air, the shock wave decouples from the source and propagates through the medium, progressively transforming into a sound wave.
The multitude of underlying processes render the complete and precise modeling of laser breakdown in gases an extremely challenging task. In addition, the various processes evolve at significantly different time scales. It is indicative that, for a laser excitation of a few tens of femtoseconds, the electrons relax within several nanoseconds, while the complete cooldown of the thermalized medium typically lasts for at least a few microseconds, depending on the amount of the deposited energy [18]. This means that the complete phenomenon spans over approximately nine orders of magnitude in the time domain, a fact that strongly hinders the development of unified models describing laser breakdown of the air or other gases. Moreover, under the described conditions, the pulse propagation is governed, among others, by nonlinear phenomena, such as plasma defocusing, self-focusing and self-phase modulation. Plasma defocusing takes place due to the change in the linear refractive index induced by the generated free electrons, while self-focusing and self-phase modulation are related to the optical Kerr effect, namely the quadratic dependence of the refractive index on the electric field of the pulse. Self-focusing causes energy from the pulse's spatial wings to shift towards the pulse core, effectively reducing the beam radius and acting as a positive lens. Self-phase modulation is the timedomain analogous of self-focusing; it causes a nonlinear phase shift that is proportional to the optical intensity, ultimately leading to chirping of the laser pulse.
Several works present models of the individual processes involved [5,17,[19][20][21]. Most of these works are based on complex systems of differential equations that are difficult to solve. Many of them focus on the early stages of the plasma formation, describing the time evolution of the densities and temperatures of electrons and heavy particles. Also, there are works describing laserinduced breakdown by ultrashort laser pulses in the multiphoton regime, where the optical intensities are relatively low and the linear and nonlinear propagation effects are negligible [22]. However, to the best of our knowledge, until now there is no complete and detailed description of LIB and its accompanying radiative phenomena in gases in the full spatiotemporal extent.
This work introduces an alternative computational approach and presents results from the modeling of the optical pulse propagation, focusing and optical energy deposition in the air from femtosecond laser pulses. The simulations are based on a Particle-In-Cell (PIC) approach that allows for the control of multiple parameters of both the optical radiation and the targeted medium in a broad range of values. Moreover, such a model enables the study of the optical energy deposition in spatial and temporal domains that cover the phenomenon in its full extent. The official release version of the PIC code "EPOCH" used here (v4.17.1) considers laser-plasma pulse propagation effects related to the generation and redistribution of free charges and currents in the plasma, while it does not account for nonlinear propagation phenomena related to the modulation of the refractive index, since displacement currents are not included in any release version. Note that, there are in-house modifications of the open-source code including the dielectric nature of the material, such as [23]; however, such modules are not included in the official release versions. Moreover, higher-order ionizations, such as double or triple ionization, are efficiently modeled by the used PIC code, which are very complicated to model via analytical ionization rate equations. The advantage of PIC codes for the simulation of such multiparametric and complex phenomena lies in that they provide multiphysics computational environments where all fundamental equations governing the phenomena under investigation are implemented, and the appropriate solvers are assigned. In this manner, only the parameters of the problem need to be calibrated and the computational domains to be defined.
The PIC code EPOCH [24] is highly suitable as an alternative approach for the description of the interaction of laser pulses with gases, due to the kinetic description of the plasma. In common simplistic fluid approximations of the plasma, the dependent variables are functions of four independent variables (x, y, z, t) by the assumption that the velocity distribution of the species is Maxwellian and, thus, only one temperature is considered in such approximations. Conversely, the kinetic description of the plasma adopted by the EPOCH code is based upon the dynamics of the phase space distribution function of the particles in which, the velocity and the position of the particles are independent variables f n (x, y, z, u x , u y , u z , t) and, therefore, may accurately describe the plasma state under the influence of the laser pulse. These numerical advantages for the modeling and simulation of laser-gas interactions are balanced by the high number of independent variables to be determined. The seven aforementioned independent unknowns demand the solution of the Vlasov, the Maxwell, the charge density and the electric current equations. The solution of this system of equations is computationally demanding, especially when large spatiotemporal domains, such as the LIBS domains, are simulated.
In this work, simulation results for three different laser pulse energies are shown and compared to experimental measurements, showing good agreement. This validates the model and demonstrates that, in contrast to the common practice of using analytical systems of Partial Differential Equations (PDEs), PIC codes can be effectively used for the simulation of multiscale phenomena, such as the laser pulse propagation and energy deposition in LIB, within large spatiotemporal domains. Also, the study shows that such an approach provides very good insight in the spatiotemporal evolution of significant parameters, such as the electron density, electron and ion temperatures, pulse energy and deposited energy, as well as the spatiotemporal evolution of the optical pulse during its propagation and interaction with the medium. It also provides clear indications about the geometry and dimensions of the plasma channel and the time scales of the described phenomena. Therefore, the first step for the development of a complete, modular and/or unified multi-scale computational approach to predict the secondary radiative phenomena following LIB in terms of the characteristics of the optical radiation and the properties of the medium, is accomplished. The basis for further development of an advanced microstructural Finite Element Method (FEM) model, with reference to the presented PIC model, that will deliver detailed results for all the sub-domains and timescales of the LIB problem, is set.

Deposition of optical energy in air
The deposition of optical energy from short or ultrashort laser pulses in air strongly depends from the parameters of the optical radiation, namely the pulse energy, duration, wavelength, and focusing conditions [16,18,20]. This effectively signifies that the breakdown phenomenon and the following secondary light and sound radiation can be controlled, to a certain extent, from the characteristics of the laser pulses. In this study, experimental measurements of the pulse propagation and energy absorption in air breakdown from tightly focused femtosecond laser pulses are presented. Laser pulses of τ FWHM = 23 fs duration, λ L = 805 nm central wavelength and energy up to 2.7 mJ are used for the plasma generation. The PIC code EPOCH [2] is used to simulate the pulse-medium interaction. This alternative numerical approach is chosen due to the numerical capabilities of the PIC codes to describe the temperature and density evolution of the plasma and the ionization processes of the neutral medium. Additionally, the Maxwell description of the laser pulse inside a plasma channel accounts for linear and nonlinear pulse propagation phenomena and, particularly, Gaussian focusing and nonlinear plasma defocusing due to spatially and temporally resolved generation and redistribution of charges and current densities, which, in turn, influence the Electric and Magnetic fields of the optical pulse and, thus also, the focusing conditions. Nevertheless, PIC simulations are computationally demanding and require high parallelization in multi-core architecture High-Performance Computer-HPC systems [25,26], especially for the modeling of the large spatiotemporal domains in which LIB takes place. The need for such large scales stems from the fact that the laser pulse has to be initialized at many Rayleigh lengths far from the focal spot to ensure that the pulse propagation starts at a distance where there is no plasma generation. Due to the excessive requirements in computational resources for the simulation of 3D domains, it is common practice in PIC models to use 2D domains. 2D PIC models have been successfully used to simulate 3D problems in other research areas, e.g., laser wakefield [26]. Here, the predictions of a 2D and a 3D PIC model were compared in a restricted spatiotemporal domain of the problem, showing perfect agreement and thus validating the 2D global approximation. Ultimately, the agreement between the model's predictions and the experimental results provides further evidence about the validity of the presented 2D approximation.

Experimental measurements
The experiments presented here were carried out in the facilities of the Institute of Plasma Physics and Lasers (IPPL). The experimental setup is shown in Fig. 1. A Ti: Sapphire femtosecond laser (Amplitude Technologies, Pulsar PW, probe beam) emitting ∼ 23 fs pulses with central wavelength λ = 805 nm and energies up to 2.7 mJ was used for the plasma generation. The pulses were focused in the air by an optically thin lens with a focal length f = 10 cm, corresponding to f/8 focusing conditions; the transmittance of the lens was more than 99%. The energy of the laser pulses was measured before and after the breakdown spot using an energy meter (Gentec, Maestro). Measurements were carried out for three pulse energies, .9, 1.8 and 2.7 mJ, corresponding to peak intensities at the theoretical focal spot in vacuum of 6.85 × 10 16 , 1.37 × 10 17 and 2.06 × 10 17 W/cm 2 . From the measurements, the absorbed energy was determined as the difference between the energy E b before and the energy E a after the breakdown: The values of E b and E a were derived by averaging over 50 laser pulses.

Multiphysics PIC model
The PIC method follows the kinetic description for collisionless, weakly coupled plasmas, based on the Vlasov-Maxwell system of equations. In particular, the propagation of the electric and magnetic field of the laser pulse is calculated via the Maxwell equations: where E and B are the magnitudes of the electric and magnetic fields, ρ τ , J T are the charge and current densities, respectively, and ε 0 , μ 0 the electric permittivity and magnetic permeability of the free space. Moreover, the description of the spatiotemporal evolution of the charged particles distribution function in the plasma due to the induced Lorenz force is calculated via the Vlasov equation: , t) is the distribution function for species n (protons, electrons, etc.), q n and m n are the charge and mass of species n, c is the speed of light in vacuum and v is the particle's velocity. The change in the particle species and motion influences the electric and magnetic fields, and the respective changes are introduced in the model via the charge and current densities of the different particle species: The generation of free electrons in the EPOCH code is modeled through electron rate equations depending on the ionization regime. Particularly, multiphoton ionization (MPI) is taken into account by a Wentzel-Kramers-Brillouin approximation based on Ammosov et al. [24,27]. Tunneling ionization (TI) is modeled via the Ammosov-Delone-Krainov (ADK) [28] ionization rate, while barrier suppression ionization (BSI) is modeled by a correction on the ADK rates, introduced by Posthumus et al. [29]. The dominant ionization regime is determined using the Keldysh parameter: where ω is the laser frequency, m e is the electron mass, ε is the ionization energy for the electron and e is the electron charge. In the EPOCH code, transition from MPI to TI takes place when the Keldysh parameter becomes smaller than 0.5; thus, for MPI γ > 0.5 while for TI γ < 0.5. Computationally, each time the conditions for ionization of a macroparticle are met, the macroparticle is removed from the system and is replaced by two new macroparticles representing the ion and electron species. Moreover, ionization occurs when: where U I is a uniform random number and W the ionization rate. In EPOCH, optical energy absorption by two distinct mechanisms is taken into account: a. photoionization (multiphoton, tunneling, barrier suppression ionization) b. plasma absorption In particular, energy conservation is accounted for by a current density correction through Poynting's theorem, while in tunneling ionization and BSI, the energy loss from the field is the ionization energy of the electron [24]. In multiphoton ionization, the energy loss is calculated as the total energy of the K absorbed photons, E ω = K ω. For a detailed description of the optical energy losses in EPOCH, see [24]. As a final remark, self-phase modulation is not taken into consideration, as the code does not support the generation of new frequencies.

Numerical model
In the EPOCH code, the standard PIC method is implemented where the particles are pushed by the Boris pusher algorithm, while the electromagnetic field is advanced by the Yee solver and the current density is calculated by the Villasenor and Buneman scheme [24,[30][31][32]. Particularly, the 2nd-order Yee solver and macro-particles with a 1st-order b-spline shape function corresponding to a 3rd-order weighting, are used [24]. Photoionization, and particularly multiphoton, tunneling and barrier suppression ionization, is included in the model via the Field Ionization module [23]. Collisional ionization and binary collisions are ignored since, during the interaction of the femtosecond laser pulse with the self-induced plasma, a very low percentage of plasma electrons undergo collisions. For typical electron densities ∼ 10 25 m −3 , which are the case in the presented simulations, the characteristic electron collision time is estimated [21] to be about 350 fs, which is at least 10 times longer than the laser pulse duration [24]. The 2D models are simulated within a moving window in the laboratory frame. The computational domain size is 800 × 45 μm and is discretized by 4000 × 900 cells. A non-ionized, 20% Oxygen, 80% Nitrogen atomic gas mixture is considered with a density of ρ = 2.6 × 10 25 m −3 . The laser pulse is modeled as a linearly p-polarized Gaussian beam with a Gaussian temporal profile with a beam waist at focus in vacuum of w o = 6.03 μm. The resulting Rayleigh Length (RL) is z R = 143 μm. The laser pulse is initialized at 40z R before the focal spot, where ionization has not yet initiated due to the low laser intensity. The laser pulse propagates at a distance of 80z R , namely 40z R before and 40z R after the theoretical focal spot. The total simulated time window is 40ps and the timestep is Δt = 180 as. In this manner, the simulation spans over

Results and discussion
In this section, experimental results and results from PIC simulations are presented from the interaction of tightly focused 0.9, 1.8 and 2.7 mJ laser pulses with the ambient air. From the PIC simulations, the pulse propagation and the plasma generation are studied in detail, providing information regarding the dimensions of the plasma channel, the generated free electron density and the pulse intensity along the propagation path. Also, the total absorbed optical energy is estimated and compared to the experimental measurements. Figure 2 presents the evolution of the laser pulse during the simulation. Figure 2a-c show the electric field strength (blue/red) of the 0.9 mJ laser pulse and the free electron density (brown/yellow) at three different positions along the propagation path. Figure 2a is before the focus, where there are no generated free electrons, while Fig. 2b is on focus where all the gas atoms are 1 s 1 ionized. Figure 2c shows the simulated pulse after the focus, where the electron density has again dropped to zero. Figure 2d, e show the 1.8 and 2.7 mJ pulses, respectively, on the focal spot. By comparison of Fig. 2b, d and e it occurs that the plasma channel starts forming at different positions in space and has a different geometry for each simulated laser energy. It is important to note that the simulations showed only first ionization taking place for all the laser pulse energies under consideration. Moreover, first ionization is saturated along a significant part of the plasma channel in all three cases. These aspects are clearly shown in the next figure. Figure 3a presents the absorbed optical energy as a function of the pulse propagation distance, expressed in Rayleigh lengths, for the three laser pulse energies. The focal spot in vacuum corresponds to 40z R = 5880 μm. From the figure, the generated plasma channels can be approximately determined by the domains where the absorbed energy increases. Indicatively, for the .9 mJ pulse, the absorption starts at 25z R and terminates at 40z R , corresponding to a plasma channel of 2.15 mm length, while for the 1.8 mJ pulse, the absorption starts at 10z R and terminates at 40z R , corresponding to a plasma channel of 4.3 mm length. For the 2.7 mJ pulse, absorption takes place almost along the whole simulation domain, indicating a plasma channel of more than 11.4 mm long. Note that, in all three curves of Fig. 3a, there is a region where the absorption rate is constant and having also its maximum value. For example, for the curve of the 1.8 mJ pulse, this region starts approximately at 25z R and ends at 35z R . It can be deduced that, in these regions the first ionization of the air particles is saturated and the free electron density equals the particle density of the medium, namely ρ e = ρ = 2.6 x 10 25 m −3 . Preliminary evaluations also showed that the peak electron temperature for the .9 mJ at the focal spot is ∼ 10 5 K, which is in good agreement with the findings in [33]. However, an extensive analysis of the electron and ion temperatures is beyond the scope of this work. Figure 3b shows the peak value of the laser E-field, which is influenced by the optical focusing and the plasma-related refocusing, as well as the energy depletion of the pulse due to the ionization and plasma absorption processes, as a function of the pulse propagation distance. As expected, the regions of maximum absorption of Fig. 3a coincide with the regions in Fig. 3b where the E-field takes its maximum values.
Moreover, by observation of the tightest beam waist and the Gouy phase of the pulse, the focal spot of the .9, 1.8 and 2.7 mJ pulses was found to be shifted with respect to the theoretical focal spot, by −160, −310 and −460 μm, respectively. This can be attributed to the combined effect of the optical absorption, which leads to a progressive energy depletion of the pulse as it propagates in the plasma, and the plasma defocusing, that leads to a wider beam waist and, hence, lower optical intensities. The focal shift estimated by the PIC model demonstrates the capability of such an approach to capture complex phenomena that significantly affect the optical energy deposition in LIB. The actual beam waists were measured to be w 0 = 11.4 μm, 11.7 μm and w 0 = 12.6 μm for the .9, 1.8 and 2.7 mJ pulses, respectively. The respective domains where the free electron density is nonzero were found to be 36z R , 58z R and ∼ 80z R .
Finally, the comparative plot of the experimentally and computationally evaluated absorbed energy in terms of the initial pulse energy is shown in Fig. 4 for the three laser pulses. The error in the measurement of the optical pulse energy before breakdown has been measured around 4%, while after breakdown the error is estimated to be 10%. For such errors, the size of the error bars is smaller than that of the markers and for this reason, they are not shown. The good agreement between the measurements and the simulations proves the validity of the presented computational approach and demonstrates that, in contrast to the common practice of using analytical systems of PDEs, PIC codes can be effectively used for the simulation of multiscale phenomena within large spatiotemporal domains. Note that, the parameter values used for the simulations correspond exactly to the actual values of the experimental parameters and no additional adjustments were made in order to match the model's predictions to the experimental measurements. Moreover, the percentages of the absorbed energy are 27%, 45% and 50% for the three pulses, respectively. It becomes evident that, in the energy range under consideration, the deposition of optical energy in the medium becomes more efficient with increasing laser pulse energy. Preliminary experimental measurements carried out by some of the authors have shown that this increase in the absorption efficiency reaches a plateau for pulse energies above 4mJ. This saturation can be potentially attributed to the strong defocusing induced by the increased plasma density, which leads the pulse intensity to significantly drop before reaching its theoretical peak value. How-

Conclusions
This work presented a computational study, validated via experiments, of the optical energy deposition in highly ionized ambient air by strong femtosecond laser pulses. A computational model was developed that allowed for the simulation of the interaction of the laser pulse with the atmospheric air across a wide range of pulse-and medium-related parameters in a multi-tens of picoseconds time scale. The model was based on the EPOCH simulation code, which implements a PIC approach, and was validated by comparison to experimental results. Specifically, the optical energy absorption in air breakdown from laser pulses with different energies was simulated computationally and measured experimentally, showing very good agreement.
For the PIC simulations, large spatial and temporal domains were used in order to capture the energy deposition process to its full extent. The pulse propagation was initiated sufficiently far from the theoretical focal spot, so that the intensity of the yet unfocused pulse to be below the ionization threshold at the beginning of the simulation. The simulation spanned over more than 4 orders of magnitude in the time domain. The model was capable of capturing potential 2nd-, 3rd-or higher-order ionizations, as well as nonlinear propagation effects such as plasma defocusing, which are very complicated to model directly through systems of energy transport and electron yield equations. The presented work demonstrated that, given the availability of the required computational power, PIC simulations are appropriate for modeling the optical energy deposition in gaseous media over large spatiotemporal domains. Moreover, they can provide valuable insight in the involved physical processes leading to the plasma formation, by estimating the spatiotemporal dimensions of the phenomenon, the degree of ionization of the medium, the time-evolution of the free electron density and the temperatures of free electrons and heavy particles in the interaction region. Such information extracted by the PIC model can be of significant importance for the development of complete equation-based models of laser-induced breakdown.
In the future, such a model based on an optical energy transport equation including a sink accompanied by electron rate equations will be developed, with reference to the results obtained from the PIC simulations. The model will be based on the transport equation with a sink term, similarly to the model developed by Bulgakov, Bulgakova et al. [34][35][36]. Such a model will account for the spatiotemporal laser pulse propagation under the influence of a Gaussian focusing lens, while the losses due to optical energy deposition in the medium by photoionization in the MPI and TI regimes, as well as due to plasma absorption, will be included in the sink term. The transport equation will be solved in combination with a spatiotemporal free electron rate equation, while an additional spatiotemporal equation will be potentially used to describe the modulation of the refractive index due to plasma defocusing and Kerr self-focusing.
The ultimate objective of such a complete model is the computational evaluation of the secondary emission of light and sound radiation following laser-induced breakdown and the determination of the dependence of the secondary radiation from the parameters of the laser pulse and the medium. This will allow for the precise modulation of the laser radiation in order to achieve secondary light or sound shaping for LIB applications such as optical spectroscopy, optical high-harmonic generation, X-ray generation, electron and ion sources, remote acoustic measurements and sound reproduction, as well as remote sensing, material diagnostics and medical applications.
ulations and analyzed the results. NP and MB had the overview of the experiments, NP and MT advised regarding the physical aspects and VD had the overview of the computational work. KK, IT and NP authored the manuscript, while all the co-authors reviewed and contributed to the final formulation of the manuscript.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: Due to the very large size of the simulation data and the specialized software needed for their processing, it is not beneficial to provide data at this stage. However, if there is interest for accessing the data, please contact the corresponding author.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.