Nonlinear Bessel vortex beams for applications

Figure 3 shows a comparison of propagation features for Bessel beams and Bessel-vortices in air, for the same cone angle θ = 2.2 deg. Very regular beam profiles were observed along the focal region, both for the Bessel beam and the Bessel vortex. From the spacing between the rings ($2{{r}_{0}}\sim 15.7\;\mu$m), the Bessel cone angle can be calculated: ${\rm sin} \theta =2.405/({{k}_{0}}{{r}_{0}})$, leading to $\theta \sim 2.24\;{\rm deg}$, Focusing a Gaussian beam with an axicon is well known to result in a Bessel–Gauss beam that remains quasi propagation invariant along a focal region, the Bessel zone, of length ${{L}_{B}}={{w}_{0}}/{\rm tan} \theta$. Adding helicity (here m = 1) to the Gaussian beam results in a higher-order Bessel–Gauss beam that remains quasi-propagation invariant over the same focal region. The highest intensity is reached on the surface of a hollow cylinder corresponding to the maximum of the ${{J}_{m}}({{k}_{0}}{\rm sin} \theta r)$ function. More precisely, in the linear propagation regime, the intensity profile is determined by the relation [16]

$$\begin{eqnarray}&&I(r,z)=\frac{4P{{k}_{0}}{\rm sin} \theta }{{{w}_{0}}}\frac{z}{{{z}_{B}}}J_{m}^{2}({{k}_{0}}{\rm sin} \theta r){\rm exp} (-2\frac{{{z}^{2}}}{L_{B}^{2}}),\end{eqnarray} \tag{ 1 }$$

where P denotes the input power of the Gaussian beam and z denotes the propagation distance from the tip of the axicon. With w₀ = 4.4 mm, we obtain a Bessel zone L_B = 11.3 cm. The highest intensity is reached at $\tilde{z}={{L}_{B}}/2$ and equals ${{I}_{{\rm max} }}=2(P/{{w}_{0}}){{k}_{0}}{\rm sin} \theta {\rm exp} (-1/2){{M}_{m}}$, where M_m = (1, 0.34, 0.24, 0.19, ...) denotes the maximum of the squared Bessel function of order m, for $m=(0,1,2,3,\ldots )$. We estimate that the beam enters a nonlinear propagation regime when the maximum intensity is sufficient to ionize the medium. With a threshold of 10¹² W cm⁻² in glass and 10¹³ W cm⁻² in air, we need a minimum beam power of ∼240 MW in fused silica (corresponding to an energy of 13 μJ for a 45 fs pulse) and 3 GW in air (0.2 mJ). With pulse energies smaller than 100 μJ, we could explore the nonlinear propagation regimes in glass.

Zoom In Zoom Out Reset image size

Propagation in the nonlinear regime was investigated by placing samples of fused silica with thicknesses 10, 20, and 40 mm close to the exit surface of the liquid immersion axicon. Figure 4 shows the fluence patterns at the exit of the 20 mm long fused silica sample when the input energy of the Gaussian beam is increased, for the Bessel beam (m = 0) and the Bessel vortex (m = 1). Again in both cases, clean Bessel-like beams are obtained, which look like the ideal Bessel beams $J_{m}^{2}({{k}_{0}}{\rm sin} \theta r)$ (with m = 0 or m = 1). We checked that these fluence patterns remain nearly identical for all sample lengths, i.e., nonlinear propagation is also featured by propagation invariance (data not shown). However, it is known that propagation invariant nonlinear Bessel beams exhibit ring compression (decrease of the ring spacing close to the main lobe with respect to the ring spacing of Bessel beams) and a loss of ring contrast, as the optical Kerr effect and multiphoton absorption (MPA) increase, respectively [17, 18]. Moreover, nonlinear propagation of Bessel beams with peak intensity above a certain threshold no longer leads to propagation invariant beams [19, 20]. As can be seen from the profiles in figure 4, ring compression and ring contrast loss is observed for the Bessel vortices as well, though higher pulse energy is required to observe the same ring spacing or contrast decrease as for Bessel beams with the same cone angle. Experimentally, white-light generation, another indicator for nonlinear propagation, was observed at 22.2 μJ for the Bessel beam and 36 μJ for the Bessel vortex, respectively.

Zoom In Zoom Out Reset image size

These observations are reminiscent of similar features previously reported for nonlinear Bessel beams: Porras and coworkers showed in 2004 the existence of Bessel-like propagation-invariant solutions to the nonlinear Schrödinger equation in the presence of Kerr effect and MPA [19]:

$$\begin{eqnarray}&&\frac{\partial E}{\partial z}=\frac{{\rm i}}{2{{k}_{0}}}\nabla _{\bot }^{2}E+{\rm i}{{k}_{0}}\frac{{{n}_{2}}}{{{n}_{0}}}|E{{|}^{2}}E-\frac{{{\beta }_{K}}}{2}|E{{|}^{2K-2}}E,\end{eqnarray} \tag{ 2 }$$

where in polar coordinates $\nabla _{\bot }^{2}E=\frac{1}{r}\frac{\partial E}{\partial r}+\frac{{{\partial }^{2}}E}{\partial {{r}^{2}}}+\frac{1}{{{r}^{2}}}\frac{{{\partial }^{2}}E}{\partial {{\vartheta }^{2}}}$. The parameters K₀, n₀, n₂ and β_K are provided in table 1. Recently, weakly localized propagation invariant solutions to equation (2) were found in the form of nonlinear high-order Bessel beams [14, 21]. These solutions have a stationary (z-independent) intensity profile and a linear dependence of the phase profile upon z:

$$\begin{eqnarray}&&E(r,\vartheta ,z)=a(r){{{\rm e}}^{{\rm i}(-\delta z+\int q(r){\rm d}r+m\vartheta )}},\end{eqnarray} \tag{ 3 }$$

where the wave vector shift is positive δ > 0, and the amplitude a(r) and phase gradient q(r) only depend on the transverse variable r. By introducing this amplitude and phase decomposition into equation (2), we find that a(r) and q(r) are solutions to Newton-like equations for a pseudo-particle with radial and angular positions a, ϕ, (with $\phi (r)=\dot{q}$), angular momentum J_r ≡ k₀ ⁻¹q(r)a²(r) and pseudo-time r:

$$\begin{eqnarray}&&\ddot{a}-{{q}^{2}}a=-\frac{{\dot{a}}}{r}+\frac{{{m}^{2}}}{{{r}^{2}}}a-\frac{\partial V}{\partial a},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\dot{{{J}_{r}}}=-\frac{{{J}_{r}}}{r}-{{\beta }_{K}}{{a}^{2K}},\end{eqnarray} \tag{ 5 }$$

where dots indicate differentiation with respect to r and the pseudo-potential is defined as $V(a)\equiv {{k}_{0}}\delta {{a}^{2}}+k_{0}^{2}{{n}_{2}}{{a}^{4}}/(2{{n}_{0}}).$ Vortex solutions must satisfy the boundary conditions a(0) = 0, q(0) = 0. The last boundary condition is governed by the decay rate as $r\to \infty$ which corresponds to that of a Bessel beam carrying helicity m: as $a(r)\to 0$, nonlinearity is negligible and V(a) ∼ k₀δa², therefore, the solution tail can be viewed as a superposition of the two linearly independent Hankel functions that are solutions to the linearized equation (4):

$$\begin{eqnarray}&&a(r)=\frac{{{a}_{0}}}{2}[{{\alpha }_{{\rm out}}}H_{m}^{(1)}(\sqrt{2{{k}_{0}}\delta }r)+{{\alpha }_{{\rm in}}}H_{m}^{(2)}(\sqrt{2{{k}_{0}}\delta }r)].\end{eqnarray} \tag{ 6 }$$

$H_{m}^{(1,2)}(\sqrt{2{{k}_{0}}\delta }r)$ represent solutions carrying an outward or inward power flux, respectively. Their respective amplitudes α_out,in differ, $|{{\alpha }_{{\rm in}}}|\gt |{{\alpha }_{{\rm out}}}|$, leading to a fundamental feature of nonlinear Bessel vortices: a helical and inward power flux ${\bf J}={{J}_{r}}{{{\bf u}}_{r}}+{{J}_{\vartheta }}{{{\bf u}}_{\vartheta }}=k_{0}^{-1}$ $[{{a}^{2}}q{{{\bf u}}_{r}}+(m{{a}^{2}}/r){{{\bf u}}_{\vartheta }}]$ transports energy from the tail of the beam to the most intense lobe where it is dissipated by MPA. In this respect, integration of equation (5) along the radial direction is equivalent to the energy conservation equation

$$\begin{eqnarray}&&2\pi r{{J}_{r}}(r)=-{{\beta }_{K}}\int _{0}^{r}{{a}^{2K}}2\pi r{\rm d}r,\end{eqnarray} \tag{ 7 }$$

which states that power losses due to MPA within a disk of radius r are compensated by the power flux flowing through its boundary (circle of radius r). This ensures propagation invariance of the nonlinear Bessel vortex [14, 21].

Examples of nonlinear Bessel beams m = 0 and vortices m = 1 integrated numerically are shown in figure 5 for a given cone angle and given arbitrary nonlinear parameter chosen to illustrate the main features of nonlinear Bessel vortices. For a given medium (n₀, n₂, β_K), cone angle (δ or θ) and helicity m, a continuum of solutions of this type exists for a range of intensities $[0,{{I}_{{\rm max} }}]$, where ${{I}_{{\rm max} }}$ can be found only numerically and is usually two orders of magnitude larger than the breakdown threshold of the material. For example for Corning0211 glass, ${{I}_{{\rm max} }}\gt 1$ PW cm⁻² [14]. As the helicity is increasing, the position of the most intense lobe is found at a larger radius. This radius is however slightly smaller than the radius of the linear solution with the same cone angle, i.e., the Kerr effect manifests itself in a compression of the most intense rings and a decrease of the ring spacing close to the center. Another manifestation of nonlinearity is the ring contrast which is decreasing as MPA increases. Asymptotically, this contrast can be readily obtained from an expansion of equation (6) as $r\to \infty$ and reads: $C=2|{{\alpha }_{{\rm in}}}{{\alpha }_{{\rm out}}}|/(|{{\alpha }_{{\rm in}}}{{|}^{2}}+|{{\alpha }_{{\rm out}}}{{|}^{2}})$. Therefore the highest contrast C = 1 is reached for linear propagation, i.e., when $|{{\alpha }_{{\rm in}}}|=|{{\alpha }_{{\rm out}}}|$, the inward and outward flux are balanced and there is no absorption (${{\beta }_{K}}=0$). In contrast, the ring disappear, leading to zero contrast C = 0 if the power flux is directed only inward and the outward component is zero $|{{\alpha }_{{\rm out}}}|=0$, corresponding to the highest possible MPA in a propagation invariant regime. The link between losses and the contrast $C\equiv 2R/(1+{{R}^{2}})$, where $R\equiv |{{\alpha }_{{\rm out}}}/{{\alpha }_{{\rm in}}}|$, can be inferred from an introduction of the asymptotic expansion of equation (6) at large r into the energy conservation equation (7): ${{a}_{0}}|{{\alpha }_{{\rm in}}}{{|}^{2}}(1-{{R}^{2}})={{k}_{0}}{{\beta }_{K}}\int _{0}^{\infty }{{a}^{2K}}r{\rm d}r$. Since Bessel vortex beams are computed numerically, the link between losses and contrast can only be entirely specified numerically. The reader is referred to [21] for further details on this point.

Zoom In Zoom Out Reset image size

We interpret our observations as a manifestation of the reshaping of the beam into quasi-stationary nonlinear Bessel vortices. Propagation invariant nonlinear Bessel vortices are ideal beams carrying infinite power as shown by the $1/\sqrt{r}$ decay rate of their tail. Without this feature, MPA would unavoidably prevent stationarity. Rather than exact stationarity, quasi-stationarity arise over a finite length, the Bessel zone, because our Bessel vortices are shaped as the ideal propagation invariant solutions and apodized by a Gaussian beam, of width much larger than the most intense lobe of the Bessel-like solutions where MPA occurs. All features of nonlinear Bessel vortices are therefore observed within the Bessel zone. To complete characterization of the beam reshaping process during propagation, we investigated the decay of the contrast of the rings of the Bessel vortex beam during nonlinear propagation in fused silica samples of different thickness (10, 20, and 40 mm). We take line-cuts (in x-direction) through the beam center, being imaged at the exit of the samples. The line-cuts are normalized and Fourier transformed. We registered the relative amplitude of the Fourier peak associated with the transverse wave number of the Bessel vortex. Results are shown in figure 6 for a range of pulse energies up to 90 μJ. The relative amplitude of the Fourier peak decreases as nonlinear propagation becomes more prominent. In each case, we reported the onset of supercontinuum generation, which is a standard signature of nonlinear propagation. Supercontinuum generation was experimentally observed, when the Fourier contrast had decreased to about half the initial value. We mark the energies, when experimentally supercontinuum generation was observed by dashed lines in figure 6 for the different sample lengths.

Zoom In Zoom Out Reset image size

By using equation (1), we can give estimations of intensity levels at the exit face of the samples. For the three energies marked in figure 6, we obtain $I(z=10\;{\rm mm})=1.2\times {{10}^{12}}$ W cm⁻² for E_in = 66 μJ, $I(z=20\;{\rm mm})=1.4\times {{10}^{12}}$ W cm⁻² for E_in = 50 μJ, and $I(z=40\;{\rm mm})=1.3\times {{10}^{12}}$ W cm⁻² for E_in = 32 μJ. In other words, the point when the contrast is lost due to nonlinear propagation can be well identified with the same intensity levels on the back surfaces for the three different sample lengths. We can therefore conclude that the nonlinear propagation we observe, when the contrast is lost, happens at about the same intensity and in the vicinity of the respective back surface and not within the sample. As shown in the next section, this would be possibly different for longer samples which would include the peak of the pulse. The onset of nonlinear propagation arises for quite low intensity levels of 10¹² W cm⁻². At this intensity, no significant plasma density can be formed for a 45 fs pulse. For this cone angle (θ = 1.5 deg in fused silica), the Bessel vortex beam does not allow for the generation of hollow plasma cylinders.

The maximum plasma density in the tubular plasma channel can be evaluated from the peak intensity, which in turn can be evaluated by energy conservation arguments similar to those applied for a Bessel beam [22]. The nonlinear Bessel vortex is assumed to have an invariant peak intensity I_p and shape over the Bessel zone. A fraction x of the input beam energy E_in is lost nearly exclusively in the main lobe to generate the plasma channel, by multiphoton ionization. This energy balance reads

$$\begin{eqnarray}&&\frac{x{{E}_{{\rm in}}}}{{{L}_{B}}}={{\beta }_{K}}\int \int {{I}^{K}}(r,t)2\pi r{\rm d}r{\rm d}t,\end{eqnarray} \tag{ 8 }$$

where the intensity distribution is approximated by the product of a high-order Bessel beam profile by a Gaussian pulse:

$$\begin{eqnarray}&&I(r,t)\sim {{I}_{p}}{\rm exp} (-4{\rm ln} 2{{t}^{2}}/{{T}^{2}})J_{m}^{2}({{k}_{0}}{\rm sin} \theta r)/{{M}_{m}}.\end{eqnarray} \tag{ 9 }$$

By introducing this profile in equation (8), the peak intensity is explicitly determined by

$$\begin{eqnarray}&&{\rm tan} \theta {{{\rm sin} }^{2}}\theta =\frac{{{\beta }_{K}}g(K,m)}{n_{0}^{2}}\frac{T{{w}_{0}}\lambda _{0}^{2}I_{p}^{K}}{x{{E}_{{\rm in}}}}\end{eqnarray} \tag{ 10 }$$

where $g(K,m)\equiv M_{m}^{-K}{{(\pi K{\rm ln} 2)}^{-1/2}}\int _{0}^{u_{0}^{(m)}}J_{m}^{2K}(u)u{\rm d}u$.

Figure 7 shows the peak intensity and electron density evaluated as functions of the cone angle from equation (10) for different helicities of the Bessel vortex. Note that electron density is evaluated by considering only multiphoton ionization processes. When it approaches the density of neutral atoms, depletion of the valence band should lead to saturation. Depletion remains however negligible for electron densities smaller than 10²¹ cm⁻³. This simple scaling law predicts that a cone angle of 1.5 deg would lead to a peak intensity of 4 × 10¹² W cm⁻² and an electron density of 1.3 × 10¹⁸ cm⁻³ for the Bessel vortex of lowest helicity m = 1. Peak intensity and electron density decrease with increasing helicity. An increase of the cone angle in the range 10–15 deg, is sufficient to reach peak electron densities in the range 4 × 10²⁰–10²¹ cm⁻³, for which energy transfer to the material can be deemed sufficiently efficient for inducing refractive index change or damage [23, 24].

Zoom In Zoom Out Reset image size

The model used for numerical simulations consists in a (2+1)-dimensional unidirectional propagation equation for the nonlinear envelope E(x, y, z), which takes a canonical form [25]. As a nonparaxial equation, it is more easily written for the Fourier components of the electric field envelope $\mathcal{E}({{k}_{x}},{{k}_{y}},z)\equiv \mathcal{F}[E(x,y,z)]$, where $\mathcal{F}$ denotes Fourier transform in the transverse plane:

$$\begin{eqnarray}&&\frac{\partial \mathcal{E}}{\partial z}={\rm i}\mathcal{K}({{k}_{x}},{{k}_{y}})\mathcal{E}+{\rm i}\frac{{{\omega }_{0}}}{c}\frac{\mathcal{P}}{2{{n}_{0}}{{\epsilon }_{0}}}-\frac{\mathcal{J}}{{{n}_{0}}2{{\epsilon }_{0}}c}.\end{eqnarray} \tag{ 11 }$$

Nonparaxial diffraction is described in the propagation constant $\mathcal{K}({{k}_{x}},{{k}_{y}})\equiv \sqrt{k_{0}^{2}-k_{x}^{2}-k_{y}^{2}}-{{k}_{0}}$. The Kerr effect is accounted for in the nonlinear polarization $\mathcal{P}({{k}_{x}},{{k}_{y}},z)\equiv \mathcal{F}[P(x,y,z)]$, with $P\equiv 2{{\epsilon }_{0}}{{n}_{0}}{{n}_{2}}IE$, where $I\equiv (1/2)$ ${{\epsilon }_{0}}c{{n}_{0}}|E{{|}^{2}}$. The current $\mathcal{J}({{k}_{x}},{{k}_{y}},z)\equiv$ $\mathcal{F}[J(x,y,z)],$ where $J\equiv {{J}_{{\rm MPA}}}+{{J}_{{\rm PL}}}$ comprises two contributions due to MPA and plasma-induced effects (PL). The contribution of MPA reads ${{J}_{{\rm MPA}}}={{\epsilon }_{0}}c{{n}_{0}}{{\beta }_{K}}{{I}^{K-1}}E$, with coefficient β_K and number of photons K.

Plasma-induced effects, i.e., plasma absorption and plasma defocusing are included in a single expression for the current, derived from the Drude model:

$$\begin{eqnarray}&&{{J}_{{\rm PL}}}={{\epsilon }_{0}}\frac{\omega _{0}^{2}{{\tau }_{c}}}{1-{\rm i}\omega {{\tau }_{c}}}\frac{\rho }{{{\rho }_{c}}}E,\end{eqnarray} \tag{ 12 }$$

where τ_c denotes a phenomenological collision time in the medium and ρ_c, the critical plasma density beyond which the medium is no longer transparent.

The evolution of the electron plasma density is modeled by a rate equation modeling multiphoton and avalanche ionization with cross section σ, as well as electron plasma relaxation with typical time τ_r:

$$\begin{eqnarray}&&\frac{\partial \rho }{\partial t}=\left( \frac{{{\beta }_{K}}}{K\hbar {{\omega }_{0}}}{{I}^{K}}+\frac{\sigma \rho I}{{{U}_{i}}} \right)\left( 1-\frac{\rho }{{{\rho }_{{\rm nt}}}} \right)-\frac{\rho }{{{\tau }_{r}}}.\end{eqnarray} \tag{ 13 }$$

We consider large cone angles, giving rise to a short Bessel zone over which dispersion is negligible. The plasma density is therefore evaluated in the frozen time approximation that consists in using a Gaussian pulse with the nominal duration of the laser.

In table 1, we provide the model parameters used for fused silica and BK7 at the pulse central wavelength of 800 nm, together with references from which the values were obtained.

Table 1.  Parameters uses for modeling the response of fused silica and BK7. Refraction index: n₀. Kerr index coefficient: n₂ [23, 26]. Number of photons: $K\equiv \langle \frac{{{U}_{i}}}{\hbar {{\omega }_{0}}}+1\rangle$. Coefficient for multiphoton absorption: ${{\beta }_{K}}$ [27]. Bandgap: U_i collision time: ${{\tau }_{c}}$. Plasma relaxation time: ${{\tau }_{r}}$. Neutral density: ${{\rho }_{{\rm nt}}}$. Critical plasma density: ${{\rho }_{c}}\equiv \omega _{0}^{2}{{\epsilon }_{0}}{{m}_{e}}/q_{e}^{2}$.

	Fused silica	BK7
n0	1.45	1.51
n2 (cm2 W−1)	3.54 × 10−16	3.45 × 10−16
K	5	3
${{\beta }_{K}}$ (cm2K−3 W−(K − 1))	6.8 × 10−50	3.9 × 10−28
Ui (eV)	7.1	4.2
${{\tau }_{c}}$ (fs)	1.3	1
${{\tau }_{r}}$ (fs)	150	2000
${{\rho }_{{\rm nt}}}$ (cm−3)	2.1 × 1022	2.1 × 1022
${{\rho }_{c}}$ (cm−3)	1.8 × 1021	1.8 × 1021

4.2.1. Nonlinear Bessel vortex in fused silica at small cone angles

Figure 8 shows simulation results for the experimental scheme A with a cone angle of θ = 1.5 deg in fused silica. First, we note that the sample thickness was approximately half the length of the Bessel zone. For a pulse energy of 7 μJ, the beam forms a Bessel vortex with its typical main ring surrounded by lower intensity rings. The highest intensity reaches 400 GW cm⁻² and does not lead to a significant plasma density. The highest calculated plasma density is about 10¹⁰ cm⁻³ (not shown). This regime therefore shows the formation of a quasi linear Bessel vortex. If the input energy is increased, there is however no abrupt threshold but a continuous transition between this regime and the nonlinear regime featured by the formation of propagation invariant nonlinear Bessel vortices, as found semi-analytically in section 3.3. By further increasing the input energy above a certain threshold, the propagation invariant regime is no longer found. Figure 8 shows that at 66 μJ, the uniform ring formed in the beginning of the Bessel zone undergoes instability, forming a high intensity peak that rotates around the propagation axis at a fairly regular rate, reflecting the angular momentum carried by the beam. The highest intensity of 5 TW cm⁻² is sufficient to generate a plasma of density of 10¹⁹ cm⁻³. However, this density might be still too low to induce damage in the material. These experimental conditions might be just below threshold for refractive index change in fused silica. As a final comparison with the experimental patterns, we note that the rotating peak appears above 32 μJ, similarly to the results in [14]. At higher energy the contrast gets smoothed out, and the rotating peak is not so pronounced as in the simulation. Dispersion of the 45 fs pulse was neglected in the simulation but could stretch the pulse by a factor of ∼2.3 over 6 cm in the experiment and therefore decrease the peak intensity and shift the instability threshold towards larger pulse energies.

Zoom In Zoom Out Reset image size

4.2.2. Nonlinear Bessel vortex in BK7 at large cone angles

Figure 9 shows simulation results for experimental parameters in scheme B, performed at higher cone angle of 11 deg in BK7 (16.95 deg in air). The experimental scheme shown in figure 2 is described in more detail in section 5. The simulations allow us to anticipate a high electron density (10²¹ cm⁻³) distributed uniformly on a hollow cylinder extending over the entire thickness of the BK7 sample (100 μm). Uniformity appears after a regularization stage in the first microns after the entrance face of the sample where the intensity is sufficiently high for secondary lobes of the Bessel vortex to also generate plasma. However, it is clear that a propagation invariant regime takes place. The peak intensity is a factor of 2 smaller than in the simulation for the small-cone-angle, 45 fs pulse. Among the reasons leading to a higher electron density, not only has BK7 a slightly smaller bandgap than fused silica but also the picosecond pulse used for BK7 facilitates the generation of electrons by avalanche while requiring a lower intensity to reach the same level as that obtained with a fs pulse. Finally, propagation invariant nonlinear Bessel vortices at a given peak intensity are more stable for large cone angles and can therefore attract the beam dynamics [14].

Zoom In Zoom Out Reset image size