1. Introduction

Femtosecond laser pulses have become versatile tools for micromachining glasses and transparent materials, such as quartz, sapphire or diamond [1, 2]. Accuracy of the order of one micron can be reached, however, multiple laser shots are usually necessary for most applications. Numerous challenging applications in nanophotonics and nanofluidics require the fabrication of high aspect ratio nano-channels, a key technological issue for laser ablation with Gaussian beams since their spatiotemporal distortion during propagation in dielectrics prevents fine control of energy deposition on the longitudinal dimension. For example, the propagation of Gaussian pulses lead to self-focusing due to the optical Kerr effect or defocusing due to their interaction with the self-generated electron plasma [3]. The need of high number of laser pulses for generation of high aspect ratio features consumes most of time in bulk micro machining. Recently, intense Bessel beams have been proposed as excellent candidates for solving this issue [4–7] and successfully used to produce high aspect ratio micro- and nano-channels in glass [8,9]. Nanochannels were drilled in glass by a single femtosecond Bessel pulse with high conical angle (17° in glass) [9].

In this paper we investigate theoretically and numerically a new type of beams, nonlinear Bessel vortices (NBVs) that started to be used recently for material processing [10, 11] or investigations of propagation invariant beams [12] and are especially suitable for uniform energy deposition on a novel tubular geometry. These laser beams fulfill in particular the following conditions: (i) Their nonlinear propagation in transparent solids is quasi-stationary. Even at high intensity, these beams do not spread and undergo only minor spatial distortions. (ii) Several degrees of freedom allow for an independent adjustment of the length and diameter of the high intensity region within the material to be modified by laser illumination. Condition (i) is fullfilled by Bessel beams, or their appodized version called Bessel-Gauss beams, which are easily produced in laboratory and known for more than quarter of a century as the simplest non-diffractive waves [13–15]. Bessel-Gauss beams propagate quasi-stationarily through a Kerr medium even if they undergo nonlinear absorption, i.e., nearly no change occurs in their intensity profiles [6, 16, 17]. Similar properties were recently demonstrated for nonlinear Airy beams [18, 19]. It may appear as a contradiction that a wave undergoes nonlinear absorption while propagating stationarily since nonlinear absorption could be expected to induce attenuation of the main lobe. However, this attenuation does not occur due to a general self-healing property shared by nonlinear conical waves, and interpreted as the result of a beam reshaping into a nonlinear Bessel beam (NBB) [20]. These beams are indeed featured by a high intensity region where nonlinear absorbtion occurs and surrounded by a large reservoir (the rings of the Bessel profile) that behaves linearly. An energy flux from the reservoir to the high intensity region balances energy losses and ensures quasi-stationarity over a certain propagation distance, the Bessel zone, that only depend on the size of the energy reservoir and the cone angle characterizing the Bessel beam. Since NBB exists as a two-parameter family of stationary solutions to the Nonlinear Schrödinger equation, it is possible to tune independently the peak intensity and the width of the high intensity core, via the cone angle of the Bessel beam. High energy and intense Bessel beams may produce long and narrow plasma channels, which later produce a material damage string with very high aspect ratio [21]. However, the Bessel zone cannot be tuned independently.

The recent development of Spatial Light Modulators (SLM), phase shift masks, and amplitude gratings opens up the opportunity to design new wave-packets so as to introduce tunability of the Bessel zone independently and fulfill condition (ii). This can be realized by using higher-order Bessel beams carrying topological charge or a phase singularity. If a vortex phase is added to a Bessel beam phase, the reconstructed beam will be a Bessel Vortex (BV), which also belongs to the class of non-diffractive solutions to the Maxwell equations. Here we show that NBVs, the nonlinear counterpart of BVs, share the properties of BVs and NBBs. The length of the quasi-stationary nonlinear propagation and the lobe width can be controlled independently by tuning three parameters: the topological charge m of the vortex, the waist of the input Gaussian beam entering the spatial light modulator, and the cone angle of the Bessel beam. We note that conventional vortex beams do not fulfill both conditions (i) and (ii) since high power optical vortices were shown to undergo self focusing. Under certain conditions the vortex in nonlinear medium splits into multiple filaments, and the filament pattern rotates along the propagation axis [22]. In contrast, during the propagation of NBV, the cylindrical volume surrounded by the main ring-lobe is not affected by nonlinear effects and therefore can be preserved.

Herein we identify NBVs and investigate their propagation in transparent solids (Kerr media with nonlinear losses). We show that the propagation of NBVs is stable in a large region of the parameter space. An unstable propagation regime is also identified where the propagation of NBV becomes oscillatory.

2. Direct numerical simulations

We first investigate numerically the propagation of Bessel-Gauss vortices, i.e., BVs with finite energy generated by propagating a large Gaussian beam with a topological charge through an axicon. Details on the initial condition are found below the presentation of the numerical model.

2.1. Numerical model

We study numerically the propagation of the intense monochromatic field envelope in Fourier space Ê(k_x, k_y, z), by means of the unidirectional nonparaxial propagation equation [23]:

2.2. Material parameters

We investigate nonlinear propagation of BV through the glass and Table 1 in the appendix lists the parameters used in the model, unless otherwise stated. Material parameters are adapted for Corning 0211 glass (with nonlinear coefficients of BK7). In particular, we used the Drude model with a very short collision time with respect to standard values in the femtosecond range found in the literature. This choice arises because we view τ_c only as a fitting parameter allowing to model the material response by the Drude model with similar quantitative permittivity as obtained by time-dependent density functional theory calculations [25, 26] However, our results are not specific to these parameter values. Qualitatively similar results are found for a wide range of material and laser parameters.

2.3. Initial condition

Since ideal Bessel beams carry infinite energy, the common way to generate them in the laboratory is to send a large Gaussian beam through an axicon, forming a Bessel-Gauss beam, which can be viewed as a Bessel beam with a narrow main lobe apodized by a large Gaussian apodizer, with quasi-stationary linear propagation over a finite distance called the Bessel zone.

This property remains valid for BV: By sending a large Gaussian beam carrying topological charge through an axicon, a Bessel-Gauss vortex is obtained which carries finite energy and propagates quasi-stationarily. We will therefore always consider the finite energy realization of Bessel beams or BV by starting our simulations from a finite energy Gaussian beam with the spatial phase corresponding to that of a BV.

2.4. Numerical simulation results

Figure 1(a)–1(c) shows the result of the propagation of a Bessel-Gauss beam in glass. For reference, Fig. 1(a) depicts the case of linear propagation of the zero-order Bessel beam. A narrow and intense peak is visible in the Bessel zone, marked by the continuous white lines showing the interference region for the inward and outward Hankel components of the Bessel beam. The peak intensity of the beam increases to its maximum value in the center of the Bessel zone and then decreases. Figure 1(b) shows the nonlinear propagation of the zero-order Bessel beam in the Corning 0211 glass and reveals a new behavior: the energy reservoir of the Bessel beam transfers energy to the central core thus enabling the intense central lobe to propagate and even maintain an almost constant intensity over an extended propagation distance, in spite of nonlinear absorption. The propagation of the intense lobe remains stable over an elongated region (delimited by white dashed lines) with respect to the linear case and the surrounding lobes propagate with minor changes as well. This observation emphasizes that the beam profile remains quasi-stationary (z-independent) and Bessel-shaped. The stationary NBBs that play the role of attractors for the beam dynamics are well identified (nonlinear unbalanced Bessel beams in [16, 17]). These are z-independent solutions to Eq. (1) in which only the most important physical effects are present, namely, diffraction, the optical Kerr effect and nonlinear losses due to ionization. The central lobe is sufficiently intense to generate a narrow plasma channel with high and nearly constant density over the Bessel zone (Fig. 1(c)). The plasma channel does not affect significantly the properties of NBB particular, it does not perturb the stationarity of nonlinear propagation. Stationary nonlinear solutions to Eq. (1) indeed not only exist for pure Kerr media as shown in Refs. [16, 17], but also for media where the refractive index varies arbitrarily with intensity, including Kerr media undergoing ionization. These solutions play the role of attractors in our simulations.

 

Fig. 1 Zero-order (a-c) and first order (d-f) Bessel beam propagation through Corning 0211 glass. Intensity profiles for (a,d) linear propagation, (b,e) nonlinear propagation. (c,f) Plasma density profile for nonlinear propagation. The slope of the white lines is equal to the cone half-angle. The Bessel zone lies between the continuous lines. The extended Bessel zone in the nonlinear regime is marked by dashed lines. The input beam power is 64.1 MW and the corresponding pulse energy is 7.8 μJ.

Download Full Size | PDF

Fig. 1(d)–1(f) show the propagation of a BV (m = 1) with the same energy as in Fig. 1(a)–1(c). The case of linear propagation (Fig. 1(d)) corresponds well to the theory developed in Refs [27]. A pipe like light channel (Fig. 1(e)) and plasma track (Fig. 1(f)) are generated. The main feature of interest is the smoothness of the variation of intensity profiles along the propagation distance. The results suggest that the simulations converge to a quasi-stationary NBV solution when the input topological charge is nonzero. This indicates the possible existence of attractors in the form of stationary NBVs, which we identify in section 3. On a more quantitative level, the ratio between energies in the first lobe of the Bessel beam (m = 0) and BV (m = 1) is around 2.1. The energy contained in the first lobe of BV is higher, but it is also distributed over a larger area, therefore the maximum intensity in the linear case is only third of the maximum intensity corresponding to the Bessel beam. BVs with the same energy but higher topological charge propagate in a qualitatively similar way but the ring diameter of the main lobe increases with m and the peak intensity and electron density correspondingly decrease. By comparing Figs. 1(a)–1(c) and 1(d)–1(f), we note that the plasma density and intensity in the case of the NBV with m = 1 is almost twice lower than in the case of the NBB (m = 0). Similar values for the plasma density could be reached by increasing the initial beam power in the case of the NBV.

Two different regimes for the nonlinear propagation of BV are presented in Fig. 2: 2(a) a stable propagation regime in the sense that the beam approaches a stationary solution, and 2(b) an unstable propagation regime featured by oscillations along the propagation distance. The unstable regime was reached when the plasma defocusing and absorption terms were switched off while keeping the same initial energy. The maximum intensity reached during propagation is much higher, leading to an oscillatory propagation. By lowering the power the stationary propagation regime is achieved Fig. 2(c). The results suggest that the plasma terms help to stabilize nonlinear propagation by lowering the maximum intensity. Note that the beam in the unstable propagation regime is not radially symmetric. It forms a speckle which rotates with respect to the beam center (see Fig. 3).

 

Fig. 2 Nonlinear propagation regimes for BV. Figures show intensity maps of the propagation of the beam sliced in one transverse direction. (a) Stable propagation when plasma defocusing and absorption terms are active, (b) unstable propagation when the plasma terms are switched off keeping the same initial beam power, and (c) stable propagation when plasma terms are switched off but initial power lowered from 64.1 MW to 38.5 MW.

Download Full Size | PDF

 

Fig. 3 (a) Stable nonlinear propagation regime for a first order nonlinear Bessel beam was reached with 38.5 MW power. (b, c) Unstable propagation for a NBV beam with the same cone angle but power 57.7 MW and 76.9 MW. (d) Unstable propagation for a third order NBV 153.8 MW. Each intensity slice was normalized to its maximum intensity. Different isosurface colors represents intensity levels. Blue color: 0.5×I_max. Green color: 0.7×I_max. Red color: 0.9×I_max. Zone 1: stable nonlinear propagation (beam smooth shrinkage and relaxation). Zone 2: break-up into rotating peaks: Zone 3: irregular rotation of the peaks.

Download Full Size | PDF

The isosurface plots shown in Fig. 3 illustrate the complex behavior of an unstable BV, suggesting a classification of different propagation zones according to observed beam dynamics. The stable nonlinear propagation of the BV (Fig. 3(a)) shows a smooth beam shrinkage and relaxation zone (Zone1) extending over and beyond the Bessel zone. By increasing the initial beam power, the propagation becomes unstable (Fig. 3(b)). The smooth beam shrinkage is still visible, but the beam quickly breaks up into a rotating peak. The peak rotates periodically in zone 2, maintaining a smooth propagation even beyond the Bessel zone. The oscillatory behavior of the peak intensity proceeds from the conjunction of two effects: on the one hand, axial intensity modulation occurs in analogy with the instability of first order Bessel beams whose central lobe oscillates in propagation [17, 28]. In our simulations, the beam starts to oscillate before the appearance of the rotating peak. On the other hand, the beam carries helicity (topological charge) that triggers symmetry breaking for the breakup and leads to an helical beam pattern. Although the helicity seems to be preserved, propagation does not preserve the intensity profile of a high-order Bessel beam. By increasing power even more, a third zone appears where the rotation of the peak becomes irregular and no longer rotates around the propagation axis (Fig. 3(c)). The unstable propagation of a higher order bessel beam (m = 3) exhibits the three types of behaviors in the three zones identified above, however, not only one but multiple rotating peaks are formed in zone 2, corresponding to the topological charge of the input beam (Fig. 3(d)).

Our simulations suggest the existence of stationary NBVs which play the role of attractors of the propagation dynamics for high-order Bessel beams, at least in zone 1. An unstable propagation can be obtained either because a given stationary NBV becomes unstable (e.g. it undergoes modulational instability with respect to specific type of perturbations), or because no stationary solution exists at all for the laser or material parameters under investigation. This calls for a detailed investigation of stationary BV.

3. Stationary solutions

In order to interpret our numerical simulations, we look for solutions which may play the role of attractors for the nonlinear dynamics, i.e., we look for nonlinear stationary solutions in the form of BV. We will then investigate the possible convergence of simulations to a stationary BV. Guided by previous works on stationary NBBs [16,17], we consider a model in the form of a nonlinear Schrödinger (NLS) equation that accounts for the main physical effects playing a role during propagation, namely diffraction, the optical Kerr effect and nonlinear losses (NLL). The propagation model reads:

Stationary BV are defined as solutions with a stationary (z-independent) intensity profile and a linear dependence of the phase profile upon the propagation variable z. The electric field envelope is therefore written as E(r, ϑ, z) = a(r)exp(−iδz + iϕ(r) + imϑ), where the wave vector shift is positive δ > 0. We further note that the phase ϕ(r) is entirely determined by the phase gradient q(r) ≡ dϕ/dr, up to a constant. By introducing this amplitude and phase decomposition into Eq. (6), we obtain the set of ordinary differential equations:

A stationary NBV is any solution of the set of equations (7,8) satisfying three boundary conditions corresponding to the requirements that a(r) → 0 as r → ∞ and a(0) = 0, q(0) = 0. Details on boundary conditions are given in Appendix 2.

3.1. Linear Bessel vortices

We first look for stationary solutions to Eqs. (7,8) in the regime of linear propagation (n₂ = 0; β_K = 0). The second equation (8) ensures that the quantity rqa² is preserved, which is satisfied by choosing q = 0. The first equation (7) then reads:

3.2. Power conservation and stationarity principle

Radial integration of Eq. (8) (after multiplication by ra²) results in the energy conservation law for weakly localized beams carrying infinite power:

3.3. Nonlinear Bessel vortices

The solutions to Eqs. (7,8) with their boundary conditions were evaluated numerically for a given medium characterized by its Kerr index n₂ and nonlinear absorption β_K coefficient. By using a fixed cone angle, (i.e. fixed axial wave-number shift, δ), we determine a family of solutions with shapes close to the BVs $J_{\pm m}\left(\sqrt{2k_0\delta }r\right)$ which exhibit a different maximum intensity. The maximum intensity I_max can take any arbitrary value up to a certain limit that only depends on the parameters n₂, β_K, and δ. This result is similar to the one known for nonlinear Bessel beams (m = 0): In Ref. [16] the region where stationary nonlinear Bessel beams exist was shown to be bounded by a frontier:

We determined the region of existence for NBVs with m = 1 (a similar region also exists for higher values of m) in Corning 0211 glass. (material parameters: multiphoton absorption β_K = 4.7 × 10⁻²⁷ cm³/W², nonlinear index n₂ = 3.45 × 10⁻¹⁶ cm²/W, refractive index n = 1.53). The boundary was determined by the requirement, that the NBV amplitude will tend to zero, while increasing the radius to infinity. Fig. 4 shows the frontier bounding the region of allowed cone angles where stationary NBVs exist in the plane spanned by cone angles and maximum intensity. A two parameter family of stationary NBVs exist, parametrized by the cone angle and the peak (maximum) intensity. A region of forbidden cone angles where numerical solutions to Eqs. (7,8) do not decay at r → ∞ is always found at large intensity. These remarks do not depend on the choice of the medium. Only the location of the frontier between allowed and forbidden cone angle regions depends on material parameters (n₂ and β_K). For example, for the angle of 16.4° the maximum intensity for which a NBV still exists in Corning 0211 is around 2.8 × 10³ TW/cm².

 

Fig. 4 Estimated intensity threshold for various cone angles defining the frontier between the decaying and non decaying NBV amplitude profiles for glass. NLL coefficient β_K = 4.7 × 10⁻²⁷ cm³/W², nonlinear refractive index n₂ = 3.45 × 10⁻¹⁶ cm²/W. The light gray zone depicts the area where the NBV amplitude does not decay as $1/\sqrt{r}$, but the increase of amplitude is not observed in computation window. Horizontal black line depicts 16.4° angle and black dots depict maximum intensity of NBV shown in Fig. 5(c).

Download Full Size | PDF

A short presentation of the properties of the tails of NBVs is particularly useful to illustrate the stationarity principle and the error bars associated with the transition between allowed and forbidden cone angles. When r is large, the amplitude of NBVs tends to zero and multiphoton absorption becomes negligible. NBVs behave as solutions to Eq. (9) and can therefore be seen as a supperposition of two linearly independent Bessel (Hänkel) functions of third kind [29], of same cone angle $\theta ={\text{sin}}^{-1}\left(\sqrt{2\delta }{k}_0\right)$, but with different weights α_out and α_in:

 

Fig. 5 Comparison of the stationary NBV when the NLL coefficient is changed. (a)β_K is changed (from 4.7 × 10⁻²⁷ to 42.3 × 10⁻²⁷ cm³/W²). The black arrow represents increasing β_K values while n₂ is fixed at 3.45 × 10⁻¹⁶. (b,d) The nonlinear refractive index n₂ is changed (from 3.45 × 10⁻¹⁶ to 17.25 × 10⁻¹⁶ cm²/W) and β_K = 4.7 × 10⁻²⁷ cm³/W². The black arrow represents increasing nonlinear index coefficients. The higher maximum intensity in case (d) induces a stronger influence of NLL compared to case (b). (c) β_K and n₂ coefficients are fixed, while leaving free the maximum intensity. For all figures, the Bessel beam angle is θ = 16.4°.

Download Full Size | PDF

3.4. Effect of material parameters

We know that for NBBs, the nonlinear index has an influence on the compression of the rings, while multiphoton absorption mostly influences ring contrast [5,17]. Nonlinear Airy beams exhibit this behavior as well [18,19]. The same behavior is found for NBV but we also highlight a new feature: both parameters have an influence on ring compression and attenuation of contrast.

Figure 5 compares NBVs showing the individual parameter effect to its shape. By increasing the NLL coefficient, it is shown in Fig. 5(a) that not only the ring contrast is drastically attenuated, but also the period of the rings is slightly altered. The period change is mostly evident for the rings which are further away from the center. The increase of the nonlinear refractive index leads to a drastic ring compression (see Figs. 5(b) and 5(d)).

Stationary BV shown in Fig. 5(b) are mostly influenced by the Kerr nonlinearity but those shown in Fig. 5(d) were obtained for parameters near the frontier of the allowed cone angles. Their higher peak intensity results in a higher sensitivity to NLL. By increasing the Kerr coefficient n₂, not only rings compress, but also the contrast increases.

3.5. Interpretation of numerical simulations

We propose an interpretation of our results of numerical simulations on the propagation of BV, based on the existence of families of stationary NBV characterized by a fixed cone angle (θ = 16.4°) and a peak intensity varying continuously up to a certain limit.

3.5.1. Stable propagation regime

In the case of a stable nonlinear propagation regime, at each propagation distance in the Bessel zone, the beam intensity exhibits a profile similar to a NBV. The beam profile from the direct numerical simulation undergoes only slight changes along the propagation distance, suggesting that it can be identified with a specific stationary NBV. Over the Bessel zone, the peak intensity varies slowly due to the fact that direct numerical simulations were performed with finite energy beams whereas stationary NBV represent ideal weakly localized waves with infinite energy. Therefore NBV constitute ideal attractors for the dynamics but the solution from direct numerical simulation never really reaches a specific attractor. Rather, we show that over the Bessel zone, the solution from our direct numerical simulations follows a well defined family of stationary NBV characterized by a specific cone angle determined by the input condition.

To identify the closest stationary NBV to our direct numerical solution, we need a practical criterion to extract the cone angle from a high-order Bessel-like beam profile resulting from nonlinear propagation of an axicon-focused large Gaussian beam with a topological charge. In a future dedicated publication, we will show that it is possible to quantify the cone angle change induced by nonlinear propagation. In the present paper, we characterize each beam at a given propagation distance by the position r_max of its peak intensity and we propose a simple link between peak intensity position and cone angle. We monitor r_max as a function of the peak intensity I_max of the beam when z varies. Fig. 6(a) and 6(b) show the dependance of r_max on peak intensity I_max. Each curve represents a direct numerical simulation for BV propagation with a given initial power. In other words, the propagation distance z increases along each curve. We first comment Figure 6(a) corresponding to the case of a stable nonlinear propagation. The small values of I_max are not relevant for our interpretation as they are out of the Bessel zone but each curve exhibits a very interesting part with growing peak intensity and smooth compression (decreasing r_max) followed by decreasing peak intensity and relaxation (increasing r_max). The point with the highest intensity on each curve is the counterpart of the center of the Bessel zone in linear propagation regime. To complete this analogy, we remind that a finite power Bessel-Gauss beam is formed at this point with the same cone angle as the corresponding ideal Bessel beam. In the case of nonlinear propagation, we therefore obtain the NBV of finite power that is the closest to the ideal stationary NBV characterized by the same cone angle. We view the stable nonlinear propagation within the Bessel zone as a smooth beam reshaping that consists in following the family of stationary NBV featured by the same cone angle.

 

Fig. 6 Position r_max dependence on intensity I_max for the peak intensity obtained by direct numerical simulation (a) with and (b) without plasma terms. The colored curves represent direct numerical simulations with different initial peak powers (6.4, 19.2, 38.5, 57.7 MW). The thick blue line represents the position-intensity dependence of the peak for stationary NBV. Position vs intensity for the peak intensity of stationary NBV profiles calculated for different cone angles is shown in (c) (blue solid curve). The black dotted curves are obtained by the analytical representation for the r_max vs I_max characterization proposed in text.

Download Full Size | PDF

To confirm this assertion, we plotted the thick blue curves in Fig. 6 representing the family of stationary NBV featured by the same cone angle (constant δ). Since stationary NBV profiles with a fixed cone angle are only known numerically, we propose an analytical representation for the position of the maximum intensity r_max as a function of the peak intensity I_max. It was obtained by assuming that in the vicinity of the beam center, the beam profile is very similar to the first-order Bessel beam with the same effective cone angle - $J_1\left(\sqrt{2k_0(\delta +{\delta }_n)}r\right)$, where δ_n ≡ 0.8 × k₀I_maxn₂/n₀ accounts for the axial phase shift induced by the Kerr effect. The coefficient 0.8 was obtained by a fitting procedure and is due to the fact that the peak intensity is located at finite distance from the propagation axis. We checked that the position for maximum intensity (thick blue curve) is accurately reproduced by $r_{\text{max}}=1.8412/\sqrt{2k_0(\delta +{\delta }_n)}$, where 1.8412 denotes the position for the maximum of the first-order Bessel function. As shown in Fig. 6(c), the analytical representation of the r_max vs I_max characterization of stationary BVs (plotted in black dots) nearly overlaps its counterpart obtained from direct calculation based on the NBV profiles (blue curve).

Similar to the results for NBBs (Fig. 3 in Ref. [17]), it is found that all curves representing the results of direct numerical solutions tend to follow the thick blue curve and cross it for the highest intensity reached in the direct numerical simulation. In other words, for a given cone angle, the nonlinear propagation of BVs with finite energy is governed by optimal stationarity, manifested in the form of a continuous beam reshaping so as to follow the family of stationary NBVs with the same cone angle.

3.5.2. Unstable propagation regime

We now comment Fig. 6(b) which incorporated the case of unstable nonlinear propagation of BV. It is obtained by switching off plasma terms (ρ = 0) in Eq. (3). Note that a similar trend is obtained from our simulation results including plasma terms but at a lower cone angle (data not shown). Again, the curves representing the results of direct numerical simulation show a regular increase of the peak intensity I_max in conjunction with a beam shrinkage (decrease of r_max) but above a certain intensity, it is followed by a second stage with drastic variations of these quantities, i.e., an oscillatory and irregular decrease of I_max and increase of r_max. Although the thick blue curve still characterizes the general trend of the direct numerical simulation, the symmetry of the beam profile is actually ruined by the oscillatory behavior, as shown in Fig. 2, indicating that stationary NBVs characterized by the Bessel cone angle of 16.4° no longer attract the dynamics, or not sufficiently to prevent large amplitude oscillations (see Fig. 6(b) black curve). An increase in input power results in an unstable propagation of the beam characterized by larger oscillation amplitudes. Note that the maximum intensity reached during the direct numerical simulation lies well below the border of forbidden cone angles in Fig. 4. A stationary NBV therefore exist in all the domain covered by the direct numerical simulation. Hence, the oscillatory behavior of direct numerical simulations is explained by an instability rather than by the loss of existence of stationary NBV. Above a certain peak intensity, stationary NBV become modulationally unstable to azimuthal perturbations. In this case, they can not play the role of attractors for the nonlinear propagation dynamics.

We finally highlight the general character of stationary NBVs identified in our semi-analytical approach. Figure 6(a) showed that for a large range of input powers, these solutions attract the nonlinear beam dynamics observed in direct numerical simulations even in the presence of plasma effects, although plasma effects were neglected in our semi-analytical approach. Their inclusion in Eqs. (7,8) can be seen as a change of the pure Kerr nonlinearity and multi-photon absorption into a more general intensity dependence for the nonlinear refractive index and nonlinear losses, respectively, which would not change the nature of our results. Only the peak intensity of stationary NBVs would be affected. In Fig. 6(a), we observed that plasma defocusing and absorption limited the growth of the peak intensity and prevented the beam from reaching the unstable propagation regime. However we underline that this is not solely due to plasma effects since nonlinear propagation of BV with the same parameters but a smaller cone angle (e.g. half, data not shown) leads to an unstable propagation as well.

4. Conclusion

We have identified families of stationary NBVs, which share several properties with NBB and nonlinear Airy beams. In particular, for a given material and a given cone angle, these solutions exist as a continuous family up to a maximum peak intensity. Direct numerical simulations show that a stable propagation regime exists for nonlinear Bessel vortices obtained by focusing with an axicon with large cone angle a Gaussian beam carrying topological charge. A BV is formed and propagates by following closely a family of stationary BVs featured by the Bessel cone angle corresponding to the axicon. High cone angles allow maintaining a stationary regime. By decreasing the Bessel cone angle, we reduce the highest intensity bounding existence and stability regions for stationary NBVs. The unstable propagation of NBV is featured by rotating peaks. In this regime stationary NBVs exist but do not play the role of attractors for the dynamics. Reaching the unstable propagation regime at lower intensity thresholds has consequences on potential applications, since it makes impossible to produce a smooth damage track in a transparent material. Our study was performed in glass but we stress all results are qualitatively valid whatever the order or type of nonlinearity. This is very promising for future novel applications of filamentation in a wide range of fields.

Appendix A: Material and beam parameters

Table 1. Material and beam parameters.

View Table

Appendix B: Boundary conditions for stationary Bessel vortices

The boundary conditions for a stationary Bessel vortex are slightly different from those for a Nonlinear Bessel beam (m = 0 in Equations (7,8)). For a Nonlinear Bessel beam, Eqs (7,8) are solved as a set of ordinary differential equations by adding a boundary condition at the origin r = 0. The boundary conditions are then expressed in terms of the unknown peak amplitude a₀ > 0 and read [16]:

(17)$$a(r=0)=a_0,$$

(18)$$\dot{a}(r=0)=0,$$

(19)$$q_0=0.$$

The requirement that a(r) → 0 as r → +∞ would require that Eqs. (7, 8) and the set of boundary conditions (17,18,19) be solved as a nonlinear eigenvalue problem where a₀ plays the role of the eigenvalue. However, a continuum of solutions exists, i.e., the requirement that a(r) → 0 as r → +∞ is automatically satisfied for all values a₀ up to a maximum [16] allowing us to avoid solving a nonlinear eigenvalue problem.

Boundary conditions for NBV are adapted from those for nonlinear Bessel beams. Equations (7,8) are still solved as a set of ordinary differential equations but the initial amplitude a₀ is equal to zero. Therefore, boundary conditions at the origin are set by requiring that for small amplitude, NBV behave as their linear counterpart. Their amplitude should therefore be given by the main order a(r) ∼ a_mr^m when r → 0 and the boundary conditions read

(20)$$a(r\to 0)~a_m{r}^m,$$

(21)$$\dot{a}(r\to 0)~m{a}_m{r}^{m-1},$$

(22)$$q(r\to 0)~-\frac{k_0{\beta }_K}{2(mK+1)}{a}_m^{2K-2}{r}^{2m(K-1)+1}.$$

As for nonlinear Bessel beams, the unknown a_m serves as an eigenvalue. Different values of a_m correspond to different solutions with different maximum amplitude. The requirement that a(r) → 0 as r → +∞ is satisfied for a continuum of solutions, allowing us again to avoid the resolution of a nonlinear eigenvalue problem. Each solution is therefore obtained by choosing a value for a_m and shooting along the unstable manifold of the fixed point (a = 0, ȧ = 0, q = 0) for the dynamical system (7, 8).

Acknowledgments

The authors acknowledge helpful discussions with Razvan Stoian. This work was supported by Region Franche-Comte and the French ANR, contract 2011-BS04-010-01 NANOFLAM.

References and links

1. R. R. Gattass and E. Mazur, “Femtosecond laser micromachining in transparent materials,” Nat. Photonics 2, 219–225 (2008). [CrossRef]  

2. G. Della Valle, R. Osellame, and P. Laporta, “Micromachining of photonic devices by femtosecond laser pulses,” J. Opt. A: Pure Appl. Opt. 11, 013001 (2009). [CrossRef]  

3. A. Couairon and A. Mysyrowicz, “Femtosecond filamentation in transparent media,” Phys. Rep. 441, 47–189 (2007). [CrossRef]  

4. P. Polesana, D. Faccio, P. Di Trapani, A. Dubietis, A. Piskarkas, A. Couairon, and M. A. Porras, “High localization, focal depth and contrast by means of nonlinear Bessel Beams,” Opt. Express 13, 6160–6167 (2005). [CrossRef]   [PubMed]  

5. P. Polesana, A. Dubietis, M. A. Porras, E. Kucinskas, D. Faccio, A. Couairon, and P. Di Trapani, “Near-field dynamics of ultrashort pulsed Bessel beams in media with Kerr nonlinearity,” Phys. Rev. E 73, 056612 (2006). [CrossRef]  

6. P. Polesana, M. Franco, A. Couairon, D. Faccio, and P. Di Trapani, “Filamentation in Kerr media from pulsed Bessel beams,” Phys. Rev. A 77, 043814 (2008). [CrossRef]  

7. A. Marcinkevicius, S. Juodkazis, S. Matsuo, V. Mizeikis, and H. Misawa, “Application of Bessel beams for microfabrication of dielectrics by femtosecond laser,” Jpn. J. Appl. Phys. 40, 1197–1199 (2001). [CrossRef]  

8. M. K. Bhuyan, F. Courvoisier, P.-A. Lacourt, M. Jacquot, L. Furfaro, M. J. Withford, and J. Dudley, “High aspect ratio taper-free microchannel fabrication using femtosecond Bessel beams,” Opt. Express , 18, 566–574 (2010). [CrossRef]   [PubMed]  

9. M. K. Bhuyan, F. Courvoisier, P. A. Lacourt, M. Jacquot, R. Salut, L. Furfaro, and J. M. Dudley, “High aspect ratio nanochannel machining using single shot femtosecond Bessel beams,” Appl. Phys. Lett. , 97, 081102 (2010). [CrossRef]  

10. B. Wetzel, C. Xie, P.-A. Lacourt, J. M. Dudley, and F. Courvoisier, “Femtosecond laser fabrication of micro and nano-disks in single layer graphene using vortex Bessel beams,” Appl. Phys. Lett. 103, 241111 (2013). [CrossRef]  

11. X. Long, W. Zhao, R. Stoian, R. Hui, and G. Cheng, “Writing of stressed waveguides with tubular depressed cladding using femtosecond hollow beams,” Opt. Lett. 37, 3138–3140 (2012). [CrossRef]   [PubMed]  

12. S. Shiffler, P. Polynkin, and J. Moloney, “Self-focusing of femtosecond diffraction-resistant vortex beams in water,” Opt. Lett. , 36, 3834–3836 (2011). [CrossRef]   [PubMed]  

13. J. Durnin, “Exact solutions for nondiffracting beams. I. The Scalar theory,” J. Opt. Soc. Am. A 4, 651–654 (1987). [CrossRef]  

14. J. Durnin, J. J. Miceli Jr., and J. H. Eberly, “Diffraction-free beams,” Phys. Rev. Lett. 58, 1499 (1987). [CrossRef]   [PubMed]  

15. F. Gori, G. Guattari, and C. Padovani, “Bessel-Gauss beams,” Optics Commun. 64, 491–495 (1987). [CrossRef]  

16. M.A. Porras, A. Parola, D. Faccio, A. Dubietis, and P. Di Trapani, “Nonlinear Unbalanced Bessel Beams: Stationary Conical Waves Supported by Nonlinear Losses,” Phys. Rev. Lett. 93 (15), 153902 (2004). [CrossRef]   [PubMed]  

17. P. Polesana, A. Couairon, D. Faccio, A. Parola, M. A. Porras, A. Dubietis, A. Piskarskas, and P. Di Trapani, “Observation of Conical Waves in Focusing, Dispersive, and Dissipative Kerr Media,” Phys. Rev. Lett. 99 (22), 223902 (2007). [CrossRef]  

18. A. Lotti, D. Faccio, A. Couairon, D. G. Papazoglou, P. Panagiotopoulos, D. Abdollahpour, and S. Tzortzakis, “Stationary nonlinear Airy beams,” Phys. Rev. A 84021807 (2011). [CrossRef]  

19. P. Panagiotopoulos, D. Abdollahpour, A. Lotti, A. Couairon, D. Faccio, D. G. Papazoglou, and S. Tzortzakis, “Nonlinear propagation dynamics of finite-energy Airy beams,” Phys. Rev. A , 86, 013842 (2012). [CrossRef]  

20. A. Dubietis, E. Gaižauskas, G. Tamošauskas, and P. Di Trapani, “Light Filaments without Self-Channeling,” Phys. Rev. Lett. 92, 253903 (2004). [CrossRef]   [PubMed]  

21. F. Courvoisier, J. Zhang, M.K. Bhuyan, M. Jacquot, and J.M. Dudley, “Applications of femtosecond Bessel beams to laser ablation,” Appl. Phys. A 112 (1), 29–34 (2013). [CrossRef]  

22. P. Polynkin, C. Ament, and J.V. Moloney, “Self-Focusing of Ultraintense Femtosecond Optical Vortices in Air,” Phys. Rev. Lett. 111, 023901 (2013). [CrossRef]   [PubMed]  

23. A. Couairon, E. Brambilla, T. Corti, D. Majus, O. de J. Ramirez-Gongora, and M. Kolesik, “Practitioner’s guide to laser pulse propagation models and simulation,” Eur. Phys. J. Special Topics 199, 5–76 (2011). [CrossRef]  

24. G. Point, Y. Brelet, A. Houard, V. Jukna, C. Milián, J. Carbonnel, Y. Liu, A. Couairon, and A. Mysyrowicz, “Superfilamentation in Air,” Phys. Rev. Lett. 112, 223902 (2014). [CrossRef]   [PubMed]  

25. K.-M. Lee, C. M. Kim, S. A. Sato, T. Otobe, Y. Shinohara, K. Yabana, and T. M. Jeong, “First-principles simulation of the optical response of bulk and thin-film a-quartz irradiated with an ultrashort intense laser pulse,” J. Appl. Phys. 115, 053519 (2014). [CrossRef]  

26. N. Shcheblanov, “Numerical study of femtosecond laser interactions with dielectric materials: application to the definition of damage threshold of optical components,” Ph.D. Thesis, Jean Monnet University (2013).

27. V. Jarutis, R. Paskauskas, and A. Stabinis, “Focusing of Laguerre-Gaussian beams by axicon,” Optics Commun. 184, 105–112 (2000). [CrossRef]  

28. E. Gaizauskas, E. Vanagas, V. Jarutis, S. Juodkazis, V. Mizeikis, and H. Misawa, “Discrete damage traces from filamentation of Gauss-Bessel pulses,” Opt. Lett. 31, 80–82 (2006). [CrossRef]   [PubMed]  

29. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover PublicationsNew York, 1972).