1. Introduction

Femtosecond laser processing is widely used in both industry and academia, especially for micro-processing applications. Compared to traditional laser processing, femtosecond laser processing is characterized by minimal thermal damage to the surrounding non-processed material; thus, it is believed to be particularly favorable for processing materials sensitive to thermal effects [1,2]. This "non-thermal" nature of femtosecond laser processing is attributed to the short pulse duration of such lasers, which is much shorter than typical heat diffusion timescales [3,4]. Recent progress in laser technologies has enabled the stable use of femtosecond lasers with even higher intensity and shorter pulse durations, typically less than 100 fs. Utilizing such lasers has begun attracting attention as a way to achieve even higher quality and precision in processing results [5–15].

A major challenge in applying sub-100 fs laser pulses to laser processing is the difficulty in controlling the laser light itself due to the many nonlinear phenomena caused by a high peak intensity. One such nonlinear effect is the distortion of the beam caused by nonlinear propagation in air, which involves various phenomena, such as self-focusing due to third-order nonlinear effects, air ionization, and plasma absorption and defocusing. Due to its complexity, it is generally difficult to quantitatively determine this propagation process. Consequently, nonlinear propagation effects remain a major factor limiting the pulse energy used for industrial laser processing in air, where high controllability is required [16]. More specifically, for sub-100 fs pulses, the temporal and spatial characteristics of the beam change significantly due to the aforementioned nonlinear propagation effects when the fluence of the laser pulse is set to several times the material ablation threshold (several to several tens of J/cm$^2$), as is typical in traditional laser processing. This leads to, for example, changes in the pulse energy dependence of the squared diameter of the ablated crater. The diameter-squared ($D^2$) method (or Liu’s method) is widely used for predicting the squared diameter of an ablated crater created with a Gaussian pulse in the linear propagation regime [17]. The experimental squared diameter agrees well with the behavior predicted by the $D^2$ method at low fluences (i.e., in the linear propagation regime), but fails at higher fluences (i.e., in the nonlinear propagation regime) [12–15].

In light of the above, it is clear that correctly simulating the nonlinear propagation of laser pulses in air and properly relating the incident fluence distribution to the ablation results are both important to achieve quantitative predictions of laser processing results in the nonlinear propagation regime. For the former, numerous theoretical studies have been performed on the nonlinear propagation dynamics of femtosecond laser pulses when focused in air, mainly to study filamentation phenomena [18–23]. For focal lengths typically used in laser processing, ranging from tens to hundreds of millimeters, it is known that plasma defocusing is the main factor that initially distorts the beam during propagation [20]. The effect of plasma defocusing is in turn dependent on the spatial density distribution of free electrons produced by air ionization. Consequently, the air ionization model employed in simulations is a crucial factor affecting model accuracy. Several models have been proposed to calculate the intensity dependence of the ionization rate of gases, including the multiphoton ionization (MPI) model [24], the Perelomov-Popov-Terent’ev (PPT) model [25,26], and the Ammosov-Delone-Krainov (ADK) model based on tunnel ionization (TI) [27]. In previous studies, each model has been used to simulate the nonlinear propagation of laser pulses within a range of intensities for which they are valid [18–21]. For the latter, it has been reported in the linear propagation regime that the ablated area created by laser ablation is determined by the area where the local fluence of the incident laser is higher than the material laser ablation threshold [28–31]. Furthermore, it has recently been shown that this property might hold even at pulse energies where nonlinear propagation occurs [12]. This suggests that if the fluence distribution of the beam could be determined, it may be possible to quantitatively predict the ablated area. As it is difficult to directly measure the beam profile in the nonlinear propagation regime [11,12], such determination would be better suited for numerical simulations.

In fact, some studies have compared the nonlinear propagation simulations results of femtosecond laser pulses in air with the ablation results of metals processed by such nonlinear propagating laser pulses [32–34]. However, these studies experimented with only a single laser pulse energy; they have not established a framework with which they could quantitatively explain ablation results over a wider range of pulse energies typically used in laser processing.

In this study, we present a method to quantitatively predict the pulse-energy dependence of ablation crater shapes for various metals when ablated by sub-100 fs pulses propagating in ambient air, applicable over a two-orders-of-magnitude range in the pulse energy. First, nonlinear propagation simulations, using the PPT model [26] as the gas ionization model, were performed at various pulse energies to calculate the fluence distribution at each condition. Then, we performed laser ablation experiments on copper under the same conditions as in the numerical simulations, followed by a comparison of the numerical simulation results with the measured ablation diameter. We found that the ablation diameter predicted by the numerical simulation agreed very well with the experimental results over a two-orders-of-magnitude range in the pulse energy. We furthermore repeated the method using gold and aluminum samples, and confirmed the method’s effectiveness regardless of the specific metal. Finally, we analyzed the pulse energy dependence of the fluence at the center of the laser pulse as well as its relationship to the ablation depth, and found quantitative agreement in their tendencies. Overall, these results indicate that our method can quantitatively predict the processing results over a wide pulse energy range and should significantly improve the controllability of laser processing, especially in the nonlinear propagation regime.

This paper is organized as follows: the nonlinear propagation model used in the numerical simulation is described in Sec. 2.1, the experimental method is shown in Sec. 2.2, then the numerical simulation results are compared with the experimental results and followed by a discussion in Sec. 3.

2. Methods

2.1 Numerical simulation

To investigate the beam propagation of a femtosecond laser pulse in air, we performed numerical simulations using the extended nonlinear Schrödinger (NLS) equation. The equation describes the propagation of the envelope function of a linearly polarized electric field. Furthermore, the NLS equation was coupled to a rate equation describing the free electron density created by ionization of nitrogen and oxygen molecules in air. The two equations were as follows [22,23,35]:

Here, $z$ represents the propagation distance, $t$ represents the time, $E$ represents the electric field, $\omega$ represents the angular frequency of the laser at the center wavelength, $\lambda$, $n_{\omega }$ represents the refractive index of air, $\epsilon _0$ represents the dielectric constant of vacuum, $c$ represents the speed of light, $k_{\omega }=2\pi n_{\omega }$/$\lambda$ represents the wave number, $k_{\omega }'=1/v_g$ (comprising $v_g$ as the group velocity) represents the first-order dispersion, $k_{\omega }''$ represents the group velocity dispersion, $N_{\mathrm {mol}}$ represents the neutral molecule density, $N_p$ and $N_p^{e}$ represent the number density of molecules and electrons, respectively ($N_{\mathrm {N_2}}=0.8 \times N_{\mathrm {mol}}$, $N_{\mathrm {O_2}}=0.2 \times N_{\mathrm {mol}}$), $\chi ^{(3)}$ represents the third-order nonlinear susceptibility, $\tau _c$ represents the electron collision time, $\sigma _{\omega }=e^2/(\epsilon _0 k_{\omega } m_e c^2 \omega \tau _c)$ represents the cross section for inverse Bremsstrahlung (comprising $e$ as the elementary charge, $m_e$ as the electron mass), $W_{\mathrm {PPT}} (E)$ represents the photoionization rate of gas molecules, and $U_p$ represents the ionization potential of the gas molecules. Table 1 summarizes the values for each parameter. Here, we used experimental values of $\chi ^{(3)}$ reported under similar conditions to our simulation conditions [36]. We note that each term on the right-hand side of Eq. (1) corresponds to the radial diffraction, first-order dispersion, group velocity dispersion, self-focusing, nonlinear absorption by photoionization of the gas molecules, plasma absorption, and plasma defocusing processes, while the first and second terms on the right-hand side of Eq. (2) correspond to photoionization and avalanche ionization processes. To calculate the photoionization rate, we used a semi-empirical model [26] which extends the original PPT model [25]. The original PPT model was derived for monatomic gases, and is not readily applicable to other gases. However, for diatomic molecules such as N$_2$ and O$_2$, it was shown that the PPT model could be applied by experimentally determining an effective charge $Z_{\mathrm {eff}}$ and correcting for the effective Coulomb potential felt by the excited electron [26]. We take the same approach here, and calculate the ionization rate from the PPT model (see Appendix A for details) using the values of $U_p$ and $Z_{\mathrm {eff}}$ given in Ref. [26] (N$_2$: $U_{\mathrm {N_2}}$ =15.576 eV, $Z_{\mathrm {eff}}$=0.9, O$_2$: $U_{\mathrm {O_2}}$ =12.55 eV, $Z_{\mathrm {eff}}$=0.53). We note that these values have previously been validated in several papers [32,39,40]. Importantly, the functional form of $W_{\mathrm {PPT}}(E)$ was directly incorporated into the NLS equations without approximation, in contrast to previous works [18,21,40]. We note that all parameters used in Eqs. (1) and (2) were taken from literature values, and there are no free parameters.

Table 1. The parameter values used in Eq. (1) and (2).

View Table

We calculated the above equations using the Split-Step-Fourier method [41] in a radial coordinate system. All simulations were performed for an initially (i.e., at low energies in the linear propagation regime) spatio-temporal Gaussian beam with a central wavelength of 800 nm, pulse duration of 43 fs, and a beam waist of 10.8 $\mathrm{\mu}$m at the focus (Rayleigh length $z_{\mathrm {R}}$=458 $\mathrm{\mu}$m). The focal point was taken as $z$=0, and the laser pulse was propagated in the range of $z$=$-$20 to 6 mm. The step size was taken as 0.55 fs in time, 9.2 $\mathrm{\mu}$m in the propagation direction ($z$-direction), and 0.76 $\mathrm{\mu}$m in the radial direction. Under these conditions, simulations were performed for varying pulse energies in the range of 2.0 to 35.0 $\mathrm{\mu}$J, or peak fluence values between 1.1 to 19.1 J/cm$^2$ at $z$=0, assuming linear propagation. During analysis, the local fluence at each position $(r,z)$ was calculated by integrating the local intensity at $r$ with respect to time. Furthermore, the center fluence was defined as the local fluence at position $r$=0, while the beam radius was defined as the distance $r$ at which the local fluence decreases to 1/$e^2$ the value of the center fluence.

We have conducted the above calculations on a computer with a single CPU (Core i9-10900K, Intel Corp.) and 32 GB memory. It took about 40 minutes to calculate one condition on a single core with this computer.

2.2 Experiment

Figure 1 shows a schematic diagram of the experimental setup. The light source was a Ti:Sapphire regenerative amplifier (Astrella, Coherent Inc.) with a center wavelength of 800 nm, repetition rate of 1 kHz, and average power of 6 W. After the laser power was reduced by a beam sampler, the beam radius was reduced to approximately 1.5 mm by a reduction optics system comprising a plano-convex lens of $f$ = 200 mm and a plano-concave lens of $f$ = $-$50 mm. The laser repetition rate was then divided into 100 Hz using an optical chopper. From this 100 Hz train, an arbitrary number of pulses could be picked by a fast mechanical shutter (Uniblitz LS6SZ1, Vincent Associates Inc.). The pulse energy could be varied from 2.0 to 34.0 $\mathrm{\mu}$J with an attenuator comprising a half-wave plate and a reflective polarizer. The actual pulse energy of the laser pulses was monitored in real-time by a photodetector pre-calibrated to the readings from a commercially available power meter (XLP12-3S-VP-D0, Gentec-EO) placed just before the laser ablation position. A thin plano-convex lens with a focal length of $f$ = 75 mm was used to focus the laser pulse onto the sample. Copper (purity: 99.96%, Nilaco Corp.), gold (purity: 99.95%, Nilaco Corp.), and aluminum (purity: 99%, Nilaco Corp.) plates of 500 $\mathrm{\mu}$m thickness were used as samples for laser ablation experiments. The sample was mounted on a 3-axis automated stage together with a beam profiler, which, with further beam attenuation by ND filters, could be used to directly measure the beam profile at the focal point [31]. This was important for correct simulations because it allowed for highly accurate measurements of the beam radius at the focal point in the linear propagation regime. The beam profile at the focus was measured to have a 1/$e^2$ beam radius of 10.8 $\pm$ 0.1 $\mathrm{\mu}$m. We adjusted the compressor inside the regenerative amplifier to give a negative pre-chirp so that the pulse duration was minimum before the sample. The pulse duration (Full Width at Half Maximum, assuming a Gaussian pulse) before the final focusing lens was measured by autocorrelation (GECO, Light Conversion), and considering the slight dispersion due to transmission through the final focusing lens, the pulse duration before the sample was 43 fs.

 

Fig. 1. Schematic of the experimental setup. The solid red line is the optical path for the laser ablation. After controlling the pulse energy with an attenuator, the laser pulse is focused onto a sample on an XYZ stage. BD: beam damper, PD: photodetector.

Download Full Size | PDF

In the experiments, five craters were ablated at each pulse energy. To avoid cross-contamination, each crater was fabricated on a new surface of the sample separated 100 $\mathrm{\mu}$m apart. The ablated craters were observed with a laser scanning microscope (VK-X1000, Keyence Corp.), which is capable of simultaneously measuring both a laser scanning microscope and optical microscope image. From the obtained laser scanning images, the ablation diameter and depth were determined as detailed in Appendix B.

3. Results and discussion

First, we discuss the effects of nonlinear propagation on the fluence distribution at $z$=0, and its qualitative effects on the ablated crater profile. Figure 2(a) shows the calculated fluence distributions for beam propagations with pulse energies of 8.0, 18.0, 26.0, and 32.0 $\mathrm{\mu}$J. It can be seen that the position of the smallest beam radius shifts towards the direction of the incident beam as the pulse energy increases. When the simulations were repeated with the terms in the NLS equations related to self-focusing and the interaction between the photoelectric field and the plasma systematically eliminated, it was found that the interaction between the electric field and the plasma had a more significant role in altering the propagation behavior (see Appendix C for details). Therefore, we concluded that the aforementioned shift in the focus position was mainly due to the effect of plasma defocusing caused by air ionization during the focusing of the laser beam. Figure 2(b) shows the fluence distributions at $z$=0, corresponding to the position of the dashed line in Fig. 2(a) for each pulse energy. While the center fluence initially increased as expected, at above 18.0 $\mathrm{\mu}$J, the fluence is seen to saturate and, eventually, decrease counterintuitively. This decrease was attributed to the increased beam divergence caused by the stronger plasma defocusing effect at higher intensities. Correspondingly, the beam radius is also larger. We also observed that while the fluence profile is close to a Gaussian at 8.0 and 18.0 $\mathrm{\mu}$J, nonlinear perturbances changed the distribution to a top-hat-like shape at 26.0 and 32.0 $\mathrm{\mu}$J. This top-hat-like shape at higher pulse energies is thought to be the result of a balance between the stronger plasma defocusing effects near the center of the beam and the continued focusing effects at the edge of the beam.

 

Fig. 2. (a) Simulated fluence distributions cropped from $-$1100 to 750 $\mathrm{\mu}$m for an 800 nm, 43 fs pulse at four different pulse energies. The white-dashed line at $z$=0 represents the position of the focal plane, which corresponds to the Gaussian beam waist position at low pulse energy. (b) The fluence distributions of the laser pulses at $z$=0 in (a).

Download Full Size | PDF

Figure 3 shows the laser ablation results of copper irradiated with 75 pulses, each at pulse energies of 7.7, 17.9, 25.8, and 31.9 $\mathrm{\mu}$J, similar to the above calculations. The optical microscope images of the ablated craters are shown in Figs. 3(a)-(d). The dark central region corresponds to the ablated crater, and it can be seen to become increasingly larger at higher pulse energies. Figure 3(e) shows the cross-sectional height profiles of the ablated craters obtained by the laser scanning microscope, corresponding to the cross-section along the dashed line in Figs. 3(a)-(d). Importantly, while the depth of the crater decreased above 17.9 $\mathrm{\mu}$J, the diameter of the crater continued to increase with increasing pulse energy. This trend of the crater diameter and depth change with pulse energy is consistent with the trend expected from the pulse energy dependence of the fluence distribution, as previously discussed with Fig. 2(b). Furthermore, the shapes of the cross-section of the ablated craters showed a parabolic shape for pulse energies of 7.7 and 17.9 $\mathrm{\mu}$J, while a flat bottom shape was observed for pulse energies of 25.8 and 31.9 $\mathrm{\mu}$J, indicating that the shape of the ablated crater corresponded well with the simulated beam profile.

 

Fig. 3. The ablation results of copper when irradiated with 75 pulses at 800 nm, 43 fs, and 100 Hz. Optical microscope images with (a) 7.7 $\mathrm{\mu}$J, (b) 17.9 $\mathrm{\mu}$J, (c) 25.8 $\mathrm{\mu}$J, and (d) 31.9 $\mathrm{\mu}$J are shown. (e) Cross-sectional height profiles obtained by the laser scanning microscope, corresponding to the cross-section along the dashed line in (a)-(d).

Download Full Size | PDF

Next, we focus on the pulse energy dependence of the ablation diameter for copper. The pulse energy dependence of the square of the measured ablation diameter is shown by the red points in Figs. 4(a),(b). Each data point represents the geometric mean of the major and minor axes lengths of an ellipse fitted to the ablated area (see Appendix B for details).

 

Fig. 4. Comparison of the simulation results with the experimental results of copper (red points) for the pulse energy dependence of the ablation diameter squared on (a) a linear scale and (b) a logarithmic scale. Five experimental data are plotted at each pulse energy. The blue line represents the results with the PPT model, the orange line represents the results with the MPI model, the green line represents the results with the PL approx. model, and the purple line represents the results with the TI model. The black-dashed line represents the fitting result of the experimental data (<10 $\mathrm{\mu}$J) using the $D^2$ method. See main text for the meanings of the abbreviations of the ionization models used in the simulations.

Download Full Size | PDF

In the linear propagation regime where the fluence distribution of the laser pulse is Gaussian, the ablation diameter $D$ is known to satisfy the following relation [17]:

On the other hand, as shown in Figs. 4(a),(b), when the pulse energy exceeds approximately 18.0 $\mathrm{\mu}$J, the measured ablation diameter deviates from the prediction of the $D^2$ method due to nonlinear propagation effects. To determine the ablation diameter from the fluence distribution obtained at these higher fluences, we assumed, as in the $D^2$ method, that ablation would occur at a diameter where the local fluence of the incident laser pulse is equal to the threshold fluence (see Appendix D for details). Here, two parameters need to be determined to correlate the simulated ablation diameter to the measured ablation diameter: one is the value of the material-specific threshold fluence $F_{\mathrm {th}}$, and the other is the $z$-position at which the sample was placed in the experiment. Adopting these values as fitting parameters, we performed least-squares fitting of the measured ablation diameters over the whole pulse energy range of the experiment. The resulting values of the threshold fluence and $z$-position were $F_{\mathrm {th}}$ = 0.47 $\pm$ 0.01 J/cm$^2$ and $z$ = 4 $\pm$ 4 $\mathrm{\mu}$m, respectively. The value of the threshold fluence agreed with the aforementioned results obtained with the $D^2$-method within the range of fitting error. We estimated the allowable experimental error for the $z$-position as $\pm$60 $\mathrm{\mu}$m, based on measurements of the tilt of the beam profiler’s sensor surface, the sample, and the sample holder relative to the $z$ axis, and find the results within this range as well. Altogether, we believed that the values of the threshold fluence and $z$-position obtained by fitting were reasonable. Figures 4(a),(b) shows the pulse energy dependence of the square of the predicted ablation diameter calculated using these values as blue lines. In all pulse energy regions, the predicted ablation diameter is in excellent agreement with the measured ablation diameter. Moreover, statistical analysis of the goodness of fit between predicted ablation diameter and measured ablation diameter also yielded good agreement (see Appendix E for details). Overall, these results show that the present method can quantitatively predict the pulse energy dependence of the ablation diameter over a two-orders-of-magnitude range in the pulse energy, including the nonlinear propagation regime.

In previous studies, other models approximating the PPT model have been used for the photoionization model of molecular gases due to their ease of incorporation into the nonlinear propagation equations. To see how sensitive prediction are to such approximations, we compared the predictive performance of the approximate models with that of the PPT model. We performed nonlinear propagation simulations employing two different approximate ionization models: the MPI and TI models. In addition, previous studies have also incorporated the intensity dependence of ionization rates into the NLS equations by a power-law approximation of the PPT model over an appropriate intensity range (PL approx.) [18,21,40]. To facilitate comparison, we conducted simulations employing the PL approx. model as well. Details of these three models are given in Appendix F. Using the aforementioned values of threshold fluence and $z$-position obtained from the fitting results ($F_{\mathrm {th}}$ = 0.47 J/cm$^2$ and $z$ = 4 $\mathrm{\mu}$m), we predicted the ablation diameter from the simulation results with each approximate ionization model. The pulse energy dependence of the squared ablation diameter from the simulation results using MPI, PL approx., and TI models are shown as orange, green, and purple lines in Figs. 4(a),(b), respectively. The simulation results predicted the experimental results well for all pulse energies only when the PPT model was employed, and the simulation results for models other than the PPT model deviate significantly from experimental values in the high pulse energy region. Compared with the PPT model, the ablation diameter was estimated to be larger in the nonlinear propagation regime for the MPI and PL approx. models. This is believed to be due to the higher ionization rates calculated at high intensities, and correspondingly, an overestimation of plasma defocusing effects in the nonlinear propagation regime. Conversely, in the TI model, the ablation diameter was estimated to be smaller in the nonlinear propagation regime. This is believed to be due to an underestimation of the ionization rate at lower intensities. Overall, the ionization models other than PPT model could not predict the pulse energy dependence of the ablation diameter because the ionization rate around 10$^{14}$ W/cm$^2$ was either underestimated or overestimated. This intensity range corresponds to the peak intensity range of sub-100 fs pulses at pulse energies typically used for material ablation (10$^{13}$-10$^{15}$ W/cm$^2$). Hence, from the above results, we concluded that it is important to incorporate the full PPT model to accurately predict diameters over a wide range of pulse energies.

To confirm the effectiveness of this method for other metals, the same laser ablation experiments were performed on gold and aluminum. Here, the pulse energy was varied from 2.0 to 34.0 $\mathrm{\mu}$J, and each sample was irradiated with 35 laser pulses. The optical microscope images of the gold and aluminum craters are shown in Figs. 5(a)-(c) and (d)-(f), and the pulse energy dependence of the square of the measured ablation diameter for gold and aluminum are shown as red points in Figs. 5(g),(i) and (h),(j), respectively. As was the case of copper, the ablated craters of gold and aluminum also became increasingly larger at higher pulse energies. The results of a $D^2$ fitting for each of the experimental data in the low pulse energy region (<10 $\mathrm{\mu}$J) are shown as black-dashed lines in Figs. 5(g),(i) and (h),(j). The threshold fluence determined here was Au: $F_{\mathrm {th}}$=0.45 $\pm$ 0.02 J/cm$^2$ and Al: $F_{\mathrm {th}}$=0.40 $\pm$ 0.01 J/cm$^2$. As was the case of copper, when the pulse energy exceeds approximately 18.0 $\mathrm{\mu}$J, the experimental data deviated from the $D^2$ method prediction. The fact that the deviation did not depend on the material further supports that this effect is not a material-specific effect. The results of the numerical simulations fitted to the experimental results for gold and aluminum, using the least-squares method, are shown as blue lines in Figs. 5(g),(i) and (h),(j), respectively. The threshold fluences and $z$-positions for gold and aluminum were Au: $F_{\mathrm {th}}$=0.47 $\pm$ 0.01 J/cm$^2$, $z$=−25 $\pm$ 6 $\mathrm{\mu}$m and Al: $F_{\mathrm {th}}$=0.41 $\pm$ 0.01 J/cm$^2$, $z$=−25 $\pm$ 5 $\mathrm{\mu}$m. All values agreed reasonably with the expected errors, as in the case of the experiment on copper. As can be seen from the figure, the predicted and measured ablation diameters were in very good agreement here as well, indicating that our developed method is applicable to a variety of metals, not just copper. These results also show that the ablated area is determined by a local threshold fluence even in the nonlinear propagation regime for these metals; the locality of laser ablation [28–31], which has been shown in the linear propagation regime, appears to hold for this material and laser combination even in the nonlinear propagation regime.

 

Fig. 5. Ablation results for gold and aluminum irradiated with 800 nm, 43 fs pulses at a repetition rate of 100 Hz. Each were irradiated with 35 pulses. Optical microscope images of gold irradiated with (a) 7.9 $\mathrm{\mu}$J, (b) 17.9 $\mathrm{\mu}$J, and (c) 32.0 $\mathrm{\mu}$J pulses, and aluminum with (d) 7.6 $\mathrm{\mu}$J, (e) 17.7 $\mathrm{\mu}$J, and (f) 31.8 $\mathrm{\mu}$J pulses are shown. Comparison of simulation results (blue line) and experimental results (red points) of the pulse energy dependence of the ablation diameter squared on a linear scale for (g) gold and (h) aluminum, and logarithmic scale for (i) gold and (j) aluminum. Five experimental data points are plotted for each pulse energy. The black-dashed line represents the fitting results of the experimental data (<10 $\mathrm{\mu}$J) to the $D^2$ method.

Download Full Size | PDF

We finally explore the correspondence between the fluence at the center of the nonlinear propagated laser pulse and the ablation results. We focus on the pulse energy dependence of the simulated center fluence and the measured ablation depth of copper. The blue line in Fig. 6(a) shows the pulse energy dependence of the center fluence at $z$=4 $\mathrm{\mu}$m determined by numerical simulation using the PPT model, while the black-dashed line shows the pulse energy dependence of the center fluence for the linear propagation case (Gaussian beam propagation with a beam radius of 10.8 $\mathrm{\mu}$m at the focal position). In the case of linear propagation, the center fluence increases linearly with pulse energy. Conversely, in the nonlinear propagation simulation, it decreases after reaching a maximum at about 18.0 $\mathrm{\mu}$J. Figure 6(b) shows the pulse energy dependence of ablation depth determined from laser ablation experiments on copper. The ablation depth decreased when the pulse energy exceeded about 18.0 $\mathrm{\mu}$J, as was previously discussed in Figs. 2 and 3. The pulse energy at which the ablation depth and the center fluence were maximized agreed well. Thus, we believe that it is possible to use such calculations to predict the optimal pulse energy to maximize ablation depth, determined by the tradeoff between a higher overall pulse energy and a lowered central fluence due to plasma defocusing. While the shapes of the two curves are similar, the relationship between center fluence and ablation depth is not a simple proportional relationship, and clarifying their relationship is a future work. We believe that these results are important basic data for investigating this relationship.

 

Fig. 6. (a) The pulse energy dependence of the center fluence. The blue line represents the simulation results at $z$=4 $\mathrm{\mu}$m, and the black-dashed line represents the center fluence for the linear propagation case. (b) Experimental results of the pulse energy dependence of the ablation depth of copper. Five experimental data points were plotted for each pulse energy. The experimental conditions are the same as that in Fig. 3.

Download Full Size | PDF

4. Conclusion

In this study, we have developed a method to quantitatively predict the pulse energy dependence of the crater shape ablated by sub-100 fs pulsed lasers over a two-orders-of-magnitude range in the pulse energy, including conditions where nonlinear propagation effects could not be ignored. Investigations revealed that the pulse energy dependence of the predicted ablation diameter agreed excellently with the measured ablation diameter. Furthermore, we found that using the modified PPT model as an air ionization model was important for quantitative predictions over a wide range of pulse energies. Additionally, we confirmed that the method was valid for several types of metals (copper, gold, and aluminum). Lastly, the pulse energy dependence of the simulated center fluence and measured ablation depth also showed good quantitative agreement between their tendencies.

These results assure that combining our simulation method with simulations modeling the interaction between the intense electric fields and materials will lead to quantitative prediction of the ablation depth of metals and, perhaps, even the ablation results of non-metallic materials, such as dielectrics, expectedly more sensitive to temporal changes in the laser pulses due to their nonlinear absorption characteristics. This study also revealed that the ablated crater shape showed a good correspondence with the calculated fluence profile, even in the nonlinear propagation regime. Such correspondence may make possible the development of reliable methods for controlling the crater shape by taking advantage of the characteristic beam profiles produced by nonlinear propagation, as was also discussed in part in [34].

On the other hand, the method presented in this study has some limitations. For example, this study did not consider the effects of heat accumulation during multi-pulse laser irradiation or the effects of the interaction between the ablated material and the laser pulses. Although such effects do not significantly affect the ablation results within the experimental conditions of this study, they may become non-negligible under conditions with higher laser repetition rates. However, even when such complex effects need to be considered, we believe that the ability of this simulation to correctly determine the intensity distribution of the incident laser pulse is an important prerequisite to predict ablation results quantitatively. Furthermore, the current framework depends on various nonlinear parameters of air incorporated within the NLS equation. Due to a lack of experimental data on the nonlinear parameters of air, the effectiveness of our method is confined to a limited range of laser conditions. We hope that the current results provide motivation to systematically determine such parameters in the future. Despite the above, even in the current form, the presented method can already quantitatively predict ablation results over a wide range of pulse energies, including the nonlinear propagation regime, and is expected to significantly improve the controllability of laser processing using sub-100 fs ultrashort laser pulses.

Appendix A. Intensity dependence of the ionization rates calculated with the PPT model

In this section, we describe the specific equations used to calculate the photoionization rates of oxygen and nitrogen molecules in our work. The photoionization rate (s$^{-1}$) of gas molecules for a linearly polarized laser in an electric field $E$ (V/m) can be written as follows [26,39]:

(4)$$W_{\mathrm{PPT}}(E)=\omega_{\mathrm{a.u.}}\sqrt{\frac{6}{\pi}}|c_{n^*,l^*}|^2 f(l,m)\frac{U_{\mathrm{i}}}{2U_{\mathrm{H}}}A_m(\gamma) \left(\frac{2E_0}{E\sqrt{1+\gamma^2}}\right)^{2n^*-|m|-3/2} \exp{\left(-\frac{2E_0}{3E}g(\gamma)\right)},$$

where $\omega _{\mathrm {a.u.}}$=4.13$\times$10$^{16}$ s$^{-1}$/a.u. represents the unit conversion coefficient (from atomic unit system to SI unit system) and $\gamma =\omega \sqrt {2m_e U_{\mathrm {i}}}/(eE)$ represents the Keldysh parameter. $\omega$ represents the angular frequency of the laser pulse, $m_e$ represents the electron mass, $U_{\mathrm {i}}$ represents the ionization energy of the molecules (oxygen molecules: 12.55 eV, nitrogen molecules: 15.576 eV [26]), and $U_\mathrm {H}$=13.6 eV represents the ionization energy of hydrogen atoms. $E_0$ is calculated as follows:

(5)$$E_0=E_\mathrm{H}\left(\frac{U_{\mathrm{i}}}{U_{\mathrm{H}}}\right)^{3/2},$$

where $E_\mathrm {H}$=514 GV/m represents the electric field in the hydrogen atom. $n^*=Z_{\mathrm {eff}} (U_{\mathrm {H}}/U_{\mathrm {i}})^{1/2}$ represents the effective principal quantum number, and $Z_{\mathrm {eff}}$ represents the experimentally determined effective charge, where $Z_{\mathrm {eff}}$=0.53 for oxygen molecules and $Z_{\mathrm {eff}}$=0.90 for nitrogen molecules [26]. $l^*=n^*-1$ is the effective orbital angular momentum. $l$ and $m$ are the orbital angular momentum and magnetic quantum number, respectively. $|c_{n^*,l^*}|^2$ is a dimensionless constant written as

(6)$$|c_{n^*,l^*}|^2=\frac{2^{2n^*}}{n^*\Gamma (n^*+l^*+1)\Gamma(n^*-l^*)},$$

where $\Gamma$ is the gamma function. $f(l,m)$ is defined as follows:

(7)$$f(l,m)=\frac{(2l+1)(l+|m|)!}{2^{|m|}(|m|)!(l-|m|)}, \ \mathrm{with} \ f(0,0)=1.$$

The remaining function $A_m(\gamma )$ is written as

(8)$$A_m(\gamma)=\frac{4}{\sqrt{3\pi}}\frac{1}{|m|!}\frac{\gamma^2}{1+\gamma^2} \sum_{\kappa>\nu}^{+\infty} \exp{\{-(\kappa-\nu)\alpha(\gamma)\}} \Phi_m\left(\sqrt{(\kappa-\nu)\beta(\gamma)}\right),$$

where $\nu =[U_{\mathrm {i}}/(\hbar \omega )][1+1/(2\gamma ^2)]$ and $\kappa$ is the next highest integer after $\nu$. The function $\Phi _m$ is defined as

(9)$$\Phi_m(x)=\frac{x^{2|m|+1}}{2}\int_0^1 \frac{e^{{-}x^2 t}t^{|m|}}{\sqrt{1-t}}dt=e^{{-}x^2}\int_0^x(x^2-y^2)^{|m|}e^{y^2}dy.$$

When $m$=0, Eq. (9) can be written as

(10)$$\Phi_m(x)=e^{{-}x^2}\int_0e^{y^2}dy=e^{{-}x^2}\times\frac{\sqrt{\pi}}{2}\times \mathrm{erfi}(x)$$

where $\mathrm {erfi}(x)=\int _0^x e^{y^2}dy$ is the imaginary error function. The other functions of $\gamma$ are as follows

(11)$$\alpha(\gamma)=2\times \left\{\left(\sinh^{{-}1}\gamma\right)-\frac{\gamma}{\sqrt{1+\gamma^2}}\right\},$$

(12)$$\beta(\gamma)=\frac{2\gamma}{\sqrt{1+\gamma^2}},$$

(13)$$g(\gamma)=\frac{3}{2\gamma}\left[\left(1+\frac{1}{2\gamma^2}\right) \left(\sinh^{{-}1}\gamma\right)-\frac{\sqrt{1+\gamma^2}}{2\gamma}\right].$$

To express $W_{\mathrm {PPT}}(E)$ as a function of the laser intensity $I$, we can use the relation $E=[2I/(cn\epsilon _0)]^{1/2}$, where $c$ is the speed of light, $n$ is the refractive index of the medium, and $\epsilon _0$ is the dielectric constant of vacuum. When considering the photoionization of oxygen and nitrogen molecules, $l=m=0$ [26,39]. The calculated intensity dependence of the ionization rates of oxygen and nitrogen molecules for a laser pulse of 800 nm wavelength is shown in Fig. 7.

 

Fig. 7. Intensity dependence of ionization rates for oxygen and nitrogen molecules calculated using the PPT model for a laser pulse of 800 nm wavelength.

Download Full Size | PDF

Appendix B. Analysis methods for the ablation diameter and depth

We detail here the methods used for determining the ablation diameter from the height image of the ablated crater measured by a laser scanning microscope. Figure 8(a) shows an optical microscope image of a typical ablation crater of copper. The dark central region corresponds to the ablated crater. The height image of the same crater measured by laser scanning microscope is shown in Fig. 8(b). Here, the base plane of the height (the plane where the height is 0) was determined from the average value in the area where the laser pulse was not irradiated in Fig. 8(b). To determine the ablation diameter, we performed a contour detection of the ablated crater by binarizing the height image with a height threshold $H_{\mathrm {th}}$. An example of a binarized height image with $H_{\mathrm {th}}$=−1.4 $\mathrm{\mu}$m is shown in Fig. 8(c). The contour of the crater can be seen as an ellipse. Figure 8(d) shows the result of a fit of the contour with an ellipse. We determined the ablation diameter at height $H$=$H_{\mathrm {th}}$ as the geometric mean of the lengths of the major and minor axes determined by the elliptical fit. In our experiments, we would like to determine the ablation diameter at $H_{\mathrm {th}}$=0 $\mathrm{\mu}$m. However, conducting the aforementioned analysis with $H_{\mathrm {th}}$=0 $\mathrm{\mu}$m was difficult due to the random roughness that existed before the laser irradiation (such as scratches on the sample surface) obstructing the contour detection step. To avoid this problem, we determined opted to the following procedure. First ablation diameters at several threshold values of $H_{\mathrm {th}}$ (−1.2 to −1.6 $\mathrm{\mu}$m) such that the random roughness of the sample surface did not affect the contour detection were determined. Next, the ablation diameter at $H_{\mathrm {th}}$=0 $\mathrm{\mu}$m was calculated by linear extrapolation from those ablation diameters. An example of such fitting is shown in Fig. 9.

 

Fig. 8. Ablation results of copper irradiated with 75 pulses at 800 nm, 43 fs pulses at 100 Hz repetition rate and 23.9 $\mathrm{\mu}$J pulse energy. (a) The optical microscope image of the ablated crater. (b) The height profile of the ablated crater measured by laser scanning microscope. (c) A binarized height image when the height threshold is −1.4 $\mathrm{\mu}$m. (d) The image of a crater contour (black line) obtained by the elliptical fitting.

Download Full Size | PDF

 

Fig. 9. The result of the linear fitting (black-dashed line) to the ablation diameters (red point) at different height thresholds for the ablated crater with 23.9 $\mathrm{\mu}$J. The height threshold $H_{\mathrm {th}}$ is taken in absolute value for ease of reading the figure.

Download Full Size | PDF

In determining the ablation depth, first, the center of the ellipse obtained from elliptical fitting to the contour of the crater was determined as the center of the crater. Subsequently, the average value of the heights within a 1.0 $\mathrm{\mu}$m radius from the center was defined as the ablation depth.

Appendix C. Comparison of the contribution of self-focusing and plasma defocusing to beam propagation

In this section, we discuss the relative effects of certain nonlinear effects to the overall propagation by comparing the results of several numerical simulations where we systematically eliminated the respective contributions of the terms related to self-focusing and the interaction between the photoelectric field and the plasma in the main NLS equation (Eq. (1)). Figure 10 shows the $z$ axis dependence of the beam radius for pulse energies of 2.0, 18.0, 26.0, and 32.0 $\mathrm{\mu}$J. The black-dashed line in Figs. 10(a)-(d) represents the results for beam propagation in vacuum (i.e., in the linear propagation regime), the red line represents results for beam propagation in air without the interaction between the electric field and the plasma by setting $\sigma _\omega$=0 in Eqs. (1) and (2), and the blue line represents results for beam propagation in air with interaction between the electric field and the plasma intact (i.e., simulated without eliminating any terms in Eqs. (1) and (2)). It is clear that when there is no interaction between the electric field and the plasma (red line), the beam propagation is similar to that in the linear propagation (black-dashed line), even when the pulse energy is increased. However, when the interaction between the electric field and the plasma is included in the simulation, the beam propagation changes significantly, especially at high pulse energies. From this result, we believed that plasma defocusing was the main factor affecting beam propagation for the considered experimental conditions. The effect of plasma defocusing on beam propagation depends on the free electron density. In our model, the generation of free electrons is described by photoionization and avalanche ionization. Under the present conditions, we have confirmed that the effect of avalanche ionization is small. Therefore, the free electron density obtained by the photoionization of gas molecules determines the threshold for the nonlinearity onset in our model.

 

Fig. 10. Simulation results of the $z$ axis dependence of the beam radius with (a) 2.0 $\mathrm{\mu}$J, (b) 18.0 $\mathrm{\mu}$J, (c) 26.0 $\mathrm{\mu}$J, and (d) 32.0 $\mathrm{\mu}$J pulses. The black-dashed line represents the beam propagation in vacuum, the red line represents the beam propagation in air without interaction between the electric field and the plasma, and the blue line represents the beam propagation in air with interaction between the electric field and the plasma.

Download Full Size | PDF

Appendix D. Method for determining the simulated ablation diameter from a calculated fluence distribution

Here, we graphically supplement our process for the extraction of the simulated ablation diameter from the simulated fluence distribution. An example of a simulated fluence distribution is shown in Fig. 11. In this study, the distance between two points where the simulated fluence and a material specific threshold fluence $F_{\mathrm {th}}$ coincide is defined as the simulated ablation diameter (the length of the blue double arrow in Fig. 11).

 

Fig. 11. Representation of the method for determining the ablation diameter from the simulated fluence distribution. The black-dashed line represents the ablation threshold fluence, and the blue double arrow represents the corresponding simulated ablation diameter.

Download Full Size | PDF

Appendix E. Statistical analysis of the goodness of fit

In this section, we discuss the statistical soundness regarding the fitting procedure of the ablation crater diameter to the simulation results. In our results, we fitted the simulated ablation diameter to the average value of the measured ablation diameter using the least-squares method with the threshold fluence and $z$-position of the sample as fitting parameters. A maximum likelihood estimate for each fitting parameter can be obtained by finding parameter values which minimize the value of $\chi ^2$ expressed by the following equation [42]:

(14)$$\chi^2=\sum_i\left\{\frac{1}{\sigma_i^2}[\bar{y}_i-y(x_i)]^2\right\},$$

where $\bar {y}_i$ is the measured value, $\sigma _i^2$ is the variance of the measured value, and $y(x_i)$ is the fitting function. In this study, $\bar {y}_i$ is the average value of the ablation diameter obtained in five experiments, $\sigma _i^2$ is the variance of the mean, and $y(x_i)$ is a function of the pulse energy containing the values of $F_{\mathrm {th}}$ and $z$ as constants. We note that $\sigma _i^2$ is different from the variance of the data points $s^2$. The variance $s^2$ can be calculated from the well-known formula for unbiased variance and written as

(15)$$s^2=\frac{\sum_j(y_j-\bar{y})^2}{n-1},$$

where $n$ is the number of experimental data points under the same condition. The relationship between the $\sigma _i^2$ and $s^2$ is defined by the following equation from the central limit theorem,

(16)$$\sigma_i^2=\frac{s^2}{n}.$$

In discussing the statistical significance of the simulation model, it is necessary to use $\sigma _i^2$ rather than $s^2$.

The average value of the squared ablation diameter measured for copper as well as its $\sigma _i$ are shown in Figs. 12(a),(b) as red points and error bars, respectively. The blue line in Figs. 12(a),(b) is the predicted ablation diameter obtained after a maximum likelihood estimate of the fitting parameters (same as the blue line in Figs. 4(a),(b)). In all pulse energy regions, the predicted ablation diameter and measured ablation diameter are in good agreement.

 

Fig. 12. Comparison of simulation results (blue line) with experimental results of copper (red points) of the pulse energy dependence of the ablation diameter squared on a (a) linear scale and (b) logarithmic scale. The experimental data of diameter squared are averaged over five data and error bar represents the standard deviation of this mean value.

Download Full Size | PDF

The goodness of fit between the predicted ablation diameter and measured ablation diameter was evaluated by a $\chi ^2$ test using the above $\chi ^2$ values. Using $\chi ^2$ and the number of degrees of freedom $\nu =N-m$ in the data ($N$: the number of data points, $m$: the number of fitting parameters), the reduced chi-square $\chi _\nu ^2$ can be calculated using the following formula [42]:

(17)$$\chi_\nu^2=\frac{\chi^2}{\nu}.$$

In this study, the number of degrees of freedom was $\nu$=32-2=30, and thus the reduced chi-square was $\chi _\nu ^2$=1.402. Generally, a value of $\chi _\nu ^2$ close to 1 implies a statistically likely fit; a value significantly smaller than 1 implies overfitting, while a value significantly larger than 1 implies an inappropriate model. Here, the present result is not rejected at a 5% risk rate [42], and we concluded that this simulation model is reasonable.

Appendix F. Comparison of the intensity dependence of ionization rates calculated by the PPT, MPI, PL approx., and TI models

In this section, we detail the intensity dependence calculations for each of the photoionization models employed in the nonlinear propagation simulations.

In previous studies, the intensity dependence of the ionization rate was often approximated by a power-law approximation in the functional form $W_{\mathrm {PL}}$=$\sigma ^{\mathrm {eff}}I^{K^{\mathrm {eff}}}$, interpreted as an effective multiphoton ionization rate at intermediate intensities between the MPI and TI regimes [18,21,40]. Here, $K^{\mathrm {eff}}$ represents the effective number of photons required for ionization and $\sigma ^{\mathrm {eff}}$ represents the effective cross-section of $K^{\mathrm {eff}}$-photon ionization. In the power-law approximation, the values of the parameters $\sigma ^{\mathrm {eff}}$ and $K^{\mathrm {eff}}$ are determined by fitting the $W_{\mathrm {PPT}}(I)$ calculated from the PPT model (blue line in Figs. 13(a),(b)) over a limited intensity range. The green line in Figs. 13(a),(b) shows the result of fitting the PPT model in the intensity range of 1$\times$10$^{13}$ to 1$\times$10$^{14}$ W/cm$^2$, which is the approximate range investigated in this experiment. Here, the values of $\sigma ^{\mathrm {eff}}$ and $K^{\mathrm {eff}}$ for nitrogen and oxygen molecules were $\sigma ^{\mathrm {eff}}_{\mathrm {N}_2}$=9.2$\times$10$^{-90}$ s$^{-1}$ (cm$^2$/W)$^{K^{\mathrm {eff}}_{\mathrm {N}_2}}$ with $K^{\mathrm {eff}}_{\mathrm {N}_2}$=7.16 and $\sigma ^{\mathrm {eff}}_{\mathrm {O}_2}$=5.5$\times$10$^{-65}$ s$^{-1}$(cm$^2$/W)$^{K^{\mathrm {eff}}_{\mathrm {O}_2}}$ with $K^{\mathrm {eff}}_{\mathrm {O}_2}$=5.44, respectively. From the figure, it can be seen that the ionization rate of the PL approx. model is larger than that of the PPT model outside of the intensity range used for the approximation.

 

Fig. 13. The intensity dependence of the ionization rate of (a) N$_2$, (b) O$_2$. The blue line represents the ionization rate of the PPT model, the orange line represents the ionization rate of the MPI model, the green line represents the ionization rate of the PL approx. model, and the purple line represents the ionization rate of the TI model.

Download Full Size | PDF

Additionally, the intensity dependence of the ionization rates of MPI and TI models can be derived by calculating the limit of the Keldysh parameter within $W_{\mathrm {PPT}}$ (in Eq. (4)) at $\gamma \gg 1$ and $\gamma \ll 1$ respectively [24]. The functional form of the ionization rate in the MPI model is $W_{\mathrm {MPI}}=\sigma I^K$. $K$ is the number of photons required for ionization and $\sigma$ is the cross-section of $K$-photon ionization. The values of $\sigma$ and $K$ for MPI of nitrogen and oxygen molecules were $\sigma _{\mathrm {N}_2}$=3.2$\times$10$^{-140}$ s$^{-1}$(cm$^2$/W)$^{K_{\mathrm {N}_2}}$ with $K_{\mathrm {N}_2}$=11 and $\sigma _{\mathrm {O}_2}$=8.2$\times$10$^{-112}$ s$^{-1}$(cm$^{2}$/W)$^{K_{\mathrm {O}_2}}$ with $K_{\mathrm {O}_2}$=9, respectively (N$_2$: $U_{\mathrm {N}_2}$ =15.576 eV, $Z_{\mathrm {eff}}$=0.9, O$_2$: $U_{\mathrm {O}_2}$ =12.55 eV, $Z_{\mathrm {eff}}$=0.53). The intensity dependence of the ionization rate of the MPI model is shown by the orange line in Figs. 13(a),(b). The ionization rate of the MPI model is consistent with that of the PPT model in the low-intensity region, whereas it overestimates the ionization rate in the high-intensity region. Finally, the functional form of the ionization rate in the TI model is as follows [24],

(18)$$W_{\mathrm{TI}}(E)=\omega_{a.u.}\sqrt{\frac{6}{\pi}}|c_{n^*,l^*}|^2 f(l,m)\frac{U_{\mathrm{i}}}{2U_{\mathrm{H}}} \left(\frac{2E_0}{E}\right)^{2n^*-|m|-3/2} \exp{\left(-\frac{2E_0}{3E}\right)}.$$

See Appendix A for the meaning and value of each parameter. The intensity dependent ionization rate calculated with the TI model is shown by the purple line in Figs. 13(a),(b). The ionization rate is consistent with that of the PPT model in the high-intensity region, whereas it underestimates the ionization rate in the low-intensity region.

Funding

Cabinet Office, Government of Japan (SIP); Ministry of Education, Culture, Sports, Science and Technology (Q-LEAP, JPMXS0118067246); New Energy and Industrial Technology Development Organization (JPNP16011); Center of Innovation Program (JPMJCE1313); Japan Science and Technology Agency (SPRING, JPMJSP2108).

Acknowledgments

This work was supported in parts by JST SPRING Grant Number JPMJSP2108, JST COI Grant Number JPMJCE1313, MEXT Quantum Leap Flagship Program (MEXT Q-LEAP) Grant Number JPMXS0118067246, Cross-ministerial Strategic Innovation Promotion Program (SIP), "Photonics and Quantum Technology for Society 5.0", and the New Energy and Industrial Technology Development Organization (NEDO) project "Development of advanced laser processing with intelligence based on high brightness and high efficiency laser technologies (TACMI project)".

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. F. Sima, K. Sugioka, R. M. Vázquez, R. Osellame, L. Kelemen, and P. Ormos, “Three-dimensional femtosecond laser processing for lab-on-a-chip applications,” Nanophotonics 7(3), 613–634 (2018). [CrossRef]  

2. X. Wang, H. Yu, P. Li, Y. Zhang, Y. Wen, Y. Qiu, Z. Liu, Y. P. Li, and L. Liu, “Femtosecond laser-based processing methods and their applications in optical device manufacturing: A review,” Opt. Laser Technol. 135, 106687 (2021). [CrossRef]  

3. B. N. Chichkov, C. Momma, S. Nolte, F. von Alvensleben, and A. Tünnermann, “Femtosecond, picosecond and nanosecond laser ablation of solids,” Appl. Phys. A 63(2), 109–115 (1996). [CrossRef]  

4. D. Nieto, J. Arines, G. M. O’Connor, and M. T. Flores-Arias, “Single-pulse laser ablation threshold of borosilicate, fused silica, sapphire, and soda-lime glass for pulse widths of 500 fs, 10 ps, 20 ns,” Appl. Opt. 54(29), 8596 (2015). [CrossRef]  

5. M. Lenzner, J. Krüger, S. Sartania, Z. Cheng, C. Spielmann, G. Mourou, W. Kautek, and F. Krausz, “Femtosecond optical breakdown in dielectrics,” Phys. Rev. Lett. 80(18), 4076–4079 (1998). [CrossRef]  

6. M. Lenzner, J. Krüger, W. Kautek, and F. Krausz, “Precision laser ablation of dielectrics in the 10-fs regime,” Appl. Phys. A 68(3), 369–371 (1999). [CrossRef]  

7. N. Sanner, O. Utéza, B. Chimier, M. Sentis, P. Lassonde, F. Légaré, and J. C. Kieffer, “Toward determinism in surface damaging of dielectrics using few-cycle laser pulses,” Appl. Phys. Lett. 96(7), 071111 (2010). [CrossRef]  

8. B. Chimier, O. Utéza, N. Sanner, M. Sentis, T. Itina, P. Lassonde, F. Légaré, F. Vidal, and J. C. Kieffer, “Damage and ablation thresholds of fused-silica in femtosecond regime,” Phys. Rev. B 84(9), 094104 (2011). [CrossRef]  

9. O. Utéza, N. Sanner, B. Chimier, A. Brocas, N. Varkentina, M. Sentis, P. Lassonde, F. Légaré, and J. C. Kieffer, “Control of material removal of fused silica with single pulses of few optical cycles to sub-picosecond duration,” Appl. Phys. A 105(1), 131–141 (2011). [CrossRef]  

10. M. Lebugle, N. Sanner, O. Utéza, and M. Sentis, “Guidelines for efficient direct ablation of dielectrics with single femtosecond pulses,” Appl. Phys. A 114(1), 129–142 (2014). [CrossRef]  

11. C. Pasquier, P. Blandin, R. Clady, N. Sanner, M. Sentis, O. Utéza, Y. Li, and S. Y. Long, “Handling beam propagation in air for nearly 10-fs laser damage experiments,” Opt. Commun. 355, 230–238 (2015). [CrossRef]  

12. C. Pasquier, M. Sentis, O. Utéza, and N. Sanner, “Predictable surface ablation of dielectrics with few-cycle laser pulse even beyond air ionization,” Appl. Phys. Lett. 109(5), 051102 (2016). [CrossRef]  

13. T. Genieys, M. Sentis, and O. Utéza, “Investigation of ultrashort laser excitation of aluminum and tungsten by reflectivity measurements,” Appl. Phys. A 126(4), 263 (2020). [CrossRef]  

14. T. Genieys, M. Sentis, and O. Utéza, “Measurement of ultrashort laser ablation of four metals (Al, Cu, Ni, W) in the single-pulse regime,” Adv. Opt. Technol. 9(3), 131–143 (2020). [CrossRef]  

15. T. Genieys, M. N. Petrakakis, G. D. Tsibidis, M. Sentis, and O. Utéza, “Unravelling ultrashort laser excitation of nickel at 800 nm wavelength,” J. Phys. D: Appl. Phys. 54(49), 495302 (2021). [CrossRef]  

16. D. Brinkmeier, D. Holder, A. Loescher, C. Röcker, D. J. Förster, V. Onuseit, R. Weber, M. Abdou Ahmed, and T. Graf, “Process limits for percussion drilling of stainless steel with ultrashort laser pulses at high average powers,” Appl. Phys. A 128(1), 35 (2022). [CrossRef]  

17. J. M. Liu, “Simple technique for measurements of pulsed Gaussian-beam spot sizes,” Opt. Lett. 7(5), 196 (1982). [CrossRef]  

18. N. Aközbek, C. M. Bowden, A. Talebpour, and S. L. Chin, “Femtosecond pulse propagation in air: Variational analysis,” Phys. Rev. E 61(4), 4540–4549 (2000). [CrossRef]  

19. R. Nuter and L. Bergé, “Pulse chirping and ionization of O₂ molecules for the filamentation of femtosecond laser pulses in air,” J. Opt. Soc. Am. B 23(5), 874 (2006). [CrossRef]  

20. Y. T. Li, T. T. Xi, Z. Q. Hao, Z. Zhang, X. Y. Peng, K. Li, Z. Jin, Z. Y. Zheng, Q. Z. Yu, X. Lu, and J. Zhang, “Oval-like hollow intensity distribution of tightly focused femtosecond laser pulses in air,” Opt. Express 15(26), 17973–17979 (2007). [CrossRef]  

21. F. Théberge, W. Liu, P. T. Simard, A. Becker, and S. L. Chin, “Plasma density inside a femtosecond laser filament in air: Strong dependence on external focusing,” Phys. Rev. E 74(3), 036406 (2006). [CrossRef]  

22. L. Bergé, S. Skupin, R. Nuter, J. Kasparian, and J. P. Wolf, “Ultrashort filaments of light in weakly ionized, optically transparent media,” Rep. Prog. Phys. 70(10), 1633–1713 (2007). [CrossRef]  

23. A. Couairon, E. Brambilla, T. Corti, D. Majus, O. De, J. Ramírez-Góngora, and M. Kolesik, “Practitioner’s guide to laser pulse propagation models and simulation,” Eur. Phys. J. Spec. Top. 199(1), 5–76 (2011). [CrossRef]  

24. A. Couairon, S. Tzortzakis, L. Bergé, M. Franco, B. Prade, and A. Mysyrowicz, “Infrared femtosecond light filaments in air: simulations and experiments,” J. Opt. Soc. Am. B 19(5), 1117–1131 (2002). [CrossRef]  

25. A. M. Perelomov, V. S. Popov, and M. V. Terent’ev, “Ionization of atoms in an alternating electrical field,” Sov. Phys. JETP 23, 924 (1966).

26. A. Talebpour, J. Yang, and S. L. Chin, “Semi-empirical model for the rate of tunnel ionization of N₂ and O₂ molecule in an intense Ti:sapphire laser pulse,” Opt. Commun. 163(1-3), 29–32 (1999). [CrossRef]  

27. M. V. Ammosov, N. B. Delone, and V. P. Krainov, “Tunnel ionization of complex atoms and atomic ions in electromagnetic field,” Sov. Phys. JETP 64, 1191–1194 (1986). [CrossRef]  

28. M. Garcia-Lechuga, O. Utéza, N. Sanner, and D. Grojo, “Evidencing the nonlinearity independence of resolution in femtosecond laser ablation,” Opt. Lett. 45(4), 952 (2020). [CrossRef]  

29. M. Garcia-Lechuga and D. Grojo, “Simple and robust method for determination of laser fluence thresholds for material modifications: An extension of Liu’s approach to imperfect beams,” Open Res Eur. 1, 7 (2021). [CrossRef]  

30. B. Zhou, A. Kar, M. J. Soileau, and X. Yu, “Invariance of the r²-ln(F) relationship and attainable precision in ultrafast laser ablation experiments,” Opt. Express 29(4), 5635–5643 (2021). [CrossRef]  

31. H. Sakurai, K. Konishi, H. Tamaru, J. Yumoto, and M. Kuwata-Gonokami, “Direct correlation of local fluence to single-pulse ultrashort laser ablated morphology,” Commun. Mater. 2(1), 38 (2021). [CrossRef]  

32. S. R. Vatsya, C. Li, and S. K. Nikumb, “Surface profile of material ablated with high-power lasers in ambient air medium,” J. Appl. Phys. 97(3), 034912 (2005). [CrossRef]  

33. J. Sun and J. P. Longtin, “Effects of a gas medium on ultrafast laser beam delivery and materials processing,” J. Opt. Soc. Am. B 21(5), 1081–1088 (2004). [CrossRef]  

34. C. Li, S. R. Vatsya, and S. K. Nikumb, “Effect of plasma on ultrashort pulse laser material processing,” J. Laser Appl. 19(1), 26–31 (2007). [CrossRef]  

35. W. Komatsubara, K. Konishi, J. Yumoto, and M. Kuwata-Gonokami, “Critical Power for Temporal-Pulse Collapse in Third-Harmonic Generation,” arXiv, arXiv:2208.08123v1 [physics.optics] (2022). [CrossRef]  

36. J. K. Wahlstrand, Y.-H. Cheng, and H. M. Milchberg, “Absolute measurement of the transient optical nonlinearity in N₂, O₂, N₂O, and Ar,” Phys. Rev. A 85(4), 043820 (2012). [CrossRef]  

37. P. E. Ciddor, “Refractive index of air: new equations for the visible and near infrared,” Appl. Opt. 35(9), 1566–1573 (1996). [CrossRef]  

38. Y. Liu, M. Durand, A. Houard, B. Forestier, A. Couairon, and A. Mysyrowicz, “Efficient generation of third harmonic radiation in air filaments: A revisit,” Opt. Commun. 284(19), 4706–4713 (2011). [CrossRef]  

39. J. Schwarz, P. Rambo, M. Kimmel, and B. Atherton, “Measurement of nonlinear refractive index and ionization rates in air using a wavefront sensor,” Opt. Express 20(8), 8791–8803 (2012). [CrossRef]  

40. S. L. Chin, N. Aközbek, A. Proulx, S. Petit, and C. M. Bowden, “Transverse ring formation of a focused femtosecond laser pulse propagating in air,” Opt. Commun. 188(1-4), 181–186 (2001). [CrossRef]  

41. G. P. Agrawal, Nonlinear Fiber Optics (Academic Press, 2019), 6th ed.

42. R. R. Bevington and D. K. Robinson, Data Reduction and Error Analysis for the Physical Sciences (McGraw-Hill Science Engineering, 2002), 3rd ed.