Direct correlation of local fluence to single-pulse ultrashort laser ablated morphology

Basic studies on ultrafast laser ablation processes are important for expanding their utility. In particular, understanding the ablated morphology in relation to the incident pulse is critical for micromachining, and an important benchmark for simulations. However, current morphological analyses rely on vast simplifications of experimental conditions, such as a singular fluence value to reduce a unique beam profile, or the maximum crater depth or diameter to describe the ablated morphology. Here, we develop a morphology analysis method in which we take the full two-dimensional information of both the input beam profile and the ablated morphology, and spatially correlate the two without data reduction. We show, using sapphire as a benchmark material, that this serves as a robust way to extract well-studied values and dependencies, such as the ablation threshold, and also as a way to probe the spatial independence of the process. We anticipate that our findings will modernize current study techniques to meet the demand for increased, high-quality data such as that required for artificial intelligence-based analysis. Describing the laser ablation process with reduction-free data is important for furthering its use in modern manufacturing. Here, fluence maps, correlating laser beam intensity and ablated depth at each point in a full two-dimensional space, provide a method to probe ablation morphology in cases of arbitrary beam and crater profiles.

U ltrafast laser ablation has allowed for the manipulation of materials with unprecedented freedom, owed in part to the vast selection of parameters available for tuning the physical processes involved 1,2 . While previously confined to laboratory environments, developments in solid-state laser technologies have enabled robust operation and lower-cost of introduction, both of which have aided its spread into a wider industrial setting. Ultrafast laser ablation has thus received increased attention as a realistic solution for faster, cleaner, and customized fabrication methods required to support the diversified needs of modern manufacturing.
While the potential of ultrafast laser ablation is undoubtedly high, it remains that the vast selection of available parameters to control its process, such as the wavelength 3 , polarization 4 , and pulse duration 5 , also complicates efforts to optimize it. Various factors compound upon each other to make trial-and-error approaches difficult to design and understand. Furthermore, the problem of laser-induced breakdown is a difficult problem in physics itself, and still under vigorous investigation [6][7][8] , with novel phenomena still being reported 9,10 . This situation is also apparent in the fact that the description of seemingly fundamental dependencies, such as understanding the ablated hole depths in relation to the incident laser fluence, is still a difficult task [11][12][13] .
Of particular value in fundamental experimental studies is the study of morphological data, or the shape of the ablated crater. Morphological data is not only of qualitative importance in cataloging the results of laser ablation, but also for quantified data, where results can readily be compared to results of numerical models. However, it remains that the vast number of parameters available makes a comprehensive and standardized study on morphology difficult, where the combinations of possible irradiation schemes complicate most analyses. In contrast, the demand for standardized data acquisition has only increased with the rise of big-data based analysis techniques, such as machinelearning, where not only a large, but a complete and reliable data set is an important prerequisite to accurate analysis 14,15 .
A central problem when standardizing morphology analysis is compensating the effect of different irradiation beam profiles. To standardize the beam profile, an approach that has become ingrained in the community is to represent all data in terms of a single scalar value: the irradiated peak fluence. The beam is often assumed to have an ideal Gaussian shape, with irregularities renormalized for by introducing effective M 2 values. From this and the pulse energy, the peak fluence is analytically calculated. In turn, the morphology is approximated by a similar scalar value, with representative values in the form of the crater depth or diameter used in lieu of a more complicated height profile. Altogether, typical morphological data in the literature consist of such data documenting the dependencies of such representative values in terms of the peak fluence.
While the progress in ablation studies in the past decades is an indication of the effectiveness of this approach, it remains far from complete. Clearly, a vast amount of data is discarded to obtain the representative values. In addition, the approximations are usually far from ideal, and errors introduced from such generalized treatment will propagate as uncertainties in the final obtained values. This may be a factor in the large deviations seen in the physical values reported in the literature, such as the ablation threshold of materials 3,16 . Another more subtle assumption that such analyses impose on the ablation phenomena is locality. By standardizing ablation in terms of the peak fluence, it is often implicitly assumed that the ablation phenomenon at a point is independent of the overall fluence distribution, and only dependent on the local fluence value at that point. This treatment should depend on the spatio-temporal timescales of various diffusion processes caused by non-uniform material excitation in relation to the laser's spatio-temporal characteristics. For example, these approximations would break down in cases with strong lateral diffusion effects, or mechanically induced ablation features 6 . However, such factors are usually ignored, and many key analysis techniques which rely on such assumptions, such as the diameter regression technique to derive the threshold fluence 17 , are applied without particular verification.
Such discussion regarding the robustness of the local fluence dependence of the ablation threshold was a central result of recent studies conducted by Garcia-Lechuga et al., where the authors found that the ablation threshold of materials was reliably determined by the local fluence of the irradiated beam independent of the total pulse energy 18 . They find these results regardless of the highly non-linear nature of the energy absorption process. Such results raise a further question: if the ablation threshold is singularly determined by a local fluence value, is the ablated depth at an arbitrary point in the beam determined likewise. Such locality would provide the prospect of mass parallelization of data acquisition, where each local fluence-depth relationship would yield information equivalent to peak-fluence-to-crater-depth relationships. The detailed and densely acquired fluence-to-depth relationship would serve to facilitate further discussions on the physical mechanisms of laser ablation and nonlinear light-matter interactions. They would also allow for the reconstruction of high-power laser profiles as was done in the X-Ray regime by Chalupský et al. 19 . While the above authors worked to recover the beam profiles of their light source by utilizing well-known Lambert-Beer fluence-to-depth relationships, determining arbitrary fluence-to-depth relationships experimentally should allow for similar procedures. In all cases, locality is an important prerequisite; however, such studies on the full locality have not been conducted.
In this study, to assess these questions and concerns, we focus on the intrinsic two-dimensional nature of imaging sensors, where the spatial arrangement of the imaging sensor should reveal critical information regarding spatial correlation. By imaging both the incident beam profile and the ablated crater profile, and subsequently linking the two, we should be able to accurately document the ablation procedure without approximation, and also reveal the presence or lack thereof of spatial correlations within imaging resolution. To implement this idea, we spatially correlate the beam profile obtained by a spatially in situ CMOS camera directly to an ablated crater profile image by numerical processing. We show the results of our analysis with the laser ablation of sapphire, a key industrial ablation target as well as an often-used material for fundamental studies. We show the first, to our knowledge, pixel-by-pixel observation of morphology depending on the local fluence of the ablating pulse, and demonstrate the advantages of this approach in detecting intra-crater features. Lastly, we explore the general consistency of values obtained by our method with that of the literature. We compare the ablation threshold determined by our method to that measured by diameter regression analysis, as well as with other values reported in the literature. We find general consistency with traditional morphology studies, especially when taking into account typical experimental accuracy 3,16 . This methodology could be of great use in improving the overall accuracy of fundamental morphology studies in the field of laser ablation, where it has often been difficult to obtain systematic and reliable data.

Results
Retrieval of the fluence map. An overview of the process, which we term fluence mapping, is shown in Fig. 1a. The method centers around correlating the local incident laser fluence (F loc ) with the local crater height (h loc ). The spatial information is then removed, and the data is correlated directly as a local fluence to local height relationship into a plot which we call a "fluence map." This mapping allows us to read experimentally important values directly from the plot, such as the ablation threshold fluence (F th ), unlike in traditional studies. Such processing also presents two unique advantages. First, as no Gaussian approximations of the beam profile are introduced, the procedure yields well-defined results regardless of beam profile, and without the introduction of systematic uncertainties. Secondly, the fluence map utilizes the full two-dimensional data obtained from traditional microscopy measurements, in contrast to traditional morphology studies, where only a single representative value, such as ablated depth or diameter, is extracted. The full data made available for analysis allows for the identification of fine dependencies, especially intracrater morphologies which were lost in traditional analyses.
The experimental setup is shown in Fig. 1b. For our ablation target, we used a commercially available, 500 μm thick, c-plane cut sapphire plate (Shinkosha Co. Ltd.). This sapphire plate was mounted onto a holder adjacent to a CMOS camera. Before each ablation, the beam profile of the processing beam was measured spatially in situ by the CMOS camera. To directly image such small profiles, a typical pixel pitch far smaller than the spot size (typically a couple tens of micrometers) is required. For this, we utilized a compact, board-type CMOS camera (Camera Module V2; Raspberry Pi Foundation). The camera, originally a color camera, was modified to have monochromatic response by manually scraping the Bayer color-filter off from the sensor surface. The CMOS camera has a specified pixel pitch of 1.12 μm, which was also confirmed by microscope; with this small pixel pitch, it was possible to measure the beam intensity profile spatially in situ at the beam focal position without any magnification optics. Compared to conventional imaging by magnified imaging of the focus, the in situ camera possesses advantages in integration and spatial calibration. We note that while it may be easier to magnify the profile to greater pixel resolution by magnified imaging, optically resolving the profile to the same spatial resolution would require using an imaging lens with numerical aperture greater than 0.5 for wavelengths near 1 μm. Ablation was performed with pulses of light from a regenerative amplifier or an optical parametric amplifier pumped by this source. Single pulses were used at each spot to avoid complications induced by multiple-pulse effects, such as damage incubation. The specifics of the light source and ablation procedure were otherwise standard, and are further described in the Methods section.
An example of the raw measurement is shown in Fig. 2. An optical microscope image of a crater ablated by a single pulse of 190 fs pulse duration, central wavelength of 1028 nm, and 30 μJ pulse energy is shown in Fig. 2a. A laser scanning microscope image of the crater is shown in Fig. 2b. The results of measurements of the beam profile used to ablate this crater is shown in Fig. 2c. For our method, each data point in the height image must be spatially correlated to an incident beam fluence value. In order to achieve this, (1) the pixel pitch discrepancy between the microscope image and the beam profile must be calibrated, and (2) the spatial position of the two images must be aligned. To achieve (1), as the beam profile was sampled with spatial density much smaller than the spot size, interpolation is possible. Specifically, we perform a smoothing B-spline interpolation of the beam profile image, and up-sample the beam intensity values with the same pixel pitch as the microscope image (approximately 46 nm). We also integrated the total signal intensity of the spline curve to calibrate the pixel count of the camera to the pulse energy. The local fluence values were calculated by dividing the local incident energy value by the pixel area.
Step (2) The intensity of the light was adjusted, then subsequently focused onto a sample holder with an integrated CMOS camera on a XYZ stage. By laterally translating the holder, either the measurement of an in situ beam profile at the focus or processing of the sample could be realized. was achieved by first numerically roughly aligning the beam profile and crater centers. Subsequently, the relative positions were fineadjusted to find a position which minimized the fluence discrepancy along the crater edge. The specific numeric procedures are outlined in the Methods section. The result of the alignment procedure is shown in Fig. 2d. The local fluence values are overlaid as a contour plot onto the local height profile data. In general, features of the crater can be seen to match well with the contour lines of the local fluence with sub-micrometer accuracy.
To better visualize the information retrieved from this procedure, we plot the correlation between the local fluence and local height data directly. The result of this correlation is depicted as a "fluence map", as shown in Fig. 3. Here, 250,000 red points are shown, each corresponding to a single pixel of the original laser microscope image of Fig. 2d. We note that each individual point is not strictly statistically independent, as any uncertainty in the profiling and correlation processes would lead to spatial correlations. The overall plot should be treated as one trend, as discussed in Supplementary Note 5. The binned average values for each fluence are shown in black. The vertical discrepancy in the data corresponds to the uncertainty of around 10 nm in the height measurement by the laser scanning microscope, in agreement with the instrument specification. The horizontal discrepancy mainly arises from uncertainty in the total beam correlation procedure, from profiling, interpolation, and alignment, and not the ablation morphology. The ablation itself was found to have excellent repeatability for 190 fs pulses, see Supplementary Note 3. Correspondingly, the width of the fluence map corresponds well to the spatial derivative of the beam profile over a distance of a couple hundred nanometers. From this single plot, the most important values in ablation studies can be directly determined. Below approximately 3.5 J cm −2 , highlighted by an orange arrow, the crater height is constant at around 0, signifying that no ablation has occurred. On the other hand, above this value, the height decreases dramatically and ablation has occurred. Thus, we can directly deduce the threshold for ablation, F th as the boundary value between these two regions: in this case 3.5 ± 0.25 J cm −2 , where we make a conservative estimate of the precision of a single fluence map based on its width. We note that uncertainties caused by pulse energy and beam profile fluctuations are an order of magnitude smaller than this value. Spatially, the threshold point corresponds to the height-fluence position at the crater edge. This definition of F th is equivalent to what would be used in the diameter regression technique 17 . However, whereas diameter regression postulates a Gaussian pulse shape and requires multiple crater measurements, the current method is applicable to an arbitrary beam profile and can be defined from a single crater. We show the case of an arbitrary beam profile in Supplementary Note 1. Returning to the plot, the peak fluence of the processing laser pulse is shown by a solid, dark-red arrow, and corresponds to the center of the pulse. The depth of the hole, highlighted by a blue, solid arrow, is the value of the lowest red points, here approximately 300 nm. The deepest points correspond to the positions where the local fluence is near the peak fluence, as would be expected. All the above demonstrates how the current technique is compatible with traditional morphology analysis, while eliminating the approximations and data reduction which could be a large source of systematic error.
Pulse energy dependence of the fluence map. The fluence map derived in the previous section reveals the complete local fluenceto-height relationship of a single crater. As a unique application of this method, we experimentally demonstrate a method for determining the local fluence dependence of the ablated depth. We do this by comparing the fluence map created by pulses with differing pulse energies. By changing the pulse energy, we can compare morphologies of craters with equivalent local fluences, but varying local fluence gradients. For example, the local depth ablated by a pulse with certain local fluence near the peak may not be the same as the local depth near the edge of a crater ablated by a higher-energy pulse with the same local fluence. This "locality" of the ablation phenomena has at times been suspected as potential sources of uncertainty in traditional morphology studies, which often postulates locality in ultrafast laser ablation due to its unique timescales 20 .
To investigate dependencies on the total pulse energy, we ablate five craters with varying pulse energies all above the ablation threshold of sapphire. The ablated craters are shown in Fig. 4a. As expected, the craters become correspondingly smaller as the pulse energy is decreased. We next construct a fluence map for each crater from these craters and the corresponding beam profile data. The fluence maps for the five craters are gathered in Fig. 4b. Each color corresponds to a different fluence map generated from a different crater. All the fluence maps reflect the same universal relationship, suggesting that the final ablated depth of laser ablation at these parameters is dictated by the local incident fluence, independent of the global pulse energy. This includes the elevated feature near the crater edge, where in a separate set of experiments, it was even found that fine-tuning the peak fluence to right near the threshold fluence allowed us to create nano-bumps without a crater formation. This result is the first, to our knowledge, to directly observe the local depth dependence on the local fluence within our current spatial resolution. This resolution in turn depends on a combination of the laser scanning microscope's resolution, the beam profiling accuracy, as well as the precision of the alignment procedure. Altogether, the current experiment was found to have an effective resolution of submicrometer order.
Verifying the spatial extent of locality is an important prerequisite for many experimental and theoretical studies. For example, in many theoretical studies 21,22 , the ablated depth dependence on the incident fluence is calculated with uniform irradiation of light over a small region. Calculated results are compared to experimental results with similar peak fluences, but such comparisons are valid only at scales in which the local fluence dominates the physical response. The ability to explore the extent of locality is also a valuable tool in itself to study the typical spatial scales of the ablation phenomena. At the spatial scales and focusing condition explored here, we have shown that sapphire can be treated as a local-fluence dominated material, in agreement with recent experiments regarding the locality of the ablation threshold 18 . Our results extend this verification of the local-fluence dominated dependency in a nonlinear process from single ablation threshold values to ablation depths for each local-fluence by comparing the whole crater morphology. We note that this is not a universal property of this spatial scale and focusing; we have also found other materials where the fluence maps for different pulse energies did not coincide with the same optical setup. Whether a material can be treated as local or non-local at the typical scales of laser processing setups depends on a combination of physical parameters, such as heat and electron diffusion properties. Our scheme allows the experimental verification of such properties. Finally, our current result also yields exciting prospects for future studies. At scales where the ablation may be effectively considered local, each fluence map accurately predicts the ablated depth of individual pulses with differing pulse energies. Thus, previous data sets consisting of tens of ablation points studying depth-fluence relationships can be effectively replaced by a single fluence map, while simultaneously yielding orders of magnitude increase in data density.
Identification of process thresholds. As mentioned in the introduction, the lack of approximations and data-reduction procedures in our analysis allows for precise observation of intracrater features across a range of local fluences. A direct consequence of the increased data density allowed by the fluence map is the ability to visualize finer dependencies in the ablated morphology. While the main ablation threshold, highlighted by the filled triangle in Fig. 4b is apparent, upon further inspection, a critical point in the depth to fluence relationship can be seen at around 6 J cm −2 , highlighted by an open triangle on the same plot.
The corresponding spatial profiles of craters ablated by pulses with energies above and below this second threshold is shown in Fig. 4c and d. The approximate spatial position where the local fluence crosses each threshold is also highlighted by their corresponding triangles. Crossing of the higher threshold was observed to be concomitant with an increase in surface roughness, as seen in the cross-section profiles. This increase in roughness can be observed as a dark shadow area near the center of the craters in the 30 and 26.7 μJ craters of Fig. 4a. Such a threshold in sapphire is not documented in the literature to our knowledge, but may suggest a transition to a violent ablation phase, such as boiling or phase explosion, although they are usually reported in the context of metals or semiconductors 6,23 . We believe this threshold to be different from the gentle/strong ablation threshold documented in the literature for sapphire 24 . In ultrafast laser processing of this material, it is known that above a certain energy threshold, the ablation rate and surface roughness both increase dramatically (from so-called "gentle" to "strong" ablation). It has been experimentally observed that in the femtosecond region, gentle ablation corresponds to a so-called Coulomb explosion ablation phenomenon, while strong ablation to a thermal process 25 . This Coulomb explosion form of ablation for single pulse illumination was measured to have an ablation depth of a few nanometers 11 , one or two orders of magnitude smaller than the strong ablation depth and consistent with theoretical models 26 . Here, we do not observe such shallow ablation features. This may be a result of a combination of focusing and measurement conditions, where we use tighter focusing than 11 , and our laser microscope vertical resolution was limited. The results hold for 800 nm light at 170 fs (discussed later and in Supplementary Note 4), where we rule out potential effects of a different illumination condition from previous works 11, 24 . We note that this is not an effect of the fluence mapping analysis, as such shallow ablation features are not present in the preprocessed crater profile, for example in Fig. 4c and d. Overall, we conclude that our second threshold is not the gentle/strong ablation threshold, but some other process threshold within the strong ablation regime. While it is currently difficult to quantitatively identify the threshold of such processes with traditional morphology studies due to its chaotic morphological attributes, the onset of different processes can be clearly recognized in the fluence map. This is only possible because the fluence map utilizes the full extent of the raw microscope data, allowing for the efficient spatial averaging of changes in tendencies. From this, we conclude that the current method is effective in elucidating new physics as well, owing to its more complete data utilization.

Discussion
Lastly, we comment on the accuracy of the derived F th . We first compare the accuracy compared to the diameter regression technique. With the craters of the varying pulse energies, we perform a diameter regression analysis of the threshold by measuring each crater diameter and plotting the results squared as a function of the pulse energy 17 , as shown in Fig. 5a. From the model fit, the ablation threshold pulse energy (the intercept of the fit) and beam radius (the slope of the fit) could be determined. From these values, it is possible to calculate the threshold fluence, which we find to be 5.00 ± 0.03 J cm −2 , significantly higher than the value derived by fluence mapping. The uncertainty here is calculated from the standard deviation of the best-fit parameters. The discrepancy can mainly be attributed to systematic errors in the beam diameter value retrieved from the diameter regression fit. Here, the radius of the diameter regression was determined to be 14.61 ± 0.05 μm and 13.54 ± 0.05 μm in the respective x-and y-direction, while the beam diameter obtained by simply fitting the beam profile yielded a minor and major radius of 14.43 ± 0.02 μm and 15.81 ± 0.02 μm, where again, the errors are simply the standard deviations of the least-squares fit. This smaller diameter-regression value is not physically impossible, but is close to the theoretical limit of focusing of the current laser pulse, and highly improbable considering our setup and the M 2 value of the laser, around 1.1. More importantly, re-calculating the fluence with the beam-profiler radius yields a threshold fluence of 4.07 J cm −2 , which, while still higher, is more consistent with the fluence mapping results. Thus, it is more natural to conclude that the diameter regression analysis failed to properly determine the laser radius from the experiment, resulting in a large systematic overestimation of the threshold fluence. Indeed, it has been found in previous studies that the diameter regression technique is particularly vulnerable to such errors, as the entirety of the beam profile is analyzed from the beam profile near the peak of the pulse 16 .
Additionally, in order to compare our results with those in the literature, we conduct the same experiment with 800 nm light from an optical parametric amplifier (OPA) with a pulse duration of approximately 170 fs. We determine a threshold value by ablating a total of 12 craters at two different axial positions within the Rayleigh range of the processing beam, and extract the threshold value from each fluence map as the peak position of the inflex feature near the rim, consistent with methodology of most diameter regression measurements. A representative fluence map is shown in Supplementary Note 4. We obtain a value of 3.73 ± 0.03 J cm −2 , with the uncertainty here corresponding to one standard deviation of the averaged individual results. This result, along with representative results from the literature, is shown in Fig. 5b 11,24,[27][28][29][30][31][32][33][34][35] . As stated previously, the current technique defines the threshold in an equivalent manner as the diameter regression technique, and hence, should yield the same result. The result matches reported values of diameter-regression determined values at similar conditions (blue circles) particularly well, and is reasonable when compared with thresholds defined in other schemes (red squares), as expected.
In conclusion, we demonstrated a method for the analysis of morphology data in laser ablation studies. Such analysis can greatly increase the robustness of quantitative measurements of ablation in regards to beam profile deviations, as well as greatly improve data utilization. The major advantages of this method include the ability to inspect the local fluence dependence of the ablation process in regards to change in pulse energy, as well as the ability to differentiate various thresholds by analyzing fine trends in the morphology.
We should also formally address the potential limitations, although many are not unique to this method. The method requires that the pulse be spatially similar (only changed in the fluence), meaning that situations where chromatic aberrations or other general spatiotemporal distortions would see inappropriate analysis. This method would also be inapplicable for situations where nonlinear propagation effects come into play for the main ablating pulse. It would also be limited in scope for highly stochastic ablation, such as that seen in pulses longer than a few picoseconds for dielectrics. Lastly, it is also ill-equipped to address multi-pulse ablation morphology, such as the effects of a slanted crater edge, or subwavelength periodic structures. These limitations are, however, not often significant for fundamental studies of ablation depths and thresholds, which is the main target of our study.
Several future prospects are opened with this study. Practically, the current method is important for quantitative comparison of experimental results between different laser systems, where the spatial deviations of the beam profile may play a non-negligible role. Even in the case of the same laser system, strong deviation may arise due to change in parameters, such as the wavelength of an OPA. Here, measurement accuracy is limited by the beam quality deviation between different wavelength settings, which was also discussed as a possible limiting factor in previous studies on the wavelength dependence of the ablation threshold 36 . In terms of basic physics, this method allows us to directly analyze the local fluence dependence of the whole ablated morphology, which may be important in materials where ultrafast carrier diffusion occurs faster than the electron-lattice thermalization times. Combined with time-resolved techniques 37 , this could open avenues in our understanding of the processes occurring at the sub-picosecond timescale for irregular irradiation. Lastly, in the case of local phenomena, the ability for this method to efficiently create unique ablation depth dependence plots should be important for data-intensive studies as well. The ability to obtain thousands of data points from a single crater opens the opportunity to methodologically categorize the behavior of various materials under laser excitation, thus creating an ablation "fingerprint" database. For example, the depth dependence on the local fluence should reflect on the effective penetration depth of the laser energy and should elucidate the nonlinear lightabsorption processes currently under discussion. While these functional forms may be difficult to analyze with current theories, deep-learning based analysis may prove successful in overcoming this obstacle 14 .

Methods
Light source. For our light source, we utilized a 190 fs, 1028 nm central wavelength regenerative amplifier system (PHAROS-1.5SP; Light Conversion) or an optical parametric amplifier (OPA) pumped by this laser (ORPHEUS; Light Conversion). The power of the pulse was modulated by a half-waveplate and polarizer set. The laser was focused onto the sample surface by a single plano-convex lens, with a focal distance 100 mm for the fundamental, and 75 mm for the OPA. The final focusing lens were chosen to keep the final focused e −2 radius of the beam constant at approximately 15 μm, where minor changes of a few micrometers of the radius did not affect fluence map results.
Single pulse ablation procedure. Before each ablation experiment, the average laser power was monitored by a power meter (XLP12-3S-VP-D0; Gentec-EO). Subsequently, neutral-density (ND) filters were inserted into the system, and the sample holder was moved via motorized XYZ stage so that the laser was incident upon the camera surface, and the laser beam profile was recorded. The camera had 10-bit resolution. The vertical displacement in regards to the sample surface and CMOS sensor plane was compensated for by measuring the 3D height information of the entire holder with a 3D scanner (VR-3200; Keyence Corp.), and moving the sample holder a corresponding distance. Finally, the ND filters were removed and the holder was translated so that the laser would be incident upon the sample surface. A single pulse, selected by a built-in pulse-picker of the regenerative amplifier, was then illuminated upon the sample surface.
We note that this procedure is highly reliant on the stability of the laser system. In cases where fluctuation of the pulse energy or beam profile are noticeable, the pulse characteristics measured beforehand would not be a good indication of the actual processing conditions. In such cases, it may be favorable to split the beam to create a separate optical path for diagnostics with identical optics and alignment. In the current study, the fluctuations for the regenerative amplifier were found to be negligible for our purposes, with 0.3% rms fluctuation in the power and less than 0.5% fluctuation in the beam radius (see Supplementary Note 2). Thus, we opt for the simplified in situ setup to avoid potential systematic errors induced by improper calibration of the diagnostic path.
Image processing procedures. To extract the center of the beam profile, a 2d Gaussian fit was conducted on the laser scanning microscope (VK-X260; Keyence Corp.) profile, yielding the beam center position. In order to up-sample this beam profile image at the pixel pitch of the crater image, a smoothing B-spline interpolation of 5-th order was conducted. This spline function was spatially integrated and equilibrated with the pulse energy, allowing us to rescale values in terms of J cm −2 .
Next, the center of the ablated crater was determined by first detecting the approximate crater edge position by a standard Canny edge detection algorithm. In order to reduce detection of random features created by localized dust and debris, a median filter of 7 × 7 pixels is conducted to the image beforehand. The median filter, as opposed to other filters, allows the removal of sharp outliers, while preserving the position of edges. From the crater edge, the region of the image corresponding to the crater was identified, and the pixel corresponding to the midpoint of the height and width was marked as the crater center.
To fine-tune the two relative center positions of the beam and crater, we worked under the assumption that points with similar morphological features would most likely identify with points with similar fluence values. To this end, we focus on aligning the beam relative to the crater lip feature, which was already identified in the previous steps. Four points along the crater edge on the x and y intercepts were selected as reference points. The square difference of the fluence values at the four points was determined as F sd ¼ ∑ i≠j À F i À F j Á 2 i; j ¼ 1; 2; 3; 4 f g ð Þ , where F i are the fluence values at the reference points. The F i were calculated for all values of the reference points spatially translated in a ± 3 μm range with 100 nm step, and the offset which minimized F sd was used as the final offset between the beam and crate center positions. While this method does not account for potential differences in the rotational degree of freedom, the mechanical alignment by the sample holder in the experimental setup proved sufficient for the current results.