Filamentation with nonlinear Bessel vortices

: We present a new type of ring-shaped ﬁlaments featured by stationary nonlinear high-order Bessel solutions to the laser beam propagation equation. Two different regimes are identiﬁed by direct numerical simulations of the nonlinear propagation of axicon focused Gaussian beams carrying helicity in a Kerr medium with multiphoton absorption: the stable nonlinear propagation regime corresponds to a slow beam reshaping into one of the stationary nonlinear high-order Bessel solutions, called nonlinear Bessel vortices . The region of existence of nonlinear Bessel vortices is found semi-analytically. The inﬂuence of the Kerr nonlinearity and non-linear losses on the beam shape is presented. Direct numerical simulations highlight the role of attractors played by nonlinear Bessel vortices in the stable propagation regime. Large input powers or small cone angles lead to the unstable propagation regime where nonlinear Bessel vortices break up into an helical multiple ﬁlament pattern or a more irregular structure. Nonlinear Bessel vortices are shown to be sufﬁciently intense to generate a ring-shaped ﬁlamentary ionized channel in the medium which is foreseen as opening the way to novel applications in laser material processing of transparent dielectrics.


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][5][6][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][14][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 higherorder 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.

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.

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]: The source terms are: where k 2 ⊥ = k 2 x +k 2 y denotes transverse wave vector, n 2 denotes the coefficient for nonlinear Kerr index, σ denotes the cross section for inverse Bremsstrahlung, τ c denotes the effective collision time, ρ denotes the density of electrons generated during propagation, ρ nt denotes the density of neutral molecules, W PI and U g denote optical field ionization rate and gap and P and J denote the nonlinear polarization and current sources. For the sake of simplicity, we use in this paper expression for the multiphoton absorption where W PI = σ K |E| 2K and σ K = β K /(Khω 0 ρ nt ), where β K is the multiphoton absorption coefficient. The quantity |E| 2 represents intensity and is expressed in units that makes the product n 2 |E| 2 dimensionless. The density of electrons in the conduction band is obtained by a simple rate equation in the form: where the terms on the right describe photoionization, ionization by avalanche and recombination respectively. In agreement with our approximation of a time independent electric field E(x, y, z), we consider the electron density ρ(x, y, z) as a function of intensity obtained by solving Eq. (4) for a pulse with maximum intensity |E(x, y, z)| 2 and fixed Gaussian pulse shape exp(−2t 2 /t 2 p ) with duration t p . This method has previously been shown to successfully capture plasma induced effects in multiple filamentation [24]. As will be shown below, the key physical effects leading to the generation of stable NBV are the optical Kerr effect and multiphoton absorption. Plasma induced effects play a secondary role. For this reason we used the rate Eq. (4) as a convenient way to calculate the order of magnitude of the density of electrons generated during propagation.

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.

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.
where r ≡ x 2 + y 2 , cos ϑ = x/r and θ denotes the half cone angle of the BV beam in the propagation medium. Most results in this work were obtained for a topological charge m = 1, however, we checked that the results are qualitatively similar for higher order BV (m ≥ 2). A Bessel beam with high cone half-angle of θ = 16.4 • in the medium is used, corresponding to a Bessel beam half angle of 25.6 • in air. 10 m   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(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.

Numerical simulation results
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). 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.

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: where dots means differentiation with respect to r. 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 4.

Linear Bessel vortices
We first look for stationary solutions to Eqs. (7,8) in the regime of linear propagation (n 2 = 0; β K = 0). The second equation (8) ensures that the quantity rqa 2 is preserved, which is satisfied by choosing q = 0. The first equation (7) then reads: Eq. (9) admits solutions in the form of Bessel functions a(r) = J ±m ( √ 2k 0 δ r). Linear conical vortices correspond to nonzero values of the topological charge m. The cone angle θ is defined by identifying the argument in the Bessel function: k ⊥ ≡ k 0 sin θ = √ 2k 0 δ . This defines the link between the cone angle and the axial phase shift: δ = (k 0 /2) sin 2 θ . For paraxial propagation (small cone angles), δ ≈ k 0 θ 2 /2.

Power conservation and stationarity principle
Radial integration of Eq. (8) (after multiplication by ra 2 ) results in the energy conservation law for weakly localized beams carrying infinite power: Interpretation of Eq. (10) provides direct insight into the mechanism responsible for propagation invariance of nonlinear weakly localized beams. The relation (10) indeed states that losses due to multiphoton absorption within a finite disk of radius r are compensated by the inward power flux from the peripheral part of the beam. This can be directly seen from the evolution equation for the power density |E(r, ϑ , z)| 2 , obtained by multiplying Eq. (6) by E * and by adding the complex conjugate equation: where J ≡ |E| 2 ∇ ⊥ (φ (r) + mϑ ). Propagation invariant beams satisfy ∂ z |E| 2 = 0. Therefore, the stationarity principle for nonlinear conical waves is equivalent to the condition that the divergence of the power flux density, proportional to the beam intensity and the phase gradient, is equal to the density of nonlinear losses. Application of the divergence theorem results in the global version of this principle: nonlinear power losses in any domain of the transverse plane are compensated by a power density flux through its boundaries. Equation (10) is nothing but an expression of this principle linking the radial component of the power flux J r = q(r)a 2 (r) to nonlinear losses within a disk of radius r. We note that the power density flux also comprises an azimuthal component J ϑ = ma 2 (r)/r. Nonlinear Besssel vortices are therefore stationary due to a helicoidal power flux, pointing toward the main lobe and in the azimuthal direction.

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 2 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 ±m ( √ 2k 0 δ r) 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 2 , β 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: where g K only depends on the order of multiphoton processes (g K = 1.67, 0.27, 0.19, 0.16, ... for K = 2, 3, 4, 5, ...).
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 −27 cm 3 /W 2 , nonlinear index n 2 = 3.45 × 10 −16 cm 2 /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 2 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 3 TW/cm 2 . 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 θ = sin −1 ( √ 2δ k 0 ), but with different weights α out and α in : Each term in Eq. (14) represent a wave associated with a power flux directed inward or outward.
In contrast with high-order Bessel beams, the inward power flux is larger than the outward power flux: |α in | > |α out |. As r → +∞, H ( 1 2 ) m (z) = 2/(πz)e ±i(z−mπ/2−π/4) and introduction of these asymptotics into Eq. (14) leads to expression for the tail of a NBV, at large r: where C ≡ 2|α in α out |/(|α in | 2 + |α out | 2 ) denotes the contrast of the Bessel tails. Introduction of the latter expression in the energy conservation equation (10) leads to the relation between amplitudes of inward/outward asymptotic waves and power losses: Therefore, multiphoton absorption determines both the unbalance between the inward and outward Hankel waves and the contrast of the Bessel tails varying from C = 1 for BVs with α in = α out to C = 0 for NBV on the frontier of allowed cone angles, with α out = 0. An accurate determination of the frontier of allowed cone angles, characterized by a vanishing contrast for NBV, requires an ideally infinite box size while our numerics are performed using large but finite box sizes. Error bars are therefore plotted in the transition region ( Fig. 4 light gray), where the amplitude of numerical solutions to Eqs. (7,8) decays more slowly than 1/ √ r, hence, the solutions do not qualify as NBV. The horizontal black line depicts the 16.4 • angle, while the black dots are the maximum intensity values of the NBV depicted in Fig. 5(c), note that the most intense NBV is in the transitional region.

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 2 , not only rings compress, but also the contrast increases.

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.

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.
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 ( 2k 0 (δ + δ n )r), where δ n ≡ 0.8 × k 0 I max n 2 /n 0 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 max = 1.8412/ 2k 0 (δ + δ 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.

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.   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.
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 semianalytical 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 multiphoton 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.

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.