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

Benoit Delobelle, Roland Salut, François Courvoisier, Patrick Delobelle

► To cite this version:

Benoit Delobelle, Roland Salut, François Courvoisier, Patrick Delobelle. A detailed study through the focal region of near-threshold single-shot femtosecond laser ablation nano-holes in borosilicate glass. Optics Communications, 2011, 284 (24), pp.5746-5757. 10.1016/j.optcom.2011.08.037. hal-02300197

HAL Id: hal-02300197

https://hal.science/hal-02300197

Submitted on 3 Feb 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A detailed study through the focal region of near-threshold single-shot femtosecond laser ablation nano-holes in borosilicate glass

B. Delobelle<sup>a,b,1</sup>, R. Salut<sup>c</sup>, F. Courvoisier<sup>b</sup>, P. Delobelle<sup>a,*</sup>

<sup>a</sup> Department DMA, FEMTO-ST (UMR CNRS 6174), University of Franche-Comté, 24 chemin de l'Epitaphe, 25000 Besançon, France

<sup>b</sup> Department DOPMD, FEMTO-ST (UMR CNRS 6174), University of Franche-Comté, 16 route de Gray, 25030 Besançon Cedex, France

<sup>c</sup> Department MN2S, FEMTO-ST (UMR CNRS 6174), University of Franche-Comté, 32 avenue de l'Observatoire, 25044 Besançon, France

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 influence of the numerical aperture (NA = 0.4 and 0.8), the pulse energy (16 nJ <$E_p$<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 profile 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.

## 1. 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 [1–30]. 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–9,22–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, carrier-carrier 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 [1–21]. 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–9].

The physical modeling of such complex processes is extremely difficult and has generated many theoretical developments from ultrafast dynamics of free electron generation to non-equilibrium thermodynamics involving large-scale molecular dynamics simulations [9–16].

However, according to the author's knowledge, the full surface profiling of nano-holes drilled by femtosecond laser single shot, in the focal region and in the vicinity of the ablation threshold has never been performed with a systematic approach in dielectrics. Indeed, the 3D morphology of the single-shot craters is very difficult to predict quantitatively due to the competition between the different physical processes; multiphotons [1–3], tunnel and avalanche ionization effects [1–3,9], energy diffusion effects [10,14], geometrical confinement effects [16], thermocapillary forces [17] and their dependence with the characteristics of the laser device; wavelength [4,9], pulse duration [6,11,14], numerical aperture [1,2], beam polarisation [4,18–20] and of the material; bandgap energy [1,4,8], refraction index, birefringence [21].

<footnote>* Corresponding author. Tel.: +33 3 81 66 60 13; fax: +33 3 81 66 67 00.
E-mail address: patrick.delobelle@univ-fcomte.fr (P. Delobelle).</footnote>

<footnote n="1">¹ Now is at: MAHYTEC, 210, Avenue de Verdun, 39100, Dôle, France.</footnote>
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 ambient 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–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.

## 2. Experiments

### 2.1. Experimental setup

The amplified laser source (Spitfire ProV, Spectra-physics) emits 100 fs laser pulses with a central wavelength $\lambda$ 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 infinitely-corrected microscope objectives were used: $\times$20 with 0.4 numerical aperture (NA) and $\times$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 $\times$20 MO and $\times$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 $\tau$ with a GRENOUILLE 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 $\tau = 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 \ \mu\text{W} < P < 3 \ \text{W}$ corresponding to a pulse energy $E_p$ in the range $2 \ \text{nJ} < E_p < 6 \ 10^5 \ \text{nJ}$. Single-shot illuminations of the sample were performed under atmospheric conditions, for different pulse energies in the range of $16 \ \text{nJ}$ to $600 \ \text{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.

### 2.2. 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 \ \text{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 nano-holes 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 $\gamma$ ($\gamma=35°$, $45°$ or $50°$) in the SEM, the projected images of the

<figure><img src="image_2.png" /><figcaption>Fig. 1. Description of the experimental procedure and definition of the used parameters $z, z', z^*, \Delta z$, describing the Gaussian beam in relation to the surface of the specimen.</figcaption></figure>
sectioned holes have been observed. Thus, the true heights h are given by h = h measured/sin η.

### 3. Results

#### 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<sub>p</sub> 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<sub>p</sub> and a decreasing function of NA. The value of the ablation threshold energy E<sub>p0</sub> 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 :

$$ \omega(\xi) = \omega_0 \sqrt{1 + \left(\frac{\xi}{z_0}\right)^2} \quad \text{with} \quad z_0 = \frac{\pi\omega_0^2}{\lambda M^2} \quad \text{and} \quad \omega_0 = \frac{\alpha\lambda}{NA}, \quad (1) $$

where z<sub>0</sub> is the Rayleigh range and ω<sub>0</sub> 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<1.1. α is a parameter which depends on the ratio of the Gaussian beam diameter d at 1/e<sup>2</sup> of the peak intensity of the laser and to the aperture diameter D of the lens. At present d/

<figure><img src="image_3.png" /><figcaption>Fig. 2. a) Evolution as a function of the pulse energy E<sub>p</sub> and the numerical aperture NA of the length ΔZ of the focal region where the laser-surface sample interaction is detected by direct imaging. b) Representation of NA<sup>2</sup>ΔZ as a function of NA<sup>2</sup>E<sub>p</sub> according to Eq. (3). Determination of the two parameters a<sub>1</sub> and a<sub>2</sub>.</figcaption></figure>

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<sub>o</sub>, thanks to the relation (1) it is easy to show that ΔZ measured in the air is given by the following relation:

$$ \Delta Z = 2z_0 \sqrt{\frac{F}{F_0}} - 1 = \frac{2\pi\alpha^2\lambda}{M^2NA^2} \sqrt{\frac{E_pNA^2}{\alpha^2\pi\lambda^2F_0}} - 1. \quad (2) $$

For F/F<sub>0</sub> = 2, ΔZ/2 = z<sub>0</sub>.

The relation (2) is recasted as:

$$ NA^2 \Delta Z = a_2 \sqrt{\frac{E_p NA^2}{a_1}} - 1 \quad \text{with} \quad a_1 = \alpha^2 \pi \lambda^2 F_0 \quad \text{and} \quad a_2 = \frac{2\pi\alpha^2\lambda}{M^2}. \quad (3) $$

The parameters a<sub>1</sub> and a<sub>2</sub> have been determined thanks to the Fig. 2b where NA<sup>2</sup> ΔZ is plotted as a function of NA<sup>2</sup>E<sub>p</sub> (Eq. 3). The continuous curve in Fig. 2b corresponds to a<sub>1</sub> = 11 nJ and a<sub>2</sub> = 1.21 µm, thus F<sub>o</sub>M<sup>2</sup> = 2.3 J cm<sup>-2</sup> and α/M = 0.491. As 1≤M<1.1, the values of F<sub>o</sub> and α should be in the range of 1.9-2.3 J cm<sup>-2</sup> and 0.49-0.54, respectively. If M=1, F<sub>o</sub> = 2.3 J cm<sup>-2</sup> and α = 0.491. The F<sub>o</sub> = 2.3 J cm<sup>-2</sup> value is in accordance with previously published works on fused silica; F<sub>o</sub> ≅ 2 J cm<sup>-2</sup> for λ = 825 nm and τ = 120 fs [25], F<sub>o</sub> ≅ 2.5 J cm<sup>-2</sup> for λ = 800 nm and τ = 120 fs [26], 2<F<sub>o</sub><3 J cm<sup>-2</sup> for λ = 798 nm and τ = 120 fs [27]. The values of F<sub>o</sub> and α reported in a previous analysis [24], F<sub>o</sub> = 2.75 J cm<sup>-2</sup> and α = 0.417 correspond to M = 0.88, thus are inconsistent with the inequality M≥1. Note that the value of F<sub>o</sub> depends on the wavelength λ [4,9] and on the pulse duration τ [6,11,14,25-27].

#### 3.2. Description of the general aspects of the nano-holes

As a function of the pulse energy E<sub>p</sub>, the NA and the z position in the focal region (0<z<Δ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]. Figs. 3a and b show two examples of such geometries with a single nano-crater (E<sub>p</sub> = 134 nJ, NA = 0.4, z = 3.75 µm for ΔZ = 8.75 µm) and with two craters (E<sub>p</sub> = 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 trapezoidal 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<sub>1</sub>, L<sub>2</sub>, h<sub>1</sub> for the main crater and L<sub>3</sub>, h<sub>2</sub> for the second crater.

#### 3.3. Study of the main crater dimensions L<sub>1</sub>, L<sub>2</sub> and h<sub>1</sub>

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.
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. a) SEM imaging of a crater after FIB drilling. E<sub>p</sub> = 134 nJ, NA = 0.4 and z = 3.75 µm < ΔZ/2 = 4.4 µm. Only one crater is visible (primary crater). b) SEM imaging of the two characteristic craters after FIB drilling E<sub>p</sub> = 108 nJ, NA = 0.4 and z = 6.5 µm > ΔZ/2 = 3.75 µm. As z > ΔZ/2 a long thin secondary crater appears. c) Diagram of the crater morphology. Definition of the different dimensions.</figcaption></figure>

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<sub>1</sub> and h<sub>1</sub> 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 <L<sub>1max</sub>> and <h<sub>1max</sub>> as a function of the pulse energy E<sub>p</sub> 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):

$$I(x, \xi, z') = I_0 \exp\left(-\frac{\xi}{A(\alpha^*, q, R, z')}\right) \exp\left(-\frac{x}{\omega(z')}\right)^2, \quad (4a)$$

with

$$\omega(z') = \beta_L \omega_0 \sqrt{f(z')} \quad \text{and} \quad A(\alpha^*, q, R, z') = A(\alpha^*, q, R) f(z'). \quad (4b)$$

<figure><img src="image_7.png" /><figcaption>Fig. 4. a) Evolution of the maximum mean values <L<sub>1max</sub>> of the width of the primary craters as a function of the pulse energy for NA = 0.4 and 0.8. b) Representation of NA < L<sub>1max</sub>> as a function of NA<sup>2</sup>E<sub>p</sub> according to the model (Eq. 7a). Determination of the β<sub>L</sub>M parameter.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 5.** a) Evolution of the maximum mean values of $<h_{1max}>$ of the depth of the primary craters as a function of the pulse energy for NA=0.4 and 0.8. b) Representation of NA<$h_{1max}>$ as a function of NA$^2$E<sub>p</sub> according to the model (Eq. 9). Determination of the β<sub>h</sub>M parameter.</figcaption></figure>

In a previous paper [24], $f(z')$ has been identified as:

$$f(z') = \sqrt{1 - \left(\frac{2z'}{\Delta Z}\right)^2}, \quad \text{thus } f(0) = 1 \text{ and } f\left(\pm\frac{\Delta Z}{2}\right) = 0. \quad (4c)$$

In these equations, $\xi$ is the distance along the propagation direction, $x$ is the lateral distance from the optical axis, $\omega$ the width of the beam function of the position $z'$ of the waist $\omega_o$ (Fig. 1). 1/A ($\alpha^*, q, R$...) is an absorption coefficient which depends on the linear absorption coefficient $\alpha^*$ of the material, the number $q$ of required photons to initiate ablation ($q=6$ for SiO<sub>2</sub>), the reflectivity $R$ of the surface, etc. [29]. $\beta_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:

$$L_1 = 2\beta_L \omega_o \sqrt{f(z')} \ln\left(\frac{F}{F_0}\right) \quad (5a)$$

and

$$h_1 = A(\alpha^*, q, R) f(z') \ln\left(\frac{F}{F_0}\right). \quad (5b)$$

As $f(0) = 1$, thus:

$$<L_{1max}> = 2\beta_L \omega_o \sqrt{\ln\left(\frac{F}{F_0}\right)} \quad (6a)$$

and

$$<h_{1max}> = A(\alpha^*, q, R) \ln\left(\frac{F}{F_0}\right). \quad (6b)$$

$<L_{1max}>$ and $<h_{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:

$$<L_{1max}> = \frac{2\beta_L \alpha \lambda}{NA} \sqrt{\ln\left(\frac{E_p NA^2}{a_1}\right)} \quad \text{and} \quad <h_{1max}> = A(\alpha^*, q, R) \ln\left(\frac{E_p NA^2}{a_1}\right). \\ \text{with } a_1 = \alpha^2 \pi \lambda^2 F_0 = 11 \text{ nJ} \quad (7a, b)$$

The relation (7a) has been identified in Fig. 4b where NA<$L_{1max}$> has been plotted as a function of NA$^2$E<sub>p</sub>. The fit of the experimental points, continuous curve in Fig. 4b, gives $2\beta_L\alpha\lambda = 535$ nm. As $\alpha/M = 0.491$ (paragraph 3.1), $\beta_L M = 0.68$. So, $\omega(0)$ in Eq. (4b) is equal to 68% of the waist $\omega_o$ if $M=1$.

For the parameter $<h_{1max}>$, the same kind of representation as for $<L_{1max}>$ in Fig. 4b, i.e. NA<$h_{1max}$> versus NA$^2$E<sub>p</sub>, has been plotted in Fig. 5b. Then, the experimental points have been fitted by the relation $7b$ in which $A(\alpha^*, q, R) = 230$ nm/NA; continuous curve in Fig. 5b. This last relation seems to show that, for a given wavelength:

$$A(\alpha^*, q, R, \dots) = \beta_h \omega_o \quad (8)$$

with $\beta_h M = 0.59$ for the studied material. Thus, the coefficient $A(\alpha^*, q, R, \dots)$ depends on NA and is equal to 59% of the waist $\omega_o$ if $M=1$. The relation (6b) is rewritten as follows:

$$<h_{1max}> = \beta_h \omega_o \ln\left(\frac{F}{F_0}\right) = \beta_h \frac{\alpha\lambda}{NA} \ln\left(\frac{E_p NA^2}{a_1}\right). \quad (9)$$

Note that the aspect ratio $h_1/L_1$ of the main crater is small, independent of NA and, from Eqs. (4b) (6a) and (9), is equal to:

$$\frac{h_1}{L_1} = \frac{\beta_h}{2\beta_L} \sqrt{f(z')} \ln\left(\frac{F}{F_0}\right), \quad \text{thus} \quad \frac{<h_{1max}>}{<L_{1max}>} = \frac{\beta_h}{2\beta_L} \sqrt{\ln\left(\frac{F}{F_0}\right)} \quad (10) \\ = 0.43 \sqrt{\ln\left(\frac{F}{F_0}\right)}.$$

As an example, for $F/F_0 = 1.5$, $<h_{1max}>/<L_{1max}> = 0.28$ and for $F/F_0 = 10$, $<h_{1max}>/<L_{1max}> = 0.66$.

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 $<L_{1max}>$, whatever the pulse energy E<sub>p</sub>, the NA and the position $z'$ of the waist, a linear relationship has appreciably been found between these two parameters:

$$L_2 = g(L_1) \approx \alpha_{21} L_1 - \beta_{21} \quad (11)$$

with: $\alpha_{21} \approx 0.85 \pm 0.05$ and $\beta_{21} \approx 100 \pm 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=\beta_{21}/\alpha_{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
<figure><img src="image_10.png" /><figcaption>Fig. 6. Dependence of the minimum diameters L₂ and <L₂max> of the primary craters as a function of the upper diameters L₁ and <L₁max>. Plot of the linear dependency according to the relation (11).</figcaption></figure>

L₂ in the focal region can be written as a function of ΔZ, F, NA, z' and the characteristics of the laser beam:

$$
\begin{align}
L_1 &= 2\beta_L \omega_0 \sqrt{f(z')} \ln \left[ 1 + \left( \frac{\Delta Z}{2z_0} \right)^2 \right] \nonumber \\
h_1 &= \beta_h \omega_0 f(z') \ln \left[ 1 + \left( \frac{\Delta Z}{2z_0} \right)^2 \right] \quad \text{with} \nonumber \\
f(z') &= \sqrt{1 - \left( \frac{2z'}{\Delta Z} \right)^2} \quad , \quad z' = z - \frac{\Delta Z}{2} \tag{12} \\
\Delta Z &= 2z_0 \sqrt{\frac{F}{F_0} - 1} \quad , \quad z_0 = \frac{\alpha^2 \pi \lambda}{M^2 N A^2} \quad \text{and} \quad \omega_0 = \frac{\alpha \lambda}{N A} \quad \text{with} \nonumber \\
\beta_L &= 0.68 \, \text{M} \quad , \quad \beta_h = 0.59 \, \text{M} \quad \alpha = 0.491 \, \text{M} \quad F_0 = 2.3 \, \text{M}^2 \, (\text{Jcm}^{-2}) \nonumber \\
L_2 &= \alpha_{21} \quad L_1 - \beta_{21} \quad , \quad \alpha_{21} \approx 0.85 \quad \text{and} \quad \beta_{21} \approx 100 \, \text{nm}. \nonumber
\end{align}
$$

3.4. Study of the second crater dimensions: L₃ and h₂

The long, thin conical profile shown in Fig. 3b suggests a self-focus-
ing of the laser pulse [9,30] or a nonlinear Gaussian-to-Bessel beam
transformation [31–33] as the beam propagates into the nonlinear sam-
ple with a sub-diffraction sized diameter. The first possibility, self-
channeling of the beam, is supported by the balance among diffraction,
Kerr-induced self-focusing and plasma-induced defocusing. The spon-
taneous reshaping of the beam is based on a dynamic spatial replenish-
ment 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

<figure><img src="image_11.png" /><figcaption>Fig. 7. Example of nearly conical crater for very small pulse energy: E<sub>p</sub> = 70.5 nJ, NA = 0.4, ΔZ/2 ~ 1.63 µm.</figcaption></figure>

the incoming Gaussian pulse into a Gauss-Bessel pulse and leads to-
wards 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 absorp-
tion and plasma generation over long distances [33]. These two possibil-
ities will be analyzed further, after the presentation of the experimental
results, i.e. condition of apparition and geometries of these craters.

3.4.1. 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₃ 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ₚ = 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₃ as a function of the z' distance in the focal region for var-
ious pulses energies and the two NA. It is thus shown that the dis-
tance z' = δ at which the first second crater occurs appears to
depend on the NA and to be a decreasing function of the energy Eₚ.
On the contrary, the maximum mean value <L₃max> of L₃ is a slowly
increasing function of Eₚ: 150<(<L₃max>)<300 nm (Fig. 10a). In this
figure, the L₃ value reported by White et al. [28] for NA = 0.85 and
Eₚ = 120 nJ has been drawn. To quantitatively compare our results
to those of the literature, the L₃ 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ₚ/Eₚ₀ in Fig. 10b, where Eₚ₀ is
the experimental value of the threshold energy reported in these pa-
pers. For Eₚ/Eₚ₀<3.3, all the experimental determinations are close,
100<L₃<300 nm, which shows that the L₃ diameter is appreciably in-
dependent on the numerical aperture (NA = 0.4, 0.8, 0.85) and to the
half conical angle αB of the Bessel beam (αB = 17° and 11° in the bo-
rosilicate glass [35]), contrary to Eₚ₀ which varies as NA². The NA
values corresponding to these two αB angles, NA = sinαB/n₀ with
n₀ = 1.5, are equal to NA = 0.19 and 0.127 for αB = 17° and 11°, re-
spectively. Now, if NA²Eₚ₀ = 11 nJ, as previously reported,
Eₚ₀ = 304 nJ and 682 nJ, which is in good agreement with the exper-
imental values, Eₚ₀~320 nJ and 610 nJ, respectively [35]. For Eₚ/
Eₚ₀ > 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]

<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Example of SEM imaging with two craters: E<sub>p</sub> = 161 nJ, NA = 0.4, and z = 8.5 µm > ΔZ/2 = 5.25 µm. Determination of the diameter L<sub>3</sub> of the secondary craters.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Evolution of the diameter L₃ of the secondary craters as a function of the z' distance in the focal region for different pulse energies and NA = 0.4 and 0.8. Determination of the maximum mean values <L₃max>.</figcaption></figure>

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<sub>pc</sub>, slightly greater than those of the threshold E<sub>po</sub> of the main craters. E<sub>pc</sub> = 80 nJ for NA = 0.4 and E<sub>pc</sub> = 20 nJ for NA = 0.8 (Fig. 10a). Indeed, as previously shown, NA<sup>2</sup> E<sub>po</sub> = a<sub>1</sub> = 11 nJ for the principal craters and E<sub>pc</sub>NA<sup>2</sup> = a<sub>3</sub> = 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 <L<sub>3max</sub>> as a function of NA<sup>2</sup>E<sub>p</sub>. It is very interesting to observe that when L<sub>2</sub> ~ L<sub>3</sub>, 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<sub>p</sub> = 31 nJ, NA = 0.8 thus L<sub>2</sub> = 330 nm (Eqs. 7a,b and 11) and L<sub>3</sub> = 210 nm (Fig. 10a).

From Fig. 9 and as previously mentioned, δ and ΔZ/2 are decreasing and increasing functions of E<sub>p</sub>, respectively, so when δ = ΔZ/2 no second crater has been detected. This condition occurs at about NA<sup>2</sup>E<sub>pc</sub> = 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 <L<sub>3max</sub>> for the Y-axis, the new coordinates in Fig. 12 are z* = (z - (δ + ΔZ/2))/(ΔZ/2 - δ) (Fig. 1) and L<sub>3</sub>/<L<sub>3max</sub>>. When z = δ + ΔZ/2, z* = 0 and when z = ΔZ, then z* = 1. In this figure, whatever the pulse energy E<sub>p</sub> 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<sub>2</sub> 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 <h<sub>2max</sub>> 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<sub>3</sub> is maximum (L<sub>3</sub>/<L<sub>3max</sub>> = 1, Fig. 12). The results are presented in Fig. 13a for the two NA. It is obvious that <h<sub>2max</sub>> greatly depends on the NA and very deep holes ( <h<sub>2max</sub>> ≥ 5 µm) have been observed for NA = 0.4 and E<sub>p</sub> > 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 < h<sub>2</sub> < 2.7 µm for 0.8 < E<sub>p</sub> < 1.2 µJ and NA = 0.85. Fig. 13b where NA<sup>2</sup> < h<sub>2max</sub>> has been plotted as the function of the pulse fluence (F α NA<sup>2</sup>E<sub>p</sub>) shows that the depth <h<sub>2max</sub>> appreciably varies as NA<sup>2</sup>. Contrary to the main craters (Eq. 10), the morphology of the secondary craters presents a high aspect ratio, i.e. : <h<sub>2max</sub>>/<L<sub>3max</sub>> is in the range 4 to 5 and 16 to 19 for NA = 0.8 and 0.4, respectively, when 20 < NA<sup>2</sup>E<sub>p</sub> < 50 nJ.

Before concluding on the experimental results some observations and remarks have to be done:

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 10.** a) Evolution of the maximum mean values <L<sub>3max</sub>> of the diameter of the secondary craters as a function of the pulse energy for NA = 0.4 and 0.8. b) Evolution of the diameter L<sub>3</sub> reported in the literature [28,35] and in the present study as a function of the fluence ratio E<sub>p</sub>/E<sub>po</sub>, where E<sub>po</sub> is the experimental value of the threshold. c) Representation of <L<sub>3max</sub>> as a function of NA<sup>2</sup>E<sub>p</sub> according to the model (Eq. 32a,b).</figcaption></figure>

- 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<sub>2</sub> 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<sub>p</sub> = 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<sub>p</sub> = 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
<figure><img src="image_17.png" /><figcaption>**Fig. 11.** Example of morphology where no discontinuity between the two craters is observed; $L_3-L_2$. $E_p = 31.6$ nJ, NA = 0.8, $\Delta Z$ ~ 1.75 µm, $L_3-L_2$ 250~300 nm.</figcaption></figure>

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 $<h_{3max}>$ values of the length of the compressed zone ($<h_{2max}>$ minus the distance of the top of the compressed zone from the surface) have been reported in Fig. 13a. Thus, for $E_p$<120 nJ ($<h_{2max}>$<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 $<h_{2max}>$ <2800 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^2E_p$ = 10.5 nJ which corresponds to $NA^2E_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^2E_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^2E_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

<figure><img src="image_18.png" /><figcaption>**Fig. 12.** Evolution of the normalized diameter $L_3/L_{3max}$ as a function of the new coordinate $z^*$ of the focal position, $z^* = (2(z - \delta) - \Delta Z)/(\Delta Z - 2\delta)$. Representation of the locus given by Eqs. (38) and (20).</figcaption></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 13.** a) Evolution of the maximum mean values $<h_{2max}>$ of the depth of the secondary craters as a function of the pulse energy for NA = 0.4 and 0.8. The depth $<h_{3max}>$ of the compressed zone has also been represented for NA = 0.4. b) Representation of $NA^2\cdot<h_{2max}>$ as a function of $NA^2E_p$ according to relation (28a). The prediction of the modeling (Eq. (28a)) with $\beta_{hs}$ = 0.77 has been drawn.</figcaption></figure>

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 $<h_{Tmax}>$ of the craters and the length $\Delta Z$ previously determined. Indeed, in the studied

<figure><img src="image_21.png" /><figcaption>**Fig. 14.** Example of SEM imaging showing the shock-affected zone at the bottom of the hole and the oscillation of diameter at the top of the secondary crater. $E_p$ = 130 nJ, NA = 0.4, $z = 5.75$ µm > $\Delta Z/2$ = 4.37 µm.</figcaption></figure>
range of pulse energy, in a first approximation, it has been observed that:

$$ \langle h_{T \max} \rangle = \langle h_{1 \max} \rangle + \langle h_{2 \max} \rangle \approx \Delta Z/2 \quad (13) $$

The relations $\Delta Z/2$ and $\langle h_{1\max} \rangle$ 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 \langle h_{2\max} \rangle$ can be written as:

$$ NA^2 \langle h_{2 \max} \rangle = \alpha \lambda \left[ \frac{\pi \alpha}{M^2} \sqrt{\frac{E_p NA^2}{a_1}} - 1 - \beta_h NA \ln \left( \frac{E_p NA^2}{a_1} \right) \right]. \quad (14) $$

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 $\delta$ tends to zero, the empirical relationship $\langle h_{Tmax} \rangle \sim \Delta Z/2$ is verified and can be utilized for practical applications; $\Delta Z$ is very simple to determine contrary to the total depth of the craters.

### 3.4.2. 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 $P_c$ given by [9,30]:

$$ P_c = \frac{\alpha^2 \pi \lambda^2}{8M^2 n_o n_2}. \quad (15) $$

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 \times 10^{-16}$ cm$^2$ W$^{-1}$ for the silica, $\alpha/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, $P_c$ 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]:

$$ I(x, \xi) = 4\pi^2 I_o \frac{\omega_o \sin \alpha_o}{\lambda} \left(\frac{\xi}{\xi_{\max}}\right) \exp\left(-\left(\frac{\sqrt{2}\xi}{\xi_{\max}}\right)^2\right) J_o^2\left(\frac{2\pi x \sin \alpha_o}{\lambda}\right) \quad (16) $$

where $\xi_{\max} = \frac{\omega_o}{\text{tg}\alpha_o}$

In this equation $x$ denotes the radial distance from the optical axis, $\xi$ the longitudinal position along the propagation direction and $\alpha_o$ the half cone angle of the resulting $k$ vectors after passing through an axis-con. $J_o$ is the zeroth-order Bessel function of the first kind and $\xi_{\max}$ the maximum range of the beam propagation. Note that $I(x, \xi)$ is maximum for $\xi_M = \xi_{\max}/2$.

Setting: $x_o = \frac{2.405}{2\pi} \frac{\lambda}{\sin \alpha_o}$, it is easy to show that:

$$ \frac{I(x, \xi_M)}{I(0, \xi_M)} = J_o^2 \left( 2.405 \frac{x}{x_o} \right). \quad (17) $$

$J_o^2(2.405 x/x_o) = 1$ for $x=0$ and is equal to zero for $x=x_o$. The relation (17) is not reversible. However, for $2.405 x/x_o < 1$ (first lobe of the beam) a good approximation of $J_o^2$ is given by:

$$ J_o^2 \left( 2.405 \frac{x}{x_o} \right) \approx 1 - \frac{1}{2} \left( 2.405 \frac{x}{x_o} \right)^2 + \varepsilon^2 \left( \frac{x}{x_o} \right). \quad (18) $$

Taking into account the relation (18) and the fluence threshold energy $F_{oc}$ for the second crater appearance, by setting $\langle L_{3\max} \rangle / 2 = \beta_{LS}$ $x$, the inversion of the relation (17) gives the evolution of the diameter $\langle L_{3\max} \rangle$ for an ideal axicon with infinite aperture as a function of the fluence:

$$ \langle L_{3 \max} \rangle = \beta_{LS} \frac{\sqrt{2\lambda}}{\pi \sin \alpha_o} \sqrt{1 - \frac{F_{oc}}{F}} = \beta_{LS} \frac{\sqrt{2\lambda}}{\pi \sin \alpha_o} \sqrt{1 - \frac{a_3}{NA^2 E_p}}. \quad (19) $$

In this relation $\beta_{LS}$ is a coefficient of the model and $a_3 = NA^2 E_{pc} = 13$ nJ.

Concerning the evaluation of the maximum depth $\langle h_{2\max} \rangle$, thanks to Eq. (16), it is easy to show that:

$$ \frac{I(x, \xi)}{I(x, \xi_m)} = \frac{I(0, \xi)}{I(0, \xi_m)} = \frac{2}{0.61} \left( \frac{\xi}{\xi_{\max}} \right) \exp\left(-\left(\frac{\sqrt{2}\xi}{\xi_{\max}}\right)^2\right). \quad (20) $$

As for relation (17), this expression is not reversible. However, as it will be shown further in Fig. 12, for $\xi/\xi_{\max} < 0.7$ a very good approximation of this expression is given by the following equation:

$$ \frac{I(x, \xi)}{I(x, \xi_m)} \approx 1 - \left( \frac{2\xi}{\xi_{\max}} - 1 \right)^2. \quad (21) $$
Thus, the two relations (20) and (21) give:

$$ \frac{I(x, \xi)}{I(x, \xi_m)} = 0 \text{ for } \xi = 0 \text{ and } \frac{I(x, \xi)}{I(x, \xi_m)} = 1 \text{ for } \xi = \xi_{\text{max}}/2. $$

Now, the inversion of the last relation (Eq. 21) allows the dependence of $<h_{2\text{max}}>$ = $\beta_{Ls}\xi$ for an ideal axicon as a function of the fluence:

$$ <h_{2\text{max}}> = \beta_{hs} \frac{\xi_{\text{max}}}{2} \left( 1 + \sqrt{1 - \frac{F_{oc}}{F}} \right) = \beta_{hs} \frac{\omega_o}{2tg\alpha_o} \left( 1 + \sqrt{1 - \frac{a_3}{NA^2 E_p}} \right). \quad (22) $$

$\beta_{hs}$ is a coefficient of the model. Note that when F is large then $<h_{2\text{max}}>$ = $\beta_{hs}\xi_{\text{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 \sim \alpha_2 L_1$). If $\theta_1$ is the half cone angle of the resulting k vectors of the generated Bessel beam, thus:

$$ \xi_{\text{max}} = \frac{L_1}{2tg\theta_1} \quad (23) $$

and moreover, if $\theta$ is the micro-axicon base angle, the Snell's law gives:

$$ tg\theta_1 = \left(1 - \frac{1}{n_o \cos\theta_1}\right) tg\theta, \quad (24) $$

where $n_o$ is the linear refractive index of the glass. The $\theta$ angle is assumed to be close to the beam divergence angle of the input Gaussian beam; i.e. of the order of:

$$ \theta = \frac{NA}{\pi\alpha} (\text{radian}). \quad (25) $$

Solving Eq. (24) with $\theta$ given by the relation (25) and considering the first order development of $tg\theta$, Eq. (24) can be recasted as Eq. (26) within a maximum error of 3% for the higher numerical aperture.

$$ tg\theta_1 = \left( \frac{n_o - 1}{n_o} \right) \frac{NA}{\pi\alpha} \left( 1 + \frac{1}{3} \left( \frac{NA}{\pi\alpha} \right)^2 \right). \quad (26) $$

The maximum value $<L_{3\text{max}}>$ of $L_3$ (Idem for $<h_{2\text{max}}>$) is obtained for $z^*=1/2$ (Fig. 12), thus for $z'=\Delta Z/4 + \delta/2$. As shown in Fig. 9, for the higher values of the pulse energy, the $\delta$ parameter can be neglected, thus $z' \sim \Delta Z/4$ ($z=3\Delta Z/4$) and $f(z') \sim \sqrt{3}/2$ (Eq.4c). On the contrary, for the lower values of $E_p$, $\delta$ tends towards $\Delta Z/2$ (Fig. 9) and for $E_p=E_{pc}$, $\delta=\Delta Z/2$. Therefore $z'=\Delta Z/2$ and $f(z')=0$. To take into account these two boundaries conditions without knowing the exact dependence of the $\delta$ parameter on the pulse energy, Eq. (5a) has been recasted as:

$$ L_1 = 2 \left(\frac{3}{4}\right)^{1/4} \beta_L \omega_o \sqrt{\ln \frac{F}{F_{oc}}} \quad \text{with} \quad \frac{F}{F_{oc}} = \frac{NA^2 E_p}{a_3} \quad \text{and} \quad a_3 = 13 \text{ nJ}. \quad (27) $$

Now, combining relations (22) (23) (26) and (27), the evolution of $<h_{2\text{max}}>$ with the pulse energy can be deduced:

$$ <h_{2\text{max}}> = \beta_{hs} \left(\frac{3}{4}\right)^{1/4} \frac{\beta_L}{2} \frac{\alpha^2 \pi \lambda}{NA^2} \frac{n_o}{n_o - 1} \frac{1}{\left[1 + \frac{1}{3} \left(\frac{NA}{\pi \alpha}\right)^2\right]} \sqrt{\ln \frac{F}{F_{oc}}} \left(1 + \sqrt{1 - \frac{F_{oc}}{F}}\right), \quad (28a) $$

of the order of:

$$ <h_{2\text{max}}> = \beta_{hs} \left(\frac{3}{4}\right)^{1/4} \frac{\beta_L}{2} \frac{n_o}{n_o - 1} z_o \sqrt{\ln \frac{F}{F_{oc}}} \left(1 + \sqrt{1 - \frac{F_{oc}}{F}}\right) \quad (28b) $$
with $\frac{F}{F_{oc}} = \frac{NA^2 E_p}{a_3}.$

According to Eq. (28a) $<h_{2\text{max}}>$ varies approximately as $1/NA^2$ as shown in Fig. 13b. The $\beta_{hs}$ parameter has been determined by fitting the experimental results in Fig. 13b, $NA^2<\hbar_{2\text{max}}>$ = $f(NA^2E_p)$, with the relation (28a) (solid lines). The best fit is obtained with $\beta_{hs}=0.77$, the other parameters being previously identified. $\beta_{hs}=0.91$ if the entrance diameter of the micro-axicon is equal to $L_2 \sim 0.85 L_1$. 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 $\alpha_o$ of the k vectors of an equivalent Bessel beam formed by a perfect axicon illuminated with a plane wave;

$$ \alpha_o = tg^{-1} \left[ \frac{(n_o - 1)tg \frac{NA}{\pi \alpha}}{\left(\frac{3}{4}\right)^{1/4} \beta_{hs} \beta_L n_o \sqrt{\ln \frac{F}{F_{oc}}}} \right]. \quad (29) $$

Note that the base angle of this axicon varies as $(\ln \frac{F}{F_{oc}})^{-1/2}$, thus $0 < \alpha_o < \frac{\pi}{2}$ when $1 < F/F_{oc} < \infty$. Indeed, the $\alpha_o$ angle cannot be greater than the opening angle $u$ of the objective ($u = \sin^{-1}NA$): i.e. $\alpha_o \le u$. From Eq. (29), the equality $\alpha_o=u$ allows to determine the condition for the apparition of the secondary craters:

$$ \frac{F_{oc}}{F_o} = \frac{NA^2 E_{pc}}{NA^2 E_{po}} = \frac{a_3}{a_1} \\ = \exp \left[ \frac{n_o - 1}{\left(\frac{3}{4}\right)^{1/4} \beta_{hs} \beta_L n_o} \frac{tg \frac{NA}{\pi \alpha}}{tg u} \right]^2 \approx \exp \left[ \frac{n_o - 1}{\left(\frac{3}{4}\right)^{1/4} \beta_{hs} \beta_L \pi \alpha} \frac{n_o}{n_o - 1} \right]^2. \quad (30) $$

For $NA=0.4$, the previous relation gives $NA^2E_{pc}=a_3=13.1$ nJ and for $NA=0.8$, $NA^2E_{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^2E_{pc}=a_3=13 \pm 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 $<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 = <h_{1\text{max}}> + <h_{2\text{max}}>$ 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<NA<0.85) or a true Bessel beam (Bhuyan et al. [35], 0.13<NA<0.19). The discrepancy obtained for NA=0.13 will

<figure><img src="image_22.png" /><figcaption>Fig. 15. Evolution of $NA^2h_T$ determined from the literature results [28,35] and in the present study as a function of the fluence ratio $E_p/E_{po}$. Predictions of the modelling: Eq. (31).</figcaption></figure>
be discussed below. The predictions of the modeling have been valued thanks to Eqs. (9) (<h<sub>1max</sub>>) and (28a) (<h<sub>2max</sub>>) which have been recasted as:

$$
\begin{equation}
\begin{aligned}
\mathrm{NA}^2 \mathrm{h}_T = \alpha \lambda & \left[ \beta_{\mathrm{hs}} \left(\frac{3}{4}\right)^{1/4} \beta_L \frac{\alpha \pi n_o}{(n_o-1)\left(1+\frac{1}{3}\left(\frac{\mathrm{NA}}{\pi \alpha}\right)^2\right)} \sqrt{\mathrm{Ln} \frac{\mathrm{E}_p}{\mathrm{E}_{\mathrm{p}_o}}} \left(1+\sqrt{1-\frac{\mathrm{E}_{\mathrm{p}_o}}{\mathrm{E}_p}}\right) \right. \\
& \left. + \beta_h \mathrm{NA} \mathrm{Ln} \frac{\mathrm{E}_p}{\mathrm{E}_{\mathrm{p}_o}} \right] .
\end{aligned}
\tag{31}
\end{equation}
$$

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<sub>T</sub> of the craters was fixed to 32 µm. So, solving Eq. (31) with h<sub>T</sub>=32 µm allows to determine the domain of validity E<sub>p</sub>/E<sub>po</sub>≤5.5, thus all the experimental points verify this inequality, but for α<sub>B</sub>=17° (NA=0.19), E<sub>p</sub>/E<sub>po</sub>≤5.5, thus all the experimental points verify this inequality, but for α<sub>B</sub>=11° (NA=0.13), E<sub>p</sub>/E<sub>po</sub>≤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 <L<sub>3max</sub>> corresponding to an infinite ideal axicon is given by the Eq. (19) and at present as α<sub>o</sub> (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 β<sub>Ls</sub>/sinα<sub>o</sub> = 1 in Eq. (19), (the minimum value of this ratio; β<sub>Ls</sub> = 1 and α<sub>o</sub> = π/2) thus:

$$
\langle L_{3\max} \rangle = \frac{\sqrt{2\lambda}}{\pi} \sqrt{1 - \frac{a_3}{NA^2 E_p}}. \qquad (32a)
$$

This relation perfectly match the experimental results as shown in Fig.10b-c (solid lines). The same result can be obtained if the α₀ angle in Eq. (19) is replaced by its complementary (π/2 - α₀), therefore:

$$
\langle L_{3\max} \rangle = \frac{\sqrt{2\lambda}}{\pi \cos \alpha_0} \sqrt{1 - \frac{a_3}{NA^2 E_p}}, \quad (32b)
$$

where α₀ is given by relation (29) in which F₀ has been replaced by F₀c. 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], multi-photon 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:

$$
\begin{align}
\langle L_{3\max} \rangle &= \frac{\sqrt{2\lambda}}{\pi} = 360 \text{ nm} \nonumber \\
\langle h_{2\max} \rangle &= \left(\frac{3}{4}\right)^{1/4} \beta_{\text{hs}} \beta_{\text{L}} \frac{n_o}{n_{o-1}} z_o \sqrt{\text{Ln} \frac{F}{F_{\text{oc}}}} \approx 1.46 z_o \sqrt{\text{Ln} \frac{F}{F_{\text{oc}}}}. \tag{33}
\end{align}
$$

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 diame-
ter 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<sub>3</sub> and of the depth
h<sub>2</sub> on the z' position of the beam waist, as for Eqs. (4a,b) for a Gauss-
ian beam, Eq.(16) has been rewritten as :

$$
I(x, \xi, z') = 4\pi^2 I_0 \frac{\omega_0 \sin \alpha_o}{\lambda} \left( \frac{\xi}{\xi_{\text{max}}(z', \delta)} \right) \exp\left( -\left( \frac{\sqrt{2}\xi}{\xi_{\text{max}}(z', \delta)} \right)^2 \right) J_0^2\left( 2.405 \frac{x}{x_o(z', \delta)} \right) \quad (34a)
$$

$$
\text{where } \xi_{\text{max}}(z', \delta) = \xi_{\text{max}}g(z', \delta) \text{ et } x_o(z', \delta) = x_o g(z', \delta). \quad (34\text{b})
$$

Following the same way as for the determination of <L<sub>3max</sub>> and
<hy>2max</sub>>, Eqs. (35) and (36) have been obtained:

$$
L_3 = \langle L_{3\max} \rangle g(z', \delta), \quad (35)
$$

$$
h_2 = \left(\frac{4}{3}\right)^{1/4} \langle h_{2\max} \rangle \sqrt{f(z')} g(z', \delta) \quad (36)
$$

and as shown in Fig. 12, setting z⁺ = z'/δ, the following relation is pro-
posed :

$$
g(z', \delta) = 1 - (2z^* - 1)^2 = 1 - \frac{2z' - \left(\frac{\Delta z}{2} + \delta\right)^2}{\frac{\Delta z}{2} - \delta}. \quad (37)
$$

Thus, for z' = δ, z* = 0 and g(z', δ) = 0, for z' = Δz/2, z* = 1, g(z', δ) = 0 and for z' = 1/2 ((Δz/2 + δ)), z* = 1/2 and g(z', δ) = 1. Moreover z' = (1 - 2Δz/2)z* - 2Δz/2 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:

$$
\langle L_{3\max} \rangle \approx 1 - (2z^* - 1)^2 = 4z^*(1 - z^*), \qquad (38)
$$

$$
\langle h_{2\max} \rangle \approx \left( \frac{4(1-z^{*2})}{3} \right)^{1/4} (1-(2z^{*2}-1)^2) \quad \text{with} \quad z^{*} = \frac{z' - \delta}{\frac{\Delta z}{2} - \delta}. \quad (39)
$$

Relation (38) has been drawn in Fig. 12. Considering the experi-
mental scattering, this relation is considered as in a fairly good agree-
ment 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<sub>p</sub> = 22 nJ and NA = 0.4 (Fig. 16). The evaluation of Eq.(39) with
δ = 0.75 µm and <h<sub>1max</sub>> = 5500 nm (see Fig. 13a) has been drawn
in Fig. 16 and a fairly acceptable agreement has been obtained.

In Fig. 12, relation (20) has been plotted with ξ/ξ<sub>max</sub> = z* and it is thus shown that the Eq. (21) also represented in this figure (idem Eq. (39) with z* = ξ/ξ<sub>max</sub>) constitutes a good approximation of the relationship (20) for ξ/ξ<sub>max</sub> ≤ 0.7.

Moreover, for the higher fluencies, <h<sub>2max</sub>> ~β<sub>hs</sub> ξ<sub>max</sub>, and if this
maximum corresponds to ξ<sub>max</sub>/2 ~z' ~ΔZ/4, thus <h<sub>2max</sub>> ~β<sub>hs</sub>
ΔZ/2 ~0.77 ΔZ/2. Taking the depth <h<sub>1max</sub>> of the primary crater (5
to 15% of ΔZ/2) into account could explain the empirical relationship
(13) :

$$
\langle h_T \rangle = \langle h_{1max} \rangle + \langle h_{2max} \rangle - \Delta Z/2.
$$

As a conclusion, Eqs.(32a) and (28a) combined with relations (38)
and (39) respectively, give the evolutions of the diameter L₃ and the
<figure><img src="image_23.png" /><figcaption>Fig. 16. Evolution of the length h₂ of the secondary crater as the function of the z* parameter. The evaluation of Eq. (39) has equally been drawn.</figcaption></figure>

depth h₂ of the secondary crater as a function of ΔZ/2, F (or the pulse energy Eₚ), NA, z' and the characteristic of the laser beam, λ, α/M, like Eq. (12) for the dimensions of the main crater, i.e.:

$$ 
\begin{align} 
L_3 &= \frac{\sqrt{2\lambda}}{\pi} g(z', \delta) \sqrt{1 - \frac{a_3}{NA^2 E_p}} \nonumber \\ 
h_2 &\approx \frac{\beta_{hs}}{2} \frac{\beta_L}{n_0 z_o} g(z', \delta) \sqrt{f(z') \ln \frac{NA^2 E_p}{a_3}} \left(1 + \sqrt{1 - \frac{a_3}{NA^2 E_p}}\right) \nonumber \\ 
&\quad \text{with } f(z') = \sqrt{1 - \left(\frac{2z'}{\Delta Z}\right)^2}, \quad g(z', \delta) = 4z^*(1 - z^*) \quad \text{and} \tag{40} \\ 
\frac{\Delta Z}{2} &= z_o \sqrt{\frac{NA^2 E_p}{a_1}} - 1, \quad z_o = \frac{\alpha^2 \pi \lambda}{M^2 NA^2}, \quad z' = z - \frac{\Delta Z}{2}, \nonumber \\ 
z^* &= \frac{z' - \delta}{\frac{\Delta Z}{2} - \delta} a_1 = NA^2 E_{p_o} = 11 \text{ nJ}, \quad a_3 = NA^2 E_{p_c} = 13 \text{ nJ}, \nonumber \\ 
\beta_{hs} &= 0.77, \beta_L = 0.68, \quad n_0 = 1.5. \nonumber 
\end{align} 
$$

## 4. Conclusion

A detailed study of the evolution of the morphology of nano-craters 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<Ep<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² 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.

## Acknowledgements

The authors thank the Conseil Regional de Franche-Comté, France, for funding, P.A. Lacourt and L. Furfaro for the setting of the femtosecond laser.

## References

[1] R.R. Gattass, E. Mazur, Nat. Photonics 2 (4) (2008) 219.

[2] W. Yang, P.G. Kazansky, Y.P. Svirko, Nat. Photonics 2 (2) (2008) 99.

[3] M. Lenzner, J. Kruger, S. Sartania, Z. Cheng, C. Spielmann, G. Marron, W. Kautek, F. Krausz, Phys. Rev. Lett. 80 (18) (1998) 4076.

[4] A.P. Joglekar, H. Liu, G.J. Spooner, E. Meyhofer, G. Mourou, A.J. Hunt, Appl. Phys. B-Lasers Opt. 77 (1) (2003) 25.

[5] A. Chimmalgi, T.Y. Choi, C.P. Grigoropoulos, C.P. Komvopoulos, Appl. Phys. Lett. 82 (8) (2003).

[6] P.P. Pronko, S.K. Dutta, J. Squier, J.V. Rudd, D. Du, G. Mourou, Opt. Comm. 114 (1-2) (1995) 106.

[7] K. Itho, W. Watanabe, RIKEN Review 50 (2003) 90.

[8] J. Koch, F. Korte, C. Fallnich, A. Ostendorf, B.N. Chichkov, Engineering 44 (5) (2005) 051103.

[9] S.S. Mao, F. Quere, S. Guizard, X. Mao, R.E. Russo, G. Petite, P. Martin, Appl. Phys. A. Mat. Sci. Process 79 (7) (2004) 1695.

[10] D. Ashkenasi, G. Muller, A. Rosenfeld, R. Stoian, I.V. Hertel, N.M. Bulgakova, E.E.B. Campbell, Appl. Phys. A-Mat. Sci. Process. 77 (2) (2003) 223.

[11] T.Q. Jia, Z.Z. Xu, R.X. Li, D.H. Feng, X.X. Li, C.F. Cheng, H.Y. Sun, N.S. Xu, H.Z. Wang, J. Appl. Phys. 95 (9) (2004) 5166.

[12] L.V. Zhigilei, Appl. Phys. A - Mat. Sci. Process 76 (3) (2003) 339.

[13] G.M. Petrov, J. Davis, J. Phys. B: At. Mol. Phys. 41 (2) (2008) 025601.

[14] S. Guizard, A. Semerok, J. Gaudin, A. Hashida, P. Martin, F. Quere, Appl. Surf. Sci. 186 (1-4) (2002) 364.

[15] L. Englert, B. Rethfeld, L. Haag, M. Wollenhaupt, C. Sarpe-Tudoran, T. Baumert, Opt. Express 15 (26) (2007) 17855.

[16] D. Bouilly, D. Perez, L.J. Lewis, Phys. Rev. B 76 (2007) 184119.

[17] T.D. Bennett, D.J. Krajrobich, C.P. Grigoropoulos, P. Baumgart, A.C. Tam, J. Heat Transfer 119 (3) (1997) 589.

[18] S. Nolte, C. Momma, G. Kamlage, A. Ostendorf, C. Fallnich, F. Von Alvensleben, H. Welling, Appl. Phys. A. Mat. Sci. Process. 68 (5) (1999) 563.

[19] C.S. Nielsen, P. Balling, J. Appl. Phys. 99 (9) (2006) 093101.

[20] D.J. Little, M. Ams, P. Dekker, G.D. Marshall, J.M. Dawes, M.J. Withford, Opt. Express 16 (24) (2008) 20029.

[21] P.G. Kazansky, W. Yang, E. Bricchi, J. Bovatsek, A. Arai, Y. Shimotsuma, K. Miura, K. Hirao, Appl. Phys. Lett. 90 (15) (2007) 151120.

[22] E.N. Glezer, E. Mazur, Appl. Phys. Lett. 71 (7) (1997) 882.

[23] I. Kazuyoshi, W. Watanabe, RIKEN Review 50 (2003) 90.

[24] B. Delobelle, F. Courvoisier, P. Delobelle, Opt. Lasers Eng. 48 (2010) 616.

[25] E.G. Gamaly, A.V. Rode, B. Luther-Davies, V.T. Tikhonchuk, Phys. Plasmas 9 (3) (2002) 949.

[26] S. Xu, J. Qiu, T. Jia, C. Li, H. Sun, Z. Xu, Opt. Comm. 274 (2007) 163.

[27] T.Q. Jia, Z.Z. Xu, X.X. Li, X. Li, B. Shuai, F.L. Zhao, Appl. Phys. Lett. 82 (24) (2003) 4382.

[28] Y.V. White, X. Li, Z. Sikorski, L.M. Davis, W. Hofmeister, Opt. Express 16 (19) (2008) 14411.

[29] D.J. Hwang, C.P. Grigoropoulos, T.Y. Choi, J. Appl. Phys. 99 (2006) 083101.

[30] L. Sudrie, M. Franco, B. Prade, A. Mysyrowicz, Opt. Comm. 191 (2001) 333.

[31] A. Dubietis, E. Gaizauska, G. Tamosauska, P. Trapani, Phys. Lett. 92 (2004) 253903.

[32] S. Juodkazis, V. Mieseikis, E. Gaizauskas, E. Vanagas, V. Jarutis, H. Misawa, Proc. SPIE 6053 (2005) 60530.

[33] P. Polesana, M. Franco, A. Couairon, D. Faccio, P. Di Trapani, Phys. Rev. A 77 (2008) 043814.

[34] J. Durnin, J.J. Miceli, J.H. Eberly, Phys. Rev. Lett. 58 (15) (1987) 1499.

[35] K.M. Bhuyan, F. Courvoisier, P.A. Lacourt, M. Jacquot, R. Salut, L. Furfaro, J. Dudley, Appl. Phys. Lett. 97 (2010) 081102.

[36] F. Courvoisier, P.A. Lacourt, M. Jacquot, M.K. Bhuyan, L. Furfaro, J. Dudley, Opt. Lett. 34 (20) (2009) 3163.

[37] S. Juodkazis, H. Misawa, T. Hashimoto, E.G. Gamaly, B. Luther-Davies, Appl. Phys. Lett. 88 (2006) 201909.

[38] G. Roy, R. Tremblay, Opt. Comm. 177 (2000) 297.

[39] D. Brzobohaty, T. Cizmar, P. Zemanek, Opt. Express 16 (17) (2008) 12688.