Mechanism of Laser Induced Filamentation in dielectrics

Laser filamentation in transparent material has a wide range of applications, from three dimensional manufacturing to biological technologies. Various experimental results showed that femtosecond laser pulse filamentation in fused silica strongly depends on laser focusing conditions. However, the physical mechanism governing each regime has not been fully understood. For the first time, single and multiple re-focusing of the laser pulse in interaction of femtosecond laser pulse with fused silica, and consequent single and multiple damage zones (filaments) have been observed in our extensive three-dimensional, high resolution FDTD (finite-difference time-domain) simulations. We show that Kerr nonlinearity plays a crucial role loose laser focusing regime, while it is not an important factor in tight laser focusing regime, where geometrical focusing becomes important. Our simulation results agree well with existing experimental findings. In addition, the improved analyti- cal model prediction gives a reasonable estimate of the shift from Kerr nonlinearity regime to linear geometrical focusing regime.


I. INTRODUCTION
Laser filamentation in dielectrics, gasses and liquids has been an active research area since the advent of powerful femtosecond lasers(see Refs. [1][2][3] and references therein). Filamentation in air has applications in atmospheric science, lightning control, rain making and remote measurement of electric field, microwave guidance and remote sensing [4]. Femtosecond laser interaction with water, is of interest to surgical techniques [5]. The critical power of self-focusing is in the Gigawatt range for gasses, and around Megawatt for solids due to higher nonlinear index. The reason is higher nonlinear index in solids [6]. In fact, the nonlinear refractive index is three orders of magnitude larger in solids than in gasses. Therefore filamentation in solids is not only interesting because of its applications but also because all physical effects can be studied on a reduced scale.
Although Gaussian beams are the most common beam shape used in the material processing by far, nonconventional beam shapes such as Bessel beams have been used for applications that couldn't be feasible with Gaussian beams [18,19]. Evolution of femtosecond laser pulses at different wavelengths have been addressed in Refs. [20,21]. The effect of single and multiple pulse length and energy as well as numerical aperture on laser energy deposition were compared in [22]. Although a large number of experimental and theoretical studies have been reported, the dynamics of filamentation at 800 nm is not well understood. A variety of experimental studies of filamentation in transparent material in the past two decades have shown that filamentation strongly depends on laser focusing conditions [23]. Various experiments reported that for tighter laser fo-cusing conditions, higher NA regime, one damage zone (or void structure) and for looser focusing conditions, lower NA regime, longer single and multiple damage zones can be generated [15]. Experimental results of Refs [12,24,25], with numerical aperture N A = 0.6 , reported the observation of sub-micron damages in the form of the voids. On the other hand, using an intermediate aperture N A = 0.25 − 0.5, single and multiple damage ares have been reported in Ref [8,26,27], however multiple filament formation is not explained. Due to the complex linear and nonlinear processes involved, filament properties are very difficult to predict. In this paper, we report on filament formation for different focusing conditions and mechanisms at work in each regime, as well as different laser energies for fixed laser spot sizes for 800 nm laser pulses. Similar studies for linear and nonlinear focusing regimes in air has been published by Lim et. al [28]. The numerical method in previous studies follows solving nonlinear Schrodinger equation (NLSE) [27] which is an asymptotic parabolic approximation of Maxwell's equations. Here we show the results of high resolution three-dimensional FDTD simulations, where the complete set of Maxwell's equations are solved. For the first time, we will show multiple filament formation in interaction of femtosecond laser pulses with the bulk of fused silica. We also improved the analytical model proposed by Lim et. al. [28] in filamentation regimes in air, to calculate the transition from nonlinear Kerr effect in loose laser focusing regime to geometrical focusing for tight laser focusing regime, which agrees with our numerical simulations. This paper is organized as follows: First we present our numerical model along with parameters used in our simulations. Then we investigate physical mechanism leading to filamention in fused silica for fixed laser energy, where we varied the laser spot size following by our improved analytical model and calculations corresponding to our simulation parameters. At last, we study the effect of the laser intensity (power) for fixed laser spot sizes.

II. NUMERICAL METHOD
In our 3D numerical simulation model, Maxwell's equations are solved using the FDTD method [29]. The Yee algorithm [30,31] is used in the code, using the constiutive relations (cgs units), H = B, D = (1 + 4π(ξ l + ξ k E 2 )E and current density J = J p + J P A , where E and B are the electromagnetic fields, D the displacement vector, H, the magnetic field auxiliary vector, ξ l is the linear susceptibility of the material, and ξ k is the Kerr susceptibility. The electromagnetic response of the generated plasma is represented by J p and laser depletion due to photo ionization (PI) by J P A . The evolution of the free electron density, n, is described as: dn dt = W P I (E), where W P I is the PI rate. We use Keldysh's formulation for the PI rate W P I [32]. Following Keldysh's Formulation [32], the adiabaticity parameter for solids γ = ω 0 m ef f U i /eE, where m ef f = 0.635m e denotes the reduced mass of the electron and the hole, U i is band gap energy, E is the laser electric field, and ω 0 is the laser frequency; γ becomes unity for the intensity of 3.5 × 10 13 W/cm 2 . For intensities considered in this paper, multi-photon ionization model is not a good approximation. We are using higher laser intensities than mentioned above and for few cases where lower intensities are used, the Kerr effect focuses the laser pulse, leading to higher intensities where the multiphoton rate model does not give the correct ionization rate. Therefore Keldysh ionization rate is used in all simulations discussed in this paper. We neglected avalanche ionization as we are using short (50 f s) laser pulse. Our simulation results including avalanche did not differ from simulations without avalanche ionization. In this paper we assume that a laser beam is focused by a perfectly reflecting parabolic mirror characterized by a given f #, corresponding to laser spot sizes of w 0 = 1.5, 1., 0.69 µm, and underfilled ratio of 0.5.The laser incident onto the mirror is considered to be a Gaussian beam. To describe the fields focused by the parabolic mirror, the Stratton-Chu integrals [33,34] are used, which specify the exact electromagnetic field emitted by the given parabolic surface. This field is calculated on the boundary of the simulation box and is used as the boundary conditions for the Maxwell solver of the parallel 3D FDTD code. For w 0 = 2 µm, the results of paraxial Gaussian beam and mirror focused beam were very close. Therefore paraxial Gaussian beam was used throughout this paper for w 0 = 2 µm to save computational resources and time. Fused silica parameters are as follows: refractive index 1.45, band gap energy of 9 eV , χ 3 = 1.9 × 10 −14 esu, saturation density is 10n cr , where n cr = meω 2 0 4π 2 , n cr is critical plasma density, and e is electron charge in cgs units. The laser pulse is Gaussian, linearly polarized along y direction and pulse duration is 100 f s. The geometrical laser focus is located at x = 40 µm. The laser wavelength is λ = 800 nm. Our simulation domain is 100µm × 16µm × 16µm, with grid size of ∆x = ∆y = ∆z = 0.02 µm. The laser pulse is propagating along x-direction.

III. RESULTS AND DISCUSSIONS
A. The effect of laser spot size for fixed laser energy In agreement with experimental results [8,12], we found that filamentation of femtosecond laser pulses in fused silica strongly depends on focusing conditions. First, we present our high resolution 3D FDTD simulation results. Figure 1 shows the contour plots of final electron densities in x−y plane for different laser focusing conditions: w 0 = 2.0, 1.5, 1., 0.69 µm which summarizes our results.
Here we show the simulation results when the laser energy is kept fixed (0.32µJ), as this is a usual practice in experimental studies, and initial spot size of the laser pulse is varied. This simply means, the smaller the spot size of the laser, the higher the input laser intensity. The magnitude of the change of refractive index ∆n corresponding to permanent damage in fused silica have been measured in the experiments in Ref. [27]. They found that the permanent damage happens in fused silica when n > 0.15n cr . The white areas in the contour plots of electron density corresponds to n > 0.15n cr , showing the permanent damage zone(s). The transition from one long filament for w 0 > 1.µm to void shape structure for tighter focusing conditions (w 0 < 1. µm) is clearly demonstrated in this figure. For w 0 = 2 µm ( Fig. 1-top panel), one long filament with multiple damage zones forms. Similar structure was reported by Sudrie et. al. [26]. As we decreased the spot size of the laser, the filament lengths became shorter (w 0 = 1.5, 1 µm, Fig. 1) and finally for very tight focusing w 0 = 0.69 µm one void shape structure with much higher plasma density forms, despite the fact that the intensity is higher for smaller laser spot sizes. From Figure 1, one can also observe that the first focus position happens farther from geometrical focus for looser laser focusing conditions and closer to geometrical focus as the spot size of the laser decreases. This means that self focusing and nonlinear Kerr effect play an important role in filamentation process. In contrary, for tighter focusing conditions, the position of the focus is very close to geometrical focus position. nonlinear Kerr effect is much more significant for loose focusing condition compared to tighter focusing regime.
To understand the physical mechanism, we performed simulations without nonlinear Kerr effect. Figure 2 shows the lineouts of the plasma density in the simulation domain at middle transverse plane for w 0 = 2, 1.5, 1, 0.69 µm, respectively. The solid lines, correspond to the simulations including the nonlinear Kerr effect and dashed lines correspond to the simulations with Kerr effect turned off. Figure 2 demonstrates the significant role of self-focusing and nonlinear Kerr effect for loose laser focusing conditions.
In simulation without the Kerr effect for w 0 = 2 µm (fig(2-top left)), the plasma density is well below the permanent damage threshold. We observed that a bulb  fig.2, bottom right). The peak density for w 0 = 0.69 µm seems sharp but it is really well resolved. Therefore, we conclude that the nonlinear Kerr effect plays a crucial role for the regime w 0 > 1µm. Contrary to w 0 > 1. µm, the nonlinear Kerr effect does not play an important role in the regime of tighter focusing conditions w 0 < 1 µm. The plasma shape in simulation with Kerr effect turned off for w 0 = 0.69 µm are very similar to the simulation results with Kerr effect included, the maximum plasma density drops by %50 (n e = 1.n cr ), however it is still well above the threshold for permanent damage. No other distinct differences is seen in fig. 2  To justify the role of Kerr effect, we calculated the global maximum over time of the on-axis laser intensity. Figure  3(a,b) shows the evolution of global maximum intensity of the laser pulse as it propagates through the medium for w 0 = 2., 0.69 µm. Figure 3(a,b) demonstrates a clear difference in the physical mechanism involved when moving from loose focusing conditions (w 0 > 1 µm) to tight focusing condition (w 0 < 1. µm). Figure 3 clearly illustrates the effect of Kerr nonlinear focusing of the laser pulse for w 0 = 2 µm. The initial growth of the field intensity is due to self-focusing. The original self-focused pulse loses its energy while it generates plasma, which defocuses the trailing part of the field. Once the plasma is not sufficient to counteract the effect of self-focusing, the trailing edge pulse focuses.
In contrast, the nonlinear Kerr focusing is not a primary mechanism for the tight focus condition as can be seen from Fig. 3-b. The laser focus position stays at the same location as the simulation without Kerr nonlinear terms. Also, the maximum laser intensity stays close to the maximum laser intensity when Kerr effect is turned off (Fig.3-right panel, dashed curve). In tight focusing regime, geometrical focusing is much more significant than the nonlinear Kerr effect. The Kerr effect has shorter propagation distance to build up and therefore the Kerr effect does not play an important role. In looser focusing conditions, the opposite is true. A more detailed picture of this dynamics can be obtained from the movies presented later in this report.
An extension to Marburger's formula [35] in a Kerr medium [36] stands out as a analytical expression to provide a good approximation of the position where filamen-tation begins: where z R is the Rayleigh length and P in is the laser power. We calculated the shift of nonlinear focus for the four focusing conditions in sample simulations. Figure 4 shows the shift of the nonlinear focus vs initial laser spot size calculated from Marburger's formula (black curve) and our simulation results (blue diamonds). Note that the laser energy is fixed. Marburger's formula prediction for the position of the first focus agrees well with our simulation results. Figure 4 clearly shows that the shift of the focus is small (1 − 2) µm for tight focusing conditions and as the laser spot size increases, the shift increases monoticaly. This suggests that for tight focusing conditions, the nonlinear Kerr effect diminishes and geometrical focusing is more significant. The laser focus position is very close to geometrical focus position and consequently, the plasma density structure is very similar to the simulation results with Kerr effect turned off. However, for larger laser spot sizes, the shift of focus position is significantly larger, the laser pulse focuses much farther from the geometrical focus of the laser pulse, and the intensity increases significantly, therefore the first plasma density above the threshold appears much farther from geometrical focus, therefore if the laser pulse power is above critical power, it has enough distance to re-focus and this process can happen several times.

B. Analytical Model
In this section, we demonstrate an extended analytical model previously used for filamentation in air, proposed by Lim et. al. [28]. They proposed an analytical method that defines the physical regime in which filamentation occurs in air. In this method, the wavefront sag(S), of a focusing Gaussian beam as a function of longitudinal position is calculated. For a focusing Gaussian beam the sag from geometrical focusing is [28]: where w 0 is the beam radius, f is the geometrical focus position and x R = πw 2 0 /λ 0 is the Rayleigh distance. The geometrical focus position is measured from one Rayleigh distance away from the position of geometrical focus. The sags from Kerr nonlinearity and plasma defocusing are calculated from optical path length differences between the center and the edge of the beam (1/e 2 of the intensity at the center): where for the nonlinear Kerr effect, ∆n(x) = n 2 I 0 (x). n 2 is the nonlinear refractive index and I 0 (x) = 2P 0 /πw(x) 2 is the peak laser intensity, P 0 is the peak laser power, and , where ρ is the on-axis plasma density and ρ c is the critical plasma density. Therefore, the sag from Kerr nonlinearity can be expressed as [28]: Lim [28] used multi-photon ionization model to calculate the sag from plasma defocusing. However, as we discussed in Sec. II , multi-photon ionization model is not a good approximation. Therefore, we used Keldysh ionization rate (See Sec.II ) and performed numerical integration to calculate the sag from plasma generation: Figure 5 shows the plot of the contributing sags for three cases discussed in simulation results (w 0 = 2, 1.5, 0.69 µm). The parameters of the laser pulse is the same as in simulations in previous section. We can see from the figure that for all three cases, when the beam is far from the geometrical focus, the intensity is low, and therefore the nonlinear Kerr effects are weak and therefore |S G |(blue) is much larger than |S K |(red) and |S p |(purple). The differences in the tight and loose focusing regimes, are as follows: In tighter focusing regime, the geometrical focusing and the plasma defocusing are the primary contributions.
The Kerr effect has a secondary effect. This is in agreement with our simulation results. From Fig. 2 (bottom right), the simulation results with and without Kerr effect show similar behavior. Figure 5-middle shows that as we move towards the looser focusing condition (w 0 = 1.5 µm), the effects of the geometrical focusing, plasma defocusing and Kerr effect become somewhat comparable. Simulation results with and without Kerr nonlinearity, show that the nonlinear Kerr is still the dominant mechanism. For even looser focusing conditions (w 0 = 2 µm), Fig.  5-left shows that the Kerr effect builds up faster than the plasma defocusing and geometrical focusing, showing that the Kerr effect plays a more important role as we move towards looser focusing regime. This also agrees well with our simulation results. Our simulation results with and without nonlinear Kerr effect are very different for looser focusing conditions (Fig. 2), demonstrating the dominant role of the Kerr effect.

C. The effect of laser energy for fixed laser spot sizes
In what follows we study how the input laser energy for fixed laser spot size (and fixed laser pulse length) affects the interaction with the bulk of the fused silica. Figure 6 shows the electron density contour plots in the simulation box, for w 0 = 0.69 µm. The peak initial laser intensities are: 5 × 10 13 , 1.1 × 10 14 , 3.3 × 10 14 , 4.9 × 10 14 W/cm 2 (Fig. 6 from top to bottom, respectively). The gen-eral feature of tight focusing condition, when the laser peak intensity is varied, is that the focus position of the laser does not move back considerably, with increasing the peak laser intensity. The plasma density generated for peak laser intensity of 4 × 10 13 W/cm 2 (not shown) was much below the permanent damage threshold. We found that the threshold for void formation happens for laser peak intensity of 5 × 10 13 W/cm 2 which leads to a ≈ 1 µm void-shape structure (6-top), with peak plasma density of 0.16n cr . Increasing the laser intensity to 1.1 × 10 14 W/cm 2 leads to a longer (≈ 2.9 µm) oval shape structure with maximum electron density 0.65n cr (Fig. 6). The damage area for laser peak intensity of 3.3×10 14 W/cm 2 is elongated (≈ 4.6 µm) and has a pear shape structure as shown in Fig. 6. Finally, increasing the peak laser intensity to 4.9 × 10 14 W/cm 2 leads to very similar structure as 3.3 × 10 14 W/cm 2 . We can see that saturation happens for laser intensities higher than 3.3 × 10 14 W/cm 2 for w = 0.69 µm. Figure 7 shows the on-axis lineouts of the electron densities corresponding to same laser intensities discussed which illustrates that positions of the damaged areas do not change significantly. It also shows that the damage structure stays almost the same for laser intensities larger than 1.1 × 10 14 W/cm 2 . The damaged structure shapes changes from void to elongated pear shape structure when increasing the laser intensity and the position of the damage areas are very similar in all cases as shown in Fig. 7. In addition, no distinct differences were observed in simulation results for tight focusing condition w 0 = 0.69 µm and laser intensities higher than 1.1 × 10 14 W/cm 2 , except that the plasma density grows with increasing of laser intensity. Next we study the effect of laser energy for looser focusing condition (w 0 = 2 µm, and fixed laser pulse duration). Figure 8 shows the simulation results for fixed laser spot size (w 0 = 2 µm) while we increase the laser energy (I = 2 × 10 13 , 4 × 10 13 , 5 × 10 13 , 6 × 10 14 W/cm2, from top to bottom). An area of very low electron density with maximum density of 0.06n cr is formed for laser peak intensity of 2 × 10 13 W/cm 2 . As can be seen also from on-axis electron density lineout (solid curve) in Fig.  9, the electron density is well below damage threshold density. The peak laser intensity of 4 × 10 13 W/cm 2 is the threshold for filament formation for laser duration of 50 f s and w 0 = 2 µm. The filament length (Fig. 8, Fig.  9) is about 4.2 µm and the first laser focus position happens at x = 31.8 µm. As we increase the laser intensity (energy), we can see that the first laser focus position is moving back (towards the left boundary), as expected. It is also clear that the length of the filament is increasing as we increase the laser energy. As we increase the laser energy to I = 5 × 10 13 W/cm 2 , we observe that the first focus position happens more closer to the left boundary (x = 20.5 µm) which leads to the first damage zone with a length of (≈ 1.9 µm) as can be seen from Fig 8. Two damage areas can be observed here. The electron density drops to 0.13n cr (which is lower than the damage threshold ) later due to absorption and ion- ization which causes energy loss of the front of the laser pulse before the second longer filament (9.4 µm) forms as the trailing part of the laser pulse focuses. This agrees well with experimental results of Ref. [8]. The numerical modeling of Ref. [26] could not explain the two damage zones which were shown in their experimental results.
Increasing the laser energy (I = 6 × 10 13 W/cm 2 ) leads to stronger laser focusing (Fig. 8). The first laser focus happens at x = 19.1 µm, and consequently larger plasma density is observed. The filament shape is very similar to Fig. Reffign except that the plasma density stays above the damage threshold after the head of the pulse loses energy.
Our key observation here is that, filament shape is not overly influenced by the peak input power beyond the saturation (I = 5 × 10 13 W/cm 2 ), however the focus position keeps moving back with increasing laser intensity.

IV. CONCLUSION
In this paper we studied filamentation mechanism in fused silica in details by performing extensive 3D high resolution FDTD simulations. We carried out simulations for both fixed laser energy where the laser spot size was varied and for fixed laser spot size where the laser energy was varied. In agreement to experimental results, we found that the filamentation strongly depends on initial laser focusing condition. Simulations without Kerr nonlinear terms showed that for loose focusing conditions, the Kerr nonlinearity plays a crucial role. The head of the filament moves back with increasing the laser spot size, meaning that the the laser focus is much farther than the geometrical focus position, and also the filament shape changes dramatically when moving to tight focusing conditions, where we observed void shape structure. In contrary, in tight focusing conditions, the Kerr effect does not play an import role. In this case, geometrical focusing is much more important than the Kerr effect. Simulations with Kerr effect turned off showed that the first focus position of the laser stays very close to the results of the simulations with Kerr effect included. The damage shape looks very similar with same length, only plasma density increases twice for simulations with Kerr effect. However, the plasma density in both simulations are well above the threshold for permanent damage. These results agree very well with experimental results of Refs[ [8,12]]. We also studied the effect of laser energy for fixed laser spot size. We have not observed distinct differences in simulations for tightly focused laser pulses when increasing the laser intensity (energy). The laser focus position stays very close to geometrical focus position. The void shape also looks very similar except that the plasma density increases with increasing the laser intensity. However, in al cases stayed well above permanent damage threshold. For loosely focused laser pulses, we found that the filament propagation is not terribly influenced by increasing the laser intensity (energy) beyond the threshold for filament formation. The first focus position moves back when the laser intensity was increased. However, the filament did not extend beyond 40 µm.