Ismail Ouadghiri-Idrissi

John M Dudley

Francois Courvoisier

Control of spatial four wave mixing efficiency in Bessel beams using longitudinal intensity shaping

Diffraction-free Bessel beams have attracted major interest because of their stability even in regimes of nonlinear propagation and filamentation. However, Kerr nonlinear couplings are known to induce significant longitudinal intensity modulation, detrimental to the generation of uniform plasma or for applications in the processing of transparent materials. These nonlinear instabilities arise from the generation of new spatio-spectral components through an initial stage of continuous spectral broadening followed by four wave mixing. In this paper, we investigate analytically and numerically these processes and show that nonlinear instabilities can be controlled through shaping the spatial spectral phase of the input beam. This opens new routes for suppressing the nonlinear growth of new frequencies and controlling ultrashort pulse propagation in dielectrics.

I. INTRODUCTION

Diffraction-free Bessel beams are formed from a conical energy flow and yield a near-uniform intensity distribution along a line focus [START_REF] Durnin | Exact solutions for nondiffracting beams in the scalar theory[END_REF]. For high power laser pulses injected into transparent dielectrics, this beam structure can sustain quasi propagation-invariant regimes of filamentation, which are highly advantageous in reducing nonlinear distortion and instabilities during propagation. This has been shown to yield significant improvement in controlling energy deposition, and has been the subject of intense interest for generating long and uniform plasma channels [START_REF] Porras | Nonlinear unbalanced Bessel beams: Stationary conical waves supported by nonlinear losses[END_REF][START_REF] Polesana | Filamentation in Kerr media from pulsed Bessel beams[END_REF][START_REF] Dota | Femtosecond laser filamentation in condensed media with Bessel beams[END_REF] and creating high aspect ratio structures in micro-and nano-machining applications [START_REF] Bhuyan | High aspect ratio nanochannel machining using single shot femtosecond Bessel beams[END_REF][START_REF] Bhuyan | Laser micro-and nanostructuring using femtosecond Bessel beams[END_REF][START_REF] Garzillo | Optimization of laser energy deposition for single-shot high aspect-ratio microstructuring of thick BK7 glass[END_REF][START_REF] Courvoisier | INVITED] ultrafast laser micro-and nano-processing with nondiffracting and curved beams[END_REF][START_REF] Stoian | Ultrafast Bessel beams: advanced tools for laser materials processing[END_REF]. The application of Bessel beams in laser processing and filamentation has also been demonstrated for Bessel vortices [START_REF] Xie | Tubular filamentation for laser material processing[END_REF].

However, although the level of nonlinear distortion during Bessel beam propagation is greatly reduced compared to Gaussian beams at comparable peak intensities, residual nonlinear instability effects can still occur and induce significant oscillations of the on-axis intensity [START_REF] Polesana | Filamentation in Kerr media from pulsed Bessel beams[END_REF][START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF][START_REF] Andreev | Bessel beams of electromagnetic waves: selfeffect and nonlinear structures[END_REF]. This is clearly detrimental for the creation of longitudinally-uniform structures in material processing, and thus controlling these instabilities is of central importance to extend the applicability of Bessel beams.

Such nonlinear instabilities were first investigated in 1989 in association with the optical breakdown in gases and discussed in terms of cubic nonlinearity and plasma dynamics [START_REF] Andreev | Bessel beams of electromagnetic waves: selfeffect and nonlinear structures[END_REF]. Gadonas et al. subsequently investigated the relation between nonlinear instabilities in Bessel beams and the distortion of their spatial spectrum [START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF]. Using phase-matching arguments and considering four wave mixing interactions, they showed that a Bessel beam with radial wavevector k r0 can sustain Kerr self-action which deforms its spatial spectrum to generate two additional spectral components: an axial wave com- * francois.courvoisier@femto-st.fr ponent with k r = 0 and a secondary conical wave with k r = √ 2k r0 [START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF][START_REF] Pyragaite | Self-action of Bessel light beams in medium with large nonlinearity[END_REF]. Their numerical simulations highlighted the development of an intensity modulation of Bessel beams along the propagation direction which was explained by the interference of the conical beam with the axial wave component. Experimental signatures of this spectrum distortion was demonstrated in [START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF][START_REF] Pyragaite | Self-action of Bessel light beams in medium with large nonlinearity[END_REF][START_REF] Polesana | Observation of conical waves in focusing, dispersive, and dissipative Kerr media[END_REF].

In order to achieve stable propagation of Bessel beams, it is necessary to overcome these Kerr-induced instabilities [START_REF] Johannisson | Nonlinear Bessel beams[END_REF]. Porras et al investigated the stability of Bessel beams in the presence of nonlinear losses and showed that nonlinear instabilities can be significantly suppressed if multiphoton absorption prevails over Kerr nonlinearity [START_REF] Porras | Nonlinear unbalanced Bessel beams: Stationary conical waves supported by nonlinear losses[END_REF]. This condition can be satisfied for relatively high input power and high cone angles [START_REF] Porras | Nonlinear unbalanced Bessel beams: Stationary conical waves supported by nonlinear losses[END_REF][START_REF] Couairon | Ultrashort laser pulse filamentation with Airy and Bessel beams[END_REF], and experimental demonstration was performed in [START_REF] Porras | Nonlinear unbalanced Bessel beams: Stationary conical waves supported by nonlinear losses[END_REF][START_REF] Polesana | Observation of conical waves in focusing, dispersive, and dissipative Kerr media[END_REF][START_REF] Gaizauskas | Discrete damage traces from filamentation of Gauss-Bessel pulses[END_REF]. However, since this imposes strong constraints on the geometrical and physical properties of the input Bessel beam, this approach is not suitable in all applications. In other work, Polesana et al. investigated the effect of the input injection condition to the nonlinear medium [START_REF] Polesana | Filamentation in Kerr media from pulsed Bessel beams[END_REF] and showed that the Kerr-induced instabilities can be attenuated if the Bessel beam is progressively formed inside the Kerr medium. In contrast, if the Bessel beam is formed prior entering the medium, it was shown to exhibit significant instabilities and quasi-periodic intensity modulation along the propagation direction.

Recently, we have shown that these instabilities can be significantly mitigated by appropriate control of the intensity evolution of Bessel beams along the propagation distance [START_REF] Ouadghiri-Idrissi | Controlling nonlinear instabilities in Bessel beams through longitudinal intensity shaping[END_REF]. In particular, we qualitatively identified the role of the spatial spectral phase in reducing the efficiency of Four Wave Mixing (FWM) and the growth of related nonlinear instabilities. In this paper we present an extended analysis of this problem, and through both analysis of the FWM process and numerical simulations, we obtain significant new insights into the physical origins of these instabilities, and identify particular quantitative parameter regimes in which they can be suppressed. Note that our approach is complementary to the work in Refs 2,[START_REF] Polesana | Filamentation in Kerr media from pulsed Bessel beams[END_REF][START_REF] Porras | Underlying conservation and stability laws in nonlinear propagation of axicon-generated Bessel beams[END_REF][START_REF] Porras | Stabilization of vortex beams in Kerr media by nonlinear absorption[END_REF] where propagation is analyzed in terms of a linearized stability analysis, and where nonlinear losses are shown to give rise to an attractor (nonlinear unbalanced Bessel beam). Here, in contrast, we neglect nonlinear losses to isolate the effect of input beam shaping on Four-wave mixing efficiency and growth of the intensity oscillations in the Bessel beam.

This paper is organized as follows. We first develop in Section II an analytical model of four wave mixing that allows us to describe the cascade of nonlinear effects that lead to the generation of new spatial spectral components. In section III, we examine the dynamics of the growth of new spatial spectral components and show that an initial phase of spectral broadening is the major driver for instabilities. In section IV, we develop a reduced model with only the essential terms responsible for spectral broadening which allows us to take into account the spectral width -in the spatial domain -of Bessel beams. In section V, we will use this model to discuss the dependence of nonlinear instabilities on the precise injection condition of an incident Bessel beams with respect to a nonlinear medium (referred to as the "soft" and "abrupt" transition cases) considering in particular the influence of the initial spectral phase. Finally in section VI, we consider previous numerical results studying the nonlinear propagation of Bessel beams with shaped on-axis intensity profiles [START_REF] Ouadghiri-Idrissi | Controlling nonlinear instabilities in Bessel beams through longitudinal intensity shaping[END_REF], with our modelling allowing us to understand why specific input spectral phases lead to the reduction of nonlinear instabilities.

II. FOUR WAVE MIXING IN BESSEL BEAMS

In 1996, Tewari et al introduced a theoretical model to describe Third Harmonic Generation (THG) in Bessel beams [START_REF] Tewari | Theory of thirdharmonic generation using Bessel beams, and self-phasematching[END_REF]. Our theoretical model follows the same approach. However, instead of THG nonlinear terms, we consider four wave mixing (FWM) interactions to study nonlinear spectral distortions in Bessel beams related to Kerr-induced instabilities. Importantly, our focus is on the growth of spatial frequencies k r in a monochromatic model such that we do not consider dynamics related to the growth of new temporal frequencies ω. However, we stress that the formalism developed here could be extended to describe the full picture of simultaneous nonlinear dynamics in both spatial and temporal domains.

We note the monochromatic approximation has proven successful in interpreting experimental results obtained with pulsed beams with typical pulse durations of some hundreds fs, peak powers in the order of TW.cm -2 , cone angles of a few degrees [START_REF] Porras | Nonlinear unbalanced Bessel beams: Stationary conical waves supported by nonlinear losses[END_REF][START_REF] Polesana | Filamentation in Kerr media from pulsed Bessel beams[END_REF][START_REF] Xie | Tubular filamentation for laser material processing[END_REF][START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF][START_REF] Polesana | Observation of conical waves in focusing, dispersive, and dissipative Kerr media[END_REF] as we will use here as numerical examples, for the nonlinear propagation of ultrafast Bessel beams in fused silica. This approximation is valid as long as temporal reshaping is not too strong.

We first compute the third order nonlinear polarization and select the relevant nonlinear terms which have direct impact on the generation of new spectral components. We will then include these terms in the Helmholtz equation which we study analytically and numerically.

A. Third order nonlinear polarization

Nonlinearity in Kerr media appears through the third order nonlinear polarization P N L , written as:

P N L = ε 0 χ (3) EEE (1)
where ε 0 is the dielectric permittivity in vacuum, χ (3) is the third order susceptibility and E is the total electric field. We consider monochromatic waves, oscillating at (temporal) frequency ω 0 . We consider four waves E j , of amplitude A j , linearly polarized along the same axis x and possessing different longitudinal spatial frequencies k zj . The total electric field is then written:

E = 1 2 x 4 j=1 A j exp[i(ω 0 t -k zj z)] + c.c (2)
If we substitute Eq. (2) in Eq. ( 1), we find that the nonlinear polarization consists of many terms covering all possible nonlinear interactions including terms responsible for Third Harmonic Generation (THG) and those preserving the same temporal frequency ω 0 . Neglecting THG terms because we investigate only the ω 0 components, the nonlinear polarization can be written as:

P N L, ω0 = 3 8 ε 0 χ (3) x p N L, ω0 e iω0t + c.c (3)
where

p N L, ω0 = 4 j=1 |A j | 2 A j e -ikzj z + 4 j=1,j =m 2 |A j | 2 A m e -ikzmz + 4 j=1,j =m A 2 j A * m e -i(2kzj -kzm)z (4) + 4 j=1,j<l<m 2 A j A l A * m e -i(kzj +k zl -kzm)z + A j A * l A m e -i(kzj -k zl +kzm)z + A * j A l A m e -i(-kzj +k zl +kzm)z
The first two terms of Eq. (4) describe Self-Phase-Modulation (SPM) and Cross-Phase-Modulation (XPM).

The other terms describe Four-Wave-Mixing (FWM), with the first term A 2 j A * m e -i(2kzj -kzm)z being the degenerate process where the two pump waves are identical while the other terms being non-degenerate processes.

Our aim is to study the generation of new spectral components and their evolution along the propagation direction. In the following, we consider the waves E 1 and E 2 as the high-intensity pump waves whereas E 3 and E 4 correspond to the signal and idler waves respectively.

We describe the evolution of the spatial spectrum with propagation as a cascading process where E 3 is generated first by cross interaction of the pump waves. Then, fourwave mixing will amplify E 3 and E 4 . This approach is consistent with the scaling performed by Gadonas et al [START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF].

We separate the terms of the nonlinear polarization in Eq. ( 4) according to the content of the exponents (which is equivalent to momentum conservation). This leads to:

p (1) N L, ω0 e ikz1z = |A 1 | 2 + 2 |A 2 | 2 A 1 (5a) p (2) N L, ω0 e ikz2z = |A 2 | 2 + 2 |A 1 | 2 A 2 (5b) p (3) N L, ω0 e ikz3z =2 |A 1 | 2 + |A 2 | 2 A 3 + A 2 1 A * 2 e -i(2kz1-kz2-kz3)z + A 2 2 A * 1 e -i(2kz2-kz1-kz3)z (5c) + A 2 1 A * 4 e -i(2kz1-kz3-kz4)z + A 2 2 A * 4 e -i(2kz2-kz3-kz4)z + 2 A 1 A 2 A * 4 e -i(kz1+kz2-kz3-kz4)z + A 1 A * 2 A 4 e i(-kz1+kz2+kz3-kz4)z + A * 1 A 2 A 4 e i(kz1-kz2+kz3-kz4)z p (4) N L, ω0 e ikz4z =2 |A 1 | 2 + |A 2 | 2 A 4 (5d) + A 2 1 A * 3 e -i(2kz1-kz3-kz4)z + A 2 2 A * 3 e -i(2kz2-kz3-kz4)z + 2 A 1 A 2 A * 3 e -i(kz1+kz2-kz3-kz4)z + A 1 A * 2 A 3 e -i(kz1-kz2+kz3-kz4)z + A * 1 A 2 A 3 e -i(-kz1+kz2+kz3-kz4)z
For p

N L, ω0 and p

N L, ω0 , we have neglected all contributions of A 3,4 as they are much smaller than A 1,2 . Similarly, for p

N L, ω0 and p (4) N L, ω0 , we have neglected all terms scaling with A 2 3,4 . We note that the two contributions A 2 1 A * 2 e -i(2kz1-kz2-kz3)z and A 2 2 A * 1 e -i(2kz2-kz1-kz3)z (first line of Eq. (5c)) will be particularly important for the rest of this paper. They arise from the XPM-like interaction with the pump, and scale as third power of the pump field and generate non-phase matched spectral broadening around the pump. The other 5 last terms scale with the second power of the pump and will contribute to amplification via four-wave mixing processes. With this expression for the nonlinear polarization, we can now use the wave equation to derive the evolution of the fields.

In the expression p

(3) N L, ω0 , the first two terms are crossphase modulation terms that do not contribute efficiently to new spectral frequency generation in our case. The third and fourth, just mentioned above, are in third power of the pump while the rest of the nonlinear polarization terms are in second order of the pump. In the following, we will see that third and fourth terms will generate a crucial non-phase matched spectral broadening around the pump while the other terms will contribute to amplification via four-wave mixing processes. As the terms corresponding to the spectral broadening are more efficient, these will be the drivers for the cascade of four wave mixing. Now, we can use the wave equation to derive the evolution of the fields.

B. Evolution of the spatial spectrum along the propagation

The starting point is the scalar wave equation describing the full field 4 j=1 E j . We separate the full-field wave equation into four independent equations by using the separation approach described above:

∆E j - ε r c 2 ∂ 2 E j ∂t 2 = µ 0 ∂ 2 P (j) N L ∂t 2 (6)
We now follow the approach in Ref. [START_REF] Tewari | Theory of thirdharmonic generation using Bessel beams, and self-phasematching[END_REF]: for each of the waves, we will only consider ideal Bessel beams defined by the J 0 Bessel function. To further simplify the analysis, we consider that the pump amplitude is undepleted with propagation. This, of course, implies that our analysis will only be valid for propagation distances shorter than the typical depletion scale length (in the examples shown, typically on the order of 1000 µm). The envelopes of the four interacting waves A j are then written as follows:

A j = a j (z) J 0 (k rj r) (7)
where k rj is the transverse spatial frequency of the envelope A j . Since we developed an expression of P (j)

N L with the same form as E j , we can develop Eq. ( 6) without the complex conjugate terms [START_REF] Tewari | Theory of thirdharmonic generation using Bessel beams, and self-phasematching[END_REF]. It becomes:

2ik zj ∂a j (z) ∂z J 0 (k rj ) = k 2 0 ε 0 p (j) N L,ω0 (8)
We multiply both parts of Eq. ( 8) by rJ 0 (k r4 r) and integrate over 0 to r f which denotes the upper integration limit such as r f 1/k r0 . In other words, we perform Hankel transformation and thus study the evolution of the waves in Fourier space. The upper integration boundary r f is chosen finite so as to avoid infinite integrals or to avoid the introduction of apodization functions such as in Ref. [START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF]. According to Ref. [START_REF] Olver | NIST digital library of mathematical functions[END_REF]:

r f 0 rJ 2 0 (k rj r)dr = r 2 f /2(J 2 0 (k rj r f ) + J 2 1 (k rj r f ))
Then using the asymptotic expressions of both Bessel functions [START_REF] Olver | NIST digital library of mathematical functions[END_REF], i.e.

J α (k rj r) = 2 πk rj r cos k rj r - απ 2 - π 4 , α = 0, 1
This integral can be approximated to:

r f 0 rJ 2 0 (k rj r)dr ≈ r f πk rj (9)
which is proportional to 1/k rj , in agreement with Ref. [START_REF] Tewari | Theory of thirdharmonic generation using Bessel beams, and self-phasematching[END_REF].

Using the simplifications described in the previous section:

A 1 ≈ A 2 ≈ A 0 = a 0 J 0 (k r0 r) and k z1 ≈ k z2 ≈ k z0 ,
where k z0 = k cos θ with θ the cone angle of the pump Bessel beam, our system of equations becomes:

∂a 0 (z) ∂z = -3iC 0 tan(θ 0 )I (00)
T P M a 0 (z) (10a)

∂a 3 (z) ∂z = -iC 0 tan(θ 3 ) 4I (33)
T P M a 3 (z) + 2I

(03)

T P M e -i(∆k03+Φ N L 0 )z (10b)
+ 4I

(34)

T P M a * 4 (z)e -i(∆k034+2Φ N L 0 )z + 4I (34)
T P M a 4 (z)e i∆k34z ∂a 4 (z) ∂z = -iC 0 tan(θ 4 ) 4I (44)
T P M a 4 (z) (10c)
+ 4I

(34)

T P M a * 3 (z)e -i(∆k034+2Φ N L 0 )z + 4I (34)
T P M a 3 (z)e -i∆k34z
where

C 0 = π k 2 r f n 2 n 0 I 0 .
Here we write the solution of Eq.

(10a) as: a 0 (z) = √ I 0 exp -iΦ N L 0 z where I 0 is the peak input intensity. Solving Eq. (10a) we obtain:

Φ N L 0 = 3C 0 tan(θ 0 )I (00) T P M (11)
The Transverse Phase Matching (TPM) integrals I (j) T P M are defined as:

I (00) T P M = r f 0 J 4 0 (k r0 r) r dr (12a) I ( 03
)
T P M = r f 0 J 3 0 (k r0 r) J 0 (k r3 r) r dr (12b) I ( 34
) T P M = r f 0 J 2 0 (k r0 r) J 0 (k r3 r) J 0 (k r4 r) r dr (12c) I ( 33
) T P M = r f 0 J 2 0 (k r0 r) J 2 0 (k r3 r) r dr (12d) I ( 44
)
T P M = r f 0 J 2 0 (k r0 r) J 2 0 (k r4 r) r dr (12e)
The indices relate to the waves involved in the last two Bessel functions in the integrals. Longitudinal wavevector mismatch terms are defined as: ∆k 03 = k z0 -k z3 , ∆k 034 = 2k z0 -k z3 -k z4 and ∆k 34 = k z3 -k z4 . The longitudinal phase matching conditions (when the wavevector mismatch defined above equals zero) are the same as those reported in Ref. [START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF]: the first one corresponds to FWM interaction of the proposed "first approximation" of the nonlinear Schrödinger Equation (NLSE). For this approximation, the same TPM integral as I

T P M was also defined. Similarly, the other two Longitudinal Phase Matching (LPM) conditions were defined for the proposed "second approximation" of the NLSE along with I (34)

T P M .
Hence, the signal and idler waves of our model can be assimilated to solutions of the first and second order approximations of the NLSE of Ref. [START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF]. Our target here is to build a fully explicit model where the mechanism actually driving the generation of new spectral components can be analytically identified. However, before integrating Eqs. (10b) and (10c), we first show numerical results so that we can later compare analytical results with numerical modelling results of the full NLSE.

III. DYNAMICS OF THE GROWTH OF NEW SPECTRAL COMPONENTS

A. Numerical model and results

Our numerical simulations are based on the nonlinear Schrödinger equation (NLSE) given in Refs. [START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF][START_REF] Couairon | Practitioner's guide to laser pulse propagation models and simulation[END_REF] for a monochromatic beam propagating in a Kerr medium :

∂A ∂z = i 2k ∆ ⊥ A + ik n 2 n 0 |A| 2 A ( 13
)
where A is the linearly polarized complex amplitude of the laser electric field, ∆ ⊥ = 1/r∂/∂r + ∂ 2 /∂r 2 is the transverse Laplacian operator, r and z are the radial and axial coordinates, k is the wavevector in the medium, n 0 and n 2 are the linear and nonlinear refractive indices. Parameters of our simulations are given in Table I and correspond to the realistic propagation of a high-intensity pulse in fused silica. Since nonlinear instabilities stem mainly from Kerr nonlinearities, we neglected other nonlinear effects (particularly nonlinear losses which are known attenuate nonlinear instabilities [START_REF] Porras | Nonlinear unbalanced Bessel beams: Stationary conical waves supported by nonlinear losses[END_REF][START_REF] Polesana | Filamentation in Kerr media from pulsed Bessel beams[END_REF]) so as to isolate the effect of intensity shaping on the control of nonlinear instabilities in Bessel beams.

The input field (a Bessel-Gauss (BG) beam) is modeled by a Gaussian beam with a spatial phase characterizing the axicon conical focusing: A BG (r, z = 0) = A 0 exp -r 2 /w 2 0 -k r sin(θ) where w 0 is the input Gaussian beam waist [START_REF] Polesana | Filamentation in Kerr media from pulsed Bessel beams[END_REF][START_REF] Jarutis | Focusing of Laguerre-Gaussian beams by axicon[END_REF].

The results of the integration of Eq. ( 13), based on a split-step algorithm are shown in Fig. 1. In Fig. 1(a), we plot the evolution of the intensity as a function of radial distance r and propagation distance z. We observe that the beam intensity undergoes longitudinal modulation not only along the central core but also in the peripheral rings. The on-axis intensity, i.e. the intensity I(r = 0, z), is shown in Fig. 1(b), where the intensity oscillations, with a period of ∼180.5 µm, are clearly apparent.

In Fig. 1(c), the spatial spectrum | Ã(k r , z)| 2 is plotted as a function of the propagation distance. We recall that the spatial spectrum of a Bessel beam, in the linear regime, is in the form of a ring centered around the Bessel transverse frequency k r0 = 0.8 µm -1 [START_REF] Jarutis | Focusing of Laguerre-Gaussian beams by axicon[END_REF]. We display the evolution of a cross-section along the propagation distance. After an initial stage of spectral broadening around the central frequency k r0 , we notice the generation of two particular spectral components at k r ≈ 0 and k r ≈ 1.5k r0 . These components are respectively referred to as the axial wave and outer ring (k r ∼ √ 2k r0 ) as mentioned in the introduction and reported in previous works [START_REF] Gadonas | Self-action of Bessel beam in nonlinear medium[END_REF][START_REF] Polesana | Observation of conical waves in focusing, dispersive, and dissipative Kerr media[END_REF][START_REF] Ouadghiri-Idrissi | Controlling nonlinear instabilities in Bessel beams through longitudinal intensity shaping[END_REF]. It is the interference of the input Bessel beam with the two new spectral components which generates the oscillations observed on the on-axis intensity distribution. (The interference pattern was initially interpreted only as interference between the Bessel beam and the axial wave, but we note that the secondary wave also generates interference with precisely the same period).

In more detail, Fig. 1(c) shows two regimes. From a propagation distance range from 0 to z ∼ 2600 µm, the spatial spectrum progressively broadens around the central frequency. It is only for further propagation distances that the growth of the axial wave and outer ring is efficient. We specifically show the evolution of these spec-tral components in Fig. 1(d). In Fig. 1(c), we also note parabolic-like structures for spatial frequencies around k r0 . Those were not discussed in previous literature and our analytical model will allow us to explain them.

B. Analysis using the Four Wave Mixing model

Here we show that the main characteristics of the first and second regime can be qualitatively described using the Four-Wave Mixing model developed in the previous section (Eqs. (10b,10c)).

For very short propagation distances, since the amplitude of a 4 is near 0, as discussed in section II A above, we first neglect the terms in a 4 in the expression of the evolution of a 3 (Eq. ( 10b)).

This becomes:

∂a 3 (z) ∂z = -iC 0 tan(θ 3 ) 4I (33)
T P M a 3 (z) ( 14)

+ 2 I 0 I (03) T P M e -i(kz0-kz3+Φ N L 0 )z Using Φ N L 3 = 4C 0 tan(θ 3 )I (33)
T P M and ∆k ef f

03 = ∆k 03 + Φ N L 0 -Φ N L
3 , a 3 is given by:

a 3 (z) = -2iC 0 I 0 tan(θ 3 )I (03)
T P M e -i(∆k ef f 03 z/2+Φ N L 3 z) × z sinc ∆k ef f 03 z/2 (15)
This result will be very important in the following sections of the paper. Then we can compute the evolution of the idler wave with Eq. (10c), using the expression of the signal wave a 3 computed above.

Similarly, using, Φ N L 4 = 4C 0 tan(θ 4 )I

T P M and ∆k ef f

04 = k z0 -k z4 + Φ N L 0 -Φ N L 4 , ∆k ef f 034 = ∆k 034 + 2Φ N L 0 -Φ N L 3 -Φ N L
4 , and ∆k ef f 34 = ∆k 34 + Φ N L 3 -Φ N L 4 . a 4 is given by:

a 4 (z) = -8iC 2 0 I 0 I (03) T P M I (34) T P M tan(θ 3 ) tan(θ 4 ) ∆k ef f 03 e -iΦ N L 4 z × 2 exp -i∆k ef f 04 z/2 z sinc ∆k ef f 04 z/2 -exp -i∆k ef f 034 z/2 z sinc ∆k ef f 034 z/2 -exp -i∆k ef f 34 z/2 z sinc ∆k ef f 34 z/2 (16)
1. First regime: Spatio-Spectral broadening characteristics Using Eq. ( 15), we show in Fig. 2(a) the evolution of the intensity of a 3 (z) along the propagation direction for different values of k r3 . In the figure, the nonlinear phases Φ N L 0 and Φ N L 3 were evaluated for an intensity of 18 TW/cm 2 as for the peak power in Fig. 1. We also choose the same pump cone angle θ = 4 • as in our numerical simulation of the NLSE. This regime qualitatively reproduces the first growth stage discussed earlier (propagation from 0 to z = 2600 µm). The wave (a 3 ) oscillates for frequencies that are away from the central frequency k r0 and exhibits parabolic-like fringe structure in the k r -z space, consistent with NLSE simulation results (Fig. 1(c)). These parabolic structures oscillate with a frequency-dependent period: p(k r3 ) = 2π/∆k ef f 03 where we can neglect the nonlinear phase to obtain:

p(k r3 ) 2π |k z0 -k z3 | 4πk |k 2 r0 -k 2 r3 |
This result is in very good quantitative agreement with the oscillation period observed in Fig. 1. For instance, for k r3 = 0.5k r0 , the expression above gives a period of 302 µm, while the numerical result is 295 µm.

Second regime: Growth of new spectral components and interference pattern

After this initial spectral broadening, numerical results show that both the axial wave and outer ring will be amplified for z > 2600 µm. Now, we show that in this second regime, an axial wave (kr 0) and an outer ring (kr (2)k r0 ) are amplified. This amplification is described by the expression of a 4 .

The transverse phase matching integral I

T P M peaks for k r3 0. Therefore, we plot in Fig. 2(b) the frequency-resolved evolution of |a 4 (z)| 2 for a signal wave a 3 at k r3 close to zero. We observe the amplification of an axial wave (k r4 0) and of an outer ring k r4 (2)k r0 . The axial wave arises from the last term in Eq. ( 16), for which ∆k ef f 34 = 0 and the outer ring arises from the second to last term in the same equation, for which ∆k ef f 034 = 0. We note that our description does not yet take into account the complete set of signal waves a 3 that are continuously generated in the first spectral broadening stage as the analysis would be extremely laborious.

In summary, we have shown that our model can explain detailed features of the nonlinear propagation of Bessel-Gauss beams. We have seen that it is the initial broadening stage (generating the wave a 3 ) that determines the efficiency of the FWM-induced amplification in the subsequent stage. In the next section, we will expand our theory to take into account the spectral phase of the pump. To simplify our analysis, we will restrict ourselves to the first broadening regime.

IV. REDUCED MODEL

Here we will describe the generation of a 3 when taking into account the fact that the input pump Bessel beam is spectrally extended. We will see in the next sections how the spectral phase impact the growth of a 3 .

We restart our analysis from Eq. ( 14). In this expression, the first nonlinear term, -4iC 0 I (33) T P M tan(θ 3 )a 3 , corresponds to cross-phase modulation and is much weaker than the second term. Therefore, to simplify our analysis, we neglect the first term in our reduced model. We also drop out the nonlinear phase terms. Now, we take into account a pump beam described by: A 0 = √ I 0 dk r S(k r )J 0 (k r r)e iφ (kr) where S(k rj ) stand for the amplitude distribution of the spectral components of the pump and signal waves. The complex spatial spectra are given by: S(k r ) = S j (k r )exp[iφ j (k r )]. This way, it is possible to take into account the input spectral distribution by associating with each spectral component the corresponding amplitude and phase values. Equation ( 14) becomes:

∂a 3 (z) ∂z = -2iC 0 I 0 tan(θ 3 )e ikz3z r f 0 J 0 (k rj r) × dk r S(k r )J 0 (k r r)e iφ(kr)-ikzz 2 ( 17
)
dk r S(k r )J 0 (k r r)e -iφ(kr)+ikzz r dr

Here, we have a triple integral over the transverse spatial frequency k r where k z = k 2 -k 2 r . To make this expression easier to analyze analytically, we define for each of these integrals a different parameter, i.e. k rj , k rl and k rm . Our equation can then be written as:

∂a 3 (z) ∂z = -2iC 0 I 0 tan(θ 3 ) k 0 dk rj S(k rj ) k 0 dk rl S(k rl ) (18) × k 0 dk rm S(k rm ) I (jlm) T P M exp(i∆Φ jlm -i∆k jlm z)
where I

(jlm)

T P M = r f 0 J 0 (k rj r)J 0 (k rl r)J 0 (k rm r)J 0 (k r3 r)rdr, ∆Φ jlm = φ(k rj ) + φ(k rl ) -φ(k rm ) and ∆k jlm = k zj + k zl -k zm -k z3 .
Now, this model will allow us to predict how the efficiency of the first spectral broadening stage is affected by the spectral phase of the input pump beam. In the next two sections, we will use our model to understand two different cases of the literature regarding the propagation of Bessel beams in Kerr media.

V. SOFT OR ABRUPT INPUT CONDITIONS

Previous work by other groups [START_REF] Polesana | Observation of conical waves in focusing, dispersive, and dissipative Kerr media[END_REF][START_REF] Sogomonian | Spatial distortions of a Bessel beam in a Kerr-type medium[END_REF][START_REF] Dubietis | Axial emission and spectral broadening in self-focusing of femtosecond Bessel beams[END_REF] showed experimentally and numerically that an abrupt transition between linear and nonlinear propagation of an intense Bessel beam yields efficient generation of outer ring and axial wave components. In contrast, this is much less efficient when the Bessel beam is smoothly forming into the nonlinear medium (see Fig. 1 in Ref. [START_REF] Dubietis | Axial emission and spectral broadening in self-focusing of femtosecond Bessel beams[END_REF]). These two input conditions are respectively referred to as soft and abrupt input conditions.

In Fig. 3, we show numerical results of the NLSE, that includes only Kerr effect, as described in section III A.

We use the same parameters as in table I except for the input power which was reduced to P in = 31.2 MW, corresponding to a beam peak intensity of I max = 9 TW/cm 2 . In the linear regime, this Bessel beam reaches its peak intensity at z = 2160 µm. For the abrupt input condition, the nonlinear medium starts at this point, whereas for the soft input conditions, the nonlinear propagation starts at z = 0.

We compare the evolution of on-axis intensity and spatial spectrum for soft and abrupt input conditions for the same Bessel beam. We see, in agreement with literature, that the abrupt transition generates pronounced on-axis intensity modulation, in stark contrast with soft input conditions. The evolution of spatial spectra can be compared from Figs. 3(b-c). For the abrupt input condition (c), the spectral intensity of the axial wave and outer ring components quickly grows, with an intensity two orders of magnitude higher than in the case of the soft input condition (b). We note that the case of soft input condition physically corresponds to the same case as previous sections. The oscillations of the on-axis intensity are reduced because of a smaller input power.

We will now use our reduced model to understand the different behaviors. We demonstrate that the spectral -c) the corresponding spatial spectra evolution along the propagation distance (logarithmic scale) for both cases respectively. Note that the intensity of newly-generated spectral components in the case of abrupt transition is two order of magnitude higher than the one reached in case of soft transition. Parameters are the presented in table I and Pin = 31.2 MW.

phase distribution largely differs between soft and abrupt input conditions, which impacts on the first stage of spectral broadening.

We compare in Fig. 4(a) the input spectral amplitude (green) and phase (red) corresponding to both cases. The distance z is the relative distance to the transition point between the linear and nonlinear medium. While the input amplitude is naturally the same (in the form of a Gaussian), we notice that the phase distribution is much steeper in the soft input condition case. In contrast, it is quasi-flat within the spectral range around the peak of the amplitude.

A quasi-flat spectral phase implies that spectral components composing the pump wave are nearly in-phase. Qualitatively, if each of these spectral components interacts according to the four-wave mixing process described above, then each signal wave generated from these interactions will be in phase with the others. The resulting signal wave at a given frequency will then be made of constructive interference between all these waves, which explains the very rapid growth of spectral components at about all frequencies around the central one in the abrupt input condition. In contrast, a steep spectral phase profile, which implies out-of-phase spectral components, leads to partially-destructive interferences and the axial wave will then be weaker. Now, we use Eq. ( 18) to obtain an analytical explanation. The triple integral over the spectrum is unfortunately too heavy and for sake of simplicity, we restrict the pump beam to only two spectral components defined at k ra and k rb , such that k ra ≈ k rb ≈ k r0 , with φ a and φ b being their respective input spectral phases and consider they have the same amplitude. The spectral distribution of each of interacting waves can then be written as

S(k r ) = δ(k r -k r0 ) [exp(iφ a ) + exp(iφ b )].
The signal wave intensity I 3 (z) = |a 3 (z)| 2 is then found to be proportional to:

I 3 (z) ∝ tan(θ 3 ) I (03) T P M 2 z 2 sinc 2 [∆k 03 z/2] (19)
× [1 + cos(φ a -φ b )] 3 .
The growth of the signal wave intensity is proportional to the term z 2 sinc 2 [∆k 03 z/2] which indicates the above discussed oscillating behavior. Of particular interest, I 3 (z, k z3 ) is proportional to the cube of the phasedependent term [1 + cos(φ a -φ b )] which shows that nonzero phase difference quickly reduces the peak value of the oscillations and thus decreases the magnitude of the axial wave seed.

In Fig. 4(b), we numerically solve our reduced FWM equation ( 19) over a more realistic case of a pump spectrum composed of five frequencies in the range k r /k r0 ∈ [0.96 -1.04]. We plot the evolution of the signal wave spectral intensity along propagation and compare results in soft and abrupt input conditions. For each case, the input spectral phase distribution is extracted from Fig. 4. In the soft input condition, initial spectral broadening is very weak and the generated frequencies are close to the central one. Notice that it gradually extends to more frequencies for longer propagation distances in cascadedlike fashion, in good qualitative agreement with numerical simulations of Fig. 3 conditions, the spectrum very rapidly broadens, again in agreement with Fig. 3(c). Therefore, we conclude that the spectral phase is an effective control parameter for the initial broadening regime and therefore a control parameter for the instabilities occurring in the second stage.

VI. CONTROL OF NONLINEAR INSTABILITIES USING SHAPED INTENSITY PROFILES

In this section, we will interpret previous numerical results on the control of nonlinear instabilities depending on the initial intensity rise in Bessel beams propagating inside nonlinear Kerr media. In Ref. [START_REF] Ouadghiri-Idrissi | Controlling nonlinear instabilities in Bessel beams through longitudinal intensity shaping[END_REF], we compared the nonlinear propagation of three Bessel beams with different on-axis intensity profiles in the linear propagation regime. Their peak maximal intensity is chosen to be the same in order to study the effect of the initial intensity rise on the growth of nonlinear instabilities. These three target intensity profiles are depicted in Fig. 5(a) and are described as follows: the first profile, denoted profile 1, is that of a conventional Bessel-Gauss beam (green dashed line) identical to that we have used previously. The second profile (profile 2) consists of a linear leading edge followed by a flat-top intensity and parabolic decaying trailing edge (blue dotted line). Profile 3 is identical to profile 2 except that it exhibits a parabolic intensity rise instead of a linear ramp (red solid line). Numerical parameters are the same as in table I.

We used the same nonlinear propagation equation described by Eq. ( 13). The evolution of the on-axis intensity of these Bessel beams is presented in Fig. 5(b). Compared to the case of the BG beam (profile 1), the two other beams present oscillations of the central core intensity that has the same period, but their amplitudes strongly differ. Although profiles 2 and 3 only differ in the initial intensity rise, the oscillation amplitude is significantly smaller for profile 3.

In the spectral domain (Fig. 5(c)), the observed weak on-axis intensity oscillations in case of profile 3 correspond to weak intensity growth of the axial wave , below -40 dB up to z = 4000 µm and its amplification remains around -30 dB afterwards. In addition, in the initial stage of spectral broadening, in the range z = [0 -2000] µm, we observe strong oscillations in the spectrum for profiles 1 and 2 while this oscillating behavior is initially absent in the case of profile 3 and only appears at a propagation distance around z = 3000 µm. Following the appearance of these oscillations, both the axial wave and outer ring increase in intensity which indicates that FWM processes become active past this propagation point. However, their growth remains noticeably weaker compared to the other two Bessel beams.

In Ref. [START_REF] Ouadghiri-Idrissi | Controlling nonlinear instabilities in Bessel beams through longitudinal intensity shaping[END_REF], we reported that it is not because profile 3 has initially a low-intensity zone that nonlinear instabilities are weaker. It is, as it is the case with soft and abrupt input conditions, the input spectral phase that is main factor of influence. Figure (5,d) shows the spectral phase distributions corresponding to profiles 1, 2, and 3. The BG beam (profile 1) has a spectral phase in the form a linear ramp in the range (k r /k r0 ∈ [0.98 -1.02]) and is flat outside this range. The spectral phase distribution of profile 2 is similar to that of profile 1 except that it exhibits very weak oscillations in the tails of the spectrum. In the case of profile 3, however, while the spectral phase also takes the form of a linear ramp in the range k r /k r0 ∈ [0.98 -1.02], it varies significantly outside, with a quasi-linear ramp.

The spectral domain of influence of the phase is the range k r /k r0 ∈ [0.8 -1.2] where a change in the phase implies a deviation of the on-axis intensity by more than 5 % in linear propagation regime [START_REF] Ouadghiri-Idrissi | Controlling nonlinear instabilities in Bessel beams through longitudinal intensity shaping[END_REF].

As in previous section, we will interpret the reduced nonlinear instabilities for profile 3 by the phase mismatch in the very broad pump spectrum. Using our reduced FWM model described by Eq. ( 18), we consider the contribution of 27 spectral components in the cho- sen spectral range. This number has simply been chosen as the minimum relevant number of waves necessary to obtain a quantitative agreement with NLSE simulation results. We numerically compute the spectral intensity of the wave using Eq. ( 18) and compare its evolution for spatial spectra corresponding to profiles 1 and 3. The input spectral phases are the same as in Fig. 5(d) and the relative amplitudes of the components are also described by the respective spatial spectra of the profiles.

Our results are shown in Fig. 6. Profile 1 (BG beam) shows the same parabolic structures as previously. In contrast, this feature is totally absent in the case of profile 3 where spectral broadening only occurs starting from a propagation distance of z ≈ 1700 µm. Particularly, the spectral expansion is in good qualitative agreement with NLSE simulation results of Fig. 5(c). We then conclude that the phase of the low intensity tails of the spatial spectrum also contribute to the initial spectral broadening stage.

We stress that these conclusions are theoretically only valid in case the phase distribution is preserved along propagation. According to our numerical simulation of the NLSE, the relative input phase remains approximately unchanged up to a propagation point of z = 1200 µm. As the spectral phase flattens because of linear propagation [START_REF] Jarutis | Focusing of Laguerre-Gaussian beams by axicon[END_REF], FWM interactions become more efficient and inevitably lead to significant growth of the axial wave and outer ring.

VII. CONCLUSIONS

In conclusion, we have developed a Four Wave Mixing model in order to interpret the different characteristics of the growth of Kerr-induced instabilities in Bessel beams. We have shown that nonlinear spectral distortions are established in two steps. The first step consists of spectral broadening and the generation of an axial wave seed. In a second step, this axial wave seed is amplified and an outer ring is generated via Four Wave Mixing. These new spatial frequency components interfere with the main input Bessel beam, yielding oscillations of the on-axis intensity.

We have then used a reduced model where only the dominant wave mixing terms were present so that we could expand the model to broad spatial spectra. This reduced nonlinear model allowed us to understand previous experimental and numerical results where input conditions of Bessel beams were yielding very different growth rates of nonlinear instabilities. We have demonstrated that these differences can be well explained from the weak differences in input spectral phases, even out-side the main peak of the spectral amplitude. We believe this approach will open new ways to control nonlinear instabilities, as well as to extend the applicability of Bessel beams for new applications. We also note that further work is needed to expand the model to ultrashort pulses with broadband temporal frequency spectrum.

FIG. 1 .

 1 FIG. 1. NLSE simulation results: (a) intensity distribution of a BesselGauss (BG) beam propagating in a pure nonlinear Kerr medium as a function of the radial and propagation distances r and z (Input peak intensity of the Gaussian beam 33.7 GW cm -2 ), (b) the corresponding on-axis intensity, (c) the spatial spectrum distribution | Ã(kr)| 2 /| Ã(kr0| 2 | (logarithmic scale dB) and (d) the spectral intensity of the axial wave (solid line) and outer ring in linear scale (dashed line) along the propagation distance. The intensity of both spectral components are normalized to the maximal intensity value of the central frequency.

FIG. 2 .

 2 FIG. 2. Evolution of frequency-resolved intensity of the (top)signal and (bottom) idler waves along propagation. The latter is computed for θ3 = 0.005θ0.

  FIG. 3. (a) On-axis intensity distribution of a BesselGauss beam propagating in a pure nonlinear Kerr medium for the soft (dashed line) and abrupt (dotted line) inputs; (b-c) the corresponding spatial spectra evolution along the propagation distance (logarithmic scale) for both cases respectively. Note that the intensity of newly-generated spectral components in the case of abrupt transition is two order of magnitude higher than the one reached in case of soft transition. Parameters are the presented in table I and Pin = 31.2 MW.

FIG. 4

 4 FIG. 4. a) Comparison between the input phase distributions at z = z -∆z = 0) where ∆z marks the position of the beam with sample input, where it is 0 for the soft case (red dashed line) and 2160 µm for the other case of abrupt input condition (red dotted line) b) Results using our FWM model: (spatial) frequency-resolved intensity (dB) of the signal wave resulting from the interaction of three pump waves composed each of five frequencies: kr/kr0 = [0.96 -1.04] in case their respective spectral phases correspond to the (top) soft and (bottom) abrupt input conditions respectively as described in Fig. 4(a).

1 FIG. 5 . 1 FIG. 6 .

 1516 FIG. 5. Simulation of the nonlinear propagation of three Bessel beams with different (a) target on-axis intensity profiles. Evolution of their respective (b) on-axis intensities, (c) spatial spectra (dB) along propagation. Spectral intensities are normalized to their respective maximal values. Numerical parameters are described in table I. (d) Input spectral phase profiles corresponding to the three Bessel beams.

TABLE I .

 I Numerical parameters used in simulations

	λ (µm)	0.8
	n	1.45
	n2 (m 2 /W)	2.48 10 -20
	θ ( • )	4
	w0 (µm)	300

The research leading to these results has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 682032-PULSAR). This work has been supported by the EIPHI Graduate School (contract ANR-17-EURE-0002).