Simple and robust method for determination of laser fluence thresholds for material modifications: an extension of Liu’s approach to imperfect beams

The so-called D-squared or Liu’s method is an extensively applied approach to determine the irradiation fluence thresholds for laser-induced damage or modification of materials. However, one of the assumptions behind the method is the use of an ideal spatial Gaussian beam that can lead in practice to significant errors depending on beam imperfections. In this work, we rigorously calculate the bias corrections required when applying the same method to Airy-disk like profiles. Those profiles are readily produced from any beam by insertion of an aperture in the optical path. Thus, the correction method gives a robust solution for exact threshold determination without any added technical complications as for instance advanced control or metrology of the beam. Illustrated by two case-studies, the approach holds potential to solve the strong discrepancies existing between the laser-induced damage thresholds reported in the literature. It provides also an appropriate tool for new studies with the most extreme laser radiations.


Introduction
The determination of the local laser fluence (or intensity) is critical when working with ultrashort laser pulses, since this marks the onset for risks in front of laser exposure (skin or corneal damage) [1][2][3] or damage of optical materials [4][5][6][7][8][9] .However, the fluence is not a directly accessible quantity, as for its determination it is necessary to characterize both the integrated pulse energy and its spatial distribution.Although the energy is easily measurable by calibrated photodiodes, thermal sensors or pyroelectrics sensors, the determination of the spatial beam distribution can be more complex depending on the considered radiation and the precision needed.
The most precise method for beam profiling is obviously by direct imaging.The use of 2D-sensors or cameras have been applied for this purpose.This methodology is of easy application for collimated beams with size comparable to the camera array.Moreover, this methodology is also applied for focused beams by introducing an optical system (normally a microscope objective and a tube lens) for re-imaging and magnifying the small laser spot onto the camera detector [10][11][12][13] .However, there are two limitations associated to this approach: (i) It is experimentally difficult to design and implement a perfect imaging system that will not distort the observation.(ii) The spectral response range of cameras is limited with our current technology.For high-resolution silicon technologies, this basically limits the applicability to the visible or near-infrared domain of the spectrum.Out of this range, different technologies exist (including InGaAs for the extended near-infrared) but they are not routinely available in laboratories and often costly despite more limited performances (pixel size and dynamic range).
Alternatively, there are strategies to retrieve characteristics of the beam shape without the need of imaging it.The first set of techniques are the ones using an obstacle to partially block the beam (knife-edge 14,15 or a wire 16 ), measuring the transmitted energy and retrieving by algorithmic calculations the beam waist.A second set of techniques are the so-called impact-based strategies.Those are of particular interest since the strict threshold response of material to ultrashort pulse irradiations makes that induced modifications can be taken as direct imprints of the laser profile 13 .This strategy analyses the shape of modifications (ablation, changes of reflectivity, etc) produced at different pulse energies, associating the modification borders to a local fluence corresponding to the fluence threshold for modification.As an example, this strategy has been effectively applied for characterising the spatial distribution of ultrashort X-rays pulses 17,18 , that is in a particularly challenging spectral range for direct imaging technologies.
Among the impact-based strategies, there is a technique that stands out for its simplicity: the Liu's method 19 .Assuming a threshold modification response, this method allows to the user to retrieve the waist of a Gaussian beam by a linear fit when representing the diameter (square) of the modification versus the pulse energy (in logarithmic scale).This method published in 1982 is still extensively used in the ultrafast laser community, an aspect that can be illustrated by more than 600 citations since 2015 (for a total of around 1050 citations, data extracted from "ISI Web of Science").The success of this technique is not only because the beam waist becomes easily accessible, but also because the modification fluence threshold of any material can be obtained by using it.This second potentiality was not commented on in the original paper, but rapidly had an impact according to evidence of its exploitation only a few years later 20,21 .
Even though Liu's method is very frequently used, it is not always applied correctly since it provides only accurate fluence analysis if the irradiation beam is perfectly Gaussian (the

Amendments from Version 1
This version contains minor corrections addressing the comments from the reviewers together with some additional clarifications.We have made in particular two changes on the terminology.First, we have replaced the term "Integrated Spatial Distribution (ISD)" by the most commonly used formulation "Effective Beam Area (EBA)".A second source of confusion in previous version was a concept appearing under two different terms ("energy fitting limits" and "maximum energy considered").In this new version, this is now everywhere referred as "maximum considered energy".This modification appears also on the new x-axis label of Figures 3(a) and 3(b).Looking to other figure modifications, we have introduced colour changes on Figure 2, following now the same colour ordering as in Figure 3. Section 5D has been expanded by including additional details on the precision aspect (error determination).Some recommendations to the good applicability of this methodology are also added in this section.Finally, we have introduced two new tables (Table 3 and Table � �) where the correction factors plotted on Figure 3 can be more easily read by users of this methodology.
Any further responses from the reviewers can be found at the end of the article assumption of the method).This condition cannot always be fulfilled as beam imperfections from laser systems or practical optical set-ups (e.g.aberrations) often occur.Therefore, the applicability of Liu's method with unperfect beams (asymmetries, pedestals, etc) could be one of the reasons on the large dispersion of the ablation fluence threshold values reported on test materials (as fused silica) for apparently similar conditions 5 .This raises an important issue that can be summarized with the following circular reasoning.How can we trust the fluence values obtained by Liu's method if we do not know if the beam is perfectly Gaussian?How can we certify having a Gaussian beam if what we wanted with Liu's method was to avoid beam imaging?
In this article, we present an extension to the Liu's method to make it valid for beams that clearly deviate from the Gaussian approximation.This extension relies on exact correction factors to account for Liu's method results when irradiating with a beam with Airy disk-like shape.This close-by Gaussian spot is a characteristic diffraction profile that can be directly generated by introducing a circular aperture in the beam path before the focusing element, a strategy commonly used in optical set-ups for laser material processing [22][23][24][25][26] .
This article is structured as follows.In the second section, a complete explanation of the original Liu's method is presented.In the third section, the calculations of the correction factors to be applied on the Liu's method when irradiating with a perfect Airy-disk are shown.In the fourth section, we repeat the calculation for more realistic cases, using truncated beams generated by different aperture sizes.Finally, in the fifth section, we make an experimental demonstration of the validity of our calculations.This demonstration is carried out for two different beams at two different wavelengths (1030 nm and 1550 nm).
We consider that the presented method can be of general application, helping for reliable comparisons and thus in solving some persistent discrepancies on fluence threshold determination that are due to methodology issues.Additionally, the advent of new laser sources in different parts of the spectrum, and in particular in the infrared domain that hold promises for new scientific and industrial applications, supports the timeliness of this report to set a general criterion for accurate determination of the fluence and that is not dependent to a metrology technology.

Liu's method
Liu's method (or D-square method) refers to a simple experimental approach that allows determination of the fluence ablation threshold by measurements of the sizes of induced modifications at different irradiation energies 19 , without the need for imaging the beam profile.This methodology assumes a Gaussian laser beam profile, being expressed mathematically as, Eq.1 where F(r) is the local fluence at a given radial position, r, F 0 the peak fluence value and w 0 the radial Gaussian beam waist at 1/e 2 the peak value.Liu's method, under the hypothesis of a deterministic ablation behaviour, defines the ablation fluence threshold, F th , as the local fluence corresponding to the border of the crater, exhibiting a radius equal to R. Therefore, the following expression is obtained, Eq.2 expression that can be transformed into a linear relationship by taking its logarithm, leading to: Eq.3 However, the experimental parameter that is usually measured is not the peak fluence but the integrated pulse energy, E.
The relationship between those parameters is linear, related with effective beam area (EBA).This is a general relationship (independent of the beam shape) and is obtained by the 2D-integration of the fluence distribution 27 , 0 ( , ) .
Calculating this integral for the Gaussian function (Equation 1) in polar coordinates (dS = r dr dθ), one obtains 2 0 / 2, EBA w π = leading to the well-known relationship: Eq.5 Accordingly, by defining the energy ablation threshold, E th as the minimum pulse energy at which ablation is produced, the fluence ablation threshold is calculated as, Eq.6 By substituting Equation 5 and Equation 6 in Equation 3, we establish the relationship between the observables R (or the ablated area, A = πR 2 ) and E, This relationship is the key point of the Liu's method.The latter proposes to represent on the x-axis the "ln(E)" values and on the y-axis the "R 2 " values.Therefore, the Gaussian beam waist (w 0 ) is retrieved through the slope of a linear regression.
Additionally, even if it was not mentioned on the original paper of Liu, throughout the x-axis intercept the energy threshold value (E th ) can be retrieved.Then, the laser fluence threshold value, F th , can be simply calculated by applying Equation 6with the two retrieved values.
This method is also called as the D-square-method according to the representation of graphics with diameter values (D 2 ) on the y-axis.Then, the Equation 7 turns into:

Extension of Liu's method for an Airy Disk: mathematical description and correction factors
An Airy disk is the diffraction pattern obtained at the focal position resulting from a uniformly illuminated circular aperture.It is described mathematically in polar coordinates as, with . Eq.8 J 1 (r′) is the Bessel function of the first kind of order one and w Airy the radial waist at 1/e 2 the peak value introduced with a change of variable.As an equivalence from the previous section for a Gaussian beam, the fluence ablation threshold is described as, Eq.9 with R the radius of the crater.
An obvious difference with the Airy spot is in the absence of a linear relationship as the one facilitating the analyses for rigorous Gaussian beams.However, we would like to highlight here the similarities and quantify the differences between these two radially symmetric profiles.As an example, an Airy disk with a waist of w Airy = 10 μm is shown in Figure 1 (a).
Together with this function is also represented a Gaussian beam of identical integrated pulse energy.The most important differences are visible near the pedestal of the distributions with a progressively vanishing Airy function oscillating around zero.Unfortunately, an Airy disk beam has not a simple analytical expression establishing the relationship between fluence and energy, as there is for Gaussians beams (Equation 5).This imposes to use directly the integration formula expressed on Equation 4 for the Airy disk function (Equation 8), which allows the numerical calculation of the corresponding EBA.For our represented case (w Airy = 10 μm), the effective beam area value obtained is EBA = 188.2μm 2 .
In order to evaluate the deviations of the Gaussian-based Liu's method when applied for an Airy disk, we present in Figure 1(b) the relationship between R 2 and ln(F 0 /F th ) for our Airy case.The values are retrieved from the data shown in Figure 1 (a), doing the pertinent calculations and swapping the axes.It can be observed with Figure 1(b) a deviation from a linear behavior at relatively high excitation levels, similar to the one experimentally observed by Bonse et al. 24 .Despite this trend, equivalent Gaussian functions (the closest) can be obtained by linear regressions of this graph, that is literally applying the Liu's method.Due to the non-linear character of this curve, the result depends on the range of considered energies.As examples, in Figure 1(b) two different linear regressions are presented, one with a maximum considered energy of E = 2 ⋅ E th and another of E = 10 ⋅ E th .
In Table 1 we provide the obtained parameters (w 0,Liu and E th,Liu ) applying the same procedure under different maximum considered energies.The obtained equivalent Gaussian beam waists (w 0,Liu ), allows for the calculation of their corresponding  After comparing those values with the EBA obtained after the Airy disk integration a correction factor, η EBA , is obtained.Knowing the relationship between EBA and the peak fluence (Equation 4) a fluence correction factor is calculated as η F = 1/η EBA .Therefore, the peak fluence of an Airy disk for a given measured energy (E) can be expressed as, This expression is directly equivalent to the Equation 5 but introducing the fluence correction factor, η F , a factor that depends on the range of energy considered when the Liu's method is applied.
Additionally, as observed with the inset of Figure 1(b), the linear fit applied to the Airy function can lead to errors on the reading of the real energy threshold by the x-intercept.For all discussed cases, this error is numerically evaluated introducing the correction factor of the ablation energy threshold, th E η .Overall, the exact ablation fluence threshold when applying the Liu's method for an Airy beam is obtained as, Complementing some correction factors already given in Table 1, more cases are calculated and a corresponding abacus is presented together with measurements in the following section (see in particular Figure 3).

Extension of Liu's method for Gaussian beams truncated using circular apertures
As previously mentioned, an Airy disk is the diffraction pattern obtained at the focal position of a lens of limited aperture when irradiated by a uniform plane wave.In practice, this corresponds to a perfectly top-hat beam facing the lens, or the use of a circular aperture much smaller than the size of a nearly-Gaussian beam.Therefore, in this section we show the calculations and corrections factors to account in the Liu's method for more realistic cases using truncated Gaussian beams with circular apertures at focusing lens position.
For obtaining the diffraction pattern for different truncation conditions we rely on the PSFLab software 28 for calculations based on rigorous vectorial theory.We calculate the point spread function and in particular the fluence radial profile at the focus for different focusing and illumination conditions.
The calculations account for a Gaussian profile filling parameter, β G , which is defined as the ratio between the radius of the circular aperture (a) and the collimated Gaussian beam radius (w): β G = a/w.This parameter describe also the power transfer (or aperture transmission) P T , as 29, Interestingly, this power transfer is a parameter that can be easily determined in an experiment by simply measuring the laser power with and without the chosen aperture.
In Figure 2, we show the obtained profiles for different power transfer values, ranging from 90 % to 25 %.In this particular case, the parameters used for the calculations were: λ=1030 nm, focal lens of 50-mm, an incident collimated Gaussian beam with diameter of 5-mm (2•w) and aperture diameter values (2•a) adapted to fit the chosen aperture power transfer (Equation 12).
All profiles in Figure 2 exhibit a ring structure, even if hardly visible with a linear scale in intensity.The relative intensities of the rings are represented in Table 2 and compared with those of a perfect Airy disk (β G = 0).Additionally, the energy distribution on the central region and the rings is also represented on Table 2, after a numerical calculation of the EBA for the full function (Equation 4) and the contribution of each part of the ring structure.The table provides already useful information on the characteristics of the ring functions resulting when truncating a Gaussian beam with circular apertures.
Even if performed for a specific wavelength and focusing conditions, it is important to highlight the generality of these calculations that only depend on the filling parameter β G .
With the obtained beam profiles, a representation similar to the one plotted in Figure 1 (b) is performed for each function (not shown).Applying Liu's method representation and the same linear regression analyses described on the previous section, correction factors are obtained.Again, the correction factor depends on the range of considered energies for the linear regression.The results of these calculations are represented in Figure 3 for various energy ranges, with energy maxima up to five times the energy threshold (5 ⋅ E th ).
We observe from Figure 3(a) that the introduction of a fluence corrective factor is needed even if the beam spot is the result of a moderate truncation (e.g.P T = 90%).This is important since one may intuitively expect negligible consequence of a   11) when irradiating with truncated beam by a circular aperture.In both cases, the horizontal axis corresponds to the maximum considered energy for the linear regression of the Liu's method.Straight lines joining the points serve for view guiding.

Table 2. Maximum intensity and energy distribution analysis of the different parts of the resulting beam profiles from truncated Gaussian beams (circular aperture) shown in
Figure 1 (a) and Figure 2.

Second ring Rest
.2•10 -3 87.9 % 5.5 % 2.1 % 4.5 % moderate truncation for the validity of a Gaussian approximation.By helping in symmetrizing an imperfect incoming Gaussian beam one may have considered the truncation of beam as an improvement to approach the Gaussian profile assumption needed to apply Liu's method.While this filtering may be beneficial for the validity of Liu's method representation, our analysis indicate one should not ignore that the resulting peak fluence and fluence threshold values will remain biased unless the appropriate correction factors, as those calculated in Figure 3, are applied in Equation 10 and Equation 11.
For the practical case of Gaussian beams truncated by using circular apertures, another interesting conclusion from Figure 3 (b) is that the correction factor of the energy threshold, th E η is less significant than the correction associated to the fluence distribution, η F .In order to make data more accessible for potential users of this methodology, we present the correction factors, η F and th E η in Table 3 and Table 4.

Experimental validation and practical relevance A. Experimental configuration
In this section, we show two experiments validating the use of the proposed corrected extension of the Liu's method for a beam truncated with a circular aperture.We show the superior relevance and robustness for threshold determinations for beams deviating from a perfect Gaussian profile.This is because the truncation tends to create, independently of incoming beam, a more controlled Airy-like pattern for which rigorous correction factors can be derived for the application of Liu's method (originally proposed for Gaussian beam only).
In the first case, this methodology was applied for a beam at 1030-nm wavelength directly generated by a femtosecond laser amplifier (Pharos, Light Conversion).On the second case, it was applied for a beam at 1550-nm wavelength generated through non-linear processes on an optical parametric amplifier (OPA, Orpheus-HP, Light Conversion).The pulse durations at both laser wavelengths were characterized by an autocorrelator (TiPA, Light Conversion), being of 170 fs at 1030 nm and 190 fs at 1550 nm.
The irradiation set-up was composed of a variable circular aperture (SM1D12C, Thorlabs) positioned as close as possible to an, aspheric lens of f=50 mm (117-2550, Eskma) and a XYZ-motorized sample holder.Single-shot irradiations were controlled by a pulse-picker integrated in the laser system.Pulse energies were externally adjusted by a set of broadband metallic filters.The irradiated sample of reference was a sapphire window of 1 mm thickness and c-cut orientation.The choice of this target material was motivated by the very neat craters produced with ultrashort pulses without apparent thermally affected zones 13 , being consequently considered as a reference dielectric for impact-based beam characterization methods.Additionally, a fused silica sample (UV-fused silica) was also used on the experiments, given the fact that it is probably the most studied dielectric material.The numerous damage fluence thresholds values reported in the literature are important for comparisons and validation of the methodology of this article.In all cases, our damage criterion is ablation, which is determined by measuring the profiles of irradiated areas by a confocal microscope (Leica DCM3D, 460 nm illumination, 150× objective lens).
Examples of the images obtained under this microscope can be found in our previous publications 13,30 .
Additionally, an imaging system composed of a microscope objective (Mitutoyo 100X NA-0.5, or Mitutoyo 50X of NA-0.42), a tube lens and infrared camera (Raptor OWL 640), mounted in a micrometric XYZ stage, was used to obtain the beam image at the focal position.In order not to introduce any distortion on the beam, images are recorded at low energies.
For each image a background subtraction procedure is applied on the basis of another reference image captured after blocking the laser beam.The 16-bit intensity image is normalized after dividing by the peak intensity value in the measured distribution.To account for the magnification of the imaging system, the image is re-scaled after imaging a resolution test target.With the image resulting from those operations the EBA in μm 2 is obtained by numerical 2D-integration.
Data produced for validation are available from Zenodo as underlying data 31 .

B. Ablation test experiment at 1030-nm
The results of the modifications induced in sapphire and fused silica with pulses of 170-fs at 1030-nm wavelength are shown on Figure 4 (a).The representation follows the Liu's method (see Equation 7) with an x-axis for the energy in logarithmic scale to perform linear regressions on the data.In particular, two experiments were performed on the sapphire sample.The first one is for a beam directly focused on the target surface and second one is for the same configuration but placing a circular aperture before the lens adjusted for a power transfer of 75%.
The crater measurements without the pinhole are represented with solid squares.A nearly perfect linear behavior is observed up to energies of about 40 µJ, which corresponds to excitation levels of ∼6 times the ablation energy threshold.This deviation could be associated to air ionization affecting the propagation of intense beams 32 .However, the corresponding irradiation intensity in our experiment (≈ 10 14 W/cm -2 ) does not directly support this hypothesis.The deviation is more surely explained through the analysis of the beam image at the focal position shown in Figure 4 (b).On this beam image and its horizontal cross-section, we observe a pedestal surrounding an almost perfect Gaussian profile.The influence of this pedestal becomes visible on the modifications only when irradiating well-above the ablation fluence threshold, as observed on Figure 4 (a).However, more importantly, this leads to a significantly biased fluence threshold determined by the original Liu's method due to some energy distributed outside the assumed Gaussian beam profile, as explained in the following paragraph.
When applying the Liu's method taking all the values below 40 µJ, we obtain E th,Liu = 7.1 µJ and ω Liu = 9.7 µm.The corresponding full-width at half maximum is FWHM Liu = 11.4 µm, being in excellent agreement with the value obtained by imaging as FWHM Image = 11.6 µm.Applying Equation 6, this leads to the fluence threshold determination F th,Liu = 4.8 J/cm 2 .
Alternatively, an integration of the fluence distribution can be obtained with the image shown in Figure 4 (b), leading to a value of EBA = 165 μm 2 .By using the ablation energy threshold obtained with the Liu's method and applying the relationship described in Equation 4, we obtain for a rigorous determination of the threshold: F th,Image = 4.3 J/cm 2 , showing a 11.6% difference when comparing with the value obtained with the Liu's method.This difference directly comes from the perfect Gaussian beam assumption used in Liu's method.Its validity is not strictly fulfilled in the considered case and we dare say in most of experiments as real beams always exhibit more or less imperfections.
According to the craters produced in sapphire when placing a circular aperture setting a power transfer of P T = 75%, the linear regression corresponding to Liu's method, also shown in Figure 4(a), leads to E th,Liu = 21.1 µJ and ω Liu = 18.1 μm.The maximum considered energy for the fit equals to 50.9 µJ, corresponding to 2.41 ⋅ E th .Therefore, under this experimental conditions and looking at Figure 3 (or Table 3 and Table 4), we extract the following correction factors: η F = 1.02 and th E η = 1.01.Applying those correction factors and the regression parameters obtained by Liu's method into Equation 11, a fluence threshold value of F th = 4.2 J/cm 2 was obtained, which is very close (2.3% difference) to the value obtained by complete numerical integration of the beam profile (see above).
Finally, this methodology is applied also for fused silica.The conventional Liu's method gives E th,Liu = 17.9 µJ and ω Liu = 17.6 µm.According to the maximum considered energy for linear regression, 47.4 µJ (2.64 ⋅ E th ), we extract the following correction factors: η F = 1.01 and th E η = 1.01 and we finally obtain F th = 3.8 J/cm 2 .Validation of this result by comparing with values in the literature is not simple, since the reported values exhibit large differences, as summarized by Gallais et al. 5 when comparing the fluence threshold values obtained by different authors in fused silica under irradiations at 800 nm or 1053 nm.Even though comparisons are difficult, we observe that our value is in good agreement with the value of 3.9 J/cm 2 obtained by Winkler et al. 33 .performed with pulses with characteristics close to the ones used in our experiment (200 fs pulse duration at 800 nm).
We can conclude, after demonstrating the obtention of similar values by two different means (beam imaging and our presented methodology), that controlled truncation of the beam leads to an improved reliability and confidence on the fluence threshold metrology.This experimental approach can be useful to solve some of the discrepancies observed in the literature due to works relying generally on Liu's method without detailed analyses on the real beam profile.

C. Ablation test experiment at 1550 nm
To further illustrate the benefit from the proposed simple method, we show in Figure 5 (a) the results of an ablation experiment by irradiation with single 190-fs pulses at 1550-nm wavelength.We compare the measurements obtained by using directly the beam delivered by the OPA (without any aperture) to those obtained with a truncated beam using an aperture set again for a power transfer of 75%.Before discussing on the results, it is worthy to look at Figure 5 (b-c), where the corresponding beam profiles at the focal position are represented.In Figure 5 (b), the beam image produced after focusing the OPA beam (without aperture) shows a notorious deviation from a perfect Gaussian beam.This beam shape, with the presence of a large pedestal where an important part of the energy is present, clearly makes Liu's method not applicable.This point will be confirmed just after.
Moreover, this imperfect shape at the focal position also suggests that the beam before the lens is also imperfect and not Gaussian.Accordingly, the Fourier transformed function resulting from the application of a circular aperture is not necessarily corresponding to the rigorous analysis made in section 4, unless the aperture is sufficiently closed.Figure 5 (c) shows the beam profile resulting after placing a circular aperture with P T = 75%.A central spot and a surrounding ring are observed, where the maximum value of the ring (after a circular analysis with ImageJ software) corresponds to 7.4 ⋅ 10 -3 .This is very close to the theoretical value of 7.7•10 -3 in Table 2 and confirms the applicability of the analysis made in section 4 suggesting that the cut pedestal was constituting the main difference with a Gaussian beam.In view of this case, it is interesting to retain P T = 75% as an appropriate level of truncation but even stronger truncation may have been needed depending on the beam quality.
We now return to the analysis of the results shown on the Figure 5 (a).When applying the Liu's method to the values obtained without an aperture we obtain an ablation threshold value for sapphire of 7.8 J/cm 2 .Due to the presence of a considerable energy in the pedestal of the beam profile and ignored by the method, this can be considered as an incorrect value.Applying the same methodology for the data obtained by placing the circular aperture, E th,Liu = 33.9µJ and ω Liu = 22.1 µm is obtained.According to the range of considered energies for linear regression, up to 55.4 µJ (1.63 ⋅ E th ) the correction factors to introduce in Equation 11 are : η F = 1.05 and th E η = 1.00.Finally we obtain F th = 4.6 J/cm 2 .This illustrates the enormous fluence determination error (>60 %) when directly applying the Liu's method for our OPA imperfect beam.Performing the same procedure for fused silica, we obtained: The values calculated by applying Equation 4 after numerical integration of the beam image (Figure 5(c)), giving a EBA = 680 µm 2 , are F th = 5.0 J/cm 2 for sapphire and F th = 4.2 J/cm 2 for fused silica, overestimating the values obtained by the proposed method of this article.This observation highlights the main limitation of fluence threshold determination accuracy when using beam imaging, since its precision lies on a perfect spatial calibration and on rigorous integration of the entire beam details.The latter consideration corresponds to ideal imaging conditions that are not accessible experimentally.In our particular case, the slightly overestimated threshold values can be associated to an underestimation of the EBA, since the intensity of secondary diffractive rings do not overpass the noise level of the IR-camera used for imaging.Considering this aspect, our obtained threshold values becomes very consistent and supports the appropriateness of the proposed analysis based on Liu's method by introducing correction factors.

D. Relevance, precision and recommendations
The proposed method of general application shows the ability to improve the reliability of the fluence threshold determination without the need for rigorous beam profile analyses.In this context, we anticipate its usability for fluence threshold determination at non-conventional wavelengths, in which beam imaging can become complex or, for the most extreme cases, impossible with the available sensor technologies.In particular and after the demonstration at 1550-nm, we consider that this methodology will be really useful for the fluence threshold determination on the short-wave (SWIR) and mid-wave infrared (MWIR) ranges.Those spectral ranges are of increasing interest for the laser material processing community due to the development of new sources 34,35 .However, there are challenges remaining for high quality beams with these new developments.This must make very appropriate the application of the method presented in this article.
Additionally, this method is not only useful for fluence threshold determination but also for peak fluence determination of an unknown beam.For doing that, after obtaining the fluence threshold with a truncated beam by a circular aperture, the peak fluence of an unknown beam (for example the original untruncated beam) is calculated following this expression, Eq.13 Where F 0 , E and E th are respectively the peak fluence, the pulse energy and the energy threshold of modification of the unknown beam, and F th is the fluence threshold value obtained after irradiating with a truncated beam using a circular aperture.E th would be the only parameter to be obtained by irradiating with the unknown beam.
Although the method presented has been shown to be efficient for accessing a robust metrology in threshold determination, some aspects should be commented for a proper usability and accuracy quantification: 1.The presented correction factors are only valid if the sample is situated on the focal position or close to it (Rayleigh length).Otherwise, the beam spot will not show an Airy-like beam shape and the described analysis is not applicable.
2. It remains crucial for accurate fluence threshold determination to have sufficient data at near threshold conditions (<2⋅ E th ), as on the classical Liu's method.Otherwise, it constitutes a new source of error, both on the E th and w 0 determination, that is not considered in our report.
3. The energy of irradiation in the experiments should remain under the condition for air ionization and nonlinear propagation effects (e.g.defocusing) affecting the spatial characteristics of the delivered beam 12 .This is of especial relevance when irradiating materials exhibiting relatively high fluence thresholds (e.g.dielectrics).It is mainly for these practical considerations that we have limited our calculations of the correction factors to levels below 5 E th (Figure 3).

4.
It is estimated, if the other considerations are fulfilled, that the obtained fluence value should not be expressed with a better accuracy of ±5 %.This percentage is a conservative estimation accounting for the following aspects: possible uncertainties on selecting the most appropriate correction factors (especially when extrapolations of theoretical calculations of Figure 3 are needed to describe the practical case) or particular material responses interfering on the measured ablated area (as for example the melted rims formed in some glasses 9,13 ).On this second aspect, imprecisions due to subjective ablated area determination are not considered here, since our metrology is based on user-independent analysis as referred on a previous publication 30 .
Finally, as recommendations for the level of truncation to be applied, we have observed that P T = 75% is an appropriate level when having a pseudo-Gaussian beam with a pedestal (section 5B and section 5C).In the case of more irregular beams or totally unknown beams, we propose to perform experiments under two levels of truncation, one with moderate truncation (P T ≥75% ) and another with strong truncation (P T ≤ 60%).
Obtaining similar fluence values under the two experiments would indicate that a moderate truncation level is sufficient to transform the beam into an Airy-like beam.Otherwise, the most accurate fluence threshold value is obtained under the strongest truncation, since strong spatial filtering reduces the influence of irregularities of the input beam.

Conclusions
In the present work, we have explored the validity limits of the Liu's method 19 which is widely applied for its usefulness for rapid assessment of material modification thresholds and achievable resolutions.The method has two requirements: (i) a strict threshold response of the material without surrounding affected zones and (ii) a perfectly Gaussian beam profile impinging on the target.While we have investigated the first of these requirements in recent works 13,30 , we have concentrated here on the more technical question of the importance of the beam profile.An important conclusion is that a modest deviation to the ideal Gaussian can lead to significant errors in threshold determinations.By calculating and measuring the errors associated with more or less diffracted imperfect beams, we show that errors exceeding 20% can be easily caused by beam imperfections that are undiagnosed if only the produced craters are analysed.This is because the upper part of most laser spot distributions can be advantageously compared to a Gaussian function, and so exhibit a linear dependence of the area to the logarithm of the energy in a more or less extended range of energies above threshold.The only strictly rigorous solution to this problem is to measure the beam for accurate determination of the profile and systematic numerical comparisons with the spatial characteristics of the produced modifications.However, such accurate measurement is not always possible (depending on radiations and associated measurement technologies) and its non-necessity represents actually the direct benefit and interest of the Liu's method.
For this reason, another important contribution with this report is with the introduction of a simple extension of Liu's method to solve this limitation.The quantitative determination of the needed correction from a truncated Gaussian beam (depending on the data considered) suggests that the introduction of a partially-closed aperture can always be used to produce a better-defined profile on target.While the associated Airy-disk pattern is in principle inappropriate to use the Liu's method, we have shown it leads to a superior reliability in threshold determination provided that the correction factors derived in this report are applied for compensation.
The reported findings give a comprehensive vision of the measurement limitations that can explain some of the strong discrepancies existing in the literature reporting damage thresholds.A general problem with norms and standards existing on this question 36 is on the a priori knowledge of all experimental conditions that is not all always accessible.The general applicability of error compensation on apertured beams makes it particularly interesting because it improves the measurement reliability without any change of the currently widely used experimental methodology.

Methods
Numerical calculation of the beams for different truncation conditions are obtained with PSFlab software 28 (version 3.5).
The calculation of Liu's method parameters (section 3, section 4 and section 5), the obtention of the Airy disk function (section 3) and the calculation of the integrals enabling to obtain the EBA (section 3 and section 4) are performed by programming respectively the Equation 7, Equation 8and Equation 4. In this manuscript those calculations were performed by using MATLAB (R2020a under a licence of Universidad Autónoma de Madrid).Other programming languages (e.g.Python or C++) would be also appropriate for those purposes.
Analysis of experimental beam images is performed by using ImageJ software (version 1.53a).
The complete experimental methodology is detailed on section 5A.

Data availability
Underlying data Zenodo: Raw data for manuscript "Simple and robust method for determination of laser fluence thresholds for material modifications: an extension of Liu's approach to imperfect beams".http://doi.org/10.5281/zenodo.4421003 31is project contains the following underlying data: -PSFlab_raw_1030nm-f50mm-T25.txt (PSFlab raw data of P T = 25% profile represented in Figure 2).
-ForLiu_log-and-R2-25T.dat (Treated data from PSFlab_ raw_1030nm-f50mm-T25.txt for directly applying the Liu's method.The first column represents the irradiation energy in log scale and the second column represents the squared crater radius).
-ForLiu_log-and-R2-40T.dat (Treated data from PSFlab_ raw_1030nm-f50mm-T40.txt for directly applying the Liu's method.The first column represents the irradiation energy in log scale and the second column represents the squared crater radius).
-ForLiu_log-and-R2-60T.dat (Treated data from PSFlab_ raw_1030nm-f50mm-T60.txt for directly applying the Liu's method.The first column represents the irradiation energy in log scale and the second column represents the squared crater radius).
-ForLiu_log-and-R2-75T.dat (Treated data from PSFlab_raw_1030nm-f50mm-T75.txt for directly applying the Liu's method.The first column represents the irradiation energy in log scale and the second column represents the squared crater radius).
-info-images.txt(Information for the spatial calibration of beam images).
threshold for material modification which constitutes a critical parameter in laser-based processing.The authors of the manuscript follow a rigorous approach to extend the method used for Gaussian beams to cases where other laser beam profiles are employed (Airy beams).They derive a solution based on analytical expressions for the radial distribution of fluence while theoretical results are tested against experimental observations.Both the predictions, results and approach constitutes a novel technique that can benefit researchers in the field of laser processing.Thus, I would recommend this work for indexing.My only concern that needs to be addressed a bit in more detail (to be elaborated further) is the concluding remark (paragraph just before the 'Conclusions') of the limitations of the technique and more specifically the range of laser energies in which the method works.Also whether the method is material and also laser wavelength-dependent.

Is the rationale for developing the new method (or application) clearly explained? Yes
Is the description of the method technically sound?Yes

Are sufficient details provided to allow replication of the method development and its use by others? Yes
If any results are presented, are all the source data underlying the results available to ensure full reproducibility?Partly Are the conclusions about the method and its performance adequately supported by the findings presented in the article?Yes Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Laser Matter Interaction, theoretical modelling of ultrafast dynamics and surface pattern formation.

I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
a list, including a new point commenting on the remaining uncertainties in the analysis.We

Andrius Melninkaitis
Laser Research Center, Vilnius University, Vilnius, Lithuania The paper of Mario Garcia-Lechuga and David Grojo entitled "Simple and robust method for determination of laser fluence thresholds for material modifications: an extension of Liu's approach to imperfect beams" is devoted to an old problem, namely reproducibility of determination for material's modification threshold and accuracy of beam radius characterization.
In general, I enjoyed reading it: the style of writing and clarity of investigated approach is very acceptable.
Below are my (minor) concerns related to the paper: The definition of "integrated spatial distribution (ISD)" perhaps should be replaced with the term "effective area" as it seems to have a physical meaning and has a dimension of area (square meters or square centimeters).Accordingly, ISD is not correctly used in Figure 1 by stating "both having an identical integrated pulse energy (ISD)."ISD -is not energy it is an area, however, it is easy to fix this discrepancy. 1.
On "generality of these calculations that only depend on the filling parameter β G .": in the present study it is assumed that aperture blocking would in all cases result in perfect airy beams (both software simulation and Pharos laser have relatively good M2 if not perfect), however, it is not clear if this assumption holds for more complicated initial beams that authors are targeting.For example: what would happen if the beam incident to the aperture is non-Gaussian or Gaussian, but has very poor M2 performance?Would the same theory work as well or also other corrections would be needed?Is "PSFLab software" able to help validate that assumption?If not then what are the limits of the presently discussed approach concerning the initial beam? 2.
In the present research image analysis was used for exemplification of improvement ideas and ISD calculation.In practice, ISD could be very much affected by the camera noise and software aperture used for ISD.Two questions related to this context: Could authors compare diameter calculated from ISD and estimated experimentally by using Liu's or their original approach? 1.
What recommendations authors could give to the integration area for measurement of energy when dealing with Airy functions: by how much energy/power detector area must be larger than the central peak diameter of Airy disk? 2.
What is the recommendation for the assessment of correction factors when using truncated beams (for any beam: is it possible to calculate those for unknown truncated input beams in the wavelength range where CCD censors do not work)?

5.
Otherwise, the paper is very good.I would support its indexing.
Finally, the reviewer would like to thank both authors and the editorial team for their patience while waiting for a response (it was unusually long due to a very busy schedule).

Are sufficient details provided to allow replication of the method development and its use by others? Yes
If any results are presented, are all the source data underlying the results available to ensure full reproducibility?Yes Are the conclusions about the method and its performance adequately supported by the findings presented in the article?Partly Competing Interests: No competing interests were disclosed.
Reviewer Expertise: Laser damage testing, beam characterization, optics, laser physics I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.

Author Response 17 Jun 2021
Mario Garcia-Lechuga, Aix Marseille Université, CNRS, LP3, UMR7341, Marseille, France Dear Dr. Melninkaitis, we really appreciate your overall positive evaluation and the constructive questions of your report.We address in detail your comments on the following paragraphs.We thank you for your time and these exchanges that are surely leading to an improved paper.Answer to question 1.We have converted the previous named ISD-integrated spatial distribution to the most common term (as you and the other reviewer referred) EBAefficient beam area.This notation change is introduced throughout all the document.Additionally, we have corrected the mistake you mentioned about Figure 1.Answer to question 2. We consider that we provide a demonstration of validity of this methodology for a non-Gaussian beam on the experiment performed at 1550 nm (section 5C).On that case, by performing a truncation equal to P T = 75% we observed that the truncated function at the focal position corresponds well with the theoretical truncated function expected when illuminating the aperture with a perfect Gaussian beam (paragraph after figure 5).Therefore, as discussed on the first version, we consider that P T = 75% is an appropriate level of truncation for pulses showing a pedestal.For more irregular beams we propose the use of strong truncation (P T < 60%), since the impact of the incoming beam shape on the spot size at the focal position will become less relevant as far as the aperture used becomes smaller.We acknowledge that it is not possible to qualify the level and nature of beam imperfections but this must be taken by the user of the method.Higher is the uncertainty on the beam or the imperfection smaller should be the aperture.Obviously, a limit will be given by the available laser power but very small aperture making possible a uniform plane wave transmitted must lead in theory to the most appropriate condition for applicability of the method In summary, even if we cannot calculate diverse beam shape with PSFlab software (the shape of the beam illuminating the aperture is always considered to be Gaussian), our proposed methodology relies more on the effect of introducing an aperture than on the non-Gaussian incoming beam.Obviously, this is relevant for common laser beam shapes, but hardly applicable to some exotic as higher-order Gaussian beams with donut profiles.This is the reason why we use term 'imperfect' to refer to beams with some uncertainties.Finally, we would like to point out that by introducing the aperture we are obtaining a worse M 2 , since we are transforming expressly the beam to a non-Gaussian beam (Airy-beam like).So, in this methodology the M 2 is a concept to avoid.We thank you again for giving us the possibility to clarify here these important aspects for the reader.Answer to question 3. First, we agree with your comment that experimental images of the beam should be treated correctly, otherwise the efficient beam area obtained can be totally biased by the noise.For minimizing this effect, after recording the image of the beam we recorded under the same conditions (same exposure time) another image by blocking the laser beam.The "background" image is subtracted on the beam image, and those images are the ones presented on the article and on the raw data provided.It should be noticed that no negative values appear, since we kept the resulting image as a 16-bit image (positive integer values) with intensity levels from 0 to 65535.However, we agree that this treatment, as well as the difficulties of perfectly beam recording (alignment and spatial calibration of the imaging system), makes not the process of obtaining the EBA (effective beam area, previously named ISD) free of incertitude.Therefore, and responding first to the second question (3.2), our recommendation with the proposed methodology is to avoid imaging.Our methodology aims to provide a simple solution in which imaging becomes unnecessary, since it relies on the diffractive pattern resulting from introducing a circular aperture.Regarding the first question (3.1), we apologize in advance since we are not sure if we fully understand it.The calculation of a beam waist from the calculated ISD (now named EBA) applicable for a general case (any beam shape) can be a subject of debate and we do not see clearly the relevance of comparisons on this aspect.As a comment, we consider the beam diameter can be only unambiguously calculated from the EBA if the beam profile is Gaussian, since the formula relating both of them are known and simple (EBA=πw 0 2 / 2).
However, generally speaking, we do not recommend extracting the beam diameter from the EBA.About the comparison of the "real beam" (the one image) and the equivalent gaussian beam retrieved with the Liu's method, on section 5B, we compared on the first version of the article the FWHM values obtained by imaging the beam and by performing the Liu's method.We find an excellent agreement between the two values, with a difference of less than a 2%.
Answer to question 4. The aperture should be placed as close as possible to the lens.On the calculations with PSFLab, this parameter is introduced as the pupil aperture of the focusing lens.We thank you to help us to clarify this very important practical aspect.This information is now added in the text (5A) of the new version.
Answer to question 5. On section 5D, now renamed it as "Relevance, precision and recommendations" we have included a new paragraph: Finally, as recommendations for the level of truncation to be applied, we have observed that P T = 75% is an appropriate level when having a pseudo-Gaussian beam with a pedestal.In the case of more irregular beams or totally unknown beams, we propose to perform experiments under two levels of truncation, one with moderate truncation (≥PT= 75%) and another with strong truncation ( PT≤ 60%).Obtaining similar fluence values under the two experiments would indicate that a moderate truncation level is sufficient to transform the beam into an Airy-like beam.Otherwise, the most accurate fluence threshold value is obtained under the strongest truncation, since strong spatial filtering reduces the influence of irregularities of the input beam." material modifications: an extension of Liu's approach to imperfect beams" by Mario Garcia-Lechuga and David Grojo.I really like the idea of a paper that focuses on the accuracy and the precision of the threshold measurements made.In some cases where it is difficult to know the energy distribution in the laser beams, these are supposed to be "perfectly" Gaussian at the focus of a lens.This is not always strictly the case, and the small deviation from this assumption can lead to large errors when Liu's method is used to determine the thresholds (ablation and/or damage thresholds) of a material or even the diameter of a beam.Therefore, the authors had the clever idea of rendering the beam of perfectly known shape using diaphragms to obtain an Airy spot on the component, a spot of known and mastered shape.This makes sense and the theoretical proof is supported by experimental validation.The authors propose a smart method with the keys to put it into practice and use it with relevance.The two aspects of the article, theory and experience, make it compelling and interesting.I therefore recommend this article for indexing and I also recommend that readers consider the tips and guidelines reported in this article.Some points need to be clarified and improved before indexing.

Few main comments:
Accuracy of the method -error bars

○
The first point deals with the accuracy of the method: the purpose of this paper is to obtain much more accurate thresholds thanks to the shaping of beams of perfectly known shape.Very good.It is therefore even more imperative to discuss the correctness of the method, its precision, its intrinsic precision as well as the gain compared to measurements obtained with imperfectly known beams.It seems to me that this point is not discussed.The article would be even more relevant if the reader could actually see the contribution of this experimental approach.In other words, error bars are not discussed with this approach.The paper would be improved with such a discussion.
Vocabulary and notations:

○
The second point deals with the vocabulary and notations: Wouldn't it be as simple to write that w 0 is the beam radius at 1/e², instead of 2w 0 is the beam diameter at 1/e²?Also, in Eq. 8, w Airy is the radial waist at 1/e².OK.In that case, report w 0 as the beam radius in Eq. 1. Otherwise, it is confusing.ISD (integrated spatial distribution) is also named equivalent area in the laser community.This equivalent area (ISD) is therefore the area at 1/e of the peak fluence F 0 .In that case, it would certainly be as simple to report all the equations of section 2 at 1/e instead at 1/e²! Eq. 8: I didn't catch why r'=r*2.5838/wairy , why 2.5838?Is it P/1.22?Limitation of the method:

○
The third point is about one limitation of the method that I would like the authors to make readers aware of: by definition the energy threshold of the material studied is unknown.Therefore how can the users choose with certainty the correction factors given in Table 1.These factors are given for different maximum energies above energy threshold.It is therefore not trivial to select the right one.
Comparison with literature:

○
The last point deals with the discussion about results in this paper and results from the literature: The comparison with Winkler data is probably not appropriate.The agreement between thresholds could be fortuitous considering first that fused silica samples are different but more than this considering that wavelengths are different (1030 nm in this work versus 800 nm in the Winkler's paper).The threshold dependence with wavelengths is reported in reference 5 of this paper [Gallais, 2015].It is also the case in this paper where it is reported that thresholds at 1550nm and at 1030nm are different for fused silica.Therefore, is there a reason to obtain the same thresholds at 1030 and 800 nm for fused silica?Considering again it could be fortuitous, I suggest avoiding such a comparison.
Other minor points that could be taken into account: Title: which beams are considered imperfect?The Airy beams?One might think subsequently that the Airy beams are imperfect!I am not comfortable with the term 'imperfect'.I do not know what an imperfect beam is.More it could be an imperfection in the time domain!I suggest the authors to propose a more precise title.
○ Is caption of Figure 3 right when reporting Equation 12?
○ Could the authors comment the fact on figure 3(a) that the correction factors are higher and smaller than the unity for a given P T (for example P T =0.75 or 0.60)?

○
In the abstract, the notion of 'spatial' should be added: '…an ideal Gaussian spatial profile…' to avoid confusion with 'temporal profile', which is also another topic that could be addressed.
○ Are the conclusions about the method and its performance adequately supported by the findings presented in the article?Yes it is based on a 2D integration of all the space (-infAiry, why 2.5838?Is it P/1.22?" we apologize for having introduced this apparently "weird" number without any context.This value is directly coming from the mathematical definition of the Airy function.According to a normalized Airy function (g(x) = [2•J 1 (x)/x]^2 ), the value corresponding to g(x 1 ) =1/e² (0.1353…) is found at a position equal to x 1 =2.5838... Therefore, for the easiness of definition of an Airy function with a desired waist (w Airy ) , it is needed to make a change of variable according to this number.This corresponds to the expression introduced on Eq. 8.
Small changes have been done to clarify this originates from a "change of variable".Answer to "Limitation of the method" First, we fully agree on your comment on the difficulties on objectively determining a threshold, by nature unknown with a beam that can also be only partly known.It can be also particularly difficult to evaluate whether there is or there is not damage when irradiating with pulse energies very close to the ablation threshold.However, our proposed methodology for determining the energy threshold is not based on 1-on-1 test strategy but in an extrapolation based on the deterministic material response with ultrashort pulses, which is the basis of Liu's method.So, what is referred on Table 1 are the energy threshold values obtained after applying the Liu's method to beams that are not strictly Gaussians.Therefore, when applying a linear regression onto a non-Gaussian beams (not fully linear on the representation D 2 vs ln(E) ) leads to a slightly incorrect determination of the energy threshold .This deviation becomes larger when including into the linear fitting modifications produces with high excitation energies (as represented on Table 1).So, we consider that users of this proposed methodology are not requested to make any assumption on the values of the energy threshold, the correction factors being applied afterwards.Users only need to know two things for using the corrections factors expressed on figure: the aperture transmission and the maximum excitation energy above the energy threshold of modification.This second parameter is a priori unknown, but the biased retrieval of the energy threshold by the linear fitting has no difficulty and can be safely used for this matter (maximum considered energy / retrieved energy threshold).We thank you for giving the opportunity to clarify here this practical aspect for the reader.Answer to "Comparison with literature" We agree that this comparison is not fully necessary and can be considered inappropriate for our 'method paper'.Our aim was only to confirm that our obtained results are somehow consistent with other reports.On the paper of Gallais et al. [ref.5] there is a important figure (fig.5) reviewing the fluence threshold values on fused silica (800 nm and 1030 nm) reported by different authors.We initially wanted to center our discussion based on this figure.Instead, due to the complexity of comparisons, we decided to focus only on the result obtained by Winkler et al. [ref.32], since the irradiation conditions were the closest to the ones we have used (in particular pulse duration of 180 fs).In order to clarify this aspect, we have introduced some changes in the text referring to this figure but also highlighting the importance to be very cautious with comparisons.We can also comment here on the wavelength dependence.We can already confirm that the difference on fluence threshold values in dielectrics on the nearinfrared range (>800 nm) are practically wavelength independent, at least at the pulse duration we performed the experiments (around 200 fs).This interesting response analysis made possible with the method that we introduce in the present paper will be the subject of another publication under preparation.Answer to "Other minor points that could be taken into account" Answer point 1.We have carefully considered the possible change in the title of "imperfect beams" into "non-Gaussian beams".We acknowledge that the term "imperfect" may sound less precise but we believe it better conveys our message.This is because our method does not only treat deviation from a Gaussian hypothesis but can be in principle applied for any imperfect or imperfectly known beam.With other words, a perfectly known beam even non-Gaussian can be considered as ideal beam for a threshold determination without our method.Answer point 2. "Table 2 should arrive before Figure 3".We agree (and it was like that on our submitted word document).We will do our best to ask to the editorial team to introduce this change on the new published version.Answer point 3.You are totally right.There is a displacement on the equation numbering that is corrected in the new version.For (a) should be eq.10 and for (b) should be eq.11.
Answer point 4. We agree that those correction factors are not totally intuitive.The explanation of having a correction factor higher than 1 means that the equivalent Gaussian beam (what the Liu's plot provides us) exhibits a larger EBA than the truncated beam.This occurs when selecting only irradiations at moderate excitation levels (e.g. 2 times the E_th), which is equivalent to thresholding only the top part of the truncated beam (similar shape than a Gaussian beam) but not the lower part (narrower than a Gaussian beam).However, as shown in table 1 on the particular case of doing this procedure on a perfect Airy beam (P T ideally equal to 0) fluence correction factors are always smaller than 1.On that case, even if the area of the central lobe is overestimated when accounting the equivalent Gaussian beam (Liu's plot), the total EBA is underestimated because of the contribution of the rings of the Airy beam (which are inexistant on the retrieved equivalent Gaussian beam).So, finally, addressing the particular truncation you mention (P T =0.75 or 0.60), and given the fact that in these cases the contribution of the surrounding rings to the total EBA is less pronounced than for a perfect Airy beam, it remains totally possible to obtain fluence correction factors above or below 1 depending on the maximum considered energy accounted for the linear fitting process.We thank you again for giving us the possibility to clarify also this important aspect for the reader.Answer point 5.We included on the abstract the notion of spatial.("spatial Gaussian beam").Answer point 6.The colors are corrected on the new version.We have followed the same color criteria than the one in Figure 3. Answer point 7.You are totally right.On the last paragraph of Page 9 we should have referred to Figure 4(a).We introduce this change on the new version.
Competing Interests: No competing interests were disclosed.

Figure 1 .
Figure 1.(a) Numerically calculated radial profile of an Airy disk (w Airy = 10 μm) and a Gaussian beam both having an identical effective beam area (EBA).The inset is the 2D profile of the Airy disk.(b) Representation according to Liu's method for the Airy disk represented in (a) and two linear regressions for different ranges of considered fluences: one with fluences up to 2 times above the ablation threshold (2 ⋅ F th ) and another up to 10 times (10 ⋅ F th ).(F 0 is the pulse peak fluence, F th is the fluence threshold for ablation, E is the pulse energy, E th is the energy threshold for ablation and R 2 is the radius square of the induced crater).

Figure 2 .
Figure 2. Numerically calculated radial beam profiles at focus position of a lens illuminated by a Gaussian beam truncated by circular apertures.The calculations assume the following parameters: λ=1030 nm, focal lens of 50-mm, an incident collimated Gaussian beam diameter of 5-mm (2•w) and varying aperture diameter adapted to correspond to the mentioned aperture transmissions.Data are obtained by using PSFLab software 28 .

Figure 3 .
Figure 3. (a) Correction factor, η F , to apply to the fluence calculation (Equation 10) when irradiating with truncated beams by a circular aperture.(b) Correction factor, th E η , to apply to the fluence threshold calculation (Equation11) when irradiating with truncated beam by a circular aperture.In both cases, the horizontal axis corresponds to the maximum considered energy for the linear regression of the Liu's method.Straight lines joining the points serve for view guiding.

Table 3 .
Correction factor, η F , to apply to the fluence calculation (Equation 10) when irradiating with truncated beam by a circular aperture.The values are represented graphically on Figure 3(a).

Figure 4 .
Figure 4. (a) Ablated areas in sapphire and fused silica as a function of the pulse energy.Craters are produced by single pulse irradiation (170-fs) at 1030 nm with two different beam profiles depending on the presence or not of a circular aperture (P T = 75 %).(b) (left) Image of the beam without any aperture at the best focal position (f=50 mm) as captured by an imaging system equipped with InGaAs array detector.(right) Horizontal beam profile at the central position of the beam image.

Figure 5 .
Figure 5. (a) Ablated areas in sapphire and fused silica as a function of the pulse energy.Craters are produced by single pulse irradiation (190-fs) at 1550 nm with two different beam profiles depending on the presence or not of a circular aperture.(b) Beam image produced at the focal position after directly focusing the OPA beam.(c) Beam image produced at the focal position after focusing the same beam but inserting a circular aperture (P T = 75%) before the lens.

Figure 2 :
Figure 2: it is difficult to distinguish the color of the different beam profiles.○