A detailed study through the focal region of near-threshold single-shot femtosecond laser ablation nano-holes in borosilicate glass

detailed study through the focal region of near-threshold single-shot femtosecond laser ablation nano-holes in A detailed study of the morphology of nano-craters drilled in borosilicate glass by single shot femtosecond laser ablation near the ablation threshold has been performed by scanning electron microscopy, atomic force microscopy and scanning electron microscopy imaging after focused ion beam sectioning. The in ﬂ uence of the numerical aperture (NA=0.4 and 0.8), the pulse energy (16 nJ b E p b 600 nJ) and the position of the specimen surface into the focal region were systematically investigated, leading to nanometric or micro- metric scales in every spatial dimension. The nanocrater ’ s size is not restricted by the diffraction limit but determined by the laser pulse stability and the material properties. If the beam is focused inside the glass, two craters are drilled, shaping very distinct morphologies. Their dimensions have been studied in details and different relationships have been proposed for the evolutions of the depths and of the various diameters of these craters as functions of the pulse energy, the numerical aperture and the position of specimen surface in the beam-material interaction region. It is suggested that the long, thin conical pro ﬁ le with very high aspect ratio of the secondary craters is due to a spontaneous reshaping of the beam which transforms the incoming Gaussian pulse into a Gaussian – Bessel pulse. As proposed in the developed model the geometry of the second craters seems to be connected with the one of the main craters.


Introduction
Femtosecond laser micromachining is a versatile materials processing technology for fabricating a wide range of micro and nano-structures in transparent media [1,2]. Moreover, direct laser writing of patterns enables rapid prototyping without mask. So, ultrafast laser processing continues to be an intense subject of research from fundamentals and modelling to a large field of applications . The deterministic machining at submicrometer scale has caught much attention through applications to nano-structure for memories, micro-optical (wave guides, optoelectronic systems, couplers,…) and fluidic devices (micro fluidic channels,…), mechanical modification, sensors, wettability modification,… [3][4][5][6][7][8][9][22][23][24].
For dielectric materials such as fused-silica or glass, non-linear absorption becomes significant when ultrashort laser pulses are tightly focused by lens with high numerical aperture onto the surface or inside the bulk of the specimen. The intensity in the focal volume becomes sufficient to initiate a series of non-linear physical processes : multiphotons and tunnel ionization, avalanche ionization, carriercarrier scattering, carrier-phonon scattering, material phase change, shock wave emission which permanently changes the refractive index, thermal diffusion, material ejection followed by condensation and resolidification of clusters or melt material . If the laser pulse energy is adjusted so that only a small part of the focused beam has irradiance above the ablation threshold of the considered material, sub-diffraction nano-structuration can be machined [3][4][5][6][7][8][9].
The purpose of this paper is to extract relevant and detailed information on the morphology through the focal volume of nano-holes drilled by near-threshold circularly polarized single-shot femtosecond laser ablation in borosilicate glass under ambiant temperature. Note that circular polarization leads to quasi-symmetrical crater, as shown further in this paper, contrary to linearly polarized femtosecond laser shot which leads to strongly anisotropic nano-craters [14,[18][19][20]. Moreover, due to light reflection on the crater walls, the energy deposition in the crater area is maximum for circularly polarized beam [19]. The dependence of the machined craters geometry on the laser pulse energy, the numerical aperture and the position of the specimen surface through the focal volume is qualitatively and quantitatively reported.

Experimental setup
The amplified laser source (Spitfire ProV, Spectra-physics) emits 100 fs laser pulses with a central wavelength λ of 800 nm, at a repetition rate of f = 5 kHz. An independent Pockels cell system with a thin-film polarizer plays the role of an optical shutter that enables single-shot illuminations. An ensemble of neutral densities allows adjusting the pulse energy. A polarizing cube and a zero-order quarter wave plate allow the production of a circular polarization before the laser beam passes through a microscope objective (MO), focusing onto a silica plate. Two types of Olympus Plan-fluor infinitelycorrected microscope objectives were used: ×20 with 0.4 numerical aperture (NA) and ×50 with 0.8 NA. As the beam diameter is far larger than the entrance aperture of the MO, the Airy spot sizes are 2.4 μm and 1.2 μm.
The glass slides were mounted over a 3D positioning motorized stage (Newport ILS M-VP25) with bidirectional repeatability better than 200 nm and the sample orthogonality with respect to the beam propagation was ensured to be less 1 mrad. The positioning of the samples was achieved by imaging on a CCD camera with depths of field of 2 μm and 0.5 μm with the ×20 MO and ×50 MO, respectively.
Specific care has been devoted to the cleanliness of the beam. The dispersion of all optics was pre-compensated with the compressor of the chirped-pulse amplifier of the laser chain. The precompensation was carefully adjusted by measuring the pulse duration τ with a GRE-NOUILLE after the laser beam passed through the microscope objective and was collimated by a thin lens with negligible dispersion. The measured pulse duration is about τ = 120 fs.
The transmitted laser power P, just after the microscope objective, has been measured with a calibrated power-meter (GENTEC, XLP 12) whose measure range is 10 μw b P b 3 W corresponding to a pulse energy E p in the range 2 nJ b E p b 6 10 5 nJ. Single-shot illuminations of the sample were performed under atmospheric conditions, for different pulse energies in the range of 16 nJ to 600 nJ, near the ablation threshold of silica. For each power, the sample was translated through the focal region, in the vertical direction Z, by steps of 250 nm over a range of 15 μm. After each laser shot, corresponding to a fixed z value, the sample was translated in plane by 5 μm in the X direction. After eight shoots on the same line in the X direction, a new line of eight shoots was realized with the same origin of the first one, but translated of 15 μm in the Y direction. Then, this sequence has been repeated ten times, covering the Z focal region over 15 μm. This procedure allows a compact area of the silica plate to be nano-structured (Fig. 1). To test the reproducibility of this method, six specimens (noted Spec. 1-6) with the same or different nano-structurations have been characterized.

AFM and SEM measurements
Atomic force microscopy measurements have been performed on each illuminated zone in non-contact mode, with a spatial resolution of 20 nm in the sample plane and with 5 nm in depth. The microscope is a PSIA XE-150 with a cantilever tip with minimal spring force constant of k = 51 N/m and whose resonant frequency is about 130 kHz. The radius of curvature of the tip is 8 nm and half-width angle of 18°. Scanning speed was typically 0.16 Hz. All measurements (depths, angles, diameters) were extracted from the AFM image by averaging four different profiles through the crater center.
While it is straight forward to measure the in-plane morphology of the nano-holes, it is not so easy to characterize the depth with high aspect ratio thanks to AFM measurements. So, in order to compare with the AFM measurements and to definitively characterize the exact shape of the crater, after metallization, selected nanoholes have been measured by direct scanning electron microscopy (SEM) imaging after focused ion beam (FIB) sectioning. A dual SEM/ FIB device (LED 4400 (SEM)/Orsay-Physics 31 (FIB)) has been used to section the nano-holes. The imaging has been performed with a high resolution SEM (Raith-Eline). After FIB milling and a sample tilt of γ (γ = 35°, 45°or 50°) in the SEM, the projected images of the sectioned holes have been observed. Thus, the true heights h are given by h = h measured/sinγ.

3.1.
Investigation of the focal region as a function of the pulse energy Ep and NA Fig. 2a shows, as a function of the pulse energy E p for NA = 0.4 and 0.8, the length ΔZ of the focal region where a visible laser-surface sample interaction has been detected by SEM imaging. Indeed, ΔZ (μm) = 0.25 N, where N is the number of observable impacts. Note that the reproducibility of these determinations from one specimen to another one is fairly good as shown in Fig. 2a. Of course, ΔZ is an increasing function of E p and a decreasing function of NA. The value of the ablation threshold energy E po is about 20 nJ and 70 nJ for NA = 0.8 and 0.4, respectively.
Considering a Gaussian beam propagating in free space, the spot size ω(ξ) along the optical axis Z is given by : where z o is the Rayleigh range and ω o is the beam waist radius in air. M ≥ 1, is a parameter which characterizes the divergence of the beam induced by the diffraction. For a perfect beam M = 1 and for the present device M b 1.1. α is a parameter which depends on the ratio of the Gaussian beam diameter d at 1/e 2 of the peak intensity of the laser and to the aperture diameter D of the lens. At present d/ D ≅ 1.9, thus α ≅ 0.43. Note that for a classical Airy disk α = 0.61. Now, considering the fluence F of the beam and the fluence threshold F o , thanks to the relation (1) it is easy to show that ΔZ measured in the air is given by the following relation: For F/F 0 = 2, ΔZ/2 = z 0 . The relation (2) is recasted as: The parameters a 1 and a 2 have been determined thanks to the  [25], F o ? ≅ 2.5 J cm − 2 for λ = 800 nm and τ = 120 fs [26], 2 b F 0 b 3 J cm − 2 for λ = 798 nm and τ = 120 fs [27]. The values of F o and α reported in a previous analysis [24], F o = 2.75 J cm − 2 and α = 0.417 correspond to M = 0.88, thus are inconsistent with the inequality M ≥ 1. Note that the value of F o depends on the wavelength λ [4,9] and on the pulse duration τ [6,11,14,[25][26][27].

Description of the general aspects of the nano-holes
As a function of the pulse energy E p , the NA and the z position in the focal region (0 b z b ΔZ), two kinds of nano-holes morphologies have been observed. The nano-holes were composed of a single crater or of two characteristic craters with very distinct geometries [24,28]. and with two craters (E p = 108 nJ, NA = 0.4, z = 6.75 μm for ΔZ = 7.5 μm), imaged by SEM after SEM/FIB cross-sectioning. The conditions of the appearance of the second crater will be described later.
The two craters are quasi-axisymmetric along the laser beam axis and surrounded by an hemi-torus rim. The profile of the main crater is trapezoïdal or more exactly looks like a diabolo with a convex bottom (Fig. 3a) and the second crater is conical with a high aspect ratio (Fig. 3b). The morphology of the main crater is typical of thermal ablation, i.e.: first, during phase explosion, small clusters are ejected and condensed, especially when the ablation is performed under atmospheric pressure, and secondly, thermocapillary forces can extract liquid material out of the center of the crater due to the thermal gradient (Fig. 3a). This effect is known as Marangoni effect [17]. On the contrary, the second crater is very clean and no rim have been observed in its vicinity (Figs. 3b or 8). This observation is typical of Coulomb explosion, soft ablation with clean craters and no protusion [10]. Fig. 3c gives the scheme of a typical profile as well as the characteristic dimensions investigated in this study : L 1 , L 2 , h 1 for the main crater and L 3 , h 2 for the second crater. Thanks to AFM and SEM imaging the dimensions of the primary craters have been evaluated through the focal zone, along the Z direction. The raw results have been presented in a previous paper [24] but are re-analyzed in this paper. The origin z = 0 corresponds to the first observable impact as the waist of the beam is outside of the glass and z = ΔZ, to the last visible impact with the focal point inside the material. To study the evolution of the parameters L 1 and h 1 with respect to z and to obtain a symmetric abscissa whatever ΔZ, those parameters have been plotted as a function of z′ = z − ΔZ/2, thus z′ = − ΔZ/2 for z = 0 and z′ = ΔZ/2 for z = ΔZ (Fig. 1) [24]. All of the obtained profiles are symmetrical with respect to the origin z′ = 0 (z = ΔZ/2) and maximum for this value when the fluence is maximum, i.e. when the beam is focused at the sample surface. At the neighbourhood of z′ = 0 and for the two NA, the maximum mean values bL 1max > and bh 1max > as a function of the pulse energy E p have been represented in Figs. 4a and 5a, respectively.
Due to non-linear nature of the interaction of the femtosecond laser pulses with transparent material, simultaneous absorption of several photons is required to initiate ablation. Multi-photon absorption produces free electrons which are accelerated by the electric field of the laser. These electrons induce avalanche ionization and optical breakdown and thus generate a micro-plasma. The expansion of this plasma creates nano-structure on the surface of the sample. If each discrete photon of the beam is described by a Gaussian distribution in both time and space and moreover if each photon is absorbed according to a look like Beer-Lambert law for Gaussian pulse, it is possible to write the relations (4a, b, c):  In a previous paper [24], f(z′) has been identified as: In these equations, ξ is the distance along the propagation direction, x is the lateral distance from the optical axis, ω the width of the beam function of the position z′ of the waist ω o (Fig. 1). 1/A (α*, q, R,..) is an absorption coefficient which depends on the linear absorption coefficient α* of the material, the number q of required photons to initiate ablation (q = 6 for Si0 2 ), the reflectivity R of the surface, etc. [29]. β L is a coefficient of the model which will be identified and discussed further. Hence, from the inversion of relations (4a) and as shown in the literature [8,16,25,29], a logarithmic dependence of the crater's width L 1 and of the depth h 1 on the fluence is deduced as follows: and As f(0) = 1, thus: bL 1max > and bh 1max > are the maximum mean values of the width and of the depth of the craters.
Taking into account the relations (1) and (2), Eqs. (6a), (6b) are recasted as: The relation (7a) has been identified in Fig. 4b where NA b L 1max > has been plotted as a function of NA 2 E p . The fit of the experimental points, continuous curve in Fig. 4b, gives 2β L αλ = 535 nm. As α/ For the parameter bh 1max >, the same kind of representation as for bL 1max > in Fig. 4b, i.e. NA bh 1max > versus NA 2 E p , has been plotted in Fig. 5b. Then, the experimental points have been fitted by the relation 7b in which A(α*,q,R) = 230 nm/NA; continuous curve in Fig. 5b. This last relation seems to show that, for a given wavelength : with β h M = 0.59 for the studied material. Thus, the coefficient A(α*, q, R,..) depends on NA and is equal to 59% of the waist ω o if M = 1.
Considering the minimum diameter L 2 of the craters (see Fig. 3c) and taking into account the half opening angle of the AFM tip (~18°), the true value of L 2 can be deduced from the measured one. From the SEM/FIB imaging, the L 2 values are directly estimated. As shown in Fig. 6 where the true values of L 2 have been plotted as a function of L 1 or bL 1max >, whatever the pulse energy E p , the NA and the position z′ of the waist, a linear relationship has appreciably been found between these two parameters: with: α 21 ≌ 0.85 ± 0.05 and β 21 ≌ 100 ± 30 nm for NA = 0.4 and 0.8. The relations (6a), (9) and (11) involve a decrease of the walls inclination as 1/L 1 with the increase of the pulse energy. For very small impacts, near the fluence threshold, L 2 = 0 and L 1 = β 21 /α 21 , the shape is nearly conical with a large opening angle (>45°) as shown in Fig. 7.
As a conclusion of the morphology study of the main craters, combining Eqs. (2), (4c), (5a), (5b), (9), (11) the evolutions of h 1 , L 1 and L 2 in the focal region can be written as a function of ΔZ, F, NA, z′ and the characteristics of the laser beam:  The long, thin conical profile shown in Fig. 3b suggests a self-focusing of the laser pulse [9,30] or a nonlinear Gaussian-to-Bessel beam transformation [31][32][33] as the beam propagates into the nonlinear sample with a sub-diffraction sized diameter. The first possibility, selfchanneling of the beam, is supported by the balance among diffraction, Kerr-induced self-focusing and plasma-induced defocusing. The spontaneous reshaping of the beam is based on a dynamic spatial replenishment where the filament appears continuously regenerated by subsequent focusing of different temporal slices of the pulse. So, once the nonlinear absorption becomes significant, it transforms the incoming Gaussian pulse into a Gauss-Bessel pulse and leads towards filamentation [31,33]. Another possibility which can generate quasi-Bessel beam is the creation of a virtual micro axicon [33,34] by the plasma confined into the primary crater. Indeed, when ultra short Gaussian laser pulse is focused by an axicon, unsteady and steady Bessel filamentation regimes appear leading to a strong multiphoton absorption and plasma generation over long distances [33]. These two possibilities will be analyzed further, after the presentation of the experimental results, i.e. condition of apparition and geometries of these craters.

Experimental results and first discussion
As a function of the pulse energy, thanks to SEM imaging of all the craters, the precise conditions of the appearance of the second craters have been observed and analyzed through the focal region. Then, its width L 3 has been measured as a function of the position of the waist of the beam. An example of such surface imaging is given in Fig. 8 (E p = 161 nJ, NA = 0.4 and z′ = 3.25 μm). Note that the second crater is perfectly centered with respect to the first. The appearance of the second crater begins in the neighborhood of z = ΔZ/2 (z′ = 0), i.e. as the beam focused at the surface or inside the material. They suddenly disappear when z′ = ΔZ/2. Fig. 9 shows the evolution of the width L 3 as a function of the z′ distance in the focal region for various pulses energies and the two NA. It is thus shown that the distance z′ = δ at which the first second crater occurs appears to depend on the NA and to be a decreasing function of the energy E p . On the contrary, the maximum mean value bL 3max > of L 3 is a slowly increasing function of E p : 150 b (bL 3max >) b 300 nm (Fig. 10a). In this figure, the L 3 value reported by White et al. [28] for NA = 0.85 and E p = 120 nJ has been drawn. To quantitatively compare our results to those of the literature, the L 3 values reported by White et al. [28] (input Gaussian beam) and Bhuyan et al. [35] (true Bessel beam) have been plotted as a function of E p /E p0 in Fig. 10b, where E p0 is the experimental value of the threshold energy reported in these papers. For E p /E p0 b 3.3, all the experimental determinations are close, 100 b L 3 b 300 nm, which shows that the L 3 diameter is appreciably independent on the numerical aperture (NA = 0.4, 0.8, 0.85) and to the half conical angle α Β of the Bessel beam (α Β = 17°and 11°in the borosilicate glass [35]), contrary to Ep 0 which varies as NA 2 . The NA values corresponding to these two α Β angles, NA = sinα Β /n 0 with n 0 = 1.5, are equal to NA = 0.19 and 0.127 for α Β = 17°and 11°, respectively. Now, if NA 2 E p0 = 11 nJ, as previously reported, E p0 = 304 nJ and 682 nJ, which is in good agreement with the experimental values, E p0~3 20 nJ and 610 nJ, respectively [35]. For E p / E p0 > 3.3, above the transition regime, due to machining by the lateral lobes, a linear relation is observed [28,35,36]. Courvoisier et al. [36]    reports that the fluence needs to stay below three times the ablation threshold to avoid such machining by the lateral lobes.
The secondary craters have only been observed above a critical energy E pc , slightly greater than those of the threshold E po of the main craters. E pc = 80 nJ for NA = 0.4 and E pc = 20 nJ for NA = 0.8 (Fig. 10a). Indeed, as previously shown, NA 2 E po = a 1 = 11 nJ for the principal craters and E pc NA 2 = a 3 = 13 ± 0.5 nJ for the secondary crater. This last value will be discussed further and compared to those of the critical energy for the self-focusing appearance. Eventually, Fig. 10c gives the evolution of b L 3max > as a function of NA 2 E p . It is very interesting to observe that when L 2~L3 , for very low pulse energies in the vicinity of the threshold, there is no discontinuity between the two craters, as shown in Fig. 11, only a smooth transition in the walls inclination, unlike general morphology (Fig. 3b). In this example, E p = 31 nJ, NA = 0.8 thus L 2 = 330 nm (Eqs. 7a,b and 11) and L 3 = 210 nm (Fig. 10a).
From Fig. 9 and as previously mentioned, δ and ΔZ/2 are decreasing and increasing functions of E p , respectively, so when δ = ΔZ/2 no second crater has been detected. This condition occurs at about NA 2 E pc = 13± 1 nJ, which agrees the previous observations. Moreover, to normalize the representations given in Fig. 9 taking into account the values of δ and ΔZ/2 for the abscissa and bL 3max > for the Y-axis, the new coordinates in Fig. 12 are z*= (z− (δ + ΔZ/2))/(ΔZ/2 − δ) (Fig. 1) and L 3 /bL 3max >. When z = δ + ΔZ/2, z*= 0 and when z = ΔZ, then z*= 1. In this figure, whatever the pulse energy E p and the NA are, all points are appreciably on the same locus. The possible analytic relationships of this locus will be discussed below.
The depths h 2 of the secondary craters have been determined from the SEM imaging of the tilted SEM/FIB sectioned holes (see Fig. 3b). For a given energy, to determine the maximum mean depths bh 2max > of the secondary craters, some measurements have been carried out in the neighborhood of z*= 1/2, (z= (3ΔZ/2 + δ)/2 ≈ 3ΔZ/4), when the hole diameter L 3 is maximum (L 3 /bL 3max >=1, Fig. 12). The results are presented in Fig. 13a for the two NA. It is obvious that bh 2max > greatly depends on the NA and very deep holes (bh 2max > ≥ 5 μm) have been observed for NA = 0.4 and E p > 200 nJ. This observation is in a fairly good agreement with the study of White et al. [28] which report craters such as : 2.3 b h 2 b 2.7 μm for 0.8 b E p b 1.2 μJ and NA=0.85. Fig. 13b where NA 2 bh 2max > has been plotted as the function of the pulse fluence (F α NA 2 E p ) shows that the depth bh 2max > appreciably varies as NA 2 . Contrary to the main craters (Eq. 10), the morphology of the secondary craters presents a high aspect ratio, i.e. : bh 2max >/b L 3max > is in the range 4 to 5 and 16 to 19 for NA= 0.8 and 0.4, respectively, when 20 b NA 2 E p b 50 nJ.
Before concluding on the experimental results some observations and remarks have to be done: -Often, after a first FIB milling, at the bottom of the secondary crater and as shown in Fig. 14 where the sample has been tilted of 5°, a shock-affected zone which was not ablated by the FIB is clearly observed. In this case, to measure the total depth h 2 of the hole it has been necessary to perform a second milling with a more energetic ions beam. Indeed, in this deep zone the material has to be expelled from the hole into a shell of very compressed material to comply with the mass conservation. In this example, E p = 130 nJ and NA = 0.4, the maximum diameter of this compressed zone is about 380 nm and its length around 1300 nm. These dimensions are increasing functions of the pulse energy; for E p = 300 nJ and NA = 0.4, values of 620 nm and 3300 nm have been measured for the diameter and the length, respectively. Another interesting  observation concerns the top of this shock-affected zone which is appreciably independent on the pulse energy and appears for a depth of about 2800 ± 300 nm. The exact bh 3max > values of the length of the compressed zone (bh 2max > minus the distance of the top of the compressed zone from the surface) have been reported in Fig. 13a. Thus, for E p b 120 nJ ((bh 2max >) b 2800 nm) there is no shell, which is the case presented in Fig. 3b for NA = 0.4 and for NA = 0.8 where no shell has been detected since bh 2max > b2800 nm. In conclusion, for depths lower than 2800 nm the ablated material is vertically ejected through the second crater. For greater depth the material seems to be radially expelled to form a shell of compressed material. These establishments agree with the work of Juodkazis et al. [37] on laser-induced confined micro-explosion in bulk silica and where the threshold for the formation of a densified region surrounding the void has been calculated as 13 nJ for NA = 0.9 in air. Thus NA 2 E p = 10.5 nJ which corresponds to NA 2 E p = 17.6 nJ (NA = 0.4, E p = 110 nJ) (Fig. 13a) determined in the present study. Moreover, the thickness of this compacted zone is about 480 nm for E p = 26 nJ (NA 2 E p = 21 nJ) [37], in agreement with the present work where this thickness is of the order of 400 nm for E p = 130 nJ (NA 2 E p = 20.8 nJ). -As shown in Fig. 14 and in all imaged second craters, a characteristic damped oscillation on the hole diameter has been observed at the top of the second craters. The period of the oscillation has been estimated at 700-800 nm, which is of the order of the wave length of the laser beam. A possible explanation will be proposed further. -The last interesting and surprising quantitative observation concerns the relation between the total depth bh Tmax > of the craters and the length ΔZ previously determined. Indeed, in the studied    range of pulse energy, in a first approximation, it has been observed that: The relations ΔZ/2 and bh 1max > as a function of the pulse energies are given by the Eqs. (2) and (9), respectively, and if the equality (13) is verified then NA 2 b h 2max > can be written as: In this relation all the parameters are known and for the two values of NA it has been plotted in Fig.13b. There is a fairly good agreement, especially if the difference between the threshold energies is considered: a 1 = NA 2 E p0 = 11 nJ for the primary craters but a 3 = NA 2 E pc = 13 nJ for the secondary craters. So, in a first approximation and particularly for the higher fluencies when δ tends to zero, the empirical relationship bh Tmax >~ΔZ/2 is verified and can be utilized for practical applications; ΔZ is very simple to determine contrary to the total depth of the craters.

Interpretation and proposed model
As previously mentioned, to analyze the experimental results two possibilities are plausible: a-self-channeling of the beam; b-spontaneous reshaping of the beam which transforms the incoming Gaussian pulse into a Bessel pulse. a. Self-channeling with a sub-diffraction sized diameter is due to the Kerr effect, the spatial variation of the refractive index causes a lens-like effect that tends to focus the beam inside the material, and to the formation of a micro-plasma which can defocus the beam and balancing the self-focusing. Self-focusing is predicted if the power peak of the laser pulse exceeds a critical power Pc given by [9,30]: where n o and n 2 are the linear and non-linear parts of the refractive index of the material. With n o = 1.48 and n 2 =3.7×10 − 16 cm 2 W − 1 for the silica, α/M = 0.491, P c is equal to 1.12 MW which corresponds to a pulse energy of E p = 140 nJ. This value is greater than the two measured values for the apparition of the secondary crater, E pc = 80 nJ and 20 nJ for NA= 0.4 and 0.8, respectively. Moreover, Pc is independent of the numerical aperture, contrary to the threshold energies since NA 2 E pc = 13 nJ. Hence, self-focusing is not a convincing scenario to explain the formation of such long secondary craters with very high aspect ratio. b. The second possible scenario which could explain the morphology of these craters is the reshaping of the Gaussian input beam into a conical wave [31,32] which is a property of Bessel beams. Following Dubietis et al. [31] this reshaping is driven by the requirement of maximum localization, maximum stationary and minimum non-linear losses. Another possibility which could generate a Bessel beam is that the complex geometry of the primary crater behaves as a virtual micro-axicon or a combination of a spherical lens with a micro-axicon. In this case, the morphology of the second craters should be linked to the one of the main crater imposed by the input Gaussian beam. Indeed, as previously mentioned, when ultra short Gaussian laser pulse is focused by an axicon, steady and unsteady Bessel filamentation regimes appear leading to a strong multiphoton absorption and plasma generation over long distances [33]. This possible scenario is quantitatively studied below.
In the case of Gaussian input beam, the intensity distribution behind an ideal axicon with infinite aperture can be calculated by resolving Fresnel diffraction integral with the stationary phase approximation [38]: where In this equation x denotes the radial distance from the optical axis, ξ the longitudinal position along the propagation direction and α o the half cone angle of the resulting k vectors after passing through an axicon. J o is the zeroth-order Bessel function of the first kind and ξ max the maximum range of the beam propagation. Note that I(x,ξ) is maximum for ξ M = ξ max /2.
it is easy to show that : Taking into account the relation (18) and the fluence threshold energy F 0c for the second crater appearance, by setting bL 3max >/2 = β LS x, the inversion of the relation (17) gives the evolution of the diameter bL 3max > for an ideal axicon with infinite aperture as a function of the fluence: In this relation β Ls is a coefficient of the model and a 3 = NA 2 E pc = 13 nJ.
Concerning the evaluation of the maximum depth bh 2max >, thanks to Eq. (16), it is easy to show that: As for relation (17), this expression is not reversible. However, as it will be shown further in Fig. 12, for ξ/ξ max b 0.7 a very good approximation of this expression is given by the following equation: Thus, the two relations (20) and (21) give: Now, the inversion of the last relation (Eq. 21) allows the dependence of bh 2max > = β Ls ξ for an ideal axicon as a function of the fluence: β hs is a coefficient of the model. Note that when F is large then bh 2max > = β hs ξ max .
In the proposed model it has been assumed that the main crater associated with the formation of a dense plasma create a kind of micro-axicon in the plasma with a tip oriented towards the surface of the impacted sample and whose diameter is L 1 (or L 2~α21 L 1 ). If θ 1 is the half cone angle of the resulting k vectors of the generated Bessel beam, thus: and moreover, if θ is the micro-axicon base angle, the Snell's law gives : where n o is the linear refractive index of the glass. The θ angle is assumed to be close to the beam divergence angle of the input Gaussian beam; i.e. of the order of: Solving Eq. (24) with θ given by the relation (25) and considering the first order development of tgθ, Eq. (24) can be recasted as Eq. (26) within a maximum error of 3% for the higher numerical aperture.
The maximum value bL 3max > of L 3 (Idem for bh 2max >) is obtained for z* = 1/2 (Fig. 12), thus for z′ = ΔZ/4 + δ/2. As shown in Fig. 9, for the higher values of the pulse energy, the δ parameter can be neglected, thus z′~ΔZ/4 (z = 3ΔZ/4) and f(z′)~√3/2 (Eq.4c). On the contrary, for the lower values of E p , δ tends towards ΔZ/2 (Fig. 9) and for E p = E pc , δ = ΔZ/2. Therefore z′ = ΔZ/2 and f(z′) = 0. To take into account these two boundaries conditions without knowing the exact dependence of the δ parameter on the pulse energy, Eq. (5a) has been recasted as: Now, combining relations (22) (23) (26) and (27), the evolution of bh 2max > with the pulse energy can be deduced: According to Eq. (28a) bh 2max > varies approximately as 1/NA 2 as shown in Fig. 13b. The β hs parameter has been determined by fitting the experimental results in Fig. 13b, NA 2 b h 2max > = f(NA 2 E p ), with the relation (28a) (solid lines). The best fit is obtained with β hs = 0.77, the other parameters being previously identified. β hs = 0.91 if the entrance diameter of the micro-axicon is equal to L 2~0 .85 L1. If no restriction is applied for the apparition of the second craters, contrary to Eq. (27), the comparison of Eq. (28a) with Eq. (22) gives the opening angle α o of the k vectors of an equivalent Bessel beam formed by a perfect axicon illuminated with a plane wave ; Note that the base angle of this axicon varies as Ln F Fo −1=2 ; thus 0<α o < π 2 when 1<F=F o <∞: Indeed, the α o angle cannot be greater than the opening angle u of the objective (u = sin − 1 NA): i.e. α o ≤ u. From Eq. (29), the equality α o = u allows to determine the condition for the apparition of the secondary craters : For NA= 0.4, the previous relation gives NA 2 E pc =a 3 =13.1 nJ and for NA = 0.8, NA 2 E pc =a 3 = 12.1 nJ. The approximation of relation (30) allows a 3 = 13.4 nJ whatever the NA. These values are close to the experimental value previously mentioned: NA 2 E pc =a 3 =13±0.5 nJ.
To estimate the validity of Eq. (28a) on a larger range of pulse energy and to compare our results to those of the literature, as for b L 3 > in Fig. 10b, the variations of the product NA 2 h T , where h T is the total length of the craters (h T = bh 1max >+bh 2max > in our study), have been plotted as a function of the pulse energy E p normalized by the threshold energy E po , in Fig. 15. All the experimental points define appreciably a same domain for an input Gaussian beam (White et al. [28] and this study 0.4 b NAb 0.85) or a true Bessel beam (Bhuyan et al. [35], 0.13b NA b 0.19). The discrepancy obtained for NA = 0.13 will Fig. 15. Evolution of NA 2 h T determined from the literature results [28,35] and in the present study as a function of the fluence ratio E p /E p0 . Predictions of the modelling: Eq. (31). be discussed below. The predictions of the modeling have been valued thanks to Eqs. (9) (bh 1max >) and (28a) (bh 2max >) which have been recasted as: In this equation, all the parameters are known. Note that for the true Bessel beam, the second term of this equation is negligible (NA is small). The curves drawn in Fig. 15 correspond to relation (31) for the different values of NA. A fairly good agreement has been obtained excepted for NA = 0.13. However, in the Bhuyan et al. [35] experiments due to the experimental set up configuration, the maximum length h T of the craters was fixed to 32 μm. So, solving Eq. (31) with h T = 32 μm allows to determine the domain of validity E p /E po where the experimental results should match relation (31). For α B = 17°(NA = 0.19), E p /E po ≤ 5.5, thus all the experimental points verify this inequality, but for α B = 11°(NA = 0.13), E p /E po ≤ 1.75 and only a small part of the experimental points have to be considered. The wrong points are noted in parentheses in Fig. 15.
The diameter bL 3max > corresponding to an infinite ideal axicon is given by the Eq. (19) and at present as α o (Eq. 29) is small, this relation predicts too large values depending of NA, which is contrary to the experimental results shown in Fig. 10b-c. However, it is interesting to remark that if β Ls /sinα o = 1 in Eq. (19), (the minimum value of this ratio; β Ls = 1 and α o = π/2) thus: This relation perfectly match the experimental results as shown in Fig.10b-c (solid lines). The same result can be obtained if the α 0 angle in Eq. (19) is replaced by its complementary (π/2 − α 0 ), therefore: where α o is given by relation (29) in which F o has been replaced by F oc . Relation (32b) is very close to (32a) as shown in Fig. 10c. This last relation tends to show that the energy along the radial direction comes from a direction appreciably perpendicular to those of the beam propagation. Indeed, as described by Polesana et al. [33], multiphoton absorption reshapes the core of the beam and tends to generate an energy flow from the periphery to the central part of the filament. As a consequence, the laser energy surrounding the core of the filament plays the principal role of refilling the intense central part of the filament [33]. For very high fluencies the two relations (28b) and (32a) can be simplified: This model tends to show that the geometry of the second crater is determined by those of the first one. The first crater combined with the formation of dense plasma behaves as a kind of micro-axicon.
Moreover, according to the work of Brzobohaty et al. [39], the shift δ observed on the apparition of the secondary crater (Fig. 9) and the damped oscillations with a period close to λ observed on the diameter of this crater (Fig. 14) could be explained by the imperfections (rounding) of the tip of the virtual micro-axicon created by this main crater. A weak unsteady Bessel filamentation regime also could explain the observed oscillations [33].
Considering the dependence of the diameter L 3 and of the depth h 2 on the z' position of the beam waist, as for Eqs. (4a,b) for a Gaussian beam, Eq. (16) has been rewritten as : Following the same way as for the determination of bL 3max > and bh 2max >, Eqs. (35) and (36) have been obtained: and as shown in Fig. 12, setting z Ã ¼ z 0 −δ ΔZ 2 −δ , the following relation is proposed : Thus, for z′ = δ, z* = 0 and g(z′,δ) = 0, for z 0 ¼ ΔZ 2 , z* = 1, g(z′,δ) = 0 and for z 0 ¼ 1 2 ΔZ 2 þ δ À Á À Á , zÃ ¼ 1 2 and g(z′, δ) = 1. Moreover z 0 ¼ 1− 2δ ΔZ À Á z Ã − 2δ ΔZ and for the higher fluencies δ → 0, z′ = z* and therefore f(z′) = f(z*). Taking into account the previous relations, Eqs. (35) and (36) are rewritten as follows: Relation (38) has been drawn in Fig. 12. Considering the experimental scattering, this relation is considered as in a fairly good agreement with the experimental results. However it is difficult to evaluate the pertinence of relation (39) due to the lack of experimental results. Only some measurements have been carried out on the z′ range for E p = 22 nJ and NA = 0.4 (Fig. 16). The evaluation of Eq.(39) with δ = 0.75 μm and b h 2max > = 5500 nm (see Fig. 13a) has been drawn in Fig. 16 and a fairly acceptable agreement has been obtained.

Conclusion
A detailed study of the evolution of the morphology of nanocraters drilled in borosilicate by a single laser shot near the ablation threshold has been performed by scanning electron microscopy, atomic force microscopy and scanning electron microscopy after focused ion beam sectioning. The influence of the numerical aperture (NA = 0.4 and 0.8), the pulse energy (16 b E p b 600 nJ) and the position with respect to the specimen surface of the waist of the Gaussian beam were systematically investigated. It has been shown that the nano-craters size is not restricted by the diffraction limit but imposed by the laser pulse stability and the material properties. Depending on the position of the waist, outside or inside the material, one or two characteristic craters with very different morphologies have been observed. The dimensions of these two kinds of craters have been studied in details. Different relationships have been proposed for the evolution of the depths and the various diameters of the main and second craters as functions of the pulse energy, the numerical aperture and the position of specimen surface in the focal region. The aspect ratio of the primary craters is small, independent on NA and lower than one. On the contrary those of the secondary crater varies as 1/NA 2 and is higher than 16 for NA = 0.4 and for pulses energies greater than twice the energy threshold. It has been assumed that these secondary craters are due to a spontaneous reshaping of the incoming Gaussian pulse into a Gaussian-Bessel pulse. Note that the geometry of the second holes seems to be linked to the one of the main crater, as proposed in the developed model.