Picosecond pulsed underwater laser ablation of silicon and stainless steel: Comparing crater analysis methods and analysing dependence of crater characteristics on water layer thickness

Liquid layer thickness dependence of 515 nm, 7 picosecond pulsed laser ablation of stainless steel 304 and silicon is analyzed. Ablated crater volume and diameter are compared to ablated craters in ambient air by means of a novel, objective numerical procedure. While silicon ablation under a water layer is found to be more efficient in terms of removed material volume per pulse than ablation in ambient air, an opposite trend is found for stainless steel 304. For both materials, the ablation efficiently drops when the liquid layer thickness is decreased to 1 milimeter. A probable reason for the ablation efficiency drop is persistent bubble formation.


Introduction
Ultra short pulsed under liquid laser ablation is a field of science which is actively studied in the context of eye surgery [1], nano particle production [2] and surface texturing [3]. In terms of material removal, ablation under a water layer is known to create deeper craters for femtosecond pulsed laser ablation of brass [4] relative to ablation in ambient air. A similar trend has been observed for nanosecond pulsed ablation of silicon [5,6] and aluminum [7]. The cause for this ablation efficiency increase in the case of nanosecond pulsed laser micromachining of silicon was attributed to an increase of plasma density created during the ablation process and the generation of a shockwave due to cavitation bubble formation [8]. The timescales at which these phenomena take place were nicely categorized by Dell'Aglio et al. [9] and play a role in under water ablation aside from photon-carrier, carrier-carrier and/or carrier-phonon interaction typical for in air ablation of metals [10,11] and semi-conductors [12]. Specifically in the context of nanoparticle generation, work has been performed on the analysis of the cavitation bubble formed in under liquid pulsed laser ablation. X-ray illumination was used to analyze bubble content [13,14] and stroboscopic shadowgraphy imaging is often employed to study bubble dynamics [15,16]. Shockwave and bubble dynamics were also studied as a function of liquid layer height over the sample [17]. The effect of liquid layer height on post-ablation crater depth has been identified for nanosecond pulsed ablation of silicon [18] as well as for aluminum [19] and Inconel 718 [20] drilling. Results in these works show a strong relation between crater depth and liquid layer thickness. In particular, the depth of ablated craters in silicon shows 0.1 mm sensitivity to layer thickness changes [18]. Ablation experiments with a layer accuracy of this liquid level has not been performed for other materials, which implies that there is room for further research on under liquid ablation using a set-up that facilitates a layer thickness with sub-milimeter precision.
Comparison of in air and under liquid experimental results would require an analysis method which allows for the unambiguous comparison of the properties (volume, diameter) of craters. Typically, crater circumferences which are drawn 'manually' in microscopy images, are used to determine a radius and circle centre that 'best' fits the crater. For craters which are not perfectly circular, this approach is far from trivial and highly subjective. The depth profile of under water ablated craters are known to be non-gaussian shaped for ultrashort pulsed ablation on silicon under certain parameter conditions [21], which hampers the effectiveness of crater diagnostics. Alternatively, crater profiles have been measured and integrated relative to a reference plane to yield a crater volume [22]. In this approach, the choice of reference plane placement is presumably based on unablated sample material surface roughness, though the exact definition of this plane is typically not defined. This would mean that local inhomogeneities in the surface roughness are not taken into account in the volume deterimnation which complicates the determination of the crater volume accurately.
The goal of our presented work is therefore twofold: First, to obtain crater data on under liquid ablation with a layer thickness that is defined with sub-milimeter accuracy and second to analyse these craters using an objective volume determination method that allows for compensation for local surface roughness inhomogeneities.

Laser setup
The laser setup is outlined in Fig. 1. A 7 ps pulsed Yb:YAG laser source (TruMicro5050 of Trumpf, Germany) with a fundamental wavelength of 1030 nm was frequency doubled to 515 nm using a second harmonic generator (SHG). The beam quality of this source equals M 2 ≤ 1.3. Hence, its fluence profile is nearly Gaussian. The pulse frequency was set to 1 kHz to avoid the laser-beam interaction with laser induced cavitation bubbles. A combination of a λ/2-plate and a polarizing beam splitter was employed to attenuate the laser beam. A galvoscanner (IntelliScan 14 of Scanlab, Germany) in combination with an Ftheta telecentric lens (F-theta-ronar lens by Linos AG, Germany) with a focal length of 100 mm was used to focus the laserbeam into an optically transparent and watertight box, see Fig. 2. The focus of the laser beam was measured outside of the watertight box to be 23 μm using a beam profiler (MicroSpotMonitor by Primes, Germany). To align the galvoscanner with respect to the box, a linear stage was used (ATS150 of Aerotech, USA). The optically transparent walls of the box consist of four 4 mm thick 50 by 50 mm square silica glass plates and a base plate of aluminum. The glass plates were coated with a visible light antireflective coating. The box was mounted to an xy-stage (two ALS20020 stages of Aerotech, USA) to allow accurate positioning of the box with respect to the incident laser beam. Two steel gauge blocks with a thickness defined with an accuracy better than 1 μm were mounted to the inside of the wall facing the incident laser beam using magnets placed on the outside of the silica glass, see Fig. 2.

Samples
Two sample materials were used. A silicon waver (thickness 1050 μm) with crystal orientation <100> was cut into samples of approximately 20 by 10 mm. Additionally, a stainless steel 304 plate was cut into samples of 20 by 20 mm, embedded into an epoxy and subsequently polished to obtain a surface rougness of Ra 0.16 μm. Prior to the experiment, samples were mounted inside the optically transparent box by pressing the sample into the gauge blocks, after which both sample and gauge block were 'locked' inside the box by means of magnets, see Fig. 2. The demineralized water was poured in the box to fully submerge the sample. During the experiments, power measurements were performed directly in front of the optically transparent box using a power meter (PM100A of Thorlabs, Germany) and a power sensor (S130VC of Thorlabs, Germany).

Theory
The height profiles of ablated craters were measured by means of confocal laser scanning microscopy (CLSM, VK-9710 of Keyence, Japan) using a 1024 times 768 pixel camera. The confocal microscope has a 1-σ repeatability error of 0.02 μm.
The height profile of under water ablated craters are known to be non-gaussian shaped for ultrashort ablation on silicon under certain parameter conditions [21], which hampers the effectiveness of crater diagnostics by means of Liu's method, also known as the D 2 -method [23]. In this section, a novel numerical method is introduced to determine both the volume and the equivalent diameter of craters.

Ablation conditions
Craters were processed using 50 consecutive laser pulses at varying levels of pulse energy on silicon and stainless steel. The ambients considered are demineralized water and air. For water, experiments with a liquid layer thicknesses of 1, 2, 3, 4 and 5 mm were performed. The effective pulse energies at the surface of the sample were determined by compensating for reflection losses [24]. At 515 nm, the refractive index of silicon is n si = 4.211 +0.0417i [25] and the refractive index of stainless steel 304 is n ss = 2.000 +3.471i which was determined by ellipsometery. The indexi of air, silica and water equals n air = 1.000 [26], n silica = 1.462 [27], n water = 1.330 [28]. Resulting transmission values are presented in Table 1 and a more elaborate procedure for the computation of the stainless steel values in Table 1 may be found in existing literature [29]. Effective pulse energies for all ambients and all samples were varied between 1 and 10 μJ. Focus conditions under liquid were determined by offsetting the focus distance in air by a distance H as a function of liquid layer thickness h l and the refractive index of water according to [30], (1)

Crater analysis method
The purpose of this section is to obtain an objective measure of the amount of removed material from the sample due to laser processing.   The region Ω covered by the confocal image is divided into two sections: a band region Ω b on the outer edge of the image, and a middle region Ω m in the centre, see Fig. 3. Ω is covered by a Cartesian array of N rectangular quadrilateral cells each with center points (x i , y i ) and cell area ΔA, where i is the sequence number of the cell and x and y are Cartesian coordinates. Altitudes of the cells are stored in an array z i with i = 1, 2, …, N. Three corrections of the altitude data z i are required: a. a correction to remove noise generated during the confocal imaging process, b. a correction to obtain altitudes relative to the unprocessed surface, c. a correction to avoid false removal contributions due to surface roughness of the unprocessed surface.
These three corrections are subsequently discussed in the following section. To remove noise generated by sharp gradients on the surface of samples, the data is smoothened as follows: where I i is the index set of the four neighboring cells of cell i and n s is the number of smoothing operations. Next, the height of the unprocessed sample surface is linearly approximated as in which the coefficients a, b 1 and b 2 are to be determined. The RMS error of the approximation over the outer region Ω b is defined as where 〈.〉 b denotes the average over Ω b , with N b the number of cells belonging to Ω b and I b the set of index values refering to grid points in Ω b , The coefficients in (3) are determined by minimizing ∊ w.r.t. a, b 1 and b 2 . With the approximation z 0 i of the unprocessed surface known, the relative altitude data z i are defined as Then points considered to be part of the crater are defined as points satisfying z i <z ☆ , where z ☆ is a threshold. This threshold is required, because without it the number of 'improper' crater points on an unprocessed sample would be equal to half the total number of points, which is evidently not useful. The number of improper crater points can of course be made zero by choosing z ☆ sufficiently small, but this would induce an unacceptable underestimation of the "real" number of crater points. The strategy chosen in this work is to derive an approximate expression for the relative error in the number of crater points N c as a function of z ☆ , and to choose an acceptable value of this error from which the corresponding value of z ☆ follows. The relative error in N c is estimated by first estimating N c itself as being approximately equal to the number of elements in the laser spot, as the diameter of the laser spot is known. The error ΔN c in N c , is equal to the number of improper crater points in the region outside the laser spot, α is defined as the fraction of improper crater points within any given set of points of unprocessed surface. The relative error β in N c is now simply The fraction α can be derived from the band region by counting the number of points in the band satisfying z i <z ☆ and dividing that number by the total number of points in the band. Once a value for β is chosen and the corresponding value of z ☆ is determined iteratively by matching the relative number of improper crater points in the band with the value of α from (10), crater area, effective crater diameter and the crater volume are computed as

Parameter validation
Suitable values of the three parameters of the numerical approach in Section 3. The value of β is chosen by observing the sensitivity of z ☆ with respect to β for the band region of samples. Fig. 5 shows the values of z ☆ , averaged over all samples, as a function of β, computed through (10) where α(z ☆ ) is determined from the band region Ω b with a fixed value of N b /N = 0.1. In addition, z ☆ is also shown seperately as determined from the band of an unprocessed sample. Remarkably, the two graphs show that the band area of ablated samples changes during the ablation process. Hence, thresholds are determined based on the band region of unablated silicon and stainless steel. An error margin of 1% is maintained for all samples considered, so β is chosen 0.01 which amounts to a threshold z ☆ of − 0.4980 μm for silicon and − 0.29637 μm for stainless steel. These thresholds are applied for all craters analysed.
To determine the number of smoothing iterations n s , N c is plotted as a function of the number of iterations in Fig. 6 for craters shot using an effective pusle energy of 3, 6 and 9 μJ. The values stabilize after about 7 iterations and therefore n s = 7 is maintained for all samples considered.

Results and discussion
In this section, the volume V c and squared diameter d 2 c as functions of pulse energy are presented and a comparison between the conventional 'D 2 -method' and the square of the numerically obtained crater diameter, d 2 c , is discussed. Finally, a selection of crater morphologies is also presented by means of a set of combined light and confocal microscopy images. Cross-sections of a selection of craters are also provided. For the complete numerical analysis, about 1100 craters were analyzed. For every pulse energy level at every ambient, 5 craters were analyzed.

Crater volume and area
The crater volume V c and diameter d c data for silicon and stainless steel are shown in Figs. 7 and 8 respectively.
In literature, it is typically assumed that the volume of a crater scales as ζln 2 E p /E th with scaling factor ζ, ablation energy threshold E th and pulse energy E p [22]. In air ablation on zinc [31], as well as various other metals [32][33][34] were analyzed using this method.
Fitting volume data in the presented work yielded unacceptably large errors for the fit coefficients, presumably because the pulse energy range used in our work is inconsistent with the chosen range in other studies. One paper for example, considered an effective pulse energy of about 1-2.5 μJ on silicon [35], other work takes into account a much larger range with less data points per energy interval [22]. For this reason, the fit was omitted in our data. Fig. 7 shows the volume data for silicon and stainless steel in different ambients. Numbers were added to the graph to indicated trend breaks. For silicon ablated under a 1 and 2 mm water layer, crater volume strongly increases for the first 2.5 μJ up to point . This increase is much steeper than the increase observed for silicon ablated in ambient air, for which volume increases up to point ②, at 5 μJ. Beyond points and ②, volume increase as a function of pulse energy seems to converge to similar values for silicon ablated under a 2 mm water layer and ambient air results. In contrast, for a 1 mm water layer, additional trend breaks at points ③ and ④ occur causing the graph as a whole to show 'oscilatory' behaviour when pulse energy is increased beyond 2.5 μJ.   -Crater volume as a function of pulse energy can be subdivided into two regimes, regardless of ambient. The pulse energy that seperates the two regimes is different for ambient air and water and also varies with ablated material. This pulse energy is 5 μJ for silicon ablated in ambient air and 2.5 μJ for silicon ablated under a water layer. For stainless steel the regime transition occurs at approximately 4.5 μJ for ambient air and 2.5 μJ for under water ablation.
-1 mm water layer results deviate significantly from results obtained under thicker water layers: on silicon, several trend breaks are observed and for stainless steel the trend break between the two regimes occurs at 3.5 rather than 2.5 μJ.
In all subsequent graphs, the two general regimes will be indicated. Interestingly, on silicon a 5 mm water layer ambient yields the largest ablated volume, far outperforming ablation in air. On stainless steel however, ablation in ambient air yields much larger craters than ablation performed under a water layer. Of all the different water layers used, a 3 mm water layer seems to yield the largest craters at a pulse energy level of 2.5 μJ. Both for silicon and stainless steel, reducing the water layer below a 2 mm water layer proves detrimental to crater size. During the ablation process, bubbles were observed in the liquid which were 'stuck' between the optically transparent box wall and the sample, possibly causing the reduction in crater volume.
Calculated diameter values are shown in Fig. 8. For in air ablated craters, stainless steel d 2 c values in Fig. 8 show two different slopes in regime I and II whereas for silicon, the second regime is identified by constant d 2 c values. This constant region will be adressed further in Section 4.3. The results in Fig. 8  Combining Figs. 7 and 8, it seems crater volume increase is accompanied by an increase in crater diameter on silicon in ambient air for regime I, whereas in regime II this increase is due to an increase in crater depth. Interestingly, this trend seems reversed for silicon ablated under a water layer. For stainless steel craters created in ambient air, volume and diameter increase occur simultaneously over both regimes. For stainless steel, the relation between volume and cross-sectional area increase under water is fairly straight forward: the initial volume increase in regime I is coupled with a slight increase in diameter and in regime II both volume and cross sectional area are mostly constant.

Numerical versus visual crater diameter determination
Conventional characterisation of a crater involves measuring the diameter d p of the perceived (by the human eye) crater. The square of this diameter plotted as a function of pulse energy forms the basis for Liu's method [23]. Results of this method applied on the data set are shown in Fig. 9. Light and confocal microscopy images for different pulse energies for air and specific water layer thicknesses were combined and are shown in Fig. 10. Note that in this figure the light microscopy image is provided in grayscale whereas regions included in the area computation and volume computation of (9) are colored, based on the colorbar indicated in Fig. 10.    To analyse the differences between the d 2 c and d 2 p graphs, crosssections of craters belonging to the different regimes are displayed in Fig. 11, in which the relative altitude of the cross-sectional areas of several craters are shown. The red dashed and green dotted line lengths are equal to d p and d c for each crater and their y-value corresponds to the relative altitude z at which they were measured or computed. The relative altitude is scaled using the threshold z ☆ . Fig. 11 shows a signifcant increase of d p from regime I to regime II in ambient air for silicon.
Conversely, d c stays nearly constant, explaining the constant d 2 c for regime II for in air ablated silicon in Fig. 8. For ablation under water, the difference between d p and d c is significant because of small portrusions on the outer edge of the crater cross-section, which have hardly any depth and are therefore not taken into account for d c while they do form part of d p . These portrusions are part of the aforementioned cauliflower structures. The relative difference between d p and d c remains largely constant from regime I to regime II for silicon ablated under a water layer, which explains why even though d c and d p values differ, the general trends in Figs. 8 and 9 are similar.
For stainless steel ablated in ambient air, the general trends of Figs. 8 and 9 for regime I and II are similar. Under water obtained d c and d p results vary a lot, mainly due to the cauliflower portrusions on the outer edge of the crater increasing the d p values relative to their d c counter parts. As cauliflower structures become more apparent for higher pulse energies, the discrepancy between d c and d p values increases from regime I to regime II. Interestingly, crater depth and width actually decreases over this pulse energy range as well, which explains the decrease in d c value as regime I changes into regime II. The decrease in crater size coupled with the severity of the cauliflower structure formation seems to suggest that a significant part of the laser deposited energy does not end up at the crater centre but rather is diverted to the perimeter of the crater where it is responsible for the creation of the cauliflower portrusions.

Crater morphology
Silicon crater evolution in air and water has been thoroughly covered in literature [21,36]. The results in Fig. 10 confirm that craters in ambient air are initially wider and shallower than their under water created counterparts. For the 8 μJ results obtained under a water layer, the crater seems to split into two excentric circles rather than a single circle shown for lower pulse energies. This behaviour was also reported in earlier work [21] and it was suggested to be caused by laser-induced effects in the water, though no specific mechanism was mentioned. The stainless steel results in ambient air show an affected zone surrounding the crater for all pulse energies and a central crater which grows as pulse energy increases. The cauliflower like structures for stainless steel ablated under a water layer cover a larger area the shallower the liquid layer as may be observed in Fig. 10. Additionally, the structures become more prominently visible for higher pulse energy levels. The speckle-like structures observed for the under liquid ablated craters were formed during the process of exposing the sample to ambient water, but do not seem to be caused by the ablation process itself as a post-process analysis of unablated material showed similar structures. A confocal data analysis shows that the structures are pits and are thus not flat regions on the sample surface.
A possible culprit for the cauliflower like crater structure under water for both sililcon and stainless steel could be bubble formation. Cavitation bubbles induced during the ablation process tend to have lifetimes of a few microseconds [9], which is much shorter than the interpulse time in our experiments and is therefore not likely to influence the process. However, persistent microbubbles are known to occur during the laser ablation process [37,38] which linger relatively long near the laser-material interaction zone and show lifetimes into the milisecond range. Such bubbles would surely cause scattering of the incident laser beam, increasing the ablated area beyond the region one would typically expect.

Threshold values
A relation between the square of the diameter of a crater and the pulse energy is typically formulated as [23] with D 2 the square of the crater diameter, ω 0 the 1/e 2 laser beam spot radius, E p the effective pulse energy and E th the energy ablation threshold. From this, a fluence ablation threshold F th can be determined as in which T is the material and ambient transmission value given in   Table 2.
For the ablation threshold of silicon, 4 mm was selected as the designated reference water layer thickness to compare the literature to.
No literature reference for under water ablation of stainless steel could be found. The literature reference for under water ablated silicon [39] refers to a 10 mm water layer thickness experiment performed used a femtosecond pulsed laser, whereas our results were obtained using a picosecond pulsed laser at smaller liquid layer thicknesses. Although much information is available on the ablation of silicon under a water layer [21,36,41], no suitable reference for the threshold of under water Fig. 11. Cross-section, obtained by confocal microscopy measurements, of craters created on silicon (left) and stainless steel (right) in air and water. For each regime identified a crater is shown. For regime I craters created using 2 μJ are shown, for regime II craters shot using 8 μJ are shown. d c and d p values are shown at their measured or computed relative altitude z. The y-axis is scaled using the altitude threshold z ☆ . The centre of gravity of each crosssection was set as the origin in each image.
ablation on silicon using a picosecond pulsed source could be found. For stainless steel ablation a reference [40] is added for the 50 pulse ablation of stainless steel using a 10 picosecond pulsed 1030 nm laser source. For under liquid ablation of stainless steel, no reference was available.
Computed threshold values as well as their literature counter parts are displayed for all ambients in Fig. 12.
For silicon ablated in air, the threshold obtained using d 2 p values corresponds well to values found in literature, whereas for stainless steel this is not the case, likely due to the difference in used wavelength. For the reference found on the ablation of silicon under water, an effective pulse energy range of about 2.5 to 22 μJ was used. Additionally, an 800 nm, 250 femtosecond laser source was used. These factors likely account for the large discrepency between existing literature and presented thresholds. Error bars for d 2 c and d 2 p obtained thresholds seem similar in size in Fig. 12, however the error in the fluence thresholds obtained by measuring the diameter of the craters is not necessarily a measure of the uncertainty with which the threshold could be determined. Rather, it is a measure of one's abillity to draw circles consistently over ablated regions. Particularly for higher pulse energies and for under water ablated craters, Fig. 10 shows crater regions may possess a form that deviates significantly from a circular one. From this perspective, the error shown in Fig. 12 for the measured thresholds is more ambiguous than the threshold determined numerically.

Conclusions
Water layer thickness dependence of silicon and stainless steel ablation was investigated and compared to ablation in ambient air. Volumes and areas of the craters were analyzed using the conventional D 2 analysis method as well as a newly created numerical objective approach. Two distinct regimes were found. Ablation data accquired using the new method agreed reasonably well with data obtained via the conventional approach, although significant deviation from literature reported values occured. Cauliflower-like structures hampered conventional D 2 analysis for under water craters, particularly for higher pulse energies. For silicon, a 5 mm water layer was found to yield optimal results in water, whereas for stainless steel this was found to be 3 mm specifically at 2.5μJ. Ablation in air yields higher volume craters for stainless steel relative to under water ablation while for silicon an opposite trend is observed. Crater areas for higher pulse energies were found to be very dissimilar to a Gaussian profile, a possible culprit for this observed phenomenom is persistent microbubbles.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  [35] and under a water layer [39]. Note that the last reference only studied under water ablation for a single layer height. For comparison purposes, this value has been extended to all water layer thicknesses in the graph. For stainless steel, a reference for ablation in ambient air is also provided [40]. A reference for stainless steel ablation under a water layer was not found.