Ionisation processes and laser induced periodic surface structures in dielectrics with mid-infrared femtosecond laser pulses

Irradiation of solids with ultrashort pulses and laser processing in the mid-Infrared (mid-IR) spectral region is a yet predominantly unexplored field with a large potential for a wide range of applications. In this work, laser driven physical phenomena associated with processes following irradiation of fused silica (SiO2) with ultrashort laser pulses in the mid-IR region are investigated in detail. A multiscale modelling approach is performed that correlates conditions for formation of perpendicular or parallel to the laser polarisation low spatial frequency periodic surface structures for low and high intensity mid-IR pulses (not previously explored in dielectrics at those wavelengths), respectively. Results demonstrate a remarkable domination of tunneling effects in the photoionisation rate and a strong influence of impact ionisation for long laser wavelengths. The methodology presented in this work is aimed to shed light on the fundamental mechanisms in a previously unexplored spectral area and allow a systematic novel surface engineering with strong mid-IR fields for advanced industrial laser applications.


A. Electron density rate equations
The absorption of light in transparent materials is a nonlinear process as a single photon does not have enough energy to excite electrons from the valence to the conduction band. Following nonlinear photoionisation (where both multiphoton and tunnelling ionisation can occur in different regimes depending on the laser intensity), impact ionisation and formation of STE lead to a variation of the electron densities. To calculate the laser absorption, the densities of the excited electron and the STE are computed by solving the following set of equations = − ( ( 1 ) + (1) ) + ( ( 2 ) + (2) ) − = − ( ( 2 ) + (1) ) (1) where N V =2.2×10 22 cm -3 1 corresponds to the atomic density of the unperturbed material while N e and N STE denote the free electron and STE electron densities, respectively . The last term in the first equation corresponds to free electron decay that is characterised by a time constant τ tr (τ tr ~150fs in fused silica) that leads to a decrease of the electron density. By contrast, ( ) + ( ) (i=1,2) correspond to the combined photo-ionisation ( ) and impact ionisation; ( ) stands for the (usually termed) avalanche ionisation rate coefficients 1,2 of the impact ionisation term that depends on which bandgap the electrons need to surpass ( (1) (=9eV) or (2) (=6eV)). The photoionisation rates are computed using the Keldysh  is the effective band gap which takes into account the oscillation energy of the free electrons in the electric field, and the conservation of energy and the momentum during the collision between the free and bound electrons. Also, c is the speed of light, ε 0 stands for the vacuum permittivity, n is the refractive index of the material, while I is the peak intensity of the laser beam, respectively, and τ c is the electron collision time.
In regard to the value of τ c , there are many reports in which the collision time has been considered to be a constant or a varying parameter. More specifically, values of τ c vary in a range between 0.1 fs and 10 fs 3 while in other studies values equal to 1.0 fs 4 , 1.5 fs 5 , 23.3 fs 6 , 10 fs 7 , 1.7 fs 8 , 0.5 fs 9 have been reported. Thus, simulations with the above choices for τ c yielded good agreement with experimental data.
On the other hand, there are other reports in which a more rigorous analysis was followed to compute the collision time: In Chimier et al 1 , the contribution of electron-phonon, electron-electron, electron-neutral and electron-ion collision frequencies were used to compute τ c . In other reports [10][11][12] , the electron density distribution function was also considered to calculate the electron collision frequency and the electron collision time. A remarkable impact of the value of the collision frequency on observable quantities such as the damage threshold has been demonstrated with simulations in other materials (i.e. Silicon 13 ) which emphasized the need for a precise evaluation. Nevertheless, in the absence of experimental results in the mid-IR, the aim of the present work was firs tly to provide a consistent multiscale modelling approach in which some approximations have been made rather than to consider a more precise evaluation of the collision frequency that could be the subject of revised model. Therefore, the collision frequenc y is taken constant in this work. In regard to the value chosen in this manuscript, the motivation was based on the work of Burakov et al 5 in which the damping term ω L τ c =3 was used that resulted into a collision time equal to1.5 fs (for λ L =800 nm). In the case of mid-IR, for laser wavelengths between 2 μm-4 μm, the collision time was calculated from the above expression yielding values for τ c in the range [3.1fs, 6.4fs].
The dependence of the avalanche ionization rate coefficient AI on the collision time to which carriercarrier or carrier-ion scattering times (and dephasing oscillations) contribute 14,15 implies that a possible future revision of the model should include a term for impact ionization in which the collision time is a varying parameter; in relevant reports 14,15 , an electron temperature dependent expression for τ c was considered.
The photoionisation rates ( 1 ) and ( 2 ) include both multiphoton and tunneling ionisation processes that are provided from the following expressions (i )  (6) The two ionisation processes become more efficient in different regimes; more specifically, tunneling ionisation becomes dominant for γ<<1 while multiphoton ionisation becomes more efficient for γ>>1 2 . In contrast to the avalanche ionisation rate coefficient and impact ionisation, multiphoton and tunneling ionisation rates do not appear to be dependent on dephasing effects.
With respect to the intensity of the laser beam I, previous studies consider that the attenuation of the local laser intensity is determined by multiphoton ionisation and inverse bremsstrahlung (Free Carrier) absorption 9 . Nevertheless, the presence of STE states and the possibility of retransfer of carriers in these states back to the conduction band through multiphoton ionisation require modification of the spatial intensity distribution; therefore a revised 3D model (in Cartesian coordinates: x =( X, Y, Z )) should include those contributions 1,16,17 where (i) ph N corresponds to the minimum number of photons necessary to be absorbed by an electron that is in the valence band (i=1) or the band where the STE states reside (i=2) to overcome the relevant energy gap and reach the conduction band. I peak (≡ 2 √ 2 √ ) corresponds to the peak value of the laser intensity, J stands for the (peak) laser fluence, R(t,x,y,z=0) stands for the reflectivity of the material while R 0 is the irradiation spot radius (R 0 =15μm in our simulations). The second equation of Eq.7 gives the spatial distribution of the intensity profile in Cartesian coordinates (and on the surface of the irradiated material).
It is also noted that, for the sake of simplicity, it is assumed that the beam shape does not change upon propagation; thus, only attenuation losses are taken into account and no Kerr or plasma distortion to the pulse phase were considered. A more accurate description of the laser energy distribution has been presented in other studies, by solving Maxwell's equations 9 or by including the shape change and Kerr effect to evaluate damage in the bulk 5 . However, we believe that a more accurate expression for the intensity profile would not lead to substantially different surface effects.

B. Dielectric constant
The refractive index of fused silica when it is in an unexcited state is denoted by n 0 18 (at I=0) ant it provided by the following expression The expression = 0 2 yields the dielectric constant of the unexcited material. Following irradiation of SiO 2 with femtosecond pulses, a temporally dependent density of excited carriers is produced that modify the dielectric constant of the material The reflectivity and free carrier absorption coefficients are given by the following expressions 22 22 2 (1 ) where k is the extinction coefficient of the material. The real part of the refractive index n is the sum of two terms, one (n 0 ) that corresponds the refractive index in an unexcited state and a second (n 2 I) that corresponds to the nonlinearities induced by Kerr effect (i.e. 02 n n n I  ). n 2 stands for the Kerr coefficient. On the other hand, the dielectric constant is given by the following expressions if the Kerr effect is taken into account 19

C. Electron and lattice heat balance
To describe the morphological properties due to laser irradiation of the material, it is important to explore , firstly, the relaxation process and heat transfer from the ionised material to the lattice system. Due to the metallic character of the excited material, a TTM model can describe the spatio -temporal dependence of the temperatures T e and T L of the electron and lattice subsystems, respectively 20 The source term has been modified properly to take into account all quantities that contribute to the total electron energy balance. Hence, the complete expression for source term S(x,t) is given by 1,11,17 The source term includes terms that are related to the dynamics of the total energy of the electron system. Processes that are taken into account, are the photoionisation of electrons (first and third terms), impact ionisation (second and fourth terms). It is noted that excitation both from VB (first and second terms) and STE (third and fourth) electrons are considered. Finally, the fifth term stands for the free electron absorption (that leads to transition to higher states in the CB) while the last two terms correspond to the energy balance due to electron density reduction (due to defect formation/trapping) and electron energy assuming carrier density variation 21 . Alternative expressions that incorporate average energies of electrons that are excited from STE states and exciton decay have been also proposed 11 , however they have not been used in this work. The reason is that on the one hand, exciton decays are supposed to last very long (in the range of some hundreds of picoseconds 22,23 ) compared to the timescales considered in the excitation process (and therefore their contribution is supposed to be insignificant); furthermore, the average energies of STE electrons contribution has been ignored in this study assuming it is not very significant compared that of the electron in CB. Nevertheless, a more precise inves tigation might be required in a future study by taking into account the role average energy of electrons in STE states in various laser conditions. We point out that a term that describes the divergence of the current of the carriers ( J  ) has not been taken into account (it turns out that the consideration of particle transport of heat diffusion does not vary significantly quantities such as the damage threshold 13 ). The temperature, electron density and temporal dependence of the thermophysical properties, C e , k e are provided from well-established expressions derived from free electron gas based on the metallic character of the excited electron system 1,17 On the other hand, C L =1.6 J/(cm 3 K) 1 while the coupling constant is estimated to be It is noted that the contribution of ambipolar diffusion of dense electron-hole plasma has not been included in the model. On the one hand, various studies have highlighted the role of ambipolar diffusion in the thermal response of the material and more specifically to studies related to damage threshold evaluation (see for example, studies by Danilov et al 24 , and Derrien et al 25 on Silicon). On the other hand, in other studies 13 , simulations did not show remarkable changes if the ambipolar diffusion (or transport) were ignored. It is evident though, that the influence might be dependent on the laser parameters (i.e. intensity, pulse duration, fluence and laser wavelength). Therefore this limitation of the present model could be the motivation for the development of a revised version.
Similarly, the temporal length of the pulse is also very important that can set limitations for the use of Eqs.1-14. One major parameter that influences the precision of the calculations is the pulse duration. The electron dynamics and relaxation processes considered in this study assumed an instantaneous thermalisation of the electron system through electron-electron scattering mechanisms. This is an important ingredient towards assuming that the electron system is in thermal equilibrium and T e is defined. By contrast, for very short pulses (for example, τ p <100 fs), this argument constitutes an overestimation and therefore alternative techniques to describe electron excitation and relaxation processes are required (see for example Ref. [24] and references therein). Alternative modelling approaches should be used in that case such as quantum mechanics -based approaches to account for the thermalisation of the electron system and dephasing effects. Similarly, appropriate corrections should be made if pulse durations smaller than the trapping time are used.

D. Fluid dynamics
To model a surface modification following irradiation with mid-IR fs laser pulses, it is assumed that the laser conditions are sufficiently high to result in a phase transition from solid to liquid phase and upon resolidification a periodic relief is induced on the surface of the material based on the series of processes described in the main manuscript. Depending on the laser intensity, mass removal is also possible if the material is heated above a critical temperature (~ >1.8 ×~4000 for SiO 2 where = 2270 K 21 ). The choice of the critical temperature is based on arguments made in previous reports 21, [26][27][28] . More specifically, a solid material subjected to ultrashort pulsed laser heating at sufficiently high fluences undergoes a phase transition to a superheated liquid whose temperature reaches values ~0.90 × ( stands for the thermodynamic critical temperature 27 ). In this work, the proposed scenario, of modeling material removal is based on a combination of evaporation of material volumes that exceed upon irradiation lattice temperatures close to ~0.90 × and evaporation due to dynamics of Knudsen layer (adjacent to the liquid-vapor interface 22,25,26 ). According to the discussion in Ref. 21 , for many materials, a typical value of is 1-2 times higher than the boiling temperature which is = 2270 K for fused silica 21 , hence ~2 × (assuming the smallest value for the attained value of ). Based on this assumption, the minimum threshold value which should be used is equal to 2 × 0.9 × = 1.8 × . Alternative scenarios for the estimation of the minimum lattice temperature that leads to mass removal could involve a lower temperature, the boiling temperature. It is evident that appropriate experimental setups could assist in a more precise evaluation, however, this is beyond the scope of the present study.
The movement of a material in the molten phase (T melting =1988 K 29,30 ) is given by the following Navier-Stokes equations (NSE) which describes the dynamics of an uncompressible fluid where 0 and stand for the density and viscosity of molten SiO 2 , while P and ⃗ are the pressure and velocity of the fluid. The fluid is considered to be an incompressible fluid (i.e. ∇ ⃗ ⃗ • ⃗ = 0) . In regard to the pressure, there are two terms that require special treatment:  the recoil pressure which is related to the lattice temperature of the surface of the material through the equation 31 where P 0 is the atmospheric pressure (i.e. equal to 10 5 Pa 33 ), L v is the latent heat of evaporation of the liquid (L v =7 kJ/gr 33 ), R is the universal gas constant, and corresponds to the surface temperature. When vapour is ejected, it creates a back (recoil) pressure on the liquid free surface which in turn pushes the melt away in the radial direction 26 which results into a depression of the surface. Furthermore, given the spatially modulated energy deposition on the material, a gradient of the lattice temperature is produced which is, in turn, transferred into the fluid and therefore a capillary fluid convection is produced.
 A precise estimate of the molten material behaviour requires a contribution from the surface tension related pressure, P σ ., which is influenced by the surface curvature and is expressed as P σ =Kσ, where K is the free surface curvature and σ=0.310 N/m 34 is the surface tension. The calculation of the pressure associated to the surface tension requires the computation of the temporal evolution of the principal radii of surface curvature R 1 and R 2 that correspond to the convex and concave contribution, respectively 35 . Hence the total curvature is computed from the expression K=(1/ R 1 +1/ R 2 ). A positive radius of the melt surface curvature corresponds to the scenario where the centre of the curvature is on the side of the melt relative to the melt surface.
Pressure equilibrium on the material surface implies that the pressure P in Eq.15 should outweigh the accumulative effect of P r +P σ . The thermocapillary boundary conditions imposed at the liquid free surface are the following where (u,v,w) are the components of ⃗ in Cartesian coordinates. Values for the thermophysical parameters that are used in the simulations are: 0 = 2.2 gr/cm 3 33 , results from a fitting procedure (Ref. 34,36 ) and σ=0.310 N/m 33 . The melting point of fused silica 1988K is taken as the threshold for a phase transition from solid to liquid while the same isothermal is considered as the same criterion for resolidification. The hydrodynamic equations are solved in regions that contain either solid or molten material. To include the 'hydrodynamic' effect of the solid domain, material in the solid phase is modelled as an extremely viscous liquid (μ solid =10 5 μ), which results into velocity fields that are infinitesimally small. An apparent viscosity is then defined with a smooth switch/step function where Δ is in the range of 10-100 0 K depending on the temperature gradient 21,26 .

E. Numerical scheme
To solve the set of the above equations, a scheme based on finite difference method is used. A common approach followed to solve similar problems is the employment of a staggered grid finite difference method which is found to be effective in suppressing numerical oscillations. Unlike the conventional finite difference method, temperatures (T c and T L ), carrier densities (N e ), pressure (P) are computed at the centre of each element while time derivatives of the displacements and first -order spatial derivative terms are evaluated at locations midway between consecutive grid points. For time-dependent flows, a common technique to solve the Navier-Stokes equations is the projection method and the velocity and pressure fields are calculated on a staggered grid using fully implicit formulations 37,38 . On the other hand, the horizontal and vertical velocities are defined in the centres of the horizontal and vertical cells faces, respectively (for a more detailed analysis of the numerical simulation conditions and the methodology towards the description of fluid dynamics, see Refs. 26,28,30,[39][40][41][42][43].
During the ultrashort period of laser heating, heat loss from the upper surface of target is assumed to be negligible. As a result, a zero heat flux boundary condition is set for the carrier and lattice systems.


For irradiation with one pulse (NP=1), Eqs.1-17 are solved assuming heating of a flat profile; a 2D numerical approach is followed by taking into account the axial symmetry of the problem.


For NP>2, the symmetry breaks and a 2D solution is no longer valid. In that case, a 3D numerical framework is developed (a finite difference methods is used again 26 ). The incident beam is no longer perpendicular to the modified profile and therefore the surface geometry influences the spatial distribution of the deposited laser energy. Hence, appropriate modification to the numerical scheme is required to compute energy absorption (for example Eq.7 needs to be corrected). Typical Fresnel equations are used to describe the reflection and transmission of the incident light. Due to multiple reflection and light entrapment, the absorption of the laser beam is modified 44 . The calculation of the pressure associated to the surface tension requires the computation of the temporal evolution of the principal radii of surface curvature R 1 and R 2 that correspond to the convex and concave contribution, respectively 35 . Hence the total curvature is computed from the expression K=(1/ R 1 +1/ R 2 ) 26 .
In regard to the material removal simulation, in each time step, lattice and carrier temperatures are computed and if lattice temperature reaches ~>1.8 , mass removal through evaporation is assumed. In that case, the associated nodes on the mesh are eliminated and new boundary conditions of the aforementioned form on the new surface are enforced. In order to preserve the smoothness of the surface that has been removed and allow an accurate and non-fluctuating value of the computed curvature and surface tension pressure, a fitting methodology is pursued 26 .
F. Impact of fluence, pulse duration, laser wavelength on excitation levels.
To evaluate the role of STE states, the electron density and its variation has been calculated as a function of the laser fluence, the pulse duration and the laser wavelength (Fig s.S1-S4) for single shot simulations. It is evident that the variation of the electron densities if STE states are formed (Eq.1 in the main manuscript) increases as the higher fluence increases. On the other hand, although simulation results (Figs.S1-S4) indicate that the role of STE can be ignored for small fluences, it is noted that experimental observations have shown that the presence of defects and incubation effects influence significantly the damage threshold of the irradiated material 16,45,46 for multipulse irradiation. Moreover, it is noted that in multiple shot experiments, the change of surface morphology (formation of crests and wells) vary the energy absorption, excited carrier densities and thermal response of the material and therefore a thorough investigation of the impact of repetitive irradiation (including the impact of STE) is required.       Simulations have been performed to compute the variation of the SP wavelength with increasing irradiation dose (NP). To compute the SP wavelength, a multiscale approach has been followed to take into account the carrier density for the profile induced as a result of the surface profile after irradiation with the NP th -1 pulse. In Fig.S5, results are shown for λ L =2.6 μm for SP wavelength as a function of NP.