1. Introduction

Ultrafast laser filamentation in gasses, liquids, and solids has been an active research area since the advent of femtosecond lasers [1–4]. Filamentation in air has been widely investigated [5], with many applications ranging from atmospheric sensing [6,7] to terahertz generation [8]. While the critical power of self-focusing is in the gigawatt range for gasses, it is in the megawatt range for condensed matter due to their much higher nonlinear index, larger by up to three orders of magnitude [2]. Therefore, filamentation in liquids [9,10] and solids [11] is interesting not only for many applications, but also because the physical mechanisms of filamentation have proven to be largely universal, and thus can be studied on a reduced scale.

Intense femtosecond laser interaction with transparent media was investigated over two decades ago in the context of material damage thresholds [12,13] and for applications including fabrication of waveguides within bulk silica [14], gratings within optical fibers [15], three dimensional data storage [16,17], tissue ablation [18] and supercontinuum generation [19]; femtosecond filamentation within bulk dielectrics was also observed [20,21]. Femtosecond laser processing via filamentation remains an active field, with ever greater precision being achieved through spatio-temporal tailoring of the input beam, including the use of Bessel and Airy beams, chirped pulses and bursts [22–27]. Further, a new regime of filamentation in transparent media has been observed in the mid-IR, within the anomalous group velocity dispersion regime [28,29], leading to, for example, unprecedented multi-octave supercontinuum generation [30].

Despite this extensive body of work, there remain aspects of intense laser dielectric interaction and filamentation that are not understood. The majority of simulations are based on an equation for nonlinear pulse propagation with plasma generation and dynamics, where approximations such as slowly varying envelopes, or slowly evolving waves, among others, are applied to obtain an evolution equation for the pulse [31]. This approach, in general, has been highly successful when the correct approximations are considered for the problem at hand [11], where the focusing is not very tight, the medium remains relatively homogeneous, and the created plasma is not very dense. However, filament properties can sometimes be difficult to predict due to rich and complex dynamics leading to the formation of complex structures such as, for example, self-assembled nanogratings [32]. Thus, a more rigorous computational approach that makes no assumptions about light propagation is sometimes required [33–35], especially for modeling filaments created by tightly focused laser sources, where overdense plasmas can be created over short lengths and with small transverse profiles.

In this paper, we revisit filamentation in silica at 800 nm under tight focusing conditions by employing high resolution three-dimensional finite-difference time domain (FDTD) simulations, where the complete set of Maxwell’s equations is solved. We report on filament formation for various energies and various focusing conditions – from what we call intermediate tight focusing ($2~\mu m$ spot size) to tight focusing (sub-micron spot size), the latter achieved via a numerical parabolic mirror [36] – in order to detangle which mechanisms are at play in each regime. Experimental studies of single-pulse filamentation in transparent material in the past two decades have shown that filamentation strongly depends on laser focusing conditions, from voids for tight focusing [37,38] to long channels of modified refractive index for loose focusing [39,40]. Multiple filamentation has also been observed and successfully explained through nonlinear envelope propagation equations for loose focusing (with spot sizes of 10s of microns) [41,42]. Through our rigorous 3D approach, we are able to address filamentation in much tighter focusing regimes than is accessible to such propagation models, and are thereby able to shed new insight on the dynamics of the intense laser-matter interaction. This includes the formation of two distinct damage zones for intermediate tight focusing, similar to what was seen but not fully understood almost two decades ago [21,43].

This paper is organized as follows: First we present our numerical model along with parameters used in our simulations. Then we investigate the physical mechanism leading to filamentation in fused silica for different focusing conditions at fixed laser energy. This is followed by the presentation of an analytical model for filamentation in solids based on that proposed by Lim et. al. [44] for filamentation in air, to visualize the transition from Kerr self-focusing to geometrical focusing as laser focusing becomes tighter. Finally, we study the effect on filamentation of the laser intensity (power) for fixed laser spot sizes.

2. Numerical method

Maxwell’s equations are solved using the FDTD method [45] via the Yee algorithm [33,46] with the constiutive relations (cgs units) $\textbf {H}=\textbf {B}, ~ \textbf {D}=(1+4\pi (\xi _l+\xi _kE^2)\textbf {E}$ and current density $\textbf {J}=\textbf {J}_p+\textbf {J}_{PA}$. $\textbf {E}$ and $\textbf {B}$ are the electromagnetic fields, $\textbf {D}$ the displacement vector, $\textbf {H}$ the magnetic field auxiliary vector, $\xi _l$ the linear susceptibility of the material, and $\xi _k$ the Kerr susceptibility, which we take to be a constant. Since the group velocity dispersion (GVD) length in fused silica at $800~nm$ for $50~ fs$ pulses (what we consider here) is around $3.5~ cm$ [2], and since we are interested here in tighter focusing with orders of magnitude shorter interaction lengths, we neglect GVD in our model. The electromagnetic response of the generated plasma is represented by $\textbf {J}_p$ and laser depletion due to photo ionization (PI) by $\textbf {J}_{PA}$. The evolution of the free electron density, $n$, is described by $\frac {d n}{d t}=W_{PI}(E)$, where $W_{PI}$ is the PI rate. We use Keldysh’s formulation for WPI [47]; while there are more modern and accurate models of strong-field ionization in solids (see, for example, Ref [48]), the computational resources required for their implementation in a full 3D Maxwell solver is prohibitive. The adiabaticity parameter for solids is $\gamma =\omega _0\sqrt {m_{eff}U_i}/eE$, where $m_{eff}=0.635m_e$ is the reduced mass of the electron and the hole, $U_i$ is the band gap energy, $E$ the laser electric field, and $\omega _0$ the laser angular frequency; for our parameters, $\gamma$ becomes unity for the intensity of $3.5 \times 10^{13}~ W/cm^2$. While the ionization rates predicted by the Keldysh model and the multiphoton ionization model coincide at lower intensities, they begin to diverge by orders of magnitude for intensities higher than $\approx 5 \times 10^{12}~ W/cm^2$ (see, e.g., Fig. 3 of Ref. [21]). While a modified MPI rate can be used to model intensities up to $\approx 3\times 10^{13}~ W/cm^2$ as we have done previously [33], we consider intensities higher than that here, so we use the Keldysh model for all cases. While we have implemented a model of avalanche ionization that follows Ref [49], we find that our simulation results for 50 fs pulses including avalanche ionization were identical to equivalent simulations that did not include avalanche ionization. Thus, to save computational resources, we did not enable avalanche ionization in our code for the simulations presented here.

We assume a laser beam focused by a perfectly reflecting parabolic mirror characterized by a given $f\#$, corresponding to laser beam waist sizes of $w_0=1.5,1.0$, and $0.69$ $\mu m$. The laser incident onto the mirror is a Gaussian beam whose beam waist is half the size of the mirror. To describe the fields focused by the parabolic mirror, the Stratton-Chu integrals [36,50] are used, which specify the exact electromagnetic field emitted by the given parabolic surface. This field is calculated on five boundaries of the 3D FDTD simulation in a total field/scattered field framework. We also considered a laser beam with $w_0=2~\mu m$, where a paraxial Gaussian beam was used to save computational resources and time; we have verified that it gives almost exactly the same results as the mirror focused laser beam. The laser pulses are Gaussian in time with a pulse duration of 50 $fs$ and a wavelength of $\lambda = 800~ nm$ and they are linearly polarized along the $y$ direction and propagating along the $x$ direction. The geometrical laser focus is located at $x = 40~ \mu m$, and the simulation domain is $100 ~\mu m \times 16~ \mu m \times 16~ \mu m$, with grid size $\Delta x = \Delta y = \Delta z = 0.02 ~\mu m$. To ensure the domain was large enough, we chose it such that, in all simulations, the laser pulse was no longer creating plasma well before it exited the domain. Further, we ran simulations where we placed the geometric focus much deeper within the simulation domain and found no difference in our results. The background refractive index of silica is $1.45$, the band gap energy is $9 ~eV$, the third order nonlinear susceptibility is $\chi^{3}=1.9 \times 10^{-4}$ esu, and the saturation density is $10n_{cr}$, where $n_{cr}=\frac {m_e \omega _0^2}{4\pi ^2}$ is the critical plasma density.

3. Results and discussions

3.1 Effect of laser spot size

Figure 1 shows the contour plots of final electron densities in the $x-y$ plane for simulations performed for different laser focusing conditions corresponding to beam waist sizes of $w_0=2.0, 1.5, 1.0, 0.69 ~ \mu m$ for a fixed laser energy of $0.32 \mu J$. Thus, the smaller the spot size of the laser, the higher the laser intensity would be if the laser propagation were in vacuum. The magnitude of the change of refractive index $\Delta n$ corresponding to permanent damage in fused silica has been measured in the experiments in Ref. [43]. They found that permanent damage happens in fused silica when $n >0.15n_{cr}$. Thus, as a rough indication of the permanent damage zones predicted by our simulations, we indicate in white in Fig. 1 (and subsequent figures) the regions in the electron density contour plots corresponding to $n>0.15n_{cr}$.

The transition from long filamentation for $w_0>1.0 \mu m$ to a compact structure for tighter focusing conditions is clearly demonstrated in Fig. 1. For $w_0=2.0~\mu m$ (Fig. 1(a)), one long filament with multiple damage zones forms, similar in shape and size to the experimental structure reported by Sudrie et. al. [21,43], proposed at the time that this could be formed from pulse focusing/defocusing/refocusing, but were not able to capture it with their simulations. In Visualization 1, we present a time-domain movie for this simulation that allows us to visualize the dynamics of the laser intensity and the free electron density. We see that plasma formation begins at about $20~\mu m$, which is $20~\mu m$ before the geometrical focus at $40~\mu m$. As the simulation progresses, the leading half of the pulse experiences a small drop in intensity due to defocusing that leads to a decrease in plasma generation. The tail half of the pulse, however, visibility deforms around this newly formed plasma due to plasma defocusing. After this, the leading half of the pulse has enough intensity that it refocusses, ultimately creating a longer filament that persists up to and beyond the geometric focus. The tail half of the pulse eventually refocusses as well, though at much lower intensity, and contributes to extending the plasma tail to $60 ~\mu m$, however, at much lower electron density levels than would be needed for permanent damage.

 

Fig. 1. Contour plots of electron density after the laser pulse (a-d) for beam waists of $w_0=2.0, 1.5,1.0$, and $0.69~\mu m$ for a fixed laser energy of $0.32~\mu J$. (see Visualization 1,Visualization 2,Visualization 3 for time-domain movies of laser intensity and free electron density simulations with $w_0=2.0 ,1.0,0.69~\mu m$, respectively.)

Download Full Size | PDF

Figure 1 shows that as we decrease the spot size of the laser, the filament length shortens, ultimately becoming a compact structure with much higher electron density for very tight focusing ($w_0 = 0.69~ \mu m$). One also observes in Fig. 1 that for looser focusing, plasma creation begins sooner, much before the geometrical focus, despite the fact that the incident laser intensity is smaller for larger spot sizes (since we keep the total energy constant). This indicates that self focusing and the nonlinear Kerr effect must play an important role in the longer filamentation regime. In contrast, for tighter focusing, the position of the focus is very close to the geometrical focus, indicating a potentially smaller contribution of the nonlinear Kerr effect.

Visualization 2,Visualization 3 show the time domain movies of filament formation for $w_0 =1.0$ and $0.69~\mu m$, respectively. In Visualization 2, we see that plasma formation does not begin until around $27~\mu m$, and that tail half of the pulse again experiences visible plasma defocusing, while the leading edge continues to create plasma until the geometric focus, after which it rapidly defocusses. The tail half does not refocus and thus does not extend the plasma tail, as it did in Visualization 1. In Visualization 3, the tightest focusing we consider, geometric focusing is clearly dominant. The plasma is confined to near the geometric focus, and is rapidly formed. Plasma defocusing is also seen, but since geometric defocusing is so strong (after the geometrical focus), it dominates over Kerr self-focusing, and thus no refocusing is observed.

To better understand the role of Kerr self-focusing, we performed simulations equivalent to those of Fig. 1, except that we have turned off the nonlinear Kerr effect by setting the Kerr susceptibility to zero. Figure 2 a-d show the final plasma density along the laser axis for $w_0=2, 1.5, 1.0, 0.69 ~\mu m$, respectively. The solid lines correspond to the simulations including the nonlinear Kerr effect (as in Fig. 1) and dashed lines correspond to the simulations with Kerr effect turned off.

We find that the Kerr effect plays a significant role in all four cases. However, in Fig. 2(a) and (b), corresponding to $w_0 = 2.0$ and $1.5 ~ \mu m$, respectively, neglecting the Kerr effect causes the plasma density to be far below the threshold for permanent damage. In Fig. 2(c), corresponding to $w_0 = 1.0~\mu m$, we find that the electron density does approach a value that could result in permanent damage, but that the damage spot would be much smaller in the absence of Kerr self-focusing. In contrast, in Fig. 2(d), for the tightest focusing with $w_0 = 0.69~ \mu m$, we see that while there is still a difference between the two curves, the effective shape and length of permanent damage area are very similar. However, the electron density is doubled when the Kerr effect is included, indicating that it plays a key role in the onset of micro-explosions and void formation.

 

Fig. 2. On-axis final electron density for $w_0=2~\mu m$ (a), $w_0=1.5~\mu m$ (b), $w_0=1~\mu m$ (c), and $w_0=0.69~\mu m$ (d). The solid lines correspond to the simulations of Fig. 1, whereas the dashed lines correspond to equivalent simulations without the Kerr effect.

Download Full Size | PDF

To further investigate the role of the nonlinear Kerr effect, we extracted from our simulations the maximum value that the light intensity reached along the laser axis at each time step. We call this the "global maximum", and we plot this as a function of time in Fig. 3 for $w_0= 2.0~\mu m$ (a) and $0.69~\mu m$ (b). As before, the solid line indicates results corresponding to the simulations of Fig. 1, and the dashed lines the equivalent simulations where the nonlinear Kerr effect is turned off. In the solid curve of Fig. 3(a), we see an initial growth of the field intensity due to self-focusing, which then oscillates around a saturation level after significant plasma is produced. Here we see an interplay between plasma defocusing and nonlinear Kerr self-focusing which ultimately creates the multiple filamentation spots we observe in Fig. 1. In the absence of the nonlinear Kerr effect (dashed curve), this saturation and oscillation behaviour is completely absent.

 

Fig. 3. Global on-axis maximum intensities as a function of time calculated from simulations for $w_0=2.0$ (a) and $0.69~\mu m$ (b). Solid and dashed curves correspond to simulations with and without Kerr nonlinearity, respectively.

Download Full Size | PDF

In contrast, we see from Fig. 3(b), that for tight focusing, the difference between the simulations with and without the Kerr effect is not as dramatic, and the the light intensity reaches its maximum at the same time for both simulations and follows the same shape. The deviation between the actual values of the maximum shows that the Kerr effect still does play an important role, however, and is thus likely relevant in determining the nature of the damaged area. The size and shape of this area, however, are primarily determined by geometrical focusing.

3.2 Analytical model

In this section, we demonstrate more visually the relative importance of Kerr self-focusing versus geometric focusing and how these interact with plasma defocusing by extending to solids, an analytical model previously developed by Lim et al [44]. They used a quantity called wavefront sag $(S)$, which is the path difference between the center and the edge of the wavefront. For a focusing Gaussian beam, the sag from geometrical focusing is [44],

The wavelength sags from the Kerr nonlinearity and plasma defocusing are obtained from calculating the optical path length differences between the wavefronts along the laser axis between the center and the edge of the beam ($1/e^2$ of the intensity at the center),

Lim et. al [44] used a multi-photon ionization model to calculate the sag from plasma defocusing. However, as we discussed in Sec. 2. , the multi-photon ionization model is not appropriate here, thus we use instead the Keldysh ionization rate to obtain,

Figure 4 shows the plot of the contributing wavelength sags for three simulations of Fig. 1 ($w_0=0.69, 1.5, 2.0~\mu m$). For all three cases, when the beam is far from the geometrical focus, the intensity is low, thus the nonlinear Kerr effect and plasma defocusing are weak and $|S_G|$(blue) is much larger than $|S_K|$(red) and $|S_p|$(purple). The position $x_k$ is defined as the position where $|S_g|=|S_K|$, where Kerr self focusing becomes non-negligible compared to geometrical focusing. Similarly, $x_p$ is the position where $|S_g|=|S_p|$, where plasma defocusing becomes non-negligible. In the tightest focusing regime (4a), we see that $x_k > x_p$. This indicates that the plasma defocusing builds up sooner than the Kerr effect, and suggests that plasma defocusing and geometrical focusing are the primary contributions. As the pulse rapidly comes into focus (due to geometric focusing), significant ionization (and thus plasma defocusing) immediately begin to occur once threshold is reached, whereas the Kerr effect takes some distance to build up once its threshold is reached. For $w_0=1.5~\mu m$ (Fig. 4(b)), $x_k$ is closer to $x_p$, therefore the Kerr nonlinearity and plasma defocusing are comparable and both play important roles. However, in the loosest focusing regime ($w_0=2~\mu m$, Fig. 4(c)), the Kerr effect builds up faster than the plasma defocusing and geometrical focusing ($x_p>x_k$), therefore the Kerr nonlinearity plays the primary role. This is in agreement with our simulation results in the previous subsection.

 

Fig. 4. Illustration of calculations of absolute values of sag contribution from Kerr self focusing ($S_K$, red), plasma defocusing ($S_P$, purple) and geometrical focusing ($S_G$, blue) as a function of normalized distance $(x/f)$ for focusing conditions $w_0=0.69,1.5,2~\mu m$, (a-c) respectively.

Download Full Size | PDF

3.3 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 duration) affects the interaction with bulk fused silica. Figure 5 shows electron density contour plots after the laser pulse for $w_0=0.69~\mu m$, for peak incident laser intensities of: $5\times 10^{13}, 1.1 \times 10^{14}, 3.3 \times 10^{14}, 4.9\times 10^{14}~ W/cm^2$ (Fig. 5(a)-d, respectively). We observe that when the laser peak intensity is varied, the focus position does not change considerably. We found that the threshold for permanent damage happens for an incident laser peak intensity of $5 \times 10^{13}~ W/cm^2$ which leads to a plasma size of $\approx 1~\mu m$, with peak plasma density $0.16 n_{cr}$. Since here geometrical focusing is very strong, there is not as much time for Kerr self-focusing to build as compared to looser focusing. Thus, a higher incident intensity is required to reach the damage threshold in the laser interaction region for very tight focusing than would be required for looser focusing. Increasing the laser intensity to $1.1 \times 10^{14}~W/cm^2$ leads to a longer ($\approx 2.9~\mu m$) oval shape structure with maximum electron density $0.65n_{cr}$. The damage area for laser peak intensity of $3.3\times 10^{14}~W/cm^2$ is elongated ($\approx 4.6~\mu m$) and has a pear shape structure the same as increasing the peak laser intensity to $4.9\times 10^{14}~W/cm^2$ leads to very similar structure as $3.3\times 10^{14}~W/cm^2$.

Figure 6 shows the on-axis values of the electron densities corresponding to the simulations of Fig. 5. While the plasma shape elongates as the intensity is increased, and the plasma density increases, we see a saturation in the electron density for $3.3 \times 10^{14}~W/cm^2$ and above.

 

Fig. 5. Contour plots of electron density for $w_0=0.69~\mu m$ for intensities of $I=5e13$ (a), $1.1e14$ (b), $3.3e14$ (c), and $4.9e14~W/cm^2$ (d).

Download Full Size | PDF

 

Fig. 6. On-axis electron density along the laser propagation direction $x(\mu m)$ from simulations for $w_0=0.69~\mu m$ for $I= 5e13$ (solid line), $1.1e14$ (dashed line), $3.3e14$ (dashed-dotted line), and $4.9e14$ (dotted line) $~W/cm^2$.

Download Full Size | PDF

Figure 7(a-d) shows the electron density distribution for a larger laser spot size ($w_0=2~\mu m$) for varying laser pulse energy, with peak intensities of $I=2\times 10^{13}, 4\times 10^{13},5\times 10^{13}$, and $6\times 10^{14}~W/cm^2$, respectively. The corresponding electron densities along the laser axis are plotted in Fig. 8.

 

Fig. 7. Contour plots of electron density for $w_0=2~\mu m$ for $I=2e13$ (a), $4e13$ (b), $5e13$ (c), and $6e13$ (d) $~W/cm^2$.

Download Full Size | PDF

 

Fig. 8. On-axis electron density along propagation direction $x(\mu m)$ from simulations for $w_0=2~\mu m$ for $I=2e13$ (solid line), $4e13$ (dashed line), $5e13$ (dashed-dotted line), and $6e13$ (dotted) $~W/cm^2$.

Download Full Size | PDF

An area of very low electron density with maximum density $0.06n_{cr}$ is formed for $I=2 \times 10^{13}~W/cm^2$, which is well below damage threshold density. $I=4\times 10^{13}~W/cm^2$ is the threshold for filament formation, with a filament length of around $4.2~\mu m$. As we increase the laser energy (and thus the laser intensity), the location where plasma is first created moves backward, as expected, and the length of total damaged area increases. Further, as intensity increases, we see the emergence of a prominent peak on the left, which gets larger with increasing intensity relative to those on the right. This trend suggests that as laser energy increases, we would expect two different type of damage zones, with the left zone being smaller and containing more structural damage, and the right zone much more elongated, perhaps reaching the threshold for permanent refractive index change but not structural damage. This agrees well with the observations described in Refs. [21,43]

4. Conclusion

Through high resolution 3D FDTD modelling, we unraveled the mechanisms underlying intense laser interaction and filamentation in fused silica at 800 nm in the tight focusing regime, bringing understanding to the transition between linear and nonlinear focusing regimes, and new insight into experiments going back two decades. Our simulations show the creation of two damage zones similar to those observed experimentally by Sudrie et al. [21] under tighter focusing conditions and find, perhaps not surprisingly, that Kerr self-focusing plays a key role in the shape of the created plasma. However, nonlinear propagation equation models had been unable to capture the creation of the two damage areas, despite their success at modelling multiple filamentation for much looser focusing. Moreover, we find that these two damage zones have different levels of carrier density for high enough laser energy, a possible explanation for the two very different types of damage zones that were observed experimentally, i.e., a small zone with obvious structural damage followed by a longer lower-damage zone. For the tightest laser focusing we considered (sub-micron spot size), we find that the compact shape of the created plasma is determined only by geometrical (i.e., linear) focusing. Kerr self-focusing, in contrast, affects only the maximum electron density attained within the plasma, signifying a possible role in the onset of structural damage and void creation through mirco-explosions. Kerr self-focusing thus plays a very different role depending on the degree of tight focusing, meaning that the concept that the physical mechanisms underlying filamentation are universal does not always apply.