Invited Article: Filamentary deposition of laser energy in glasses with Bessel beams

We investigate the nonlinear absorption of laser energy in the bulk of transparent dielectrics for femtosecond and picosecond laser pulses focused by a conical lens. We highlight the inﬂuence of the pulse duration, laser pulse energy, and cone angle on laser energy absorption in transparent dielectrics. We provide a semi-analytical model allowing the calculation of maps for the density of nonlinear absorption of energy in BK7 and in SiO 2 as a function of the pulse duration and peak ﬂuence in the focal region. The comparison of the density of nonlinear absorption of energy with the available energy density determines optimal pulse durations and Bessel beam cone angles compatible with uniform laser energy deposition in the Bessel zone. The results reproduce quantitatively the transmission measurements for experiments in with picosecond pulses and suggest that the loss of propagation invariance and uniform laser energy deposition is responsible for a previously reported transition between different types of damage morphology in SiO


I. INTRODUCTION
By focusing a femtosecond or picosecond laser pulse inside a transparent dielectric, it is possible to deposit energy in a localized region corresponding to the domain where intensity exceeds the ionization threshold. 1,2 Energy initially deposited on the electron sub-system is subsequently transferred to the lattice and may induce an alteration of the refractive index or the formation of a void for tight focusing conditions. Laser-matter interaction in this regime has been thoroughly studied with Gaussian beams focused by standard devices (lenses and microscope objectives) leading to the confinement of energy in three spatial dimensions. 3 Recently, high aspect ratio laser structuring was demonstrated by focusing a femtosecond pulse in borosilicate glass and fused silica with a conical lens, both with multiple pulse illumination and in single shot, [4][5][6][7] opening the way to high-speed deep-drilling and cutting of glass samples as well as to nanoscale volume structuring with ultrashort laser pulses. [8][9][10][11][12] As sketched in Fig. 1, focusing a Gaussian beam by a conical lens generates a Bessel-Gauss beam over an elongated focal region; i.e., a light string with a narrow core and almost a Bessel beam profile J 0 (k 0 sin θr) is formed over a distance L B = 0 /tan θ called the Bessel zone. 13,14 Here θ denotes the cone half-angle of the Bessel beam in the medium, k 0 = 2πn 0 /λ 0 is the wavenumber in the medium of refractive index n 0 at the laser wavelength λ 0 , and 0 is the e −2 -radius of the initial Gaussian beam of intensity distribution I(r) = I 0 exp(−2r 2 / 2 0 ). In the linear propagation regime [ Fig. 1(b)], the fluence profile along z was calculated by Jarutis et al. 15 and takes the form F(z) = 4πF 0 (k 0 sin θ/ 0 )z exp(−2z 2 ), wherez = z/L B and F 0 denotes the peak fluence of the Gaussian beam. The radius r 0 of the Bessel filament can be approximated by the first zero of the Bessel beam profile, which satisfies k 0 sin θ r 0 = j 0,1 , where j 0,1 ∼ 2.405 is the first zero of the J 0 Bessel function.
Several nonlinear propagation regimes were reported for Bessel beams: 16,17 (i) The regime of small cone angles is featured by oscillations of the peak intensity of the Bessel beam along the propagation axis when the beam power exceeds a certain threshold. This regime exhibits the typical pulse dynamics occurring in ultrashort laser pulse filamentation. 18 For instance, Gaižauskas et al. 19 have observed discrete, equidistant damage spots in a borosilicate glass for a cone angle of 5 • . (ii) The regime of large cone angles is characterized by a uniform light channel over the Bessel zone with a quasi-constant intensity [see Fig. 1(c)], associated with a propagation invariant nonlinear Bessel beam profile. 16,20 In the present work, this regime was observed slightly above the transition for cone angles of ∼10 • in air (6.6 • in borosilicate glass); however, the transition also depends on the pulse duration and electron dynamics. The latter regime is of high interest for applications in precision micro-and nano-channel drilling of transparent dielectrics since the laser pulse forming the Bessel filament leaves in its wake a uniform plasma track, generated by the main Bessel lobe, that is, the main support for the absorption of laser energy. After this stage, the absorbed energy is transferred to the lattice and induces thermo-mechanical stress leading to a void nanochannel. While different mechanisms from microexplosion 21 to cavitation in the liquid phase 22 were proposed to be responsible for the formation of a void nanochannel, efficiency of this technique for nano-channel drilling has already been reported. 5 Figure 2 shows typical examples of 100 µm long tracks generated with a single laser shot in borosilicate glass (Corning 0211) by means of a Bessel beam produced by a spatial light modulator (SLM) out of the beams of a Ti-sapphire laser delivering 135 fs, 800 nm pulses with ∼15 µJ per pulse. For the same focusing conditions and pulse energy, refractive index changes induced by the Bessel beam are clearly visible and significantly different for pulses of picosecond duration. The shortest pulses, with potentially higher intensity in the focal region, lead to less visible tracks. The work we present is mainly theoretical, and the main question we address is the existence of optimal laser parameters for absorption of laser energy in the bulk of transparent materials, along the focal line of a Bessel beam.
In this aim, we present maps of the density of laser energy absorption in BK7 and in SiO 2 as a function of the pulse duration and the peak fluence in the Bessel zone. These maps are qualitatively representative of nonlinear absorption in any dielectric medium, though quantitatively specific to each medium as they depend on nonlinear absorption coefficients. Other experimental control parameters, including the cone angle of the Bessel beam and the pulse energy, determine the peak fluence in the Bessel zone. Given the damage tracks, nonlinear absorption appears to be significant only in the central core of the Bessel zone, even if the secondary lobes of the Bessel beam are important for nonlinear propagation and for energy balance considerations. Mapping the density of energy absorption can be performed semi-analytically owing to the properties of nonlinear propagation of Bessel beams. Propagation invariance in the Bessel zone suggests a simplified flat-top profile for the peak intensity in the Bessel zone [ Fig. 1(d)], with the underlying idea that the beam is intense only in the central core of the Bessel zone, of volume L B πr 2 0 , and that the localized absorption of laser energy can be evaluated by limiting integration volumes to the main lobe of the Bessel beam. A comparison of the nonlinear absorption density maps with the density of available energy sets the validity limit of this evaluation, interpreted as a loss of propagation invariance. Regions where propagation invariance holds are presented for BK7 and SiO 2 as a function of the pulse duration, cone angle, and pulse energy.
More generally, the method to evaluate nonlinear absorption maps applies in the context of laser-matter interaction and laser damage measurements of optical materials and components [23][24][25] and should help us develop well-controlled micromachining processes 26,27 for applications in the photonics industry.

II. PLASMA GENERATION AND NONLINEAR ABSORPTION OF ENERGY
Nonlinear propagation is usually described by a propagation equation governing the evolution of the laser electric field coupled with a model for the medium response in the form of a constitutive relation linking the nonlinear polarization and current entering in the propagation equation to the electric field. 28 Here, we are mainly interested in the propagation regime of Bessel beams with large cone angles, which is featured by a uniform absorption of laser energy due to the quasi-invariance of beam propagation. 17 Earlier studies have shown numerical evidence that for a sufficiently large cone angle, a Gaussian beam focused by an axicon reshapes into a propagation invariant nonlinear Bessel beam, 16,17,20 which preserves the amplitude of the linear Bessel beam that the axicon would generate in linear propagation. 29,30 The nonlinear Bessel beam is featured by a balance between the incoming conical energy flux toward the center and nonlinear losses in the intense part of the beam. The beam shape is similar to that of a linear Bessel beam with slightly compressed rings due to Kerr self-focusing and an attenuation of contrast due to the effect of nonlinear losses. 31 In contrast to linearly propagating finite energy Bessel-Gauss beams, the peak intensity of nonlinear Bessel beams along the propagation axis remains fairly constant over the Bessel zone. 17 We therefore evaluate the nonlinear absorption for a perfectly propagation invariant nonlinear Bessel beam over the entire Bessel zone [ Fig. 1(d)], where the beam exhibits a profile almost similar to a J 0 Bessel function, with a sufficient peak intensity to generate a plasma. Propagation invariance is associated with nonlinear energy losses which exactly balance the conical energy flux toward the high intensity core of the Bessel beam, where energy is absorbed nonlinearly, i.e., deposited to the medium. 20 We evaluate the absorption coefficient as the ratio of the density of energy absorption to the available laser energy density in the Bessel zone. When the assumption of perfect propagation invariance of the Bessel beam is no longer valid, the absorption coefficient takes values larger than unity, defining boundaries for pulse durations, cone angles, and pulse energies where nonlinear absorption maps cannot be used. These boundaries are calculated a posteriori since we focus on the constitutive relations of the medium in a first stage, i.e., the generation of electron density which is of particular importance for the absorption of laser energy with Bessel beams. In a second stage, we will determine the validity region of the nonlinear absorption maps.

A. Ionization
Intense laser fields can lead to ionization of the medium by multiphoton ionization and avalanche. Following a standard model, the evolution of the electron plasma density ρ is governed by a simple rate equation whereρ ≡ ρ/ρ 0 denotes the ionization degree and ρ 0 denotes the electron density in the valence band. The quantityĨ ≡ I(t)/I p denotes the time-dependent pulse intensity, normalized by I p , the peak intensity of the pulse. Ionization rates are expressed as ν av ≡ σI p /U g and ν mpi ≡ β K I K p /U K ρ 0 , where β K denotes the multiphoton absorption coefficient, U K = K ω 0 denotes the energy of K photons of frequency ω 0 necessary to promote an electron from the valence to the conduction band, with gap U g , and σ denotes the cross section for inverse Bremsstrahlung entering in the avalanche ionization rate. The last term in (1) denotes the recombination law, with either a quadratic dependence of the rate upon electron density (Model Mq) or a linear dependence (Model Ml) to model processes with exponential decay (see the second column in Table I). The model Mq describes electron decay at the rate ν rec due to bimolecular recombination with positive ions, 32 whereas the model Ml describes trapping of electrons at the rate ν r . 33,34 Both processes are slow, and the results do only moderately depend on the specific choice for the recombination model.

B. Absorption of energy
Two physical effects are responsible for the absorption of laser energy by the medium: (i) Multiphoton absorption corresponds to the energy necessary to initiate the generation of the electron plasma by the intense field. (ii) Plasma absorption corresponds to the energy necessary for the laser field to accelerate the seed electrons promoted in the conduction band. Energy is initially transferred to electrons but eventually ends up heating the medium.
Energy deposition to the medium corresponds to the sum of both contributions. Let us define the volume averaged density of absorbed energy, u NL , as the ratio of the absorbed energy per length unit by the cross section of the focal volume, where the integration volume is limited to the focal region r < r 0 , defined as the domain where the beam is intense enough to generate a plasma; i.e., r 0 denotes the radius of the plasma filament. For convenience, we will use below the normalized density of nonlinear energy absorptionũ NL ≡ u NL ρ 0 U g , where ρ 0 U g represents an energy scale corresponding to the energy that would be deposited in the medium if all valence electrons underwent a transition to the conduction band.  Model

III. NONLINEAR ABSORPTION MAP FOR A FLAT-TOP PULSE
We first present a fully analytical solution to Eqs. (1) and (2). In this aim, we simplified beam and pulse shapes to easily carry out the calculations using flat-top beam and pulse shapes. In Sec. III, we will numerically integrate Eqs. (1) and (2) for the Gaussian pulse and Bessel beam shapes. The comparison will show that the analytical solution provides a good approximation for the nonlinear absorption of laser energy (order of magnitude and trends vs parameters), even if the choice of pulse and beam shapes has a quantitative effect on the results.

A. Plasma density
Equation (1) is solvable analytically under certain conditions. For instance, if the laser pulse is approximated by a flat-top pulse of duration τ p and constant intensity I = I p for 0 ≤ t ≤ τ p , I = 0 otherwise, Eq. (1) becomes a Ricatti equation in the form where the coefficients ν 1 and ν 2 are listed in Table I for each recombination model. Only the expressions for parameters ν 1 and ν 2 differ when the recombination model is changed from Mq to Ml; otherwise, the formalism is the same for both models. For this reason, we present analytical expressions valid for both models, but we will illustrate results solely for model Mq with a recombination rate given by ν rec ≡ τ −1 r , where τ r is given in Table II for BK7. For 0 ≤ t ≤ τ p , the solution to Eq. (1) can be expressed as where

B. Nonlinear absorption of energy
Equation (2) can also be evaluated analytically in the case of a flat-top beam shape of radius r 0 with uniform intensity I = I p for 0 < r < r 0 . In this case, the radial integration amounts to a multiplication by the surface of the plasma filament πr 2 0 which cancels out with the denominator. We can thus work with the time integration only to calculate the density of energy absorption, normalized to the product of the gap by the density ρ 0 , whereĨ = I/I p . For a flat-top pulse, the normalized density of energy absorption reads where the integrals over the pulse duration are calculated from the electron density equation (3), resulting in closed analytical expressions and These expressions are general for transparent materials. We illustrate results obtained for BK7 and SiO 2 , with material parameters listed in Table II.

C. Nonlinear absorption maps for BK7 and for a flat-top pulse
The pulse intensity I p in the focal region is linked to experimental parameters by considering that the pulse energy, the pulse duration τ p , and focusing conditions are fixed, which determines the pulse fluence F in the focal region. The link between these quantities I p = F/τ p allows us to express the density of energy absorption and the electron density as functions of the pulse duration and peak fluence. Figure 3(a) shows the nonlinear absorption map obtained from Eqs. (4) and (8) for model (Mq). The highest density of energy absorption is obtained in the region of short pulses and high fluence (lower right corner) where multiphoton absorption dominates over plasma absorption. The density of energy absorption decreases monotonically as the pulse duration increases for a given pulse fluence or as the fluence decreases for a given pulse duration. This result may seem to differ from the conclusion by Liu et al. that the maximum electron density and absorption are associated with avalanche ionization and breakdown for microjoule femtosecond pulses propagating in transparent condensed matter. 41 However, our analysis in Sec. IV will show that a more realistic pulse shape leads to the appearance of a maximum in the density of absorbed energy, at picosecond pulse durations, thus stressing the prevailing role of avalanche ionization around the maximum, in agreement with the conclusion of Ref. 41.
Proceeding with the link between electron density and density of energy absorption, we note that Eq. (7) can be rewritten as (11) showing that if the recombination were not taken into account, the density of nonlinear absorption would be equal to the product of the gap, U g , by the radially averaged electron density. For the flat-top pulse and model (Mq), Eq. (11) becomes The first term on the right-hand side simply represents a radially averaged energy density u NL =ρ(τ p )ρ 0 U g equal to the product of the energy gap U g by the density ρ(τ p ) of conduction band electrons at the end of the pulse. This dominant contribution represents the minimum cost for the transition from the valence to the conduction band and varies monotonically with pulse duration τ p and fluence F. Two terms actually increase the density of absorbed energy above ρ(τ p )U g : (i) the second term on the rhs of Eq. (12) represents the product of the energy gap by the density of electrons that not only underwent a transition from valence to conduction band but also recombined during the pulse. (ii) The third term on the rhs of Eq. (12) represents a residual energy absorption due to the fact that the energy of the K photons involved in multiphoton ionization of an electron exceeds the energy gap. A good order of magnitude for the density of absorbed energy is therefore given by the product of the energy gap by the density of electrons in the conduction band immediately after the pulse. The peak electron density is shown in Fig. 3 multiphoton ionization prevails over avalanche ionization. An ionization degree of one, corresponding to all valence electrons in the conduction band, is asymptotically approached in the lower right corner of Fig. 3(b). It has been established that the damage fluence threshold for dielectric materials varies as τ 1/2 p for pulse durations larger than 10 ps (see Ref. 42). However, it is for short pulse durations that the red dashed line in Fig. 3(b) coincides with isolevels of the electron density, the slope of which corresponds to the τ 1/2 p scaling of the breakdown threshold. For longer pulses, the blue dashed-dotted line F ∝ τ p fits better the curves of constant electron density. Hence, from the standard empirical definition of breakdown threshold stating that breakdown occurs when the electron density exceeds a certain threshold, the value of which ranges from 10 17 to 10 20 cm −3 (see Refs. [43][44][45], and the large fluence values reached in the lower right corner of these maps (corresponding to electron densities above 2 × 10 21 cm −3 ∼ 10 −1 ρ 0 for short pulses) must be considered as well above the breakdown threshold and cannot be expected to accurately represent nonlinear absorption of energy. Nevertheless, Du et al. have proposed a scaling law F ∝ τ −1 p for the damage threshold of dielectric media under femtosecond laser pulse radiation, 42 potentially extending the validity domain to fluence levels larger than a few J/cm 2 for subpicosecond pulses. This is especially relevant with Bessel beams since energy deposition in the focal region is only weakly affected by material damage in the region of the central lobe owing to the conical energy flow featuring Bessel beam propagation.

IV. NONLINEAR ABSORPTION MAPS FOR BK7 AND SiO 2 WITH A GAUSSIAN PULSE
In this section, we still assume the pulse to be propagation invariant in the Bessel zone, but we relax the approximation of flat-top pulse and beam shapes. We write the intensity distribution in the Bessel zone in the form of the product of a Gaussian pulse profile with a Bessel beam profile, where a = 4 ln 2 and τ p denotes the full width at half maximum (FWHM) pulse duration. Integration of Eq. (1) can be performed semi-analytically with a Gaussian pulse (see, e.g., Ref. 46); however, the spatial integration of the Bessel beam intensity profile in Eq. (2) requires a numerical integration. We therefore obtained the plasma density and the density of nonlinear absorption of energy from a fully numerical integration of Eqs. (1) and (2). To express the results as functions of the pulse duration and peak fluence as in the case of the flat-top beam, we used the expression F = I p τ p √ π/a linking I p and the peak fluence F. We solved Eqs. (1) and (2) and plotted the nonlinear absorption maps for BK7 and for SiO 2 in Fig. 4 as a function of fluence and pulse duration. The density of nonlinear absorption of energy exhibits a similar general dependence upon τ p and F as for a flat-top pulse. A difference is however found for the region where the density of absorbed energy exceeds 2ρ 0 U g (contour level marked "2") and reaches a maximum for the highest values of fluence and pulse durations of several picoseconds in Fig. 4(a), whereas the highest absorption is found for the lowest right corner in Fig. 3(a) with pulse durations smaller than 100 fs. For fused silica, the material with the higher bandgap, this maximum is obtained for a fluence that is lower by a factor of ten. For both materials, the position of the maximum at picosecond pulse durations means that longer pulses are likely to be more efficient than femtosecond pulses for laser energy deposition, if the corresponding range of fluence can be reached in the Bessel zone. For pulse durations of several tens of picosecond, these fluence levels lie at the border of the damage threshold. For instance, the isocontour labeled "1" (i.e., 30.3 kJ/cm 3 ) in Fig. 4(b) is in excellent agreement with the damage threshold curve reported by Du et al. for fused silica, suggesting another empirical definition of the damage threshold as the fluence necessary to reach a certain density of energy deposition.
We note that the radius of the plasma filament r 0 formally depends on the cone angle, but this dependence cancels out in the numerical integration of Eq. (2) due to the plasma filament surface in the denominator. Therefore, the density of energy absorption depends on material parameters but not on the cone angle of the Bessel beam. In fact, the fluence in the Bessel zone is determined by the cone angle, but the normalized quantityũ NL plotted in Figs. 3 and 4 does not directly depend on cone angle. It represents the density of energy absorption that the medium can support under irradiation by a pulse with duration τ p and peak fluence F and is typical of the material.

A. Optimal absorption of laser energy
The energy absorption maps in Fig. 4 provide information about the ability of the medium to absorb a certain amount of energy density independently of focusing conditions. An additional analysis must be performed in order to determine whether this energy density is available for given beam and pulse parameters in specific conditions of axicon-focusing. In particular, if the pulse energy and cone angle are fixed by experimental conditions, we are interested in the possibility to play on the pulse duration, by chirping the pulse, in order to optimize laser energy deposition. How can we link material and laser parameters to find this optimum?
To answer this question, we first determine the peak fluence set by the constraints on the laser pulse and focusing conditions, sketched in Fig. 1. The peak fluence in the Bessel filament is determined by uniformly distributing the energy of the beam within the column of length L B and radius r 0 . The input beam can be viewed as a superposition of rings carrying a fraction 2r 0 / 0 of the total energy. Each ring is focused with the same cone angle θ at a certain distance within the Bessel zone; i.e., its energy is distributed over the surface πr 2 0 of the Bessel filament. Hence, the corresponding fluence is F 0 = (2r 0 / 0 )E in /(πr 2 0 ) = 2E in /π 0 r 0 . The peak fluence in the focal region is thus evaluated as a function of the pulse energy and the cone angle of the Bessel beam, Equation (14) exhibits a similar dependence upon the parameters to that of the peak fluence obtained in the linear propagation regime for an axicon-focused Gaussian beam. 15 This is explained by the fact that Bessel beam propagation is tantamount to a uniform distribution of the Gaussian beam energy in the focal volume, which is our assumption to obtain Eq. (14). Thus, we only find a different numerical prefactor which reflects the difference between the peak fluence and the averaged fluence over the focal volume. From Eq. (14) and Fig. 4, we extract a cross section representing the density of energy absorptionũ NL (τ p , F 0 ) supported by the medium. For a uniform distribution of the pulse energy over the whole volume of the Bessel filament, the available energy density in the focal region reads and depends on the cone angle, as shown by introducing the θ-dependence of both L B and r 0 into Eq. (15). The available energy density must be compared to the density of energy that the medium can absorb nonlinearly.
In this aim, we define the absorption coefficient by A = u NL /u B and write it as a function of the pulse energy, duration, and cone angle, where F 0 is given by Eq. (14). Note that Eq. (16) generalizes the model proposed in Ref. 47 to the cases where nonlinear absorption is not solely determined by mutiphoton absorption (see the Appendix for details). When the condition 0 < A < 1 is satisfied, the density of available laser energy is sufficient for nonlinear absorption to be homogeneous within the Bessel zone, in keeping with the propagation invariance of the Bessel beam. This usually occurs for long pulses. For short pulses or, depending on the cone angle, for pulse durations within a certain range, the level of fluence can be such that u NL (τ p , F 0 (E in , θ)) > u B (E in , θ); i.e., the incoming energy flux at the given pulse energy, pulse duration, and cone angle cannot sustain the corresponding high nonlinear absorption of energy. It is therefore impossible to satisfy our assumption of propagation invariance. Pulse propagation at these levels of fluence is non-stationary, featured by cycles of fast nonlinear absorption rapidly exhausting the incoming flux and decreasing the intensity in the focal region and followed by a slower restoration of the incoming energy flux due to the self-healing property of Bessel beams. 17 The validity range of our assumptions can therefore be derived by combining Eqs. (14) and (16) asũ Figure 5 illustrates this condition. The condition for propagation invariance expressed by Eq. (17) is not satisfied within non-stationarity regions. For a fixed laser beam width ( 0 = 180 µm in Fig. 5), these regions take the form of tongues in the plane spanned by the pulse duration and the cone angle. The boundaries of these tongues depend on laser pulse energy. For example, nonlinear Bessel beam propagation of a 70 µJ pulse in BK7 [ Fig. 5(a)], for a cone angle of 6.6 • , should undergo a transition from a propagation invariant to a non-stationary propagation regime when the pulse duration is decreased below 1 ps (crossing of the green boundary). The non-stationarity tongues in Fig. 5 show that for a given pulse energy and duration, a range of cone angles exists for which absorption of energy cannot be uniform in the Bessel zone. For a given beam width and pulse energy, long (picosecond) pulses and large cone angles are predicted to be compatible with invariant Bessel beam propagation and a uniform laser energy deposition. The optimal pulse duration and cone angle for laser energy deposition with a Bessel beam are found outside of and close to the corresponding non-stationarity tongue, where the absorption coefficient is approaching unity.
A comparison of Figs. 5(a) and 5(b) shows that the material with the lower bandgap, BK7, allows for a uniform laser energy deposition with Bessel beams in a propagation invariant regime over a much wider range of cone angles and pulse durations compared to fused silica. Pulses with significantly higher energy can be used. In fused silica, with moderate cone angles (θ < 10 • ), pulses of a few µJ can lead to a uniform energy deposition if their duration exceeds a certain threshold. This threshold increases with pulse energy. For instance with a cone angle of 10 • , this threshold reaches 1 ps for a pulse energy of 10 µJ. At the same cone angle and pulse energy, a uniform energy deposition can be obtained in BK7 for pulses longer than ∼20 fs. These results show that it is easier to obtain a uniform laser energy deposition with a Bessel beam in BK7 compared to fused silica, in keeping with the report of single shot, high aspect ratio laser drilling by focusing a femtosecond pulse in borosilicate glass with a conical lens. 5

B. Comparison with measurements and numerical simulation results
We performed experiments to compare the results given by the semi-analytical model [Eqs. (14) and (16)] with measurements of the integrated transmission through BK7. We used an infrared Ti:sapphire laser delivering 40 fs pulses and operating at the repetition rate of 20 Hz. The grating FIG. 5. Non-stationarity tongues for a fixed beam width ( 0 = 180 µm) and a set of pulse energies (see the inset) in (a) BK7 and (b) SiO 2 . Each curve forms a tongue inside which the density of absorbed energy u NL is higher than the available laser energy density u B . The condition (17) for propagation invariance in the Bessel zone is fulfilled outside the tongue corresponding to the pulse energy, i.e., for large pulse durations or large cone angles. pairs in the laser were detuned to lengthen the pulse duration up to 10 ps. The setup is similar to that used in Ref. 48: The Bessel filament was generated in BK7 by inserting in the path of the laser beam an axicon between two telescopic systems so as to image the Bessel zone after the axicon inside the sample. The cone angle was 10 • in air (θ ∼ 6.6 • in BK7). The equivalent beam width at the entrance of the Bessel zone is 0 = 180 µm, and the corresponding Bessel zone length is L B = 1.5 mm. The transmission was estimated just after the sample by measuring the total output energy of the beam with an energy meter. The measurement was performed by averaging the energy value over 20 pulses. Figure 6 shows the measured transmittance Tr for a pulse duration of 5 ps and pulse energies E in = 20, 30, 40, 50, and 70 µJ (F ∼ 4, 6, 8, 10, and 14 J/cm 2 ) and the transmittance obtained by our model Tr = Tr 0 (1 − u NL /u B ) for these parameters, where Tr 0 ∼ 0.85 represents the transmittance measured at low energy, i.e., the effect of linear losses and reflectivity (Fresnel losses). We also performed numerical simulations resolving a unidirectional envelope propagation equation coupled with Eqs. (1) and (2). The propagation model is presented in detail in Ref. 39. In contrast to our semi-analytical model, numerical simulations do not assume a propagation invariant beam and pulse, allowing for pulse shortening along propagation.
We obtain a fair quantitative agreement between the measurements (triangles) and the results of numerical simulations (circles), even if the latter slightly overestimate the transmittance. The general trend observed in the measurements is also qualitatively well reproduced by the semi-analytical model, the results of which are in good agreement with numerical simulation results for pulse energies below 50 µJ. At larger pulse energies (E in ∼ 70 µJ), the semi-analytical model clearly underestimates the transmittance. Our assumption of propagation invariance in the semi-analytical model is justified by the good agreement obtained for pulse energies below 50 µJ and by the pulse profiles obtained by numerical simulation results, which are similar to those in Ref. 39 and indeed confirm that after an initial beam reshaping and pulse shortening, the pulse propagates through the Bessel zone with only a smooth variation of the peak fluence and small oscillations, the amplitude of which increases with the pulse energy. Figure 5(a) shows that the working point for a pulse duration of 5 ps and a cone angle of 6.6 • belongs to the stationarity region for all pulse energies up to 100 µJ, but close to the frontier. Simulations tell us that the actual working point should have a smaller pulse duration due to the initial transient pulse shortening. Figure 5(a) also shows that for increasing pulse energies, the region of non-stationarity (interior of the tongues) extends toward larger pulse durations: A range of picosecond pulse durations enters the regime of non-stationarity while a range of subpicosecond pulse durations enter the regime of stationary propagation. For these reasons, the working point in Fig. 5(a) could travel from a stationarity region to the non-stationarity region when pulse energy is increased. We therefore interpret the appearance of the small oscillations in the simulation and the discrepancy for the transmittance at 70 µJ in Fig. 5 as a progressive transition to the regime of non-stationarity.
In addition to the progressive loss of stationarity at large pulse energies, the deviation between the semi-analytical model and simulations/experiments at E in = 70 µJ may be due to the fact that (i) material parameters for BK7 (Table II) are not known with high accuracy; (ii) our model considers multiphoton and plasma absorptions as the only phenomena for laser energy deposition, with a single functional form (dependence upon pulse intensity) in the entire parameter space; (iii) a frozen (flat-top) intensity profile has been assumed along the propagation distance in the Bessel zone; this assumption could be relaxed to allow for a slow variation of the peak fluence or plasma channel radius with propagation distance; (iv) absorption in secondary lobes of the Bessel beam was neglected; (v) our model entirely neglects plasma defocusing in the evaluation of the energy concentration in the Bessel zone. In an experimental situation, the fluence could be slightly lower than that evaluated from Eq. (14). These five possibilities are ranked from the most to the less important for the origin of the discrepancy at 70 µJ. More precisely, at the frontier between (i) and (ii), the more critical parameter is the collision time used to parameterize collisional absorption in BK7, in the context of the Drude model. Future work will be devoted to an improvement of the modeling of collisional absorption.
As a final remark, we stress the good matching of the boundaries of non-stationarity tongues calculated from the semi-analytical model with the frontiers of the region of uniform void formation obtained from recent observations of different damage morphologies in fused silica with Bessel beam illumination (see Fig. 3 of Ref. 7). For a given cone angle of 8 • as in Ref. 7, the non-stationarity tongues in Fig. 5(b) show an increase of the pulse duration above which a propagation invariant regime and uniform energy deposition are obtained when the pulse energy is increased. We interpret a uniform void formation within the Bessel zone as the result of a uniform laser energy deposition associated with an invariant propagation of the Bessel beam. The frontier between the regime of uniform void formation and the regime of void with fragmentation in Fig. 3 of Ref. 7 then qualitatively agrees with our results plotted in Fig. 5(b). Measurements from Ref. 7 show a significantly higher energy deposition in the case of picosecond pulse durations and a lower energy density concentration for femtosecond pulses. This is also the case in other materials, as indicated by the weak trace in Corning 0211 in Fig. 2(a). This optimum arises because of two key ingredients: (i) the local maximum on the map for the density of nonlinear absorption of energy in Fig. 4 lies in the picosecond range, making picosecond pulses more advantageous for laser energy deposition and (ii) for a uniform energy deposition, the laser pulse energy and focusing conditions are selected out of the corresponding non-stationarity tongue in Fig. 5(b).

VI. CONCLUSION
We have proposed a semi-analytical model to map the density of absorption of laser energy in glasses (BK7 and SiO 2 ) as a function of pulse duration and peak fluence in the focal region. We have shown how to select optimal focusing conditions and laser pulse parameters for a uniform distribution and absorption of energy in the focal volume corresponding to axicon-focusing of a Gaussian beam, the Bessel zone, where a plasma is generated by optical field ionization (multiphoton ionization) and by avalanche.
From the maps of laser energy absorption and given laser parameters (pulse energy, beam width), we determined the region of parameters compatible with a uniform deposition of laser energy. We showed that in a tongue-shaped region of the plan spanned by the pulse duration and the cone angle for the Bessel beam, the nonlinear absorption of laser energy is too fast with respect to the available laser energy for energy deposition to be uniform. Optimal uniform energy deposition for a given pulse energy is predicted to occur for pulse durations and cone angles lying close to the frontier of the non-stationarity tongue.
Results reproduce qualitatively the transmittance in BK7 measured in a dedicated experiment, provided conditions of propagation invariance remain fulfilled. Numerical simulation results, in agreement with the measurements, helped us to confirm the loss of propagation invariance when the pulse energy was increased, leading to an overestimation of nonlinear absorption. Earlier observations of different types of damage morphology obtained by Bessel beam propagation in fused silica 7 were also interpreted in light of our findings. The transition between the regime of uniform void formation and the regime of void with fragmentation occurs for picosecond pulse durations in fair agreement with the frontier of the non-stationarity tongues predicted from our model.
A uniform laser energy deposition with Bessel beams opens the way to applications such as highaspect-ratio structuring of glasses. The conjunction of maps of nonlinear absorption of laser energy and the non-stationarity tongues help us select optimal parameters for these applications. Since a nonlinear absorption map depends on material parameters only, it can also be used for different focusing conditions, e.g., tight focusing of Gaussian beams with lenses or parabolic mirrors, in the aim of laser energy deposition in the bulk of dielectric media.

ACKNOWLEDGMENTS
We acknowledge support from the PICS program (CNRS), Projet International de Coopération Scientifique for the support of collaboration, and exchange visits between the University of Insubria and CNRS, France.

APPENDIX: LIMIT OF EQ. (16) IN THE ABSENCE OF AVALANCHE AND RECOMBINATION
Equation (16)