where G is a constant that depends on the transition probability. The parameter R is determined by the optical transmittance type and is R = 0.5 for directly allowed optical transmissions [29,30]. Plotting αEP over EP and using a linear fit of the straight portion of the absorption edge, the intersection at αEP = 0 determines EBG (Figure 2b). The results demonstrate optimal conditions for a strong fundamental single photon absorption of the pump beam (λpump = 355 nm). The optical penetration depth is determined by dp = 1/α with the absorption coefficient α = ln(T−1)/dL. dp = 62 nm and corresponds by definition to the depth at which the beam intensity has dropped to 1/e. This value is smaller than the layer thickness dL and thus almost all of the pulse energy is used to promote electrons to the CB, which are available for the subsequent ESA-based ablation process.

By analyzing the sinusoidal interference fringes in the weak and medium absorbing region of the TZnO-spectrum by the envelope method of Swanepoel [26], the interference free transmittance TE was derived. Here, the TZnO-maxima and minima were used to fit an upper TM and lower Tm envelope function, respectively. Both envelopes were computer-generated using the program OriginPro 2020b (Origin-Lab Corp.) [30]. At shorter wavelengths in the strong absorbing region above EBG, the envelope functions converge to a single curve, merging with the measured values of TZnO. The interference free absorption coefficient was calculated from the determined TE values (Table 1) [31].

**Table 1.** Optical parameters (TE, n, α, dp) at the utilized λ calculated from the transmittance data by the method of Swanepoel (Figure 2a) [26], as well as ablation parameters for single-beam ablation (Eth, 2ωf, Fth) determined by the method of Liu (Figure 3a) [25].

**Figure 2.** (**a**) UV/VIS transmittance spectra of 1 mm thin uncoated fused silica substrate (**---**) and fused silica coated with a 200 nm thin ZnO-layer ((ŷ)). The dotted lines labeled TM, Tm, and TE are the fits of the upper and lower envelope function as well as their geometric mean, respectively [26]. In the emission spectrum of the ZnO film ((ŷ)), λ = 325 nm was used for excitation. Selected emission peaks assigned to specific transitions are numbered. (**b**) Tauc plot [28] of TZnO ((ŷ)). The linear fit ((ŷ)) at the straight portion of the band edge allows to determine EBG by the intersection at αEP = 0. (**c**) Schematic of the location of intrinsic point defect levels within the band gap by full-potential linear Muffin-tin orbital method calculations according to Xu et al. [32]. The numbers on the lines indicate ΔEP of the defect states compared to the energy of the valence band (VB)-maximum with filled ((ŷ)), partially filled (---), and empty (**---**) states. The positions of the energy levels are true to scale with respect to EBG = 3.26 eV. Selected transitions are labelled with dashed arrows and a number corresponding to the emission wavelengths in (**a**).

Furthermore, the amplitude and oscillation of the interference fringes contain information about the ZnO-layer thickness dL and the refraction index n, respectively. The refraction index of the substrate was derived from the transmittance of the uncoated substrate. The interference free refraction index in the weak absorbing region was calculated using the envelope functions as demonstrated by Swanepoel [26]. By using λ and n of adjacent fringe maxima or minima and the order number, dL was calculated to be 203 nm. The measured and calculated parameters for the utilized λ are listed in Table 1 [26,31].

Below the band gap energy EBG, n and α steadily increase and the small variation of λ below EBG indicates the comparably weak absorption of the ESA wavelengths without additional pump excitation. The difference between TE and TS is caused primarily by the deviating n, as the upper envelope function that intersects the fringe maxima almost coincides with TS [26]. With increasing absorption at EBG, α increases rapidly.

The sinusoidal interference fringes in the spectrum are caused by reflection losses due to interference effects of the radiation between the air–film and film–substrate interfaces [33]. This indicates homogeneous, smooth, and low scattering surfaces [34]. Thus, only at sufficiently high beam fluences, an independent ESA beam ablation is expected by either a multi-photon absorption or processes including various intrinsic point defects. The latter include deep level defects inside the ZnO layer or shallow donors that are located typically at the surface [20,22,35,36]. As illustrated in Figure 2c, the location of the intrinsic defect levels were calculated by Xu et al. with a full-potential linear Muffin-tin orbital method [32]. These include vacant zinc (VZn), vacant oxygen (VO), interstitial zinc (Zni), interstitial oxygen (Oi), substituted oxygen at zinc positions (OZn), and complexes of V0 and Zni (V0Zni).

The emission spectrum in Figure 2a shows a weaker emission in the UV-range from 360–500 nm (close to EBG), compared to λ = 550–800 nm. The sharp band emission at λ1 = 380–390 nm overlaps with the band edge and is caused by recombination of photoinduced charge carriers through exciton–exciton collision processes (Figure 2c) [21,22,37]. The specific emissions at λ2 = 430 nm and λ3 = 480 nm are attributed to Zni → VB and

Zni → VZn transitions, respectively [22]. Their weak emission indicates low roughness and homogeneous surface quality, in line with the findings in the TZnO spectrum as these transitions occur primarily at the surface [22]. The strongest emission in the visible region corresponds to deep level defects in the bulk of the layer [22,32,36]. Here, the broad emission is comprised of several transitions at λ4 = 625 nm, λ5 = 650 nm, λ6 = 700 nm, and λ7 = 750 nm, induced by Zni → Oi, Zni → OZn, CB → OZn, and Zni → V0Zni [22].

#### *3.2. Single-Beam Ablation*

It becomes evident from Table 1 that Fth of the ESA beams is about one order of magnitude larger than the pump beam. The pump beam with λpump = 355 nm is absorbed fundamentally by single photon absorption, indicated by the linear increase in D2 (at the semilog scale in Figure 3a) and the transmittance in Figure 2a (Epump > EBG) [38,39]. In contrast, the ESA beams require either populated defect levels that absorb the laser wavelengths and promote electrons to CB [35], or a multi photon process (EESA < EBG) [7]. However, such nonlinear processes would cause an exponential deviation from the linear trend [38] and is therefore unlikely. Thus, the ESA at defect levels is expected as the primary absorption mechanism. The steeper slope angle of the pump beam indicates a larger beam diameter 2ωf = 28 μm [40]. It is comparable or larger than the ESA beams (Table 1) and therefore it provides ideal conditions for optimal utilization of the pulse energy for the ESA process.

**Figure 3.** (**a**) Semilog plots of D2 over E with the corresponding linear fits, to determine the ablation parameters Eth, 2ωf, and Fth of single-beam absorption by the method of Liu [25]. The exact values at the utilized λ are listed in Table 1. (**b**–**e**) Semilog plots of D2 over FESA for ESA based ablation at different λESA with (**b**) λ = 450 nm, (**c**) λ = 500 nm, (**d**) λ = 550 nm, and (**e**) λ = 600 nm. For excitation, different Fpump values of 0.2 J/cm2 ((်ȣ်)), 0.1 J/cm2 (Ÿ), and 0.02 J/cm2 (ź) were used. FESA was varied below their respective Fth (Table 1), indicated by the vertical dashed line. The values for Fpump = 0.2 J/cm2 were linearly fit to D2 = 0. The diagonal dotted lines highlight the Fth fits of λ = 355 nm (---) and the respective λESA (---) plotted in (**a**).

The behavior at λESA = 450 nm deviates slightly from the other wavelengths, where Fth is smaller and 2ωf is larger. This suggests a partially different excitation process due to the high EP closest to EBG. This might be explained by the excitation of VB-electrons directly to defect levels close to CB. Moreover, defect levels could be able to absorb the beam and excite electrons directly to CB (Figure 2c), which would lower Fth [32,35]. As the defect concentration is expected to be highest at the ZnO–fused silica interface, most of the energy would be absorbed here. At sufficient FESA, the abrupt absorption at this interface can cause an explosion-like ablation of the whole layer without the necessity to heat the entire ZnO inside the laser spot. For single-beam ablation of the ESA-beams, the resulting spots are mostly inhomogeneous at low FESA (Figure 4b,g) and the further increase of FESA leads to substrate damage (Figure 4a,f). Here, an inhomogeneous ablation with ablation depths da of several microns were observed in the dark areas of the ablation spot. The intense absorption at the fused silica interface exceeds the substrate ablation threshold and limits the FESA-range significantly. Thus, for the Fth calculations, FESA was decreased and the spots at higher FESA exceeding the damage threshold of the fused silica substrate were excluded from analysis, as the measured D2 are highly inconsistent and do not represent the ablation characteristic of the ZnO-layer.

**Figure 4.** Microscopic images of single-beam ablation spots with (**a**) F450 = 5.7 J/cm2, (**b**) F500 = 6.8 J/cm2, and (**c**) Fpump = 1.1 J/cm2, as well as ESA-based dual-beam ablation spots at (**d**) Δt = 5 ns with Fpump = 0.2 J/cm2 + F450 = 0.9 J/cm2 and (**e**) Δt = 10 ns with Fpump = 0.2 J/cm2 + F500 = 1.0 J/cm2. (**f**–**j**) White light interferometric depth profiles of spots (**a**–**e**).

#### *3.3. Dual-Beam Ablation*

3.3.1. Influence of the ESA Beam Fluence

At first, FESA was varied (Figure 3b–e) and the values for optimum ablation conditions were set to a small positive temporal delay of Δt = 5 ns and Fpump = 0.2 J/cm2. In the case of an ablation process solely based on ESA, the pump beam induces the transient absorption of the spatially and temporally superimposed ESA beam (Figure 1). This requires to prevent surface modification of the individual beams by ablation [8]. Therefore, all measurements were performed using laser peak fluences below their respective Fth (Table 1). This is indicated in Figure 3b–e where all values are below the vertical dashed lines, originating at Fth of the respective ESA-beams. Additionally, to compare the dual-beam process with single-beam ablation, the diagonal dashed lines show the linear fits of the Fth calculations, determined by the fits of Eth in Figure 3a.

Decreasing FESA results in a decreased D2 as less additional energy is available to further excite the already excited electrons. This limits the ablation to an area closer to the beam center, where the intensity of the Gaussian beam profile is still sufficient for ablation. Here, the utilizable fluence range for ESA-based ablation extends about one order of magnitude down to FESA ≈ 0.3 J/cm2 for λ = 450 nm. Slightly below these values, the ablation spots are characterized by their non-circular shape. By further decreasing FESA, inhomogeneously distributed substructures were observed within the irradiated area. These surface modifications at low laser fluences are caused by locally enhanced carrier densities that are generated by either an inhomogeneous laser beam profile or locally enhanced ablation at dust, scratches, or crystal defects [41]. It has to be noted that the analysis only includes ablated structures that resulted in distinct and well pronounced ablation spots, excluding the aforementioned value-ranges.

The comparison of the different ESA wavelengths reveals a similar behavior of the beams with the highest photon energy at λ = 450 and 500 nm showing highly increased D2 when compared to λ = 550 and 600 nm at similar fluences (Figure 3b–e). Higher EP allow an increased number of populated defect levels to absorb these ESA beams and to populate the CB by an ESA (Figure 2c). It has to be noted that a single photon absorption of VB-electrons directly to the CB is not possible at λESA with EP < EBG. In addition, a multiphoton absorption is negligible at the low applied fluences. Furthermore, a transition to another defect level is unlikely due to the specific required wavelength. Thus, besides the intra-band transition of CB-electrons, a transition from defect levels close to the VB into the CB is the most likely absorption mechanism. The existence of a considerable amount of these defect levels is confirmed by the characteristic fluorescence (Figure 2a), where the highest number of emitted photons is attributed to these defect level transitions [42–44].

Thus, with the highest EP, using λESA = 450 nm is optimal to generate the biggest ablation spots at the lowest FESA (Figure 3b). Here, at Fpump = 0.2 J/cm2 (red markers), even less additional FESA was necessary than the single pump beam requires for ablation. This is indicated by values above the purple line (linear Fth-fit of Figure 3a).

To analyze the effect of Fpump on the ESA-based ablation, Δt = 5 ns at the identical values of FESA for each wavelength were used. As seen in Figure 3d,e, below Fpump = 0.2 J/cm2 no distinct ablation spots can be observed for λ = 550 and 600 nm. The utilizable FESA range reduces drastically with decreasing Fpump. Without a sufficient excited electron concentration, the ESA-based ablation process stops abruptly. Only for λ = 450 nm at Fpump = 0.02 J/cm2 distinct ablation spots can be observed (Figure 3b,c).

#### 3.3.2. Ablation Efficiency

When compared to single-beam ablation, the ESA-based process allows to decrease the amount of laser energy required for ablation (Table 1). This improvement was evaluated based on the method of Liu [25]. For this purpose, D2 of the ESA-based ablation spots in Figure 3b–e (for Fpump = 0.2 J/cm2 at Δt = 5 ns) was plotted in a semilog plot as a function of FESA. The linear extrapolation allows to calculate the effective ESA-beam diameter 2ωf eff (slope of the fit) and the effective ESA-beam threshold peak fluence Fth eff (fit intersect at D2 = 0) (Table 2).

**Table 2.** Calculation of dual-beam ablation parameters (2ωf eff, Ftheff) based on the method of Liu [25] and dual-beam efficiency parameters (FE, FET) determined from Equations (3) and (4).

Here, the values of 2ωf eff for λ = 450 and 500 nm are close to their respective 2ωf of single-beam ablation. Consequently, almost the entire area of the fluence profile exceeds the ablation threshold and contributes to the ESA-based ablation. This suggests an effective utilization of the pulse energy. For λ = 550 and 600 nm, the values deviate significantly by a factor of two, which becomes evident by the decreased slope of the fits in Figure 3d,e.

Similarly to Fth that defines the minimum required fluence for single-beam ablation, Fth eff defines the minimum FESA at given Fpump that is required for ESA-based dual-beam ablation. These values decreased by a factor of 5–10, but are still proportional to Fth. With the highest EP, λESA = 450 nm can utilize the energy, induced by the pump, most optimally. Here, Fth eff = 0.26 J/cm is even below Fth = 0.49 J/cm2 of the single pump beam, indicated by values above the purple line, as described earlier (Figure 3b).

To evaluate the improvements of the ESA-based process, compared to single-beam ablation, fluence efficiency

$$\mathbf{F\_E = \frac{F\_{\rm th}^{\rm eff}}{F\_{\rm th}^{\rm pump}}} = \frac{\mathbf{F\_{pump}}}{\mathbf{F\_{th}^{\rm pump}}} \tag{3}$$

and total fluence efficiency as the sum of the normalized single-beam fluences

$$\mathbf{F\_{ET}} = \frac{\mathbf{F\_{th}^{eff}}}{\mathbf{F\_{th}^{ESA}}} + \frac{\mathbf{F\_{pump}}}{\mathbf{F\_{th}^{pump}}} \tag{4}$$

were defined (Table 2). For FE-values below one, the combined fluences of the ESA- and the pump-beam are still below Fthpump for single-beam ablation (FE = 1, if Fth eff + Fpump = Fthpump). As already mentioned, only λ = 450 nm at Δt = 5 ns and Fpump = 0.2 J/cm2 reached FE = 0.86. Similarly, replacing the single-beam threshold fluence in Equation (3) to Fth of the ESA-beams, allows the comparison of the efficiency improvement with the ESA-beams. In order to exclude an ordinary accumulation effect of the pump and ESA beam, FET was defined. For FET < 1, an additional induced transient absorption must be considered as part of the ablation process (Table 2).

#### 3.3.3. Ablation Quality

Figure 4a,b show selected single-beam ablation spots for the utilized ESA wavelengths. At high FESA, the absorption of defect levels at the fused silica–ZnO interface leads to substrate damage (Figure 4a,f). At FESA close to Fth, the ablation spots become increasingly inhomogeneous (Figure 4b). As described earlier, the further decrease of FESA results in only partially interconnected spots inside the beam. In comparison, the pump beam hits the ZnO–air interface first (Figure 1b) and the fundamental absorption leads to an almost circular and homogeneous single-beam ablation of the entire ZnO layer (Figure 4c). When reaching the ZnO–substrate interface, the pump beam is already highly attenuated preventing abrupt heating and damaging of the substrate (Figure 4h). The maximum da observed in the depth profiles (Figure 4g–j) are in good agreement with the calculated layer thickness dL = 203 nm.

Compared to single-beam ablation, the spot diameters are decreased by the ESA-based ablation process, indicated by decreased values of 2ωf eff when compared to 2ωf (Tables 1 and 2). The ablation area is restricted by the superposition of both pulses. This limits the ablation to the mostly uniform high-intensive center of the pulses and results in an enhanced ablation quality with circular ablation spots. In comparison to single pump beam ablation, the absorbed energy and induced heat at the ZnO–air interface is reduced and homogeneously distributed throughout the entire ZnO-layer. This results in a reduced melt formation at the edges protruding above the layer surface (Figure 4g–j), as well as decreasing the surface contamination by resolidified ZnO-droplets around the ablation spots (Figure 4d,e). By increasing Δt, these effects are reduced even further, as the heat dissipates and an ESA at defect levels dominates the ablation process (Figure 4i,j).

The optimization of the parameters allows the complete removal of the ZnO layer without damaging the fused silica substrate (Figure 4d,i) and a precise control of da (Figure 4e,j). The reduced FESA prevents the fused silica–ZnO interface from abruptly heating to the point of substrate damage (Figure 4a,f). As the fundamentally absorbed pump beam is focused on the ZnO–air surface, the initial abrupt heating together with a high concentration of excited carriers diffuse gradually and homogeneously through the ZnO-layer (dp = 62 nm). Optimal ablation conditions were found at Δt ≈ 5 ns. The ZnO layer is already excited by the pump beam, in particular at the ZnO–substrate interface, where the defects are located primarily. It thus requires only a low excess of additionally absorbed ESA photons and energy to ablate the ZnO layer from the fused silica substrate.

#### 3.3.4. Influence of Temporal Pulse Delay

To determine the importance of defect levels compared to absorption solely based on intra-band transition, the temporal delay between the pulses was varied for different FESA at the fixed optimum of Fpump = 0.2 J/cm (Figure 5). At the highest FESA, for λ = 450 and 500 nm (Figure 5a,b), delays above 300 ns still resulted in a measurable ablation. Compared to λ = 550 and 600 nm (Figure 5c,d), this delay is much larger with Δt mostly below 50 ns. As the fluorescent lifetime of ZnO band gap electrons is well below 10 ns [42,44], an ablation at these high Δt cannot be explained solely by intra-band transition-based absorption. Thus, the defect levels with their increased fluorescent lifetime serve as the main driving force [35]. First, the pump beam promotes electrons into CB (Figure 1a). From there, besides relaxation back to the VB, electrons transfer rapidly to deep defect levels, which are suitable as absorption centers for the subsequent ESA beam (Figure 2c). The decrease of D2 and the Δt-range is almost proportional to FESA. Only FESA-values close to their respective Fth deviate from this trend. Especially for λ = 550 and 600 nm, these values increase overproportional at the highest FESA, where the ESA-beams can almost ablate the ZnO solely by single-beam ablation.

**Figure 5.** D2 plotted versus Δt to determine the optimal delay between the pulses. For excitation, Fpump = 0.2 J/cm2 was used at (**a**) λESA = 450 nm, (**b**) λESA = 500 nm, (**c**) λESA = 550 nm, (**d**) λESA = 600 nm. FESA was varied in the identical increments as for Figure 3b (color coded with FESA and the FESA/Fth-ratio). In the top-right corners, the section of short Δt is extended.

Figure 5c,d reveals deviations from this trend, especially for λ = 550 and 600 nm, where the values increase overproportionately at the highest values of FESA. These fluence values are close the respective ablation thresholds, which almost allow ablation of the ZnO layer solely by the ESA beam.

At low temporal delays (Δt < 10 ns), the impact of λESA is less pronounced (Figure 5). Here, the effect of an intra-band transition and sample heating dominates as the beams still overlap. For λ = 550 and 600 nm, an increasing Δt leads to a rapidly decreasing D2 because the defect level absorption emerging as the dominant process is only utilized effectively at λ = 450 and 500 nm. As shown in Figure 1c, at Δt = +5 and −5 ns, the temporal overlap is similar, but the resulting D2 differ significantly, as prior excitation of electrons by the pump beam is required for ESA.

Exemplary, for λ = 500 nm at FESA = 4.51 J/cm2, Δt was varied in one ns-steps between Δt = −5 to +10 ns to analyze the area of temporal beam overlap in more detail (Figure 5b). At the optimum conditions around Δt ≈ 5 ns, the highest number of prior excited electrons is available for an ESA (Figure 1c). At Δt = 0 ns, where both beams overlap optimally, the excited electrons are only partially used for the ESA, as parts of the radiation excite carriers after the highest intensity of the ESA beam reaches the sample. At negative temporal delays and further decreasing beam overlap, only minor ablation is observed for FESA close to Fth. Here, the number of excited electrons is insufficient to reach the required ablation threshold. Increasing Δt above the optimum also leads to a decreasing D2-values. However, the rate is much lower than for negative delays.

At Δt < 50 ns, the heat generated by the pump excitation must be considered, as it effects the ablation process considerably, particularly close to the temporal beam overlap (Figure 1c). Here, increased sample temperatures act as an offset for the ESA, as less additionally absorbed energy is required for ablation. Thus, a wider area of the ESA beam surpasses Fth, resulting in increased D2. Fpump = 0.2 J/cm2, well below Fth = 0.49 J/cm2 was used to reduce sample heating.

For polycrystalline ZnO-layers, thermal conductivities kth < 10 W/m·K are found in literature, which is well below kth ~ 30–100 W/m·K for epitaxially grown or bulk ZnO [45–49]. The thermal diffusion length is determined by μ = 2(αth·Δt)0.5 with αth = kth/cV where cV is the volumetric heat capacity with about 2.4–2.8·106 J/m3·K [47,49]. With dp = 62 nm, most of the heat is created by the pump beam at the ZnO–air interface. It takes about 5 ns to reach the ZnO–fused silica interface, where the highest defect concentration and thus, the strongest ESA is expected. This is in good agreement with the ablation maximum measured at a delay of about 5 ns.

Furthermore, the rapid D2-decrease in the delay range 20–50 ns (Figure 5) can be attributed to heat dissipation and correlates with literature values [45,50]. As 2ωf is about 100 times larger than dL and μ for the time scale in consideration, one-dimensional heat conduction can be assumed [47,51,52]. Here, the heat should be mostly dissipated. Consequently, the defect level absorption is now the dominant factor for ESA, and it is no longer superimposed by pump beam heating. At higher Δt, the ablation slowly decreases due to the finite fluorescent lifetime of the defect level electrons that slowly lowers the available excited electron density for the ESA [43,53].

#### **4. Conclusions**

In the present study, we investigated a dual-beam laser process based on ESA by applying several wavelengths (λpump = 355 nm, λESA = 450, 500, 550, 600 nm) at various Fpump and FESA. The ESA-beams with the highest photon energies produced the largest ablation spots at the lowest beam fluences. A small temporal delay of +5 ns is optimal for the dual-beam ablation process, as it utilizes a combination of heat and excited carriers for increased absorption. Longer temporal delays result in heat dissipation and quickly lower the ablation diameters. Now, the ESA beam absorption is determined solely by the excited carrier concentration. The quality of the ablation spots was improved by reducing the debris and melt formation and the required laser energy was decreased. Furthermore, a model for the evaluation of the ESA-based efficiency improvement was developed. Based on these findings, future investigations on a combination of STED and ESA with regards to STED direct writing processes were advanced.

**Author Contributions:** Conceptualization, C.W. and S.E.; data curation, C.W.; formal analysis, S.G. and F.A.M.; funding acquisition, S.E., S.G., and F.A.M.; investigation, C.W. and S.E.; methodology, C.W. and S.E.; project administration, S.G. and F.A.M.; resources, F.A.M.; software, C.W.; supervision, S.G. and F.A.M.; validation, C.W., S.E., and F.A.M.; visualization, C.W.; writing—original draft, C.W.; writing—review and editing, S.E., S.G., and F.A.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the German Research Foundation (DFG), grant number GR4710/4-1.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Backward Flux Re-Deposition Patterns during Multi-Spot Laser Ablation of Stainless Steel with Picosecond and Femtosecond Pulses in Air**

**Tong Zhou 1, Sebastian Kraft 2, Walter Perrie 1,\*, Jörg Schille 2, Udo Löschner 2, Stuart Edwardson 1 and Geoff Dearden 1**

**Abstract:** We report on novel observations of directed re-deposition of ablation debris during the ultrafast laser micro-structuring of stainless steel in the air with multi-beams in close proximity on the surface. This interesting phenomenon is observed with both 10 ps and 600 fs NIR laser pulses at 5 kHz repetition rate. Ablation spot geometries could be altered with the use of beam splitting optics or a phase-only Spatial Light modulator. At low fluence (F ~ 1.0 J cm−2) and pulse exposure of a few hundred pulses, the debris appears as concentrated narrow "filaments" connecting the ablation spots, while at higher fluence, (F ~ 5.0 J cm−2) energetic jets of material emanated symmetrically along the axes of symmetry, depositing debris well beyond the typical re-deposition radius with a single spot. Patterns of backward re-deposition of debris to the surface are likely connected with the colliding shock waves and plasma plumes with the ambient air causing stagnation when the spots are in close proximity. The 2D surface debris patterns are indicative of the complex 3D interactions involved over wide timescales during ablation from picoseconds to microseconds.

**Keywords:** laser ablation; ultrafast; re-deposition

#### **1. Introduction**

The interaction of intense ultrafast laser radiation with atoms, gases and solid materials has led to remarkable scientific advances such as High Harmonic Generation for attosecond spectroscopy [1,2], the use of filamentation for atmospheric probing [3,4], two-photon microscopy [5,6] and intraocular fs-LASIK in eye surgery [7]. In particular, laser ablation (LA) of materials with ps and sub-ps temporal pulse length has advanced materials analysis in areas such as Laser Ablation Inductively Coupled Mass Spectrometry (LA-ICP-MS) [8–10] and Laser-Induced Breakdown Spectroscopy (LIBS) [11,12]. With temporal pulse lengths <10 ps in ultrafast laser ablation, the instantaneous energy deposition reduces plasma absorption while minimizing melt and thermal diffusion during the pulse and yield high precision material removal [13,14]. This allows for easy generation of plasmonic, interferometric and hierarchical microstructures on, for example, stainless steel [15–17]. In metals, which have absorption coefficients typically α ~ 106 cm−1, a thin layer with thickness d ~ 10–30 nm is converted to a plasma at a solid density which expands well after the pulse is absorbed [18,19].

In vacuo, laser ablation can be described as a free expansion [20–23] while the laser plume which develops in the presence of an environmental gas is physically much more complex involving timescales from picosecond to nanosecond and microsecond. After electron heating during the pulse, the lattice heats typically over a few picoseconds due to e-phonon coupling, raising the surface temperature well above the evaporation point.

**Citation:** Zhou, T.; Kraft, S.; Perrie, W.; Schille, J.; Löschner, U.; Edwardson, S.; Dearden, G. Backward Flux Re-Deposition Patterns during Multi-Spot Laser Ablation of Stainless Steel with Picosecond and Femtosecond Pulses in Air. *Materials* **2021**, *14*, 2243. https://doi.org/10.3390/ ma14092243

Academic Editors: Jörg Krüger and Jörn Bonse

Received: 20 March 2021 Accepted: 23 April 2021 Published: 27 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

This results in fast surface electron emission due to thermionic [24] or multi-photon absorption [25,26] and a remarkable early stage "air plasma" can develop during the first 150 ps [27,28], with electron density ne ~ 1020 cm−3 and a transient electric field developing above the surface [29].

On the nanosecond timescale in air, the plasma plume expands rapidly and initially longitudinally (at speed v ~ 106 cm s−1), colliding with the high density of surrounding gas molecules, (N1atm ~ 2.6 × 1019 cm−3), confining the plasma to the surface region while creating a blast (shock) wave which moves ahead of the plume [30–34]. As the longitudinal expansion decelerates rapidly, the plasma expands laterally after a few ×100 ns (vlat ~ 105 cm s−1) and can form a vortex structure with plasma plume concentrated at the ablation spot edge [35]. Transient electron temperatures in the plasma can reach Te >> 104 K [24], then plasma cooling (over a few hundreds of ns to μs) occurs due to plasma expansion, electron-ion recombination and collisions with the background gas molecules. When thermal equilibrium is achieved in the plasma (~100 ns), this is the regime for LIBS and material analysis with fluorescence from highly excited atoms/ions [33,34]. As the plasma cools further, molecular species are formed at later times (~2–50 μs) due to atomic collisions, recombination and oxidation occur after the shock wave collapses [36]. Aerosols and nanoclusters are generated during a nucleation-condensation process with nanoparticle agglomerates re-deposited symmetrically by backward flux around a single spot [35,37]. While time-resolved spectral information using gated ICCDs can help elucidate the plasma dynamics, time-resolved scattering can detect nanoparticle formation and re-deposition on the microsecond timescale [35].

Debris re-deposition was first investigated during UV, excimer micro-structuring of polymers (e.g., Polyimide, PI), responsible for the developing cone structure with debris field dependent on spot shape, inferring strong horizontal forces close to the substrate surface [38]. Apparent field rotations with non-spherical spots were numerically modelled as a purely gas dynamic effect [39]. The debris radius was also shown to follow blast wave theory during excimer ablation of polyimide in air and inert gases [40]. Timedelayed (20–80 μs delay) weak emission detected at the plume periphery consisted mainly of continuum radiation consistent with thermal radiation from solid particles [20]. The difference in ablation geometries observed with a linear array of closely spaced ablation spots inferred stronger plume interactions in the middle [41]. Monte Carlo simulation of the expansion of a copper plasma in the presence of a background gas (Ar) shows compression of the ambient gas atoms by the plume (snowplough effect) for pressures >50 Pa and narrowing of the plume with increasing pressure while predicting backward plume motion and vortical flow at the plume periphery [42]. In the case of excimer ablation of stainless steel in air, (F = 10 J cm−2/308 nm/25 ns) time-resolved ICCD images of the ablation plume up to 10 μs delay show this developing vortical structure after 1 μs and fully developed at 10 μs delay with a high density of nanoparticles deposited at the periphery [35]. Time-resolved light scattering during ablation also demonstrated that debris re-deposition occurred over a time delay of 1–60 μs after the pulse.

Polymers (such as PI, a photoresist) have strong absorption coefficients α > 104 cm−1 in the UV leading to electronic excitation and "cold" ablation via direct bond breaking [43]. In the air, a plasma plume, carbonization (thermal component) and re-deposition around the ablation spot occurs but can be minimized using, for example, H2 as a reactive, ambient gas [44]. Low ablation thresholds Fth < 50 m J cm−2 and precise etch rates from 25 nm/pulse are observed [45]. On steel, UV (248 nm/25 ns) ablation re-deposition generates iron oxide nanoparticles with a wide distribution ϕ = 2–20 nm while at 500 fs/248 nm, this narrows to ϕ = 3 ± 2 nm [46]. The use of NIR, ns pulses result in the incorporation of nitrogen, as well into the re-deposited layer [47].

More recently, there has been an emphasis on temporal plume expansion with bursts of ultrashort pulses [48], and with a very high repetition rate USP Lasers [49,50] driven by the requirement to increase material ablation rates through a better understanding of the plume expansion and plasma absorption by the following pulses.

The interaction dynamics of shock waves during fs laser multi-spot ablation at a water/air interface were recently observed using transient reflection creating colliding supersonic airflows [51,52]. The authors point to the potential for producing focussed ablation plumes relevant to pulsed laser thin film deposition [53,54]. This idea has been investigated here using the observed backward re-deposition patterns from multi-spot ultrafast laser ablation on stainless steel as a novel approach highlighting the resulting plasma plume and shock wave interactions. Symmetric jets, with debris ejected at right angles to the spot axis are reported for the first time, and we present a tentative explanation of the physics behind these observations. The results may well impact the crucial area of Laser-Induced Forward Transfer of nanoparticle thin films at ambient pressure [55,56].

#### **2. Materials and Methods**

Experiments were carried out with two separate laser and optical systems based at the Universities of Liverpool and Mittweida where the temporal pulse lengths available were 10 ps and 600 fs, respectively. Figure 1 shows a schematic of the optical system in Liverpool. The output from a Nd:VAN seeded Regenerative amplifier (High-Q IC-355-800 ps, 10 ps/1064 nm Photonic Solutions Ltd., Edinburgh, UK) was attenuated, expanded (×3), then directed at low AOI to a phase-only Spatial Light Modulator (SLM, Hamamatsu-10468- 03, Hamamatsu Photonics, Hamamatsu, Japan) for generating multi-spot patterns and focussed to a sample mounted on a three-axis stage able to bring the substrate surface to the focal plane. Computer Generated Holograms (CGHs) controlling spot numbers and separation on the substrate were generated using Inverse Fourier Transforms in the Hamamatsu software, while a pick-off optic (not shown) could direct spot patterns to a CCD camera prior to multi-spot ablation. This allows precise electronic control of spot separation with no mechanical movements. A 4f optical system (f1 = f2 = 400 mm) relayed the modulated beam from the SLM (addressed with CGHs) to the input aperture of a galvo system which directed the beams to the target surface, focussed by an f-theta lens (f = 100 mm). A fastmechanical shutter (Thorlabs SH05, Thorlabs Ltd., Ely, UK) allowed the pulse number on target to be varied, synchronised to the scanning software (SCAPS GmbH). The expanding plumes were imaged to a time-resolved spectrometer (Andor Shamrock, model SR303i with intensified CCD, iStar 734, 2 ns gate width, Andor Technology Ltd., Belfast, UK) and synchronised from the 5 kHz Regenerative amplifier.

**Figure 1.** Optical set-up for 10 ps/5 kHz multi-spot ablation of stainless steel in the air. An attenuated and expanded beam from a Nd:Van Regen amplifier is directed to the phase-only SLM at low AOI then imaged via a 4f optical system (f = 400 mm) to the galvo input aperture. Appropriate phase-only CGH's could create arbitrary spot geometries and separations at the steel surface. A time-resolved spectrometer with triggered ICCD (synchronised to Regen amplifier) allowed measurement of the spectral plasma dynamics by focusing the plasma emission on a fibre coupler.

Figure 2 shows a schematic of the optical system for 600 fs/5 kHz double spot ablation (Mittweida). The pump laser (FX-Series, edgewave GmbH, 1030 nm/600 fs, Würselen, Germany) output was attenuated by HWP/polarising beam splitter BS. Afterwards, a second HWP/PBS combination produce two optical lines (1,2) with adjustable power distribution. Line 1 path length was altered by translating mirrors M (1.1) and M (1.2) to synchronise the ablation spots temporally. This synchronisation is proven by the ablation interference patterns with Line 2 at the sample surface. The ablation spot separation was adjusted by slightly tilting mirror M (1.4). Both lines were focussed onto the sample surface by an *f* = 100 mm optic L2. At the focal plane, the laser spot radii are *w*86 = 15 μm. A spot monitor (CCD) on the level of the material surface to enable a precise temporal and geometrical adjustment. A probe laser (Cavilux, Fa. Cavitar Ltd., 688 nm/13 ns, Tampere, Finland) electronically synchronised to the pump laser allowed pump-probe shadowgraphy of expanding plasma plumes. The shadowgraphs were recorded with a 14-bit cooled CCD camera (pco.1600, PCO AG, Kelheim, Germany).

**Figure 2.** Scheme of the optical set-up for double spot ablation with simultaneous shadowgraphy of stainless steel in the air. The pump laser, probe laser and axis and shadowgraph capture were electrically synchronised. The attenuated and split pump beams are focussed and temporally overlapped at the sample surface by one optic lens L2. The temporal (by moving M (1.1) and M (1.2)) and geometrical adjustment (by tilting M (1.4)) can be monitored with a spot monitor in the focal plane.

#### **3. Results**

#### *3.1. 2-Spot Ablation with 10 ps Pulses at 1064 nm/5 kHz Repetition Rate*

The substrate used was ANSI 304 stainless steel and optically polished to a roughness Ra ~50 nm. The focussed single beam diameter was measured from the observed ablation crater diameters with increasing pulse energy [57] and found to be ϕ = 22.2 ± 0.2 μm. No significant variation with multi-spot geometry was found. Single-pulse ablation threshold (*N* = 1) was measured to be Fth = 0.29 ± 0.01 J cm−2 decreasing with pulse number, and incubation coefficient measured to be S = 0.85 ± 0.01 in excellent agreement with the literature [58]. Figure 3 shows the re-deposition patterns observed during two spot ablation with 10 ps laser ablation of stainless steel in air while varying separation and pulse number at fluence F = 0.90 J cm−2, (a) 200 pulses, (b) 400 pulses (c), 800 pulses, all 75 μm separation, (d)–(f) 95 μm separation, (g)–(i) and 145 μm separation, respectively. It required multihundred pulse exposure to observe debris for good optical contrast. The ablation debris is concentrated between the spots at *d* = 75 μm separation with a width comparable to the ablation spot diameter, while at 95 μm separation this narrows to around a 10 μm wide "filament". The concentration of this directed re-deposition, observed at low fluence, is an interesting phenomenon, indicating transient forces during the plume expansion and collisions with the air, resulting in some ablation debris acquiring momentum components preferentially directed along the axis between the spots. There is evidence also that at the

midplane, some debris is expanding normal to the axis. As spot separation increased to 145 μm at this fluence, there was negligible interaction between the plumes. The debris radius near spots Rd ~ 30 μm.

**Figure 3.** Optical images of observed two spot re-deposition during 10 ps laser ablation of stainless steel in air while varying spot separation and pulse number at F = 0.90 J cm−2, (**a**) 200 pulses, (**b**) 400 pulses, (**c**) 800 pulses, all 75 μm separation, (**d**–**f**) 95 μm separation, (**g**–**i**) 145 μm separation respectively. Debris is concentrated in a filament between the spots at 75 μm, thinning to a line about 10 μm wide at 95 μm spot separation. No interaction between the spots is observed at 145 μm separation at this low fluence.

The effect on debris re-deposition of increasing fluence to F = 4.51 J cm−2 on stainless steel is shown in the optical images of Figure 4a–i. At this higher fluence in Figure 4a–c, we observe diverging debris jets ejected normal to the spot axis. This extends to a radial distance of Rjet ~200 μm, well beyond that around each spot Rd ~ 30 μm, inferring that an energetic process may be involved. In Figure 4b,c, with increased exposure, removal of debris from the surface between the spots (in the form of two slightly curved lines) supports the view that strong shock wave interactions between the colliding plumes during ablation may be responsible, clear in Figure 4c,f. These effects essentially disappear at the highest, 145 μm separation with a return to the concentration of material between the spots Figure 4g–i, similar to the patterns at low fluence (Figure 3). The tiny ablation spots along the axis are due to low energy ghost beams appearing during multi-pulse exposure, while the top spot is the remaining zero order. Their presence here helps detect local physical effects of shock wave and air movements affecting the debris motion.

**Figure 4.** Optical images of observed two spot re-deposition during 10 ps laser ablation of stainless steel in air while varying spot separation and pulse number at fluence F = 4.51 J cm−2, (**a**) 200 pulses, (**b**) 400 pulses (**c**) 800 pulses, all 75 μm spot separation, (**d**–**f**) 95 μm separation, (**g**–**i**) 145 μm separation. At this higher fluence in (**a**–**c**), we observe diverging debris jets ejected normal to the spot axis while in (**c**,**f**), removal of surface debris by shock waves is apparent. These effects essentially disappear at the highest separation with a return to the concentration of material between the spots (**g**–**i**). The tiny ablation spots along the axis are due to low energy ghost beams while the top spot remains zero order.

#### *3.2. Ablation Rates and Debris Radii (10 ps Pulses)*

The measured ablation volume/pulse of single and two spot geometries with pulse number N is shown in Figure 5a for a fluence F = 2.9 J cm−2 (Ep = 6 μJ/pulse). These results confirm that ablation rates are essentially independent of spot geometry and spot separation when *N* ≥ 200, whether single or double spot, within experimental error. Hence, the proximity of the spots does not affect ablation rates significantly. The ablation volume/pulse V ~ 5.7 μm3/pulse corresponding to a mass ablation M ~ 0.05 ng/pulse. Crater volumes and geometries were measured with a white light interferometer, Wyko NT3300. The single spot debris radius and jet radius with pulse energy and exposure are shown in Figure 5b on an Ln-Ln plot. This confirms that the debris radii follow a power-law R ∝ E0.47 for a single spot while the much higher Jet debris radii, Rjet ∝ E0.41. This relationship approximates that predicted by blast wave theory [59], however, higher than the expected R ∝ E1/3 which has also been observed by other authors [60].

**Figure 5.** (**a**) Measured ablation volume versus pulse number for single and double spot ablation at *d* = 75 μm, 95 μm and 145 μm separation with fluence F = 2.9 J cm−2/spot (Ep = 6 μJ/spot). Spot geometry does not affect ablation volume for *N* ≥ 200, within experimental error, (**b**) Observed single spot debris radius and jet radius with pulse energy and exposure on an Ln-Ln plot. The radii follow a power-law R ∝ E0.47 for debris and Rjet ∝ E0.41 respectively.

#### *3.3. SEM Imaging of the Debris Fields, Two Spot Patterns*

Figure 6a–d shows a series of SEM images of ablation debris from two spot ablation of stainless steel in the air with 10 ps pulses at low fluence F = 0.9 J cm−2, *N* = 800 pulses and spot separation *d* = 75 μm. The concentration of debris between spots is clear in Figure 6a, with evidence of some material ejected normal to the spot axis in the centre. At low fluence, there are two main components to the debris—particle agglomerates appear on the collision plane (and around the spots), while solid spherical nanoparticles with diameters 50–150 nm appear at the spot edge, Figure 6b. These are likely formed during collisions of the expanding plumes and condensed from the stagnation region at the midplane. Figure 6c (2000×) shows the particle agglomerates near the centre at higher magnification while Figure 6d (35,000×) confirms that the particle agglomerate consists of both solid np's and agglomerated chains of fine np's. There is a significant change in the nature of the debris generated at higher fluence, F = 4.8 J cm−2, as shown in the SEM images of Figure 6e–h. Deposition now produces a thick deposit near the ablation spots while the shock waves lift material from the surface between the spots Figure 6e and deposit well away from the spots, Figure 6f. Figure 6g,h, with increasing magnification show that this jet debris consists of np chain agglomerates and is almost devoid of solid nanoparticles. This is consistent with higher surface temperature achieved during ablation, well above the evaporation temperature, Tev = 2861 K. The high concentration of the np chain agglomerate debris supports the assertion that this material is formed during strong stagnation of the plume collisions in the midplane.

**Figure 6.** SEM images of the ablation debris from two spot ablation of stainless steel at low fluence (F = 0.9 J cm−2 (**a**–**d**) and high fluence F = 4.8 J cm−2 (**e**–**h**), *N* = 800 pulses, spot separation *d* = 75 μm. (**a**) image of the whole debris field, (**b**) 50–150 nm diameter np's observed at spot periphery in (**a**) red ring, (**c**) particle agglomerate at the centre, (**d**) high magnification image of (**c**) showing round np's and np agglomerate, (**e**) High fluence re-deposition pattern showing material lifted directly from the surface by the energetic colliding plumes and shock waves, (**f**) jet debris ejected normal to spots, (**g**) high magnification image of jet np chain agglomerate, (**h**) highest, 100,000× magnification of jet np chain agglomerate which still shows little structure.

#### *3.4. Time-Resolved Plasma Emission Spectra*

Figure 7 shows the two spots (E = 20 μJ/spot, *d* = 75 μm) time-resolved (0–95 ns) plasma emission spectrum of stainless steel in air over Δλ = 395 nm–415 nm. Gate width was set at 5 ns, gate delay interval 2 ns and data accumulated over 50 spectra. The substrate was scanned at 2 mm/s while an f = 125 mm bi-convex lens imaged the whole plasma emission to a fibre coupler (NA ~0.2, Figure 1) then to the spectrometer (Andor Shamrock 303i, 50 μm slit, 1800 L/mm grating) and cooled ICCD camera. The ICCD was triggered externally from the Laser Regenerative amplifier. Continuum dominated the spectra at early times, likely black body radiation from the hot plasma near the surface [21], while spectral line intensities rise sharply, then decrease along with the background continuum as the ablation plasma cools. Spectral line widths also decrease with time as electron density decreases, reducing Stark broadening [61]. The spectral lines in this region have been identified as due mainly to excited neutral atoms of Fe I, Cr I and Mn I: Fe I: 395.667 nm, Cr I: 396.368 nm, Cr I: 396.974 nm, Mn I: 397.708 nm, Fe I: 398.396 nm, Mn I: 399.161 nm, Fe I/Fe II blended line: 403.130 nm, Mn I: 403.307 nm, Fe I: 404.581 nm, Fe I: 406.359 nm, Fe I: 407.581 nm (centre), Fe I: 411.854 nm, Fe I: 413.290 nm, Fe I/Fe II/Fe III blended line: 414.26 nm. Single spot (20 μJ/pulse) spectra were very similar with lower intensities [62].

**Figure 7.** (**a**) Two spots (E = 20 μJ/spot, *d* = 75 μm,) time-resolved plasma emission spectrum (Δλ = 395–415 nm) of stainless steel over time delay 0–95 ns. Continuum dominates at early times while atomic lines appear as the plasma cools. The three intense lines near the centre are due to Fe I: 404.581 nm, Fe I: 406.359 nm, Fe I: 407.581 nm, (**b**) Integrated intensity (Δλ = 395 nm–415 nm) with time for single and double spot (*d* = 75 μm) showing exponential decrease with plasma lifetimes τ1/e ~ 9.2 ± 0.5 ns (single) and 13.9 ± 0.7 ns (double) respectively. The fits are exponential.

Figure 7b shows the time-integrated plasma emission (Δλ = 395–415 nm) for single and double spot (*d* = 75 μm) confirming that the plasma lifetimes τ1/e ~ 9.2 ± 1.0 ns and 13.9 ± 0.7 ns, respectively, hence double spot plasma lifetime increasing over the single spot. These lifetimes are similar to those observed with 20 μJ, 150 fs laser ablation of Al in the air [21]. The short plasma lifetimes are a consequence of both the low pulse energies and rapid collisional cooling of the dense ambient air.

Plasma excitation temperature for single spot ablation (E = 20 μJ) has been estimated by the well-known Boltzmann method [11] from the Fe I line intensities Imn, transition probabilities (gAmn, g degeneracy) and upper energy levels, Em. A plot of Ln(λImn/gAmn) versus Em yielded a linear plot inferring Te ~7500 K for single spot near 40 ns delay time. We also estimated the electron density from the Stark broadening of the Fe I line at λ = 404.58 nm yielding Ne ~ 1018 cm−3.

By removing the ICCD from the spectrometer and placing this at the image plane of the focus lens, (f = 75 mm, M ~ 4), the time-resolved plasma plumes were observed, Figure 8. Gate width here is 5 ns, energy/spot = 20 μJ and spot separation *d* = 75 μm. With spots normal to the optic axis, we can observe plasma expansion, collision and stagnation after 15 ns. When the plumes are imaged parallel to the optic axis and 10–15 ns delay, we see some interesting structure at right angles to the spot axis which may be connected to the jets. The lateral plume expansion velocity can be estimated to be v⊥ ~3.5 × 103 ms−1, decelerating after 15 ns, while the elliptical plume shape which develops supports lateral plasma expansion [35,63].

**Figure 8.** Time-resolved imaging of single and double spot ablation plume emission (20 μJ/spot, *d* = 75 μm) of stainless steel with time delays 5–20 ns. ICCD gate width was 5 ns (**a**) 2-spot normal to optic axis showing plume collisions and stagnation at 20 ns delay (**b**) two spots parallel to optic axis with an indication of jet structure 15–20 ns.

#### *3.5. Two Spot Ablation with 600 fs/5 kHz Temporal Pulses at 1030 nm*

The effect of increasing peak intensity by over 1 order of magnitude was investigated with the experimental system of Figure 2, while maintaining the same repetition rate of 5 kHz. Results are shown in Figure 9 for pulse numbers 200, 500 and 1000 at fluence F = 1.41 J cm−2 (peak intensity I = 2.4 × 1012 W cm−2). Again, we observe symmetric debris ejection or jets at right angles out to a radius of Rjet ~ 150 μm with spot separation *d* = 75 μm. As spot separation increases, we return to the directed, filamentary re-deposition (*d* = 100 μm) between the spots, and just apparent at *d* = 150 μm, *N* = 1000), similar to patterns observed with 10 ps pulses, Figure 3a–f.

The results of high fluence two spot ablation with fluence F = 5.7 J cm−2 (20 μJ/spot, I = 9.6 × 1012 W cm−2) are shown in Figure 10. where plume interactions create strongly diverging jets at spot separation *d* = 75 μm. At *d* = 100 μm separation, jets are more collimated with re-deposition jet radius Rjet > 300 μm. There is a slight tilt in the spot axis here relative to the horizontal. With *N* = 500 and 1000 at *d* = 100 μm, shock wave and plume interactions also remove debris from the surface (compare Figure 4c,f), evidence of the quasi-stationary shock waves [51,52]. At *d* = 150 μm material is again concentrated between the spots as interactions weaken. Peak intensity therefore appears to play a minor role in the debris re-deposition during the plume interactions from 600 fs to 10 ps pulse length on stainless steel in the air.

Figure 11a,b show SEM images of the diverging symmetric jet debris from 600 fs double spot (*d* = 100 μm) ablation of stainless steel in air (20 μJ/spot, F = 4.5 J cm−2, *N* = 1000). This material likely consists of np chain and is concentrated at the jet ends, almost 400 μm from the spots.

**Figure 9.** Resulting debris re-deposition and jet formation after 600 fs laser double pulse ablation with fluence F0 = 1.41 J cm−2 (EP = 5 μJ/pulse). Optical images of stainless steel surfaces after applying for various pulse numbers (200, 500, 1000) at different spot separations (50 μm, 100 μm and 150 μm). Debris jets are visible at 75 μm at the investigated *N*. At larger separations, the debris is concentrated in a filament between the spots (*d* = 100 μm) and becomes apparent at *d* = 150 μm separation, *N* = 1000. Fluence here is 1.5 times higher than with 10 ps pulses (Figure 3).

**Figure 10.** Resulting debris re-deposition and jet formation after 600 fs laser double spot ablation with fluence F0 = 5.7 J cm−2 (EP = 20 μJ/pulse). Within the investigated pulse numbers *N* (200, 500, 1000) and pulse separations *d* of 50 μm, 100 μm and 150 μm, the optical images reveal an interaction between the two spots. Strongly diverging debris jets at 75 μm and collimated jets at *d* = 100 μm are visible. At the largest separation *d* = 150 μm, where the interactions are weaker, the material is again concentrated between spots.

**Figure 11.** Tilted SEM images of the jets after double spot femtosecond ablation of s. steel in the air (EP = 20 μJ/pulse, *N* = 1000, *d* = 100 μm). (**a**) shows the symmetric diverging debris jets (np chain agglomerates) with an accretion of material at the jet ends, almost 400 μm from the double spot centre. (**b**) higher magnification of concentrated jet debris near the jet end. The debris radii near the spots and jet radii are indicated (Rdeb and Rjet respectively).

The debris and jet radii generated with 600 fs pulses are shown in Figure 12. For comparison, the data from 10 ps two spot ablation. Logarithmic fits are included, and the exponents vary from n = 0.33 to 0.47. The fs and ps debris radii are close exponent while the fs jet radii appear to have the lowest, n = 0.33. The higher radii for given energy with 10 ps pulses is due to the smaller spot separation of 75 μm, while this was 100 μm with 600 fs pulses. There may also be differences due to uncertainties in the estimation of the radii.

**Figure 12.** Ln-Ln plots of the two spots debris and jet radii ([R] = μm) measured for 600 fs and 10 ps ablation. The logarithmic fits ([EP] = μJ) are shown to be reasonably close in exponent except for femtosecond ablation which is lower. The spot separations were 100 μm (600 fs) and 75 μm (10 ps), respectively. Error bars represent 1 σ.

The reproducibility of the results reported is supported by the excellent observed fits to the Ln-Ln plots with indicated errors. Ultrafast laser ablation minimises thermal diffusion and melt leading to a deterministic evaporation process [19]. The observed nanoparticle agglomerated jets can be deflected with a significant airflow over the substrate. However, when the surrounding ambient air is stable, the jets emanate at right angles and have the same length, Figure 11a.

#### *3.6. Time-Resolved Shadowgraphy (600 fs/5 kHz).*

Ablation plumes and their development after ablation were imaged with the pumpprobe experimental set-up shown in Figure 2. The probe beam (*λ* = 688 nm, *τ*H = 13 ns) was electronically synchronised to the pump beam (*λ* = 1030 nm). Time-resolved measurements are often carried out with a single pulse exposure–but here, after simultaneous multipulse, multi-spot exposure, we can see the developing plumes and their interactions. The connection axis of two spots is aligned perpendicular to the imaging plane. The total delay time regarding the arrival of the first pulse pair with 5 kHz repetition rate is τ = 1002 μs, hence a delay time of 2 μs after the last, *N* = 6th pulse, Figure 13. There is a strong plasma plume overlap confined to the surface at *d* = 100 μm separation with diverging plume above the spots containing solid (dark) micron size particles, strongly absorbing/scattering. The previous pulses have formed these. These also appear at 150 μm spot separation but disappear entirely at *d* = 200 μm, where the plasma plumes near the surface are distinct and plume absorption much more uniform. This particle is likely the np chain agglomerate observed on the surface, shown in the SEM images of Figure 11. which appear in the jets due to the strong plume interactions. Note in the shadowgraphs a visible expanding spherical shock wave from the last pulse pair overlaps the material previously ejected. The shock wave speed, from the time delay, *v*s ~375 m/s just above the speed of sound while the lateral plasma expansion has slowed significantly to around *v* ≤ 50 m/s after 2 μs delay.

**Figure 13.** Time-resolved shadowgraphs of the interacting plumes with two spot ablation in the air after *N* = 6 and delay time of τ = 2 μs after ablation and τ (total) = 1002 μs. These images demonstrate that large condensed absorbing particles are apparent in the debris above the spots when the separation is smallest and where the plume collisions are strongest, (**a**) *d* = 100 μm, (**b**) *d* = 150 μm, (**c**) *d* = 200 μm where these particles have all but disappeared. Note that the spherical shock waves from the last pulse overlap the material ejected beforehand.

#### **4. Discussion**

Laser ablation on metal in the ambient atmosphere can be likened to a mechanical detonation, creating an ablation plume expanding supersonically against the background gas, causing a shock wave. As the plume does work by expanding against atmospheric pressure, the expansion velocity decreases with time, and the shock wave radius is given by Taylor's blast wave theory [59],

$$R = S(\gamma) \left( E / \rho\_0 \right)^{1/5} t^{2/5} \tag{1}$$

where *S*(*γ*) ~1 is a function only of the air specific heat *γ* (~1.4), *E* is the energy released, *ρ*0 is the undisturbed background gas density and t is the time after ablation. The effect is to force most of the air within the shock front into a thin shell just inside the front, compressing and heating the air. The initial pressure driving the front *p*max *p*0, (1 atm). When *p*max ~ *p*0, the self-similar solution of Equation (1) is no longer valid. In this case, the radius for this upper limit is given by [20],

$$\mathcal{R} = (E/p\_0)^{1/3} \tag{2}$$

The ablation plasma plume follows behind the shock front and a contact discontinuity can be observed inside the shock front [33]. For the range 2 μJ ≤ *E* ≤ 20 μJ, single spot debris radii range from 40–120 μm while *R* calculated from Equation (2) yields 270 μm ≤ *R* ≤ 580 μm. The debris radius can be related to the shock radius through RD ~ fR wheref<1 and here, f ~ 0.15–0.2, similar to that observed by other authors [60]. From the time-resolved plasma emission (10 ps), backward re-deposition with single spot ablation likely starts about ~100 ns after ablation and from the pump-probe observations, (600 fs) continues over 10 μs and longer. The debris power laws yield Rdeb ∝ E0.3–0.5 for both single spot and jets but the much larger jet radii if applied to RD ~ fR yields f ~ 0.7–0.8 using Equation (2) inferring that an additional energetic process is involved during the plasma plume interactions.

Multi-spot ablation in vacuum between two independent and closely spaced seed plasmas collide due to the lateral plasma expansion. With a low density and high relative velocity, plasmas tend to interpenetrate, relevant in collisionless astrophysical plasmas, leading to ion reflection and particle acceleration [64]. However, with higher density and a low relative velocity, the plasmas rapidly decelerate at the collision plane, forming a stagnation layer (SL). Accretion and compression of the material within the SL leads to a local increase in density and temperature. The degree of stagnation can be described by a collisionality parameter *ξ* = *d*/*λii*, where *d* is the distance between the two plasmas and *λii* is the ion-other ion mean free path, given by, [65,66]

$$
\lambda\_{ii\ (1-2)} = \frac{4\pi\varepsilon\_0^2 m\_i^2 v\_{12}^4}{q^4 Z^4 n\_i \ln(\Lambda\_{1-2})} \tag{3}
$$

where *ε*0 is the permittivity of free space, mi is the ion mass, *v*12 is the relative ion flow velocity (prior to impact), *q* is the elementary charge, *Z* is the average ionisation state, *ni* is the plasma density at the collision plane and ln (Λ1–2) is the Coulomb Logarithm, a slowly varying function, with a value *O* (10–20) [67]. The parameter *ξ* is very sensitive to the relative plasma velocity term *v*124 while only linearly dependent on separation d. If used in ambient air, this description is more complex than in vacuum due to the presence of shock waves and rapid plasma deceleration due to collisions with the air. From timeresolved plasma emission at 15 ns, Figure 8. We estimated that *v*⊥ ~ 3.5 × 103 m s−1 so that *v*12 = 2*v*⊥ ~ 7 × 103 m s−1. Inserting this value into Equation (3), we find *λii* ~ 3.1 μm so that the plasma collision region is limited to the micron scale with collisionality parameter *ξ* = *d*/*λii*= 75 μm/3.1 μm ~ 24.

In a plasma, the Debye length is the characteristic length over which electrons and ions can be separated, and in an ideal plasma (e.g., astrophysical) has many particles per Debye sphere *ND* given by

$$N\_D = N\_\mathfrak{t} \left(\frac{4\pi}{3}\right) \lambda\_D^3 \gg 1 \tag{4}$$

The classical plasma parameter g = 1/*NeλD*3 << 1 so that collective effects dominate the plasma. In laser ablation, we can calculate g for each plasma plume knowing *Ne* and *λD* which is given by [67] *λD* = 743 (Te) 1/2 (*Ne*) −1/2 where the electron temperature Te is in eV while *Ne* is in cm−3. At delay time τ = 50 ns, we measured Te (10 ps, Ep = 20 μJ) from a Boltzmann plot to be ~7500 K (1 eV = 11,600 K) hence plasma temperature Te ~ 0.65 eV while *Ne* ~ *Ni* ~ 1018 cm−3 from measured Stark broadening of the Fe I line at *λ* = 404.58 nm. The Debye length in each plasma plume is then *λD* ~ 6.0 × 10−7 cm (6 nm) while the particle number in a Debye sphere is, from Equation (4), *ND* = 0.92. The collision parameter g = 1/*NeλD*3 = 4.6 hence the plasma plumes are each highly collisional as expected. We

can also estimate the Coulomb Logarithm [67], LnΛ = 9*ND*/*Z* = 8.3 assuming *Z* ~ 1 is the average ionisation state, supported by the time-resolved spectroscopy, Figure 7a.

Recently, the interaction of two colliding Al plasmas (in the air) and their shock waves were observed in two spot nanosecond laser ablation using shadowgraphy, schlieren images and interferometry [68]. From refractive index profiles, expanding shock fronts and their reflections were observed after collision along with compression and stagnation and of the air behind the shock fronts. With pulse fluence F ~ 27 J cm−2, focus intensity I ~ 1.6 GWcm−2 and 1 mm spot separation, compressed air density reached *n*0 ~ 5 × 1020 cm−3 behind the interacting shocks while plasma density *ni* ~ 5.1018 cm−3. The stagnation was described as "soft" [69] with a plasma hill developing over times scales from 0.7–3 μs.

More recently, stationary transient straight shockwaves (on a 10 ns timescale) were detected during multi-spot fs ablation of water in the air creating supersonic air flows which collided [51,52]. Fluence per spot F ~ 18 J cm−2 with peak intensities I ~ 1014 W cm−2 and spot separations from d ~ 14–20 μm. The length of the shock waves was related to the local flow velocity at a given radial position from the spots and observed when the relative speed of shock wave propagation exceeded the velocity of sound in the air. The symmetric geometry of the observed two spot transient shock waves is highly reminiscent of the observed patterns of debris removed between double spot ablation in this work at higher fluence (Figures 4, 6 and 11) with ps and fs pulse lengths. This indicates that material removal between spots during ablation at high fluence may be evidence of shock wave interactions.

#### **5. Conclusions**

The jets appearing at higher fluence (and small spot separation) could be created as follows using a simplified 2D approach. When the high-pressure shock waves meet at the collision plane, the pressure due to each shock wave *p*max = 0.155 E/R3 (*γ* = 1.4) [59]. If we set E ~ Ep = 20 μJ and R = 50 μm, *p*max = 2.48 × 107 Pa ~248 atm. The air between these shocks will be highly compressed by this "piston" as they meet and the curved expanding shocks could compress, accelerate and expel the air symmetrically in a diverging jet in both directions normal to the axis, Figure 14a. The plasmas following then stagnate and ion Coulomb repulsion off-axis at low impact parameters converts their axial momentum to transverse momentum assuming elastic collisions. The observed material lifted from the surface leaves patterns very similar to the transient stationary shock waves observed during multi-spot induced supersonic air flows at a water/air interface [51,52].

At low fluence and larger separations, shock wave pressures are much reduced, plasma density decreases, reducing stagnation to "soft", hence allowing plasma interpenetration between the spots. Weak shock waves pass through each other and likely reflect back along the axis from the plasmas. Multiple collisions can decelerate ions and a degree of ion reflection may occur [64] so that particles acquire momentum directed along the spot axis then arrive at the surface through backward flux, Figure 14b.

Two colliding air breakdown plasmas created with high energy, nanosecond laser pulses demonstrated reflection and transmission of expanding shock waves depending on pulse energy and plasma seed separations [70]. One might ask if oxidation during ultrafast laser ablation is significant here, but we have also observed the Jets in pure Nitrogen at 1 atm (10 ps, not reported here) with the aid of a gas cell. We estimate that oxidation of Fe atoms to Fe3O4 could at most add 10% additional energy, not nearly enough to explain the phenomena here.

Time-resolved plasma emission yielded important plasma parameters such as lifetime and temperatures, and plasma collisions between spots show plume stagnation, resulting in bi-directional jets. The effects of quasi-stationary shock waves are very clear, compressing the air and lifting debris from the surface between the spots. Time-resolved shadowgraphy with 600 fs pulses showed micron size particles (which constitute the jets) only when the ablation spots were in close proximity. In fact, the re-deposited jet debris consists of agglomerated nano-particle chains.

**Figure 14.** Schematic diagrams summarising current view of two spot ablation at high and low ablation (**a**) at high fluence, the air between the shocks is compressed by a "piston" as they meet and the curved shocks will accelerate and expel the air symmetrically in a diverging jet in both directions normal to the axis. The plasmas then stagnate and Coulomb interactions off-axis at low impact parameters convert their axial momentum to transverse momentum, (**b**) at low fluence and larger separations, shock wave pressures are much reduced, plasma density decreases, reducing stagnation to "soft", hence allowing plasma interpenetration and directed deposition between the spots.

Our observations of the directed deposition at low fluence and jets during highfluence, multi-spot ultrafast laser ablation of stainless steel in air constitute a potentially interesting new area of research, providing a basis for further investigations to build a thorough understanding of these complex phenomena. Our study and its findings are considered important for a deeper understanding of the interaction of multi-spot arrays in laser patterning [71,72] and enhanced laser deposition of thin films [53,55]. The phenomena of plasma collisions and shocks observed and explored here should also be relevant to astrophysics and physics research, such as in magnetic re-connection [73] and phenomena such as bi-directional jet formation and particle acceleration [74,75] observed at ultrahigh laser intensities (1014–1015 W cm−2) in two spot ablation with high transient B fields.

**Author Contributions:** Conceptualization, G.D., U.L. and W.P., methodology, W.P., J.S. and S.E.; validation, T.Z. and S.K.; formal analysis, W.P., T.Z. and S.K.; investigation, T.Z. and S.K.; resources, G.D. and U.L.; data curation, T.Z. and S.K.; writing—original draft preparation, W.P.; writing—review and editing, all authors; visualisation, J.S., S.E.; supervision, G.D. and U.L.; project administration, G.D., S.E., J.S. and U.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors would like to express their thanks to Andy Snaylam (Laser group, Liverpool University) and the SEM SRF Albert Crewe Centre for Electron Microscopy (Liverpool) for their support & assistance in this research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Squared Focal Intensity Distributions for Applications in Laser Material Processing**

**Henrike Schlutow 1, Ulrike Fuchs 1, Frank A. Müller 2 and Stephan Gräf 2,\***

**Abstract:** Tailored intensity profiles within the focal spot of the laser beam offer great potential for a well-defined control of the interaction process between laser radiation and material, and thus for improving the processing results. The present paper discusses a novel refractive beamshaping element that provides different squared intensity distributions converted from the Gaussian output beam of the utilized femtosecond (fs) laser. Using the examples of surface structuring of stainless-steel on the micro- and nano-scale, the suitability of the beam-shaping element for fs-laser material processing with a conventional f-Theta lens is demonstrated. In this context, it was shown that the experimental structuring results are in good agreement with beam profile measurements and numerical simulations of the beam-shaping unit. In addition, the experimental results reveal the improvement of laser processing in terms of a significantly reduced processing time during surface nano-structuring and the possibility to control the ablation geometry during the fabrication of micro-channels.

**Keywords:** beam shaping; squared top-hat; fs-laser; laser surface structuring; laser-induced periodic surface structures; micro-channels

#### **1. Introduction**

In modern laser material processing, beam shaping is receiving growing attention, as it allows the focal intensity distribution to be adjusted to the interaction process between the laser radiation and the material [1–3]. This makes laser processing even more flexible, precise and effective. Moreover, the range of achievable types of structures is extended [4–6] and surface properties can be adjusted in a targeted manner [7–9]. Top-hat profiles are characterized by a uniform intensity distribution over the beam cross-section in a specific working plane. They can be generated by a large variety of available beam-shaping techniques, including apertures [10], refractive [11–13] and diffractive optical elements [10,14], as well as beam integrators or homogenizers [10,15]. Finally, a spatial light modulator provides a flexible optical device that allows to create various intensity profiles based on the dynamical modification of the amplitude, phase and polarization of the incoming laser beam [16–18].

The advantage of refractive beam shapers is provided by their simple structure and a very high conversion efficiency. They are less sensitive to wavelength changes, easier to manufacture and can withstand higher laser powers when compared to diffractive elements. Therefore, the present work aims at the introduction of a highly flexible and compact refractive, freeform beam-shaping element that is suitable for laser material processing applications.

Large-area structuring of surfaces on the micro- and nano-meter scale usually requires scanning of the focused laser beam across the material surface. In the case of the commonly used Gaussian intensity distribution, a certain inhomogeneity of the deposited energy

**Citation:** Schlutow, H.; Fuchs, U.; Müller, F.A.; Gräf, S. Squared Focal Intensity Distributions for Applications in Laser Material Processing. *Materials* **2021**, *14*, 4981. https://doi.org/10.3390/ma14174981

Academic Editors: Jörg Krüger and Jörn Bonse

Received: 14 July 2021 Accepted: 27 August 2021 Published: 31 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

occurs, which in turn has an influence on the processing result. A first improvement was demonstrated in this context by using a beam shaper that allows to generate a rotationally symmetric top-hat profile [19]. However, due to the rectangular symmetry of the (unidirectional) scan procedure with a specific overlap of the adjacent scan lines, the use of a squared top-hat seems to be even more advantageous with respect to homogeneous energy coupling. This squared top-hat is derived from the collimated Gaussian intensity profile of the input beam by combining the investigated refractive beam shaper with a focusing optics. The latter can be a simple focusing lens as well as a complex f-Theta objective, as is very often used for laser-based surface structuring. Therefore, a main objective of this study was to evaluate whether the x-y-symmetry of the focal distribution can further improve material processing when working with a scanner and an f-Theta lens. This aspect plays an essential role with regard to modern laser systems and their high repetition rates, whose maximum utilization ensures high processing efficiency. For this purpose, the beam-shaping element was practically tested in different applications on stainless-steel as a substrate material. This includes surface nano-structuring in the form of laser-induced periodic surface structures (LIPSS), which allow to engineer surfaces with versatile functional properties [7,20]. Furthermore, the ablation-based micro-structuring of the sample surface using femtosecond lasers (fs-lasers) was demonstrated as a conventional processing technique [21].

#### **2. Materials and Methods**

#### *2.1. Focal Beam-Shaping and Characterization*

The most common approach for shaping focal intensity distributions is based on Fourier optics, where the focusing optics transforms the incoming beam profile into the desired shape. This is well-explained for rotational symmetry in [19,22]. However, for non-rotationally symmetric focal intensity distributions, the symmetry of the beam-shaping element needs to be broken, while the focusing optics remains the same. Therefore, the beam-shaping element becomes a freeform surface. As a phase plate, it imprints a phase modulation on the incoming plan wave of the collimated beam, which stays collimated as this modulation is in the order of λ (PV) and thus very small. The focusing optics performs the Fourier transformation of this modified wave front and therefore determines the width of the generated profiles inverse proportional to the numerical aperture (NA). Consequently, this approach is modular with regards to the focusing optics. The phase plate (a|SqAiryShape) itself is a freeform distribution based on a super Gaussian profile:

$$\mathbf{z}(\mathbf{h}) = \mathbf{H}\mathbf{e} - ^2(\frac{4\mathbf{h}^2}{\mathbf{D}^2})^{\frac{\mathbf{M}}{2}} \tag{1}$$

of order M = 8, diameter D and plateau height H in x-y-symmetry. The element is made of fused silica via CNC-based grinding and polishing. It can therefore withstand high-power laser intensities as well as ultrashort laser pulses, as employed in the experimental studies of the present work.

As a measure of how much the beam shaping will be affected by diffraction effects, the β-factor:

$$
\beta = \frac{2\pi \text{RD}}{\text{f}\_{\text{L}}\lambda} \tag{2}
$$

can be used [23], with R being the half diameter of the entrance beam, D being the half diameter of the focal distribution, fL being the focal length and λ being the wavelength. Using the specific parameters of the present study (R = 10 mm, D = 25 μm, fL = 100 mm and λ = 1030 nm) leads to β = 15.25. Having such a large value ensures stable beam-shaping profiles, which are essential for well-controlled material processing.

The squared focal intensity distribution was characterized in order to compare the experimental results with simulations. For this purpose, a supercontinuum laser (SC-400-4, Fianium, Southampton, UK) was used as a light source, from which λ = 1064 nm was

selected via a spectrometer (LLTF Contrast, Fianium, Southampton, UK), coupled to a fiber, collimated (ACM25-20-L-C-1064, asphericon, Jena, Germany) and expanded (BAM25-175- D-C-1064, asphericon, Jena, Germany) to a collimated beam with 10 mm diameter. This beam fulfilled the input conditions for the freeform beam-shaping element. All following optical components used for characterization of the beam shaping are depicted in Figure 1. The beam-shaping element (a|SqAiryShape) was placed in the collimated beam, which was subsequently focused by a plano-convex lens with 300 mm focal length (#38-631, Edmund Optics, Mainz, Germany). The focal region was then imaged onto a camera (SP928, Ophir, Germany) by using a microscope objective (Aspheric Lens Kit 5720-B, Newport, Irvine, CA, USA). The setup was calibrated to obtain the real physical size of the focal distribution.

**Figure 1.** Experimental setup for characterization of the freeform beam-shaping element.

In the simulation, only the freeform beam-shaping element and the plano-convex focusing lens were considered. The focal intensity distribution was calculated via the point spread function with CodeV (Synopsis). This calculation is based on the squared Fourier transform of the pupil function [24]. The latter includes all aberrations and thus the phase manipulation introduced by the beam-shaping element.

Close to the focal point, various intensity distributions, such as a squared top-hat or squared donut, emerge. The following laser processing experiments were carried out with the squared top-hat, which occurs closest to the focal point. This is the top-hat shown in Figure 2. Comparison with the simulation, especially through the cross-sections in the bottom of Figure 2, reveals close agreement between measurement and simulation in terms of shape and size of the focus spot.

A sequence of those cross-sections taken along the optical axis is plotted in Figure 3. The upper panel (Figure 3a) shows the beam waist of the Gaussian input beam without the freeform beam-shaping element. With the beam-shaping element (panel b), the Gaussian distribution is notably flattened to a top-hat. Note that the beam shaping only occurs close to the focal point (z = 0 mm), whereas the profile remains almost Gaussian for larger distances. The blue-shaded distribution corresponds to the one shown in Figure 2.

**Figure 2.** Color maps and cross-sections of the squared focal intensity distribution of the top-hat closest to the focal plane (λ = 1064 nm, f = 300 mm) normalized to the maximum intensity for comparison between measurement and simulation.

**Figure 3.** Profile cross-sections measured along the optical axis (z) in the focal region (**a**) without and (**b**) with freeform beam-shaping element.

#### *2.2. Material Processing Setup and Sample Preparation*

The operation of the beam-shaping element was investigated by means of surface structuring experiments using the experimental setup illustrated in Figure 4. The linearly polarized Gaussian output beam (M2 ~ 1.08) of the diode-pumped Yb:KYW thin-disc fs-laser system (JenLas D2.fs, Jenoptik, Germany) with a central wavelength λ = 1025 nm is characterized by a pulse duration τ = 300 fs (FWHM) and pulse energies Eimp ≤ 40 μJ at a repetition frequency frep = 100 kHz. Before passing the beam shaper, the output beam was widened by a beam expander (5×). The shaped beam was then focused onto the sample surface by a galvanometer scanner (IntelliScan14, Scanlab, Puchheim, Germany). The focal length of the associated f-Theta objective (JENar, Jenoptik, Jena, Germany) was fL = 100 mm.

**Figure 4.** Experimental setup used for surface structuring experiments.

The structuring experiments were carried out on commercially available austenitic stainless-steel (X2CrNiMo17-12-2, Outokumpu, Krefeld, Germany) at normal incidence and under air atmosphere. Prior to laser irradiation, the sample surface was manually grounded and polished to a mirror finish, with an average surface roughness Ra ~4 nm. This procedure involves the use of SiC abrasive paper of 800, 1200 and 2400 grit and 6, 3 and 1 μm polycrystalline diamond suspension, respectively. Besides, the sample surfaces were ultrasonically cleaned in acetone and isopropanol both before and after laser structuring. The morphology of the samples was characterized by scanning electron microscopy (SEM) (SigmaVP, Carl Zeiss, Oberkochen, Germany) at an accelerating voltage of 5 kV using a secondary electron detector. The spatial periods and the orientation of the LIPSS were quantified by 2D-Fourier transform analysis (2D-FT) of the SEM micrographs. White light interference microscopy (WLIM) (CCI HD, Taylor Hobson, Weiterstadt, Germany) was utilized to evaluate the surface topography. The microscope was equipped with a 50× objective.

#### **3. Results and Discussion**

#### *3.1. Single-Spot Experiments*

The first part of the experiments was to investigate whether the theoretical calculations for the focal intensity distribution correlate with the geometry of the ablated region within the focal spot on the material surface. For this purpose, the surface was irradiated with 10 single fs-laser pulses with a temporal pulse spacing of 10 μs (frep = 100 kHz) and a single-pulse energy of Eimp = 6.5 μJ. Without the beam-shaping element in the beam path, the utilized f-Theta objective results in a focal spot diameter of 2wf = (24 ± 0.5) μm, determined using the method of Liu [25]. This corresponds to a fs-laser peak fluence of F = 2Eimp/(πwf 2) = 2.9 J/cm2 at the pulse energy used.

Figure 5 shows SEM micrographs of the ablation spots obtained at different distances (z-positions) of the stainless-steel surface from the focusing objective using the freeform beam-shaping element. At the reference position (z = 0 mm) and at z = −0.6 mm, squared top-hat intensity profiles are available that result in square material ablation of size (30 × 30) μm2 and (40 × 40) μm2, respectively. As described later, these intensity profiles are predestined for numerous laser material processes that are based on scanning procedures with x-y-symmetry (Figure 6). Assuming a perfectly rectangular fluence distribution of both top-hat profiles, the larger focal dimensions together with the pulse energy Eimp = 6.5 μJ result in a lower fs-laser peak fluence of F = 0.72 and F = 0.41 J/cm2 respectively, when compared to the initial Gaussian focal spot. It should be noted that the peak fluence, F, is identical to the average fluence, Fav, of the fluence distribution. The

larger fluence atz=0 compared to z = −0.6 mm is indicated by the stronger melt formation in the ablation spot fabricated with the smaller top-hat profile.

**Figure 5.** SEM micrographs of the stainless-steel surface after irradiation of 10 linearly polarized single fs-laser pulses at different z-positions of the sample surface within the focal intensity distribution. For comparability, the single-pulse energy was kept constant at Eimp = 6.5 μJ.

**Figure 6.** LIPSS structuring of stainless-steel surfaces using the squared top-hat profile (z = 0 mm) of the linearly polarized fs-laser beam (2wf ~30 μm, Eimp = 3.5 μJ, v = 1 m/s, Δx = 10 μm): (**a**) SEM micrograph and (**b**) 2D-FT spectrum and its corresponding cross-section calculated from the SEM micrograph quantifying the homogeneity and orientation of the LIPSS.

In addition to the top-hat profiles, the beam-shaping element provides other intermediate stages, which are characterized by more exotic intensity distributions. At z = −1.0 mm, a donut-shaped focal intensity distribution with square beam cross-section can be observed, characterized by the zero-intensity center on the beam axis. In contrast, at z = −0.35 mm, the intensity distribution is almost completely reversed, which is why a stronger intensity can be observed at the center. Both examples are intensity distributions suitable for the generation of novel surface structures (e.g., bio-inspired surfaces) [5].

The SEM micrographs of the single-spot experiments reveal a well-defined modulated surface topography within the ablation crater, whose periodicity is in the sub-μm range and therefore much smaller when compared to the diffraction-limited focal spot diameter of a few tens of μm. These periodic patterns are termed laser-induced periodic surface structures (LIPSS), which for metals are typically oriented perpendicular to the linear beam polarization [26]. The most widely accepted theory describes the formation of LIPSS as a result of interference of the incident laser radiation with surface electromagnetic waves generated by scattering at the nano-rough surface [27,28]. These scattering processes are possible despite the polished initial surface, since the LIPSS formation process is generally a multi-pulse process.

As Figure 5 shows, the exact morphology of the LIPSS pattern within the focal spot depends on the fluence distribution in the beam cross-section, i.e., on the intensity profile. A particular contribution is provided by melt formation in the intensive regions of the beam cross-section. However, it becomes evident that a particularly homogeneous LIPSS distribution can be realized by the appropriate choice of the z-position, which suggests here a constant fluence in the beam profile. This is a key aspect for homogeneous large-area structuring with LIPSS [29]. In addition, the possibility to control the LIPSS distribution within the focal spot by choosing a certain beam profile leads to an extension of the current state-of-the-art and, thus, to a higher flexibility of the direct-writing structuring method. The fluence values calculated above for the top-hat profile are comparable with the fluence required for ablation-based LIPSS formation on stainless-steel at comparable irradiation parameters [30,31].

#### *3.2. Laser Surface Structuring on the Nanoscale*

Based on the demonstrated LIPSS formation within the focal spot, the nanostructures can be scaled to large surface areas by a relative movement between the focused fs-laser beam and the sample surface. This is particularly achievable with metals due to the linear absorption of the fs-laser radiation. In this context, the specific optical properties make stainless-steel a well-suited candidate for the production of highly regular LIPSS [31]. This concerns in particular the exact alignment and a small dispersion of the spatial frequencies of the LIPSS, which are coherently linked to each other despite the unidirectional scanning procedure [30,31].

The SEM micrograph in Figure 6a shows the stainless-steel surface obtained from scanning the squared top-hat profile (z = 0 mm) with a scan speed of v = 1 m/s and a line spacing of Δx = 10 μm across the sample surface. The single-pulse energy was reduced to Eimp = 3.5 μJ, i.e., to F = 0.39 J/cm2, in order to avoid the extensive melt formation observed in Figure 5. The calculated 2D-FT spectrum (Figure 6b) confirms the high regularity of the LIPSS aligned perpendicular to the linear beam polarization by the very narrow characteristic peaks. The spatial periods of the LIPSS were determined to be around (950 ± 20) nm.

The obtained regularity is comparable with our previous work on rotationally symmetric beam shaping [19], although the following main differences should be discussed:

**Figure 7.** Overlap and contribution of individual pulses during scanning of a single line using (**a**) the Gaussian focal intensity distribution and (**b**) an ideal top-hat beam.

#### *3.3. Surface Micro-Structuring*

In order to investigate the influence of the intensity distribution on the micro-structuring of surfaces, channels were created using the top-hat profile (z = 0 mm) and 10 over-scans with v = 0.1 m/s. Earlier studies with rotationally symmetrical top-hat and donut-shaped beam profiles already showed that steeper channel walls and a better surface quality (especially at the bottom) can be achieved in comparison to the Gaussian profile [19].

These findings can also be observed in Figure 8, which shows 3D-WLIM micrographs of channels produced with the top-hat profile using different pulse energies, Eimp. It becomes evident for moderate pulse energies that the homogeneous fluence distribution across the beam profile leads to a very high quality of both the walls and the bottom of the channel (Figure 8a,b). The corresponding arithmetical roughness values of the channel surface were measured via WLIM to be Sa = 0.1072 μm for Eimp = 2.5 μJ and Sa = 0.1844 μm for Eimp = 4.6 μJ. As indicated by the cross-sections in Figure 8d, the corresponding depth of the channels increases with increasing pulse energy, from about 6.2 μm at Eimp = 2.5 μJ (Figure 8a) to about 12.5 μm at Eimp = 7 μJ (Figure 8c). In addition, a slight increase in the channel width can be observed. However, the WLIM micrographs indicate that at very high pulse energies and thus large laser peak fluences, an inhomogeneous and rough channel bottom results from extensive ablation and the melt formation already detected in the single-spot experiments. The irregularity is expressed in particular in hole-like depressions that are not adequately displayed by the perspective of the 3D micrograph in Figure 8c. However, a qualitative measure is given by the grey-shaded area in the cross-section of Figure 8d that illustrates the remarkably larger variation, especially at the channel bottom, when compared to smaller Eimp. This is also confirmed by the larger surface roughness that was measured to be Sa = 0.5174 μm for Eimp = 7 μJ.

**Figure 8.** WLIM micrographs of the stainless-steel surface illustrating the 3D geometry of channellike structures obtained from 10 over-scans with v = 0.1 m/s using the squared top-hat profile at z = 0 mm and different single-pulse energies: (**a**) Eimp = 2.5 μJ, (**b**) Eimp = 4.6 μJ and (**c**) Eimp = 7 μJ. A comparison between the cross-sections of the channels is illustrated in (**d**). Here, the solid line corresponds to the mean profile and the shadow illustrates the variation of the profile series measured along the channel.

Structuring with the donut-shaped intensity profile at z = −0.9 mm allows to create channels with steep walls and a w-shaped cross-sectional geometry that is determined by the specific intensity distribution (see Figure 5). According to the cross-sectional height profile of the channel (Figure 9c), 10 over-scans with v = 0.1 m/s and Eimp = 4.6 μJ result in a width of about 40 μm and a maximum ablation depth of about 10 μm. The channel walls are characterized by their very homogeneous and regular surfaces in the form of superimposed nanoscale LIPSS (Figure 9a). The high quality is confirmed by the fact that the modulation depth of LIPSS is typically in the order of 200 nm [30]. Contrary to a rotationally symmetrical donut, the used squared intensity profile provides the zero-intensity center and two axes in x- and y-directions with almost zero intensity. Consequently, the areas with lower (zero) intensity are not passed over by highly intensive beam parts during scanning in these directions, resulting in a channel with a less ablated central area (Figure 9b). This could be used, for example, to periodically propagate these w-shaped structures to large surface areas based on a suitable scan line distance. In dependence on the focal length of the focusing optics, this would result in a microscale periodic grating with periods in the range of 10–20 μm. Similar surface structures with comparable structural sizes can be obtained, e.g., from Direct Interference Laser Patterning, which is used among others to generate surfaces with functional properties (optics, tribology, biology, wetting) [32].

Generally, the results show that the donut and top-hat profiles provided by the beam-shaping element are advantageous for realizing micro-channels with specific wall inclination angles as well as tailor-made cross-sectional geometry. This option is also available (in situ) during the structuring process, since the individual intensity profiles can be addressed in a relatively simple way by changing the z-position of the sample surface.

**Figure 9.** Micro-channel fabricated on stainless-steel using the squared donut-like intensity distribution at z = −0.9 mm, 10 over-scans with v = 0.1 m/s and a single-pulse energy of 4.6 μJ: (**a**) SEM micrograph, (**b**) WLIM micrograph and (**c**) cross-sectional height profile. The blue line in (**c**) corresponds to the mean profile and the grey shadow illustrates the profile series measured along the channel in (**b**).

#### **4. Conclusions**

The present study deals with a versatile refractive beam shaper that allows to convert the initial Gaussian intensity distribution of the fs-laser beam into different squared beam profiles at different z-positions. Numerical simulations for the beam-shaping unit were verified experimentally by spatial beam profile measurements and by surface structuring on the micro- and nano-scale. It was demonstrated that the squared top-hat profiles available in the focal region are advantageous to significantly reduce the processing time during large-area LIPSS structuring, while preserving their high regularity. Using the beam shaper for the generation of channel-like structures via multiple line scans, it was shown that the different focal intensity distributions resulted in precise channels with a controllable cross-sectional geometry.

**Author Contributions:** All authors contributed equally to conceptualization, methodology, validation, writing—original draft preparation, and visualization. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data underlying this article will be shared on reasonable request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **High Aspect Ratio Structuring of Glass with Ultrafast Bessel Beams**

**Christian Vetter, Remo Giust, Luca Furfaro, Cyril Billet, Luc Froehly and Francois Courvoisier \***

FEMTO-ST Institute, University Bourgogne Franche-Comté, CNRS, 15B Avenue des Montboucons, CEDEX, 25030 Besançon, France<EMAIL_ADDRESS>(C.V<EMAIL_ADDRESS>(R.G<EMAIL_ADDRESS>(L.F<EMAIL_ADDRESS>(C.B<EMAIL_ADDRESS>(L.F.) **\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Controlling the formation of high aspect ratio void channels inside glass is important for applications like the high-speed dicing of glass. Here, we investigate void formation using ultrafast Bessel beams in the single shot illumination regime. We characterize the morphology of the damages as a function of pulse energy, pulse duration, and position of the beam inside fused silica, Corning Eagle XG, and Corning Gorilla glass. While a large set of parameters allow for void formation inside fused silica, the operating window is much more restricted for Eagle XG and Gorilla glass. The transient formation of a molten layer around voids enables us interpreting the evolution of the morphology with pulse energy and duration.

**Keywords:** glass; ultrafast laser processing; Bessel beam

#### **1. Introduction**

Cutting and drilling glass is an important technological problem, particularly when the width over which the material has to be modified is much smaller than the length of the modification. This is high aspect ratio structuring. In this context, ultrafast laser processing is very attractive because the infrared wavelength and the very high intensity of the pulses enable energy deposition in three dimensions, including inside the bulk of glass [1]. The laser-induced micro/nano-structures can be index modifications, nanogratings or even voids [2]. These modifications can be combined with chemical etching, which led to a number of advances in the field of photonics, microfluidics chips, mechanics [3,4].

Because point-by-point scanning is relatively slow, a number of techniques were developed to improve the overall speed for the micro-structuring. In addition, elongated voids cannot be formed by scanning, since a nano-void produced by a first pulse would be erased by the microexplosion produced by the second pulse next to it [5]. Bessel beams were revealed to be excellent candidates to control the formation of high aspect ratio elongated structures inside transparent materials [6,7]. Bessel beams are formed by a conical superposition of plane waves, where the central hotspot maintains its diameter over the whole length of the interference region [8,9]. Importantly, when the focusing angle of the Bessel beam is sufficiently high, the conical flow of energy crosses the high-intensity region only in the central hotspot. This drastically reduces the influence of the Kerr effect, so that longitudinally-invariant structures can be processed within dielectrics [10–12].

Laser structuring with high aspect ratio is particularly useful in the context of stealth dicing, a concept in which a series of pulses under high-speed translation of the sample is used to create an entire modification plane [13]. A mechanical, chemical, or thermal stress is then used to cleave the glass along the modified plane. The main benefit of this approach is that cutting is possible at very high speed (typically 10 cm to 1 m per second). Filamentation, i.e., ultrafast pulse self-channeling in glass, was first used for stealth dicing purposes [14,15]. Bessel beams offer in addition the possibility to control very easily the onset of the material damage by simple geometrical considerations since nonlinear effects

**Citation:** Vetter, C.; Giust, R.; Furfaro, L.; Billet, C.; Froehly, L.; Courvoisier, F. High Aspect Ratio Structuring of Glass with Ultrafast Bessel Beams. *Materials* **2021**, *14*, 6749. https://doi.org/10.3390/ma14226749

Academic Editors: Jörg Krüger and Jörn Bonse

Received: 15 October 2021 Accepted: 4 November 2021 Published: 9 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

are reduced. They were used to cut glass and sapphire [16–20]. Since beam length and focusing angle can be independently adjusted, it is also possible to produce very long Bessel beams with the same angle that allows for cutting thick glass. This however requires high pulse energies since the energy is spread over the entire focal line of the beam [21]. With millijoule energy, Bessel beams were used to cleave glass up to 1 cm thick [22,23]. Recently, breaking the symmetry of the Bessel beam was used to improve the precision and ability of guiding fractures [24–27]. Controlling the formation of voids with ultrafast laser pulses is key to improve the stealth dicing of glass.

Here, we show that, depending on the material, voids do not form for the same parameters. We compare three glass types: fused silica, Corning Eagle XG, and Corning Gorilla glass. We investigate the regimes in which various structures are produced using single shot Bessel beams depending on pulse energy, pulse duration, and position in the sample. We observed that these three parameters control the formation of various types of structures inside glass. Future development of applications such as laser cutting should therefore operate in the right parameter window.

#### **2. Materials and Methods**

Our experimental setup (see Figure 1) is composed of a chirped pulse amplified (CPA) Ti-Sapphire laser source, at 800 nm central wavelength and of a reflective axicon (Cailabs) associated to a telescope arrangement to magnify the cone angle [28]. The laser pulses can be stretched from 100 fs to 3 ps using the grating compressor of the amplifier. The experimental characterization of the beam is shown in Figure 2. The Bessel beam has a cone angle of 19.5◦ in air, which corresponds to a central spot diameter of 0.9 μm Full Width at Half Maximum (FWHM). The Bessel beam length is ∼220 μm FWHM. One can notice the high degree of circular symmetry of the Bessel beam. All samples were 700 μm thick.

**Figure 1.** Experimental setup. Ultrafast laser beam is shaped by a reflective axicon (Cailabs). It is then demagnified by a telescopic arrangement in a 2f-2f configuration. Laser-written structures are observed using an incoherent white light source.

**Figure 2.** Characterization of experimental Bessel beam. (**left**) Longitudinal cross-cut of fluence distribution as a function of relative propagation distance in air. Please note aspect ratio of axes. (**right**) Cross-section of fluence distribution at center of beam.

After processing, the laser-modified regions were characterized by scanning the sample with white light, brightfield illumination, in an optical microscope arrangement with a numerical aperture of 0.8 and aberration-correction. Stitching the microscope images provided high resolution images of the whole high aspect ratio structure. Figure 3 shows the result of single shot illumination in fused silica at a pulse energy of 52 μJ, for increasing pulse durations. A periodic background pattern can be observed and simply arises from the illumination inhomogeneity associated to the stitching operation.

**Figure 3.** Microscope images of damages produced in fused silica glass for several pulse durations, for same beam position inside sample and for a constant single shot energy of 52 μJ. Please note image aspect ratios (left insets show 1:1 zooms of panels to their right). Beam propagates from left to right.

#### **3. Results**

#### *3.1. Microstructure Types*

In Figure 3, at a pulse duration of 125 fs, only a line of faint index change is observable. In the cases of 250 and 500 fs, high- aspect ratio voids are observed, they are characterized by black marks. We readily see that the void structure is shorter than the Bessel beam length. The void is present only in the moderate intensity region while in the highest intensity region, the void structure has almost vanished. From 1 ps, the black mark is split around its center into small dots: these are void bubbles (see scanning electron microscopy image in Fig. 4 of reference [29]). This process is called bubbling hereafter. Faint bubbling appears here in the central region of the Bessel beam and is slightly stronger for the case of 2 ps.

Importantly, the morphology of the material transformation can vary along the beam. An example is shown in Figure 4 where we report the results for the interaction of a 3 ps Bessel beam in Gorilla glass for different focusing positions inside the sample. In this example, the material modification undergoes minor variations with the position. We can observe three sections with different morphologies from left to right: (i) a void with faint bubbling in the relatively low intensity region of the Bessel beam (onset), (ii) a wide index modification in the highest intensity segment of the Bessel beam, and (iii) a void. The central panel allows us to understand the ordering of the material modification type with increasing intensity: index modification on a small diameter (central lobe), void formation, structure with more or less randomly distributed bubbles, wide index index modification

on a diameter much larger than the beam central lobe. The wide modification is developed longitudinally step-wise, suggesting a threshold-based mechanism. It erased all structures like void or bubbles. Note that the bubbles appear elliptical because we compressed the horizontal axis scale by a factor 5.

**Figure 4.** Brightfield microscopy images of damages produced in Gorilla glass with 3 ps pulse duration at a constant energy of 52 μJ as in Figure 4. Position of beam in sample was varied from top to bottom. Please note image aspect ratios (left insets show 1:1 zooms of panels to their right).

We further explored, in a parametric study, the evolution of the laser-induced structures as a function of pulse energy, pulse duration, and beam position. We summarized our results in Figure 5 for three positions of the Bessel beam inside the fused silica sample. In Figure 5a, the beam is positioned at the center of the sample such that none of the entrance and exit surfaces are damaged. We classified the modifications in 5 groups, ranging from no damage to strong bubbling. We note that we did not make a separate class for the wide index modification, such as the one shown in Figure 4, because those modifications are conventionally useless for applications. Also, they often occur alongside voids, which we are more interested in. For the drilling and cutting of glass applications, the narrow index modifications and nice voids, shown respectively as gray and green markers, are the most interesting features. When bubbling occurred on a substantial section of the beam, the structure was classified as bubbling.

In Figure 5a, we see that the relatively large range of parameters allow for inducing nice voids in fused silica, similarly as in reference [30]. Higher energies and longer pulse durations tend to create bubbling. When the pulse duration is increased, the energy required to open a void is reduced, but bubbling also arises earlier. The energy range over which nice voids can be opened is reduced at longer pulse durations. In Figure 5b, the beam crosses the entrance side. In this case, the shortest pulse durations provide conditions where voids can be opened on a larger range of parameters. The presence of the interface increases the ability of void opening after energy deposition, as was already noticed in a former reference [6]. In Figure 5c, the beam crosses the exit interface. In this case, the short pulse durations induce voids with very small diameters. A much larger range of energies and pulse durations enable void formation. In almost all the cases tested here, nice voids could be produced. We hypothesize that our results can be understood by the following: the presence of the interface allows the phenomenon of ablation cooling [31] where a fraction of the energy deposited inside the material can be evacuated via the open channel at the surface. This allows for a fast quenching of the material relaxation which increases the capability of "freezing" a void before the compressed melt part surrounding the void can close it. In the case where the beam crosses the entrance surface (Figure 5b), the dynamics of plasma formation is slightly different from the two other cases because a dense plasma can form on the surface which can then shield light preventing void formation in the very first micrometres from the surface. This could prevent the effect of the ablation cooling described above.

**Figure 5.** Morphology of single shot-induced damage in fused silica as a function of pulse energy and pulse duration. This is performed for a Bessel beam position: (**a**) at center of sample, (**b**) crossing the entrance surface, (**c**) crossing exit surface.

#### *3.2. Other Glasses*

In Figure 6, we show the results for Eagle XG glass (Figure 6a) and Gorilla glass (Figure 6b) when the beam is enclosed within the sample, as in Figure 5a. For these two materials, the situation highly differs from the case of fused silica. Overall, it is much more difficult to create voids, particularly for short pulses, where no modification is apparent below a threshold pulse duration of 1 ps.

The operating window for void formation is restricted to a small energy range before bubbling sets in. This window is even smaller in the case of Gorilla glass. For most of the pulse durations in Gorilla glass, above the energy for the bubbling regime, only faint structures are found, whereas this phenomenon occurs in Eagle XG only for 1.5 ps and never in fused silica. It could be attributed to a slower relaxation in Gorilla glass than in that of the other media.

We attribute the differences between the different glasses in terms of void opening to the difference in softening point. Indeed, it exceeds 1500 ◦C for fused silica, while it is only of 971 ◦C for Eagle XG, and 853 ◦C for Gorilla glass. Therefore, assuming a similar profile of laser-deposited energy, the molten region is wider in the case of Eagle XG and Gorilla glass than in the case of fused silica. This would lead to a faster quenching of the relaxation in fused silica, with stronger gradients, in contrast with the other two glasses.

In Eagle XG and Gorilla glass, beam positioning close to the entrance and exit surfaces yield to similar tendencies as in the case of fused silica. While no major differences could be observed with Figure 6 when the beam crosses the entrance surface, we observed a much larger operating window for void opening when the beam crosses the exit surface for both glass types.

**Figure 6.** Morphology of single shot-induced damage in (**a**) Eagle XG and (**b**) Gorilla glass as a function of pulse energy and pulse duration. Beam is centered inside sample.

#### *3.3. Large Heat-Affected Zone for Picosecond Pulses*

Figure 7 shows the result of a 3 ps laser shot in Eagle XG. One can notice that the heat affected zone is very wide, typically above 8 μm in diameter. Importantly, in this regime, the bubbles are not aligned: it seems that they have randomly moved from the optical axis during the relaxation stage. This well supports the interpretation of a large molten layer around the void at large pulse energies and duration.

**Figure 7.** Strong bubbling in Eagle XG glass. Pulse duration is 3 ps, pulse energy 141 μJ. Please note image aspect ratio. Bottom inset shows a 1:1 zoom of marked area in panel above.

#### **4. Conclusions**

We observed that void formation in transparent glasses strongly depends on energy and pulse duration. In fused silica, voids can be formed for a relatively large set of parameters. By increasing pulse energy and pulse duration, the high aspect ratio void channel splits into bubbles with a diameter of approximately 1 μm. This structure can even vanish for higher energies. In Eagle XG and Gorilla glass, only very long pulse durations could yield void formation, but for a very limited energy range. For a pulse duration of 2 ps and 16 μJ pulse energy, a void can be formed in all three materials, with a length of 150 ± 10 μm in Eagle XG, of 170 ± 10 μm in fused silica and the highest length is for Gorilla glass with 230 ± 10 μm. For all materials, we remark that void formation is possible over a larger range of energies when the beam crosses the exit surface.

We interpret the evolution of the morphology of the damages induced by ultrashort Bessel beams in these glasses in the following way. For sufficiently high energies, a void can form inside glass, either after a micro-explosion [32,33] or by cavitation inside a molten zone [34]. We infer that this void can be transiently surrounded by a molten layer, whose thickness depends on the glass type and pulse duration. The longest pulse durations and the highest energies tend to increase the thickness of the molten layer. A low softening point, such as in the case of Eagle XG and Gorilla also increase the molten layer. Thicker molten layers increase the occurrence of bubbling or are even closing the void. We understand the formation of bubbles as a redistribution of the material during the cooling stage by splitting or even fully closing the void transiently formed. We interpret the role of the exit surface as a cooling agent. The evacuation of a fraction of the hot material via the void channel produced at the center of the Bessel beam, allows for quenching the relaxation and maintains the void structure even for parameters where the void would be partially or totally reclosed by the molten material around it. Recent results on double pulse illumination let us anticipate that double pulse could partially reduce the amount of the molten layer, and increase the temperature contrast between the optical axis and the surrounding lobes [35].

These results shed light on the formation of high aspect ratio void channels inside different glasses. We believe they will impact on future technologies for ultrafast laser material processing of glass.

**Author Contributions:** Conceptualization, F.C., R.G., L.F. (Luc Froehly) and C.V.; methodology, C.V., L.F. (Luc Froehly), F.C.; software, C.V.; formal analysis, F.C.; investigation, C.V.; resources, L.F. (Luca Furfaro) and C.B.; data curation, C.V.; writing, F.C.; writing—review and editing, all; visualization, C.V.; supervision, F.C. and R.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Corning Incorporated, Region Franche-Comte council and Agence Nationale de la Recherche (EQUIPEX + SMARTLIGHT platform ANR-21-ESRE-0040), I-SITE BFC project (contract ANR-15-IDEX-0003), and the EIPHI Graduate School ANR-17-EURE-0002.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available upon reasonable request from the corresponding author.

**Acknowledgments:** We gratefully acknowledge the technical support of E. Dordor.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data, or in the writing of the manuscript.

#### **References**

## *Article* **Prediction of Optimum Process Parameters Fabricated by Direct Laser Interference Patterning Based on Central Composite Design**

#### **Mikhael El-Khoury 1,\*, Bogdan Voisiat 1, Tim Kunze 2 and Andrés Fabián Lasagni 1,2**

Received: 30 June 2020; Accepted: 5 September 2020; Published: 15 September 2020

**Abstract:** In this study, we report on the optimization of the direct laser interference patterning process by applying the design of experiments approach. The periodic line-like microstructures of a 8.50 μm spatial period were fabricated by a two-beam interference setup with nanosecond laser pulses, varying laser fluence, pulse overlap, and hatch distance. Central composite design with three factors and five levels was implemented to optimize the required number of experiments. The experimental and numerical results show the impact of various structuring process parameters on surface uniformity. The responses measured are the structure height, height error, and waviness of the pattern. An analysis of the microstructures on the patterned surface was conducted by confocal microscopy and scanning electron microscopy. A 3D-characterization method based on morphological filtering, which allows a holistic view of the surface properties, was applied, and a new qualification scheme for surface microstructures was introduced. Empirical models were also developed and validated for establishing relationships between process parameters and performance criteria. Multi-objective optimization was performed to achieve a minimal value of structure height errors and waviness.

**Keywords:** nanosecond laser; direct laser interference patterning; design of experiments; central composite design; morphological filtering; surface texture homogeneity; micro structuring; bearing steel

#### **1. Introduction**

The functionalization of technical surfaces by producing deterministic topographies today represents an innovation carrier of modern materials engineering. Nature has shown to be the best surface engineer, and the design of these textured surfaces often follows a biomimetic approach motivated by natural designs [1,2]. Therefore, mimicking natural designs helps in understanding the role of surface microstructures and to correlate topographies to macroscopic surface properties. For instance, these well-defined and highly-oriented structures fabricated at the micron and sub-micron scale on surfaces of modern industrial products enable a clear innovation potential to improve products performance significantly [3]. This upgrade concerns a wide range of applications, such as antifouling [4], wetting control, [5,6], tribology [7–10], electrical conductivity improvement [11], cell adhesion [12], as well as surface optical appearance alteration [13,14]. However, the reproduction of these versatile surface structures represents one of the most significant technical challenges today due to their complexity.

In this frame, non-contact manufacturing processes, such as laser-based microprocessing, arose as an extremely viable approach for mimicking natural surfaces because it cannot only provide both the required technological and economic aspects but also ensures the capability to produce high-resolution features [15]. Nowadays, the most prominent laser patterning approaches are direct laser writing, laser-induced periodic surface structures, and direct laser interference patterning [15].

In direct laser writing, the focused laser beam is scanned over the material surface employing pulse-to-pulse strategies, and the resolution is limited by optical diffraction at the focal position. In the case of laser-induced periodic surface structures, repetitive patterns are obtained based on self-organization processes with feature sizes even below the diffraction limit and in the range of the wavelength, or even much smaller [16]. On the other hand, direct laser interference patterning (DLIP) takes advantage of the physical principle of interference of coherent light waves to produce periodic structures on a surface by transferring the pattern shape directly to the material when sufficient laser energy per unit of area (fluence) is applied [17–20]. The interference patterns are formed by splitting a coherent laser beam into multiple beams and hereafter overlapping them on the samples' surface (see Figure 1a,b). Therefore, due to constructive and destructive interference, a specific intensity pattern is obtained, and the size of the periodic structures can be controlled by varying the angle between the interfering beams and the laser wavelength, according to Equation (1):

**Figure 1.** (**a**) Schematic representation of DLIP optical setup; (**b**) schematic view of the interaction zone (with a diameter *Ø)* of two beams intercepting at angle 2θ and forming an interference pattern with a particular period of Λ2; (**c**) distribution of pulses during the DLIP process: Ø, d, and *hd* denote the laser beam diameter (at level 1/e2), pulse overlap, and hatch distance, respectively.

Despite the simple optical arrangement required for DLIP for fabricating periodic structures, it is also important to assure in several cases a uniform distribution of these structures to guarantee, overall, a specific surface functionality.

As reported in several studies, the quality of the final surface morphology strongly depends on laser-material interaction and, thus, on the used laser beam profile and the structuring strategy [21,22]. During the DLIP process, laser-material interaction occurs predominantly at the positions corresponding to the interference maxima, inducing various metallurgical processes, such as melting, ablation, and recrystallization [20]. During nanosecond-pulsed laser processing of metals, structuring mechanisms are mainly based on recoil vapor pressure as well as on Marangoni convection, which have an effect on the overall picture of melt flow [22–24]. In most of the commercial laser systems, the laser beam profile has a Gaussian distribution (TEM00 profile), which necessitates the use of advanced processing strategies, such as consecutively overlapped irradiation of interfering laser pulses in order to form homogeneous structures on large surface areas. For example, Aguilera et al. [22] used this technology and varied different processing parameters in order to produce homogeneous structures on a large surface area. Experiments were carried out in such a way that one process parameter (factor) was applied and varied through different values (levels), then the response was analyzed, whereas the other factors remain unaffected. This procedure is called one variable at a time [25]. For instance, in the mentioned study, discrete values of overlap distance, hatch distance, and fluence were chosen for specific special periods, and it was found that for fluence F = 1.42 J/cm2, 98.5% overlap and hatch distance of 20 μm, for 120 μm beam diameter spot, the structured surface was uniform and, thus, obtaining a homogeneous periodic pattern [22]. However, this approach is time-consuming since several variables have to be to screened independently, resulting in a large number of experiments. Furthermore, the number of experiments increases exponentially with the increase of the number of the factors and their levels.

Additionally, in one variable at a time method, it is impossible to single out the effects of factors interaction that can only be observed when varying multiple factors at the same time. Furthermore, by applying one variable at a time approach, it is possible to miss a process window with optimal settings of factors that will give the desired response. Therefore, these factors should not be examined independently. They should be taken into account simultaneously and must be investigated together since one factor might depend on the level of the other factor. In this case, a more effective approach is a statistical design of the experiment, which aims to decrease the number of experiments and to study the effect of interactions between different factors. For this purpose, many experimental design methods, namely, Plackett and Burman, factorial, Box–Behnken, and central composite design (CCD), have been developed [26–28]. Among these experimental methods, CCD developed by Box and Wilson is a very efficient experimental design method to reduce the number of experiments in the studies with a large number of factors and levels [28]. CCD has more advantages compared to other experimental design methods. For instance, it provides high-quality predictions in studying linear and quadratic interaction effects of factors influencing a system. Whereas interactions, unobserved in Plackett and Burman's experimental methods and Box–Behnken has less coverage than in the case of CCD [28,29]. Therefore, CCD has been widely used in the fields of engineering and science [29–34]. The CCD consists of three main parts and of 2*k* + 2*k* + *m* runs.