NUMERICAL INVESTIGATION OF THE PLASMA FORMATION IN AIR GENERATED BY 355 nm Nd:YAG LASER PULSES

In the present work, a numerical analysis is performed to investigate the comparative contribution of the mechanisms responsible for electron gain and losses in laser-induced breakdown. In this regard, we adopted a simple theoretical formulation relying on the numerical solution of a rate equation that describes the growth of the electron density due to the joined effect of multiphoton and cascade ionization processes. The rate equation also includes the effect of electron loss due to diffusion, attachment and recombination processes. The analysis considered atmospheric air irradiated by a Nd:YAG laser radiation at a wavelength of 355 nm with 5 ns pulse duration full-width halfmaximum (FWHM).


INTRODUCTION
Since its first discovery in the mid-1960s (ZEL'DOVICH and RAIZER, 1967), laserinduced breakdown in the air has been extensively studied due to its potential and promising applications in the laboratories ranging from few-cycle optical pulses and high order harmonic generation (TOSA et al., 2003), light detection and lightning protection (TZORTZAKIS et al., 2001) to diagnostic techniques such as laser-induced breakdown spectroscopy (LIBS) and laser electronic excitation tagging (DOUCET et al., 2011). Parallel to the experimental progress in the field, many theoretical approaches have been developed to describe the temporal evolution of the free electron density and temperature in laser-induced plasmas (LIP) (MORGAN, 1975;MARCH et al., 2009). Nevertheless, the complete description of the mechanisms involved in the laser-induced breakdown of the air is still an open scientific filed.
Laser-induced breakdown (LIB) can be defined as the formation of an ionized gas during or by the end of the laser pulse (SINGH and THAKUR, 2007). As described in the pioneering work of Zel'dovich and Raizer (ZEL' DOVICH and RAIZER, 1967), during the rise time of the laser pulse, the seed electrons required for the cascade ionization are produced by direct ionization of the gas atoms or molecules via photoionization, also called multiphoton-mediated electron cascade (TOZER, 1965). Following the increasing phase of electron density, the abovementioned processes of cascade ionization and photoionization are in competition with loss processes by diffusion, attachment and recombination within the focal volume. Ιt is important to describe the evolution of the plasma properties after the formation of the initial LIP, especially the density of free electrons. This study aims at clarifying the correlation between the pulse length and the dominant ionization mechanism.
Accordingly, in the present work, a modified model previously developed by (THIYAGARAJAN et al., 2012) based on the solution of the rate equation is proposed to numerically investigate laser-driven plasma formation in atmospheric air induced by a Gaussian pulse with a full width at half the maximum (FWHM) of 5 ns at 355 nm. The model takes into account the generation of electrons due to the combined effect of both multiphoton and cascade ionization processes. These processes are opposed by the loss of electrons within the interaction region through diffusion out of the focal volume, attachment and recombination. This work focuses on studying the individual contributions of each process to the gain and loss of electrons. This is achieved by using Wolfram Mathematica (WOLFRAM, 1999) as the platform for numerical solving of the rate equation.
The paper is structured as follows: in Section 2, the basic concepts of LIP from nanosecond lasers in air and the gain and loss processes involved are briefly described. Numerical solution of the electron density growth equation is conducted by using the Keldysh-PPT theory and Drude model, in which, the significance and behavior of the photoionization and cascade ionization are described for the given experimental parameters, using Wolfram Mathematica scripts. Section 3 discusses the obtained results, while Section 4 provides a summary and conclusion.

Electron Density Rate Equations
The general form of the rate equation that describes the time variation of the free electrons density generated by nanosecond laser pulses in ambient air can be expressed as (MORGAN, 1975): where the first two terms describe the evolution of the electron density generated by the combined effect of multiphoton ionization, ( ( )), and cascade ionization, Photoionization in gases can be associated with both multiphoton and tunneling, which essentially are the two limiting cases of the same physical process of nonlinear photoionization. The regime in which the tunnel or the multiphoton processes take place is determined by the Keldysh parameter, (KELDYSH, 1965). This parameter is defined as a ratio of the laser frequency, , to the tunneling frequency, , characterizing the time of electron tunneling through the potential barrier: = / . In case of large laser intensities (low frequency limit, < ), the typical Keldysh parameter is smaller than 1, ≪ 1, while for low laser intensity (high-frequency limit, > ), Keldysh parameter is much greater than 1, ≫ 1. Multiphoton ionization takes place when ≫ 1, while tunneling ionization dominates when ≪ 1. After the appearance of the Keldysh work, its results were refined by Perelomov, Popov and Terent'ev (PPT) (PERELOMOV et al., 1967). Based on the generalized Keldysh-PPT formulation, the photoionization rate can be written in the following form (ETTOUMI et al., 2010): where ( ( )) denotes the peak optical amplitude ( ( ( )) = √2 ( )/ 0 , where is the velocity of light in vacuum and 0 is the electrical permittivity of vacuum), is the ionization energy of the air molecules (THIYAGARAJAN and THOMPSON prolonging the latter to low intensity multiphoton ionization regimes. By doing so, Eq. 2 is able to describe the air molecule ionization after the substitution (TONG et al., 2002): A general analytical solution of Eq. 2 is not available, and therefore, the Keldysh-PPT ionization rate has to be determined numerically. Our numerical analysis allowed us to explicitly distinguish multiphoton and tunneling contributions to the total photoionization rate using Eq. 2. Results are presented in Fig. 1. The figure describes the total photoionization rate, ( ( )), as well as the significance and behavior of the individual roles of the multiphoton, ( ( )), and tunneling ionization rates, ( ( )), versus the laser intensity, ( ), using Wolfram Mathematica scripts. The rate ( ( )) is given as Eq. 3.7 in (FERRIS, 2014) and is as follows: while for the multiphoton rate ( ( )) we used Eq. 3.9 of (FERRIS, 2014): ( ( )) = Here, the laser wavelength is assumed to be 355 nm. The solid red line in Fig. 1 corresponds to the total photoionization rate, ( ( )). The blue dashed and green dash-dotted lines represent ( ( )) and ( ( )) respectively. For the sake of simplicity all rates are normalized to unity. The ionization rates were modeled for a 355 nm laser wavelength.
In Fig. 1 it can be seen that the full rate, ( ( )), only agrees with multiphoton rate for < 5 × 10 13 W/cm 2 and solely with tunneling rate for > 5 × 10 14 W/cm 2 . Tunneling and multiphoton rates overlap with each other in the range 5 × 10 13 − 1 × 10 14 W/cm 2 , i.e. when the Keldysh parameter is ~ 1. This intermediate range of corresponds to the nonadiabatic tunneling regime (YUDIN and IVANOV, 2001). As intensity rises above ~1 × 10 14 W/cm 2 , we approach the pure tunneling regime and the tunneling approximation in Eq. 2 is valid. It is important to mention that the sudden changes in the multiphoton ionization process causes a step-like behavior in the full photoionization rate ( ( )). According to (JI et al., 2019) this happens because the effective ionization energy increases as the laser intensity increases, thereby increasing the number of photons required to ionize electrons.
To summarize the results presented in Fig. 1, the validity of Eq. 2 can be limited within the multiphoton regime and used as an ( ( )) rate in Eq. 1, only when the laser intensities are below 5 × 10 13 W/cm 2 .
As we already stressed, a free seed electron produced via multiphoton ionization gains kinetic energy by absorbing the energy from the laser radiation. Furthermore, in collisions with other atoms, new electrons are generated and a cascade effect takes over. Additionally, according to (THIYAGARAJAN and THOMPSON, 2012) from the total ionization process, 88% are due to cascade and 12% from multiphoton effects at pressure of 760 torr. In contrast to multiphoton ionization, no wavelength dependence is expected for this initiation effect. The rate of this process is commonly adopted from the Drude model which has been widely used in gases, liquids and solids and is given by the following expression (KENNEDY, 1995): Here, represents the electron-phonon momentum relaxation time, is the refractive index of air and is taken as the average mass of molecules consisted in air (SHEN, 1984). In this work, we adopted a value of estimated by (THIYAGARAJAN and SCHARER, 2008) for the atmospheric air. The cascade rate has a negligible contribution for pulses shorter than 20 ps, since electrons have insufficient time to couple their energy during the laser pulse [8]. In the present study, for LIB of air with nanosecond pulse widths the cascade term cannot be neglected.

Numerical Calculation
Considering a Gaussian temporal profile of the laser intensity within the focal volume, we solved Eq. 1 numerically to obtain the time evolution of the free electron density. Thus, the laser intensity ( ), can take the following form (GOCIĆ et al., 2009): where 0 is the peak intensity, 2 is the pulse width and is the time at which intensity is computed.
As has been discussed previously (GAABOUR et al., 2012; THIYAGARAJAN et al., 2012), the electron density equation (Eq. 1) can be solved by the fourth order Runge-Kutta algorithm. This algorithm has been developed by (MAEDER, 1991) and written into a package called RungeKutta.m in Wolfram Mathematica. Because of that, Eq. 1 can be solved accurately with routines that already exist in this package, using function NDSolve[eqns, , {t, t min , max }]. This internal function of Mathematica can handle a wide range of ordinary differential equations. In a system of ordinary differential equations (ODE), there can be any number of unknown functions , but all of these functions must depend on a single independent variable , which is the same for each function. Also, the function NDSolve finds solutions iteratively. It starts at a particular value of , covers the full range from to with a selected time step. Α time step of ∆ = 1 ps was found to be sufficient for the evaluations presented in this study.
The Runge-Kutta method is a widely used numerical method because of its high accuracy and high stability. This method is divided into two types: the explicit Runge-Kutta (ERK) method and the implicit Runge-Kutta (IRK) method. The IRK method has wider stable region and higher accuracy compared to ERK methods. Therefore, the IRK method is suitable for solving stiff differential equations (BURRAGE, 1978), and abundant results have been mentioned in the literature (JAMESON, 2017).
To study the influence of the gain and loss terms of Eq. 1 on the breakdown, we included them sequentially (step by step). As a result, we formed the system of differential equations and solved them numerically using Wolfram Language's as follows: The code presented in Eq. 7 produces the numerical solutions for every electron density equation.

RESULTS AND DISCUSSION
As mentioned in the last subsection, Eq. 1 is incorporated into a Wolfram Mathematica script together with the relations which describe the rate coefficients of the physical processes encountered into the model and solved numerically applying Runge-Kutta fourth order technique. The calculations are performed under the experimental conditions given by Tambay  and Thareja (TAMBAY and THAREJA , 1991). They reported LIB studies of laboratory air using 355 nm radiation for focal spots varying from 30 μm to 100 μm, producing laser intensities up to 0.1 -1 TW/cm 2 . Lastly, it should be stated that the generation of free electrons has to reach a critical electron density, equal to (ELIEZER and MIMA, 2008): [cm −3 ] ≅ 1.1×10 21 2 [ ] ( = 8.7 × 10 21 cm −3 for 355 nm laser) for optical breakdown to occur. The critical value of ~10 21 cm −3 is widely accepted in the field of LIB as the onset of air breakdown.
The relative contribution of cascade and multiphoton ionization to the free electron generation rate, which is characterized by the rate parameters, ( ( )) and ( ( )), is strongly influenced by the dependence upon the time-dependent laser intensity, ( ), and the respective frequency, , of the electromagnetic radiation. From the numerical analysis of the multiphoton rate, ( ( )), in the nanosecond regime it follows that when the pulse duration lies between 0.01 ns and 100 ns the rate is proportional to 0 , if ≤ 10. In addition, the Keldysh coefficient, , depends on the radiation intensity exponentially. Such observations are completely consistent with the reports of ZEL'DOVICH and RAIZER, 1967 and FERRIS, 2014. As far as the intensity dependence of the cascade rate, ( ( )), during the effect of the nanosecond pulse is concerned, it can be considerably stronger than the multiphoton rate, ( ( )).
To examine the exact contribution of the cascade and multiphoton ionization to the growth rate of the free electrons, numerical calculations are carried out to obtain the results displayed in Fig. 2 as curves ( ( )) and ( ( )). Laser intensity is fixed to the value: Figure 2. Time evolution of the free electron density in the presence of ( ( )) + ( ( ), ( )) and in the absence of ( ( )) of the cascade process. The red solid curve is the critical electron density, set as the damage criterion.
As demonstrated in Fig. 2, the multiphoton-cascade and multiphoton-only produced electron density showed similar behaviors during the early stages of the laser pulse. This may be attributed to the competition between the two ionization mechanisms during this time which hinders the electron density growth. The curve ( ( )) showed a gradual increase reaching a value of 8.8 × 10 15 cm −3 at the peak of the laser pulse, then it keeps this value until the end. Such a result is in accordance with (ELIEZER et al., 2018) who stated that in multiphoton ionization of the atmospheric air, the plasma density has to be in the range of 1 × 10 15 cm −3 − 1 × 10 16 cm −3 or less to ensure that the contribution of plasma nonlinearities to the refraction index is negligible. Also, from this figure one can observe that cascade ionization contributes effectively to the breakdown of air during the nanosecond regime and exceeds the rate of multiphoton ionization by orders of magnitude. The curve ( ( ), ( )) started with a slower increase followed by an almost constant value at about 1 × 10 12 cm −3 until its end. This gradual increase of the electron density indicates the slow rate of cascade ionization during the first half of the laser pulse. After ~20 ns the free electron density becomes large enough so that the cascade ionization turns into the dominant uncontrolled process and may cause a breakdown (ELIEZER et al., 2018). We would like to note that the critical density of 8.7 × 10 21 cm −3 is reached at 40 ns.
The main physical processes reducing the number of electrons in a laser-induced air plasma are recombination with air molecules and ions, attachment and diffusion out of the interaction volume. They are described by the corresponding rate coefficients in Eq. 1. The total rates of the processes involved in the gain and loss of electrons are given in Fig. 3. All parameters are the same as in Fig. 2.
From Fig. 3 it is clear that the four curves exhibit almost the same behavior during the first 10 ns. The fast increase shown during the early stages of the laser pulse assures that the electron growth proceeds mainly via multiphoton and cascade ionization. This result is consistent with (CHENG et al., 2016). From this figure it becomes apparent that diffusion and recombination are the most important electron loss mechanisms. The electron losses due to attachment are orders of magnitude smaller in this particular analysis domain thus can be neglected. In addition, our analysis showed that with the end of the laser pulse, the plasma gradually relaxes as a result of the cumulative contribution of diffusion, attachment, and recombination of ions and electrons, until local thermodynamic equilibrium with the surrounding gas is restored. Theoretical and experimental studies of these aspects can be found

CONCLUSION
In summary, we analyzed the time evolution of the free electron density in the absence and presence of each gain and loss process included in the rate equation describing the time evolution of the free electron density in the laser-induced breakdown. The presented results verified that in the nanosecond regime, the dominant electron gain process is cascade ionization (which however cannot be initiated without the contribution of multiphoton ionization) and the respective dominant loss processes are diffusion and recombination. Our results can be extended to any medium for which the index of refraction, ionization potential, multiphoton ionization, cascade ionization, electron-ion recombination, attachment, and diffusion rates are known.