Stationary and stable light beam propagation in Kerr media with nonlinear absorption with controllable dissipation patterns

TL;DR

This paper demonstrates the stable propagation of light beams with customizable intensity and dissipation patterns in Kerr media with nonlinear absorption, enabling the creation of tailored plasma channels.

Contribution

It introduces a method for controlling light beam dissipation patterns in Kerr media with nonlinear absorption, including the formation of plasma channels with specific geometries.

Findings

Stable propagation of arbitrary intensity patterns

Generation of tailored plasma channels

Analysis of formation laws from Bessel beam superpositions

Abstract

We report on the stationary and robust propagation of light beams with rather arbitrary and controllable intensity and dissipation transverse patterns in self-focusing Kerr media with nonlinear absorption. When nonlinear absorption is due to multi-photon ionization at high beam powers in transparent media such as glasses or air, these beams can generate multiple plasma channels with tailored geometries. Their nature and spatial characteristics are discussed in detail, as well as the laws determining their spontaneous formation from coherent superpositions of Bessel beams of different amplitudes and topological charges.

Figures

Equations

$${\partial }_{z}A=\frac{i}{2k}{\Delta }_{\perp }A+i\frac{k{n}_2}{n}|A{|}^{2}A-\frac{{\beta }^{(M)}}{2}|A{|}^{2M-2}A,$$

$$\rho \equiv k\theta r=\sqrt{2k|\delta |}r,\zeta \equiv |\delta |z,\tilde{A}\equiv {\left(\frac{{\beta }^{(M)}}{2|\delta |}\right)}^{\frac{1}{2M-2}}A$$

$${\partial }_{\zeta }\tilde{A}=i{\Delta }_{\perp }\tilde{A}+i\alpha |\tilde{A}{|}^2\tilde{A}-|\tilde{A}{|}^{2M-2}\tilde{A},$$

$$\alpha \equiv {\left(\frac{2|\delta |}{{\beta }^{(M)}}\right)}^{1/(M-1)}\frac{k{n}_2}{n|\delta |}.$$

$$\frac{1}{2}{\partial }_{\zeta }{\tilde{a}}^2+{\nabla }_{\perp }\cdot \mathbf{j}=-{\tilde{a}}^{2M},$$

$${\partial }_{\zeta }\Phi =\frac{{\Delta }_{\perp }\tilde{a}}{\tilde{a}}-{\left|{\nabla }_{\perp }\Phi \right|}^2+\alpha {\tilde{a}}^2,$$

$${\nabla }_{\perp }\cdot \mathbf{j}=-{\tilde{a}}^{2M},$$

$${\Delta }_{\perp }\tilde{a}-{\left|{\nabla }_{\perp }\varphi \right|}^2\tilde{a}+\tilde{a}+\alpha {\tilde{a}}^3=0,$$

$$-{\oint }_C\mathbf{j}\cdot \mathbf{n}dl={\int }_S{\tilde{a}}^{2M}ds,$$

$$F_{\infty }\equiv -{\int }_0^{2\pi }{j}_{\rho }(\rho ,\phi )\rho d\phi =N_{\infty },$$

$$N_{\infty }\equiv {\int }_{{\mathbb{R}}^2}{\tilde{a}}^{2M}ds,$$

$$H_s^{(1,2)}(\rho )\simeq \sqrt{\frac{2}{\pi \rho }}{e}^{\pm i(\rho -\pi /4-\pi s/2)},$$

$$\tilde{A}\simeq \sum _j\frac{1}{2}\left[b_{\mathrm{out},j}{H}_{s_j}^{(1)}(\rho )+b_{\mathrm{in},j}{H}_{s_j}^{(2)}(\rho )\right]e^{i{s}_j\phi }{e}^{-i\zeta },$$

$$\tilde{A}\simeq \frac{1}{\sqrt{2\pi \rho }}\left[b_{\mathrm{out}}(\phi )e^{i(\rho -\pi /4)}+b_{\mathrm{in}}(\phi )e^{-i(\rho -\pi /4)}\right]e^{-i\zeta },$$

$$b_{\mathrm{out}}(\phi )$$

$$b_{\mathrm{in}}(\phi )$$

$$2\pi \rho |\tilde{A}{|}^2$$

$$C(\phi )=\frac{2|b_{\mathrm{out}}(\phi )||b_{\mathrm{in}}(\phi )|}{|b_{\mathrm{out}}(\phi ){|}^2+|b_{\mathrm{in}}(\phi ){|}^2}$$

$$R(\phi )=|b_{\mathrm{out}}(\phi ){|}^2+|b_{\mathrm{in}}(\phi ){|}^2.$$

$$|b_{{\rm in},j}|=|b_{L,j}|\quad\mbox{for all $j$.}$$

$$-F_{\infty }=-{\int }_0^{2\pi }\frac{|b_{\mathrm{out}}(\phi ){|}^2-|b_{\mathrm{in}}(\phi ){|}^2}{2\pi }d\phi =N_{\infty },$$

$$\sum _j|b_{\mathrm{in},j}{|}^2-\sum _j|b_{\mathrm{out},j}{|}^2=N_{\infty }$$

$$\tilde{A}(\rho ,\phi ,0)=b_{L,s}{J}_s(\rho )e^{is\phi }+b_{L,-s}{J}_{-s}(\rho )e^{-is\phi }$$

$$\tilde{A}$$

$$b_{\mathrm{out}}(\phi )$$

$$b_{\mathrm{in}}(\phi )$$

$$\tilde{A}=\frac{1}{\sqrt{2\pi \rho }}\left[b_s(\rho )e^{is\phi }+b_{-s}(\rho )e^{-is\phi }\right]e^{-i\zeta },$$

$$b_s(\rho )$$

$$b_{-s}(\rho )$$

$$2\pi \rho |\tilde{A}{|}^2$$

$$C(\rho )=\frac{2|b_s(\rho )||b_{-s}(\rho )|}{|b_s(\rho ){|}^2+|b_{-s}(\rho ){|}^2}$$

$$\mathbf{j}\simeq \frac{[|b_{\mathrm{out}}(\phi ){|}^2-|b_{\mathrm{in}}(\phi ){|}^2]}{2\pi \rho }{\mathbf{u}}_{\rho }+\frac{s}{\rho }\frac{[|b_s(\rho ){|}^2-|b_{-s}(\rho ){|}^2]}{2\pi \rho }{\mathbf{u}}_{\phi }.$$

$$L=\text{mbox}Im\left\{{\tilde{A}}^{\star }\frac{\partial \tilde{A}}{\partial \phi }\right\}\simeq \frac{s}{2\pi \rho }[|b_s(\rho ){|}^2-|b_{-s}(\rho ){|}^2],$$

Full text

I Introduction

The achievement of localized, stationary and stable propagation of multidimensional waves in homogeneous, isotropic nonlinear media is one of the main problems in nonlinear optics, and constitutes today a broad area of research MALOMED1 ; KIVSHAR ; CHEN . Many of the elaborated solutions, some of them realized experimentally, rely on two broad models. The first one is the conservative model based on nonlinear Schrödinger equations with dispersive, Kerr-type nonlinearities and others, describing passive propagation of intense waves in unbounded transparent media. The second model is the complex Ginzburg-Landau equation in systems with gain and loss, as laser cavities AKHMEDIEV ; ACKEMANN ; MALOMED2 . More or less broad families of solitons and dissipative solitons (also called auto-solitons or cavity-solitons) with a variety of geometries have been reported to be supported by these models. Given their relation with the present article, we mention the families of vortex-solitons —solitons with a nested vortex— KRUGLOV , their azimuthally modulated versions called azimuthons DESYATNIKOVPrl1 ; MINOVICHOe and solitons clusters DESYATNIKOVPrl2 , see DESYATNIKOVReview for a review. In the conservative models, a stable balance between diffraction and the dispersive nonlinearities supports the robust stationary propagation. In the dissipative model, an additional balance is needed between gain and losses PAULAU ; CRASOVAN ; ZHANGPra .

Non-conservative models involving only losses are usually left aside because they do not support solitons, in the sense of self-trapped structures. However the nonlinear Schrödinger equation with purely loss terms and no gain terms arises naturally in situations as common as light beam propagation in transparent dielectrics media at intensities such that multi-photon absorption processes are relevant. The existence of light beams, called nonlinear unbalanced Bessel beams, with the ability to propagate without any change, including no attenuation, in media with Kerr-type nonlinearities and nonlinear absorption, was described in 2004 PORRASPrl . They provided a theoretical basis for interpreting some features of light filaments generated from intense, ultrashort pulses COUAIRONPhysRep ; BERGERep ; BERGEDis , where nonlinear absorption arises from multi-photon ionization, particularly the replenishment and self-healing mechanisms involved in their robust and quasi-stationary propagation DUBIETIS ; GAIZAUSKAS ; FACCIOOptComm ; ROSOAppPhysA ; ZEMLYANOV ; DUBIETISLithuanian ; FACCIOPra . Stationarity with nonlinear absorption is indeed possible for less strongly localized waves than solitons, as nonlinear versions of Bessel-like and Airy-like beams DURNINPrl ; SIVILOGLOUOptLett carrying, ideally, an infinite amount of power, or a power reservoir. In these beams PORRASPrl ; PORRASJosaB ; JUKNAOe ; XIESci ; LOTTIPra the power in the reservoir is permanently flowing towards the central region of high intensity, where most of nonlinear power losses take place, replenishing it. After the fundamental nonlinear unbalanced Bessel beam PORRASPrl , the same mechanism has been more recently seen to work with high-order (vortex-carrying) nonlinear Bessel beams PORRASJosaB ; JUKNAOe ; XIESci . Here we will refer to all them as nonlinear Bessel vortex beams (nonlinear BVBs), the fundamental beam being the particular vortex-less situation.

Nonlinear BVBs have proven to be particularly fruitful to explain the filamentation dynamics seeded by ultrashort (linear) Bessel beams POLESANAPra ; AKTURKOptComm ; POLESANAPre ; POLESANAPrl ; POLESANAOe ; MAJUSJosaB ; AKTURKSpie ; KAYAAip ; JUKNAOe ; XIESci . Even if the physics of filamentation involves complex light-matter interactions and rich spatiotemporal dynamics COUAIRONPhysRep ; BERGERep ; BERGEDis , the three key phenomena determining substantially the filamentation dynamics with conical beams are just those supporting the stationarity of nonlinear BVBs, namely, diffraction, Kerr self-focusing and the nonlinear absorption POLESANAPra ; JUKNAOe ; XIESci . The steady and unsteady regimes of filamentation observed in the filamentation with the lowest-order Bessel beam POLESANAPra have been successfully explained in terms of nonlinear BVB attractors and their stability properties PORRASPra1 . More recently, tubular filamentation from vortex-carrying Bessel beams has been created JUKNAOe ; XIESci , and has also been interpreted as the formation of a vortex-carrying nonlinear BVB that is stable against radial and azimuthal perturbations PORRASPra2 . The disintegration of the tubular regime into rotatory or random filaments observed in related experiments is also explained in terms of the development of the azimuthal instability of unstable nonlinear BVBs PORRASPra2 .

The intensity profiles of nonlinear BVBs are all circularly symmetric, and so they are the dissipation channels that they can generate. The tubular regime is a sophisticated example of a “tailored” or controlled transversal pattern, which is expected to open new perspectives in laser-powered material processing and other applications XIESci ; COURVOISIERApplPhys ; BHUYAN ; KRISNALaserPhysics ; ARNOLDJPhysB ; COURVOISIER2 . Complex filamentation patterns other than circular arise naturally in multiple filamentation from small-scale self-focusing of laser beams carrying many critical powers COUAIRONPhysRep ; BERGERep ; BERGEDis , or from the azimuthal instability of tubular filaments JUKNAOe ; XIESci ; PORRASPra2 , but originating from noisy perturbations, these patterns are largely uncontrollable and far from being propagation-invariant.

In this paper we investigate if propagation-invariant light beams with other geometries, preferably with controllable intensity and dissipation patterns, are supported by transparent dielectrics with self-focusing nonlinearities and nonlinear absorption. We find the affirmative answer when studying by means of the numerical simulations the nonlinear propagation dynamics seeded by illuminating the medium with coherent superpositions of Bessel beams of different topological charges and amplitudes. These (linear) Bessel beam superpositions are easily generated using spatial light modulators FROEHLYJosaA ; CLARKOe ; VALILYEUOe , illuminating an axicon with several coherent Gaussian beams with embedded vorticities, or with superpositions of Laguerre-Gauss beams ARLTOptCommun . For a range of the parameters that define the input Bessel beam superposition and the material medium, stationary and stable propagation states dissipating constantly their power along organized and narrow channels are seen to emerge spontaneously. We confine ourselves to the superposition of two Bessel beams, where the channels geometry is primarily controlled by their amplitudes and topological charges, and the variety of intensity and dissipation patterns is already very vast.

To clearly distinguish these states from other known localized structures, we first present the more general states without any particular symmetry, and then we focus on those with $n$-fold rotational symmetry formed from superpositions of two Bessel beams of opposite topological charges, resembling azimuthons DESYATNIKOVPrl1 ; MINOVICHOe ; PAULAU ; CRASOVAN ; ZHANGPra . These structures are nevertheless neither conservative nor supported by a gain-loss balance, and do not rotate either. Given the distinctive properties of having arbitrary shapes, and of not being supported by gain, i. e., of being strictly dissipative, we have called them “dissipatons”.

Outside their central region, dissipatons are always seen to have linear conical tails, and therefore a power reservoir, shaped like a superposition of unbalanced Bessel beams (i. e., with unequal inward and outward Hänkel beam amplitudes), all them with the same cone angle but different topological charges. The excess of the inward Hänkel beam amplitudes creates an inward radial current, as that in nonlinear BVBs but with an azimuthal dependence, and such that net power flux coming from the reservoir equals the total power loss rate in the dissipative core. We find simple rules that describe the dissipaton formation. Given an input superposition of Bessel beams, the dissipaton that emerges preserves the cone angle, the topological charges and the amplitudes of the inward Hänkel components of all input Bessel beams. Although the dissipative center has a more complex structure of intertwined hot spots and vortices, it reflects up to certain extent the structure of the linear tails, and therefore that of the input Bessel beam superposition, the dissipation pattern being then quite controllable. Numerical simulations also show that dissipatons are stable under large perturbations strongly distorting the dissipative center. In fact, the only condition for the whole dissipaton structure to emerge spontaneously is the presence of the inward currents carried by its inward Hänkel beam components.

II Basic relations

According to the above, we consider the propagation of a paraxial light beam $E=\mathrm{Re}\left\{A\exp[i(kz-\omega t)]\right\}$, with carrier frequency $\omega$, propagation constant $k=n\omega/c$ ($n$ is the linear refractive index and $c$ the speed of light in vacuum) and with complex envelope $A$. We assume that cubic Kerr nonlinearity and nonlinear absorption are the dominating nonlinearities of the medium. As with nonlinear BVBs, other dispersive nonlinearities (quintic, saturable) could be considered as well without substantial changes in the results. The evolution of the envelope is then ruled by the nonlinear Schrödinger equation (NLSE)

[TABLE]

where $\Delta_{\perp}$ is the transverse Laplacian given by $\Delta_{\perp}=\partial_{r}^{2}+(1/r)\partial_{r}+(1/r^{2})\partial_{\varphi}^{2}$ in cylindrical coordinates $(r,\varphi,z)$, $n_{2}>0$ is the nonlinear refractive index, and $\beta^{(M)}>0$ is the $M$-photon-absorption coefficient.

At low enough intensities (neglecting the nonlinear terms), Eq. (1) is satisfied by the linear (and paraxial) Bessel beams $A(r,\varphi,z)\propto J_{s}(k\theta r)e^{is\varphi}e^{i\delta z}$ of cone angle $\theta$, with shortened axial wave vector by $\delta=-k\theta^{2}/2<0$, and carrying a vortex of topological charge $s$. To investigate on nonlinear beams characterized also by a cone angle $\theta$ but at high intensities, we find it convenient to use the scaled and dimensionless radial coordinate, propagation distance and envelope

[TABLE]

[in many Figures below, we will also use the dimensionless cartesian coordinates $(\xi,\eta)\equiv k\theta(x,y)$], to rewrite Eq. (1) as

[TABLE]

where now $\Delta_{\perp}=\partial_{\rho}^{2}+(1/\rho)\partial_{\rho}+(1/\rho^{2})\partial_{\varphi}^{2}$, and

[TABLE]

With these variables, the linear Bessel beams reads as $\tilde{A}\propto J_{s}(\rho)e^{is\varphi}e^{-i\zeta}$. Values of $|\tilde{A}|^{2}$ of the order of unity correspond to the typical values of the intensities involved in filamentation experiments in solids or gases, e. g., at $800$ nm with cone angle $\theta=1^{\circ}$ in water ($n=1.33$, $n_{2}=2.6\times 10^{-16}$ cm2/W, $M=5$, and $\beta^{(5)}=8.3\times 10^{-50}$ cm7/W4 COUAIRONSpie ), or $0.1^{\circ}$ in air ($n\simeq 1$, $n_{2}=3.2\times 10^{-19}$ cm2/W, $M=8$, and $\beta^{(8)}=1.8\times 10^{-94}$ cm13/W7 COUAIRONSpie ), $|\tilde{A}|^{2}\simeq 1$ corresponds to $4$ TW/cm2 and $20$ TW/cm2, respectively. The value of the parameter $\alpha$ is determined by the cone angle and the medium properties at the light wavelength. At $800$ nm in water, the values of $\alpha$ range from $0.17$ to $5.3$ for cone angles $\theta$ from $10^{\circ}$ to $1^{\circ}$, or in air, $\alpha$ varies between $0.08$ and $4.2$ for $\theta$ varying between $1^{\circ}$ and $0.1^{\circ}$.

We will further write the complex envelope as $\tilde{A}=\tilde{a}e^{i\Phi}$, where $\tilde{a}>0$ and $\Phi$ are its real amplitude and phase. The NLSE (3) is then seen to be equivalent to the energy transport and eikonal equations

[TABLE]

where $\nabla_{\perp}=\mathbf{u}_{\rho}\partial_{\rho}+\mathbf{u}_{\varphi}(1/\rho)\partial_{\varphi}$, $\mathbf{u}_{\rho}$ and $\mathbf{u}_{\varphi}$ are unit radial and azimuthal vectors, and $\mathbf{j}=\tilde{a}^{2}\nabla_{\perp}\Phi=\mbox{Im}\{\tilde{A}^{\star}\nabla_{\perp}\tilde{A}\}$ is the current of the intensity $\tilde{a}^{2}$ in the transversal plane. Nonlinear BVBs were described as particular solutions to Eq. (5) and Eq. (6) that generalize Bessel beams to nonlinear media, whose intensity profiles are propagation-invariant and circularly symmetric PORRASPrl ; PORRASJosaB ; JUKNAOe ; XIESci . Here we search for conical beams of the most general form $\tilde{A}=\tilde{a}(\rho,\varphi)e^{i\phi(\rho,\varphi)}e^{-i\zeta}$, whose intensity profiles $\tilde{a}^{2}(\rho,\varphi)$ are also localized but may no have any particular symmetry. According to Eqs. (5) and (6), their amplitude and phase patterns, $a(\rho,\varphi)$ and $\phi(\rho,\varphi)$, must satisfy

[TABLE]

where $\mathbf{j}=\tilde{a}^{2}\nabla_{\perp}\phi$ is the $\zeta$-independent intensity current. Equation (7) expresses the general condition for a conical beam to propagate without any change in a medium with absorption. Integrating over an arbitrary area $S$ in the transversal plane, and using the divergence theorem, this condition is readily seen to be equivalent to

[TABLE]

or $-F_{C}=N_{S}$ for short, where $dl$ is a differential element of the contour $C$ of $S$, $\mathbf{n}$ is the outward normal unit vector, and $ds$ and elemental area of $S$. This condition thus establishes that the power loss rate $N_{S}$ in any region $S$ of the transversal plane is offset by an inward power flux $-F_{C}$ though its contour for the stationary propagation to be possible. In particular, applying Eq. (9) to a circle of large enough radius $\rho$ for the total nonlinear power losses to occur inside it (assumed they are finite), this condition requires

[TABLE]

for $\rho\rightarrow\infty$, i. e., an inward power flux from the power reservoir of the conical beam that is independent of the large radius $\rho$, and that replenishes the power losses,

[TABLE]

in the entire transversal section.

A direct search of solutions of the nonlinear partial differential equations (7) and (8) without simplifying assumptions, such as the circular symmetry, separability of the intensity and phase patterns in $\rho$ and $\varphi$, or particular boundary conditions, is complicated. An alternate procedure is suggested by a well-established connection between linear Bessel beams and nonlinear BVBs in the circularly symmetric case, as summarized below.

III Circularly symmetric case: Nonlinear Bessel vortex beams

If a powerful Bessel beam is introduced into the medium, its nonlinear dynamics is determined by the existence and properties of an attracting nonlinear BVB PORRASPra1 ; PORRASPra2 . Given the input Bessel beam $\tilde{A}(\rho,\varphi,\zeta=0)=b_{L}J_{s}(\rho)e^{is\varphi}$, or equivalently $\tilde{A}(\rho,\varphi,\zeta=0)=\frac{1}{2}\left[b_{L}H_{s}^{(1)}(\rho)+b_{L}H_{s}^{(2)}(\rho)\right]e^{is\varphi}$, the attracting nonlinear BVB, $\tilde{a}(\rho)e^{i\phi(\rho)}e^{is\varphi}e^{-i\zeta}$, has the same cone angle and topological charge as the linear Bessel beam, and behaves asymptotically at large radius $\rho$, where it remains linear, as the unbalanced Bessel beam $\tilde{A}\simeq\frac{1}{2}\left[b_{\rm out}H_{s}^{(1)}(\rho)+b_{\rm in}H_{s}^{(2)}(\rho)\right]e^{is\varphi}e^{-i\zeta}$, with the property that $|b_{\rm in}|=b_{L}$ PORRASJosaB ; PORRASPra1 ; PORRASPra2 . The conservation of the cone angle, the topological charge and the amplitude of the inward Hänkel beam component specifies unambiguously an attracting nonlinear BVB for any input Bessel beam. Using the equivalent expressions

[TABLE]

of the Hänkel functions at large $\rho$, the outward Hänkel beam, $(b_{\rm out}/2)H_{s}^{(1)}(\rho)e^{is\varphi}e^{-i\zeta}$ is seen to carry the current $\mathbf{j}=(|b_{\rm out}|^{2}/2\pi\rho)[\mathbf{u}_{\rho}+(s/\rho)\mathbf{u}_{\varphi}]$ spirally outwards, while the inward Hänkel beam $b_{\rm in}H_{s}^{(2)}(\rho)e^{is\varphi}e^{-i\zeta}$, the current $\mathbf{j}=(|b_{\rm in}|^{2}/2\pi\rho)[-\mathbf{u}_{\rho}+(s/\rho)\mathbf{u}_{\varphi}]$ spirally inwards. From the asymptotic form of nonlinear BVBs and Eq. (10), one then gets the relation $-F_{\infty}=|b_{\rm in}|^{2}-|b_{\rm out}|^{2}=N_{\infty}$ for the inward power flux restoring the total power losses in the nonlinear BVB.

The particular Bessel beam dynamics depends on the stability properties of the attracting nonlinear BVB. According to the linear-stability analysis in Refs. PORRASPra1 ; PORRASPra2 , there exist nonlinear BVBs that are stable against radial and azimuthal perturbations. In case of stability, the Bessel beam dynamics ends in the complete formation of the nonlinear BVB. The example of Fig. 1 summarizes these previous results. The input Bessel beam represented by the light gray curve in Fig. 1(a) transforms into the nonlinear BVB represented by the black curve, having a radial inward component of the intensity current that sustains the stationarity, as seen in Fig. 1(b). After a transient, the peak intensity and nonlinear power losses $N_{\infty}$ reach the constant values corresponding to the attracting nonlinear BVB [light gray curves in Fig. 1(c)].

The fact that the linearly-stable nonlinear BVB is formed starting from so distant initial condition as the linear Bessel suggests stability under large perturbations within a certain attraction basin. Additional numerical simulations indicate that the same final nonlinear BVB is formed when the input beam carries at large radius only the inward Hänkel beam component that is preserved in the propagation. The attraction basin of the nonlinear BVB then includes all beams of the type $\tilde{A}=\frac{1}{2}b_{L}H_{s}^{(2)}(\rho)e^{is\varphi}+f(\rho,\varphi)$, where $f(\rho,\varphi)$ does not carry any additional inward radial current at large radius. The dark gray curves in Figs. 1(a) and (c) show the radial intensity profile of an input beam with these characteristics that transforms into the same nonlinear BVB, and the peak intensity and power losses reaching the same constant values as with the input Bessel beam.

On the opposite side, if the attracting nonlinear BVB determined by the conservation of the cone angle, the topological charge and the amplitude of the inward Hänkel beam component turns out to be unstable (according to the linear-stability analysis), the Bessel beam dynamics consists on the growth of its unstable modes, developing into large, periodic or chaotic oscillations around the nonlinear BVB if radial instability is dominant PORRASPra1 , or leading to circular symmetry breaking and a complex dynamics of rotatory or random spots if the azimuthal instabilities dominate PORRASPra2 .

IV Dissipatons

With the above in mind, we investigate on the beam propagation dynamics that results from illuminating the medium with the coherent superposition $\tilde{A}(\rho,\varphi,\zeta=0)=\sum_{j}b_{L,j}J_{s_{j}}(\rho)e^{is_{j}\varphi}$ of several Bessel beams with same cone angle but different topological charges $s_{j}$ and amplitudes $b_{L,j}$. A large variety of input beam profiles can be obtained with different choices of the parameters of these superpositions, which are easily generated in laboratories FROEHLYJosaA ; CLARKOe ; VALILYEUOe . In this section we consider the most general situation where this input beam does not have any particular symmetry in its transversal sections. As above, we will also consider illuminating the medium with other beams carrying the same inward intensity currents, as $\tilde{A}(\rho,\varphi,\zeta=0)=\sum_{j}\frac{1}{2}b_{L,j}H_{s_{j}}^{(2)}(\rho)e^{is_{j}\varphi}+f(\rho,\varphi)$ at large radius.

In Fig. 2 we have chosen two superposed beams with $s_{1}=2$, $b_{L,1}=2$ and $s_{2}=1$, $b_{L,2}=1$ in a medium and with cone angle such that $M=5$ and $\alpha=0.5$. The intensity patterns of the input Bessel beam superposition [Fig. 2(a)] and of the input Hänkel beam superposition [Fig. 2(b)] differ substantially, so also their propagation in the medium. The input Bessel superposition is strongly absorbed at first, while input Hänkel superposition builds up [Fig. 2(c)]. With both input beams, however, the peak intensity and total power losses reach the same positive constant values at long distances. These values correspond to the same final, nonlinear stationary state that attracts the two beams, whose intensity and intensity current patterns in a transversal plane are shown in Figs. 2(d) and (e).

These numerical simulations, and others below, show that the Kerr medium with nonlinear losses supports stationary propagation of beams that do not have any specific symmetry, to which we have referred in the introduction as “dissipatons”. In the example of Fig. 2, the final state has even lost the bilateral symmetry of the input beam. The two vortices of unit charge in the final state are slightly displaced in the dissipaton from those in the input Bessel beam superposition.

In Figs 3(a) to (d) we have plotted the peak intensity and nonlinear power power losses as functions of the propagation distance for the same couple of input Bessel and Hänkel superpositions as above but several increasing values of $\alpha$ (e. g., increasing $n_{2}$, or diminishing the cone angle). Above a certain threshold, say $\alpha_{\rm max}$, ($\alpha_{\rm max}\simeq 1.5$ in this example) that depends on the particular input superposition and $M$, the dynamics does not lead to a stationary propagation regime, regardless of whether the Bessel or the Hänkel superpositions are introduced in the medium. The propagation-invariant intensity profiles of the dissipatons from $\alpha\simeq 0$ (large enough cone angle or negligible Kerr nonlinearity) and to the maximum value $\alpha_{\rm max}\simeq 1.5$ [Figs. 3 (e) and (f)] feature increasing compression and irregularity of the nonlinear core.

Similar results are seen to hold for other input beams. In the example of Fig. 4 we wanted to produce a triangular dissipaton with a a fourth hot spot in its center. For this we chose Bessel or Hänkel superpositions with topological charges $s_{1}=3$ and $s_{2}=0$. All intensity patterns in Fig. 4 correspond to the same cone angle and the same medium ($\alpha=0.5$ and $M=5$), but increasing intensities ($b_{L,1}=b_{L,2}\equiv b_{L}=0.5,1.5$ and $3$). As above, these patterns become increasingly irregular, and no stationary state is reached above a certain value of $b_{L}$ that depends on the values of other parameters (above $b_{L}=3$ in this example). Figure 4 also shows the respective current patterns $\mathbf{j}$ and their associated sinks equal to the dissipation patterns $|\tilde{A}|^{2M}$. An input power such that $b_{L}=1.5$ produces four approximately equal dissipative spots. These examples, and others considered in Sec. V, illustrate that dissipation only occurs in strongly localized spots and not in the outer rings of the intensity pattern, and that these dissipative spots can be organized almost arbitrarily with an adequate choice of the topological charges and intensities of the input beam superposition.

Asymptotic behavior of dissipatons

The structure of the non-dissipative outer part of the dissipatons is found to be simpler than that of the dissipative spots. At fixed azimuthal angle $\varphi$ and increasing radius the intensity is always seen to decay as $\rho^{-1}$ while oscillating harmonically. The same happens with the circularly symmetric nonlinear BVBs PORRASPrl ; PORRASJosaB , but now the decay rate and the oscillation contrast depend on the specific azimuthal angle. This leads us to conjecture a linear conical structure for these tails of the form

[TABLE]

i. e., a linear superposition of unbalanced Bessel beams. Using Eq. (12) for Hänkel beams at large radius, Eq. (13) can also be written as

[TABLE]

where we have defined

[TABLE]

Equation (14) then yields an intensity profile at large radius given by

[TABLE]

where $\kappa(\varphi)=\mbox{arg}b_{\rm out}(\varphi)-\mbox{arg}b_{\rm in}(\varphi)$. As a function of $\rho$, Eq. (17) represents harmonic oscillations of $\varphi$-dependent contrast

[TABLE]

about the $\varphi$-dependent mean value

[TABLE]

The linearity of these tails allows us to generalize the result in the circularly symmetric case PORRASJosaB that the cone angle, all topological charges $s_{j}$ and amplitudes of the inward Hänkel beams in the asymptotic expression (13) of a dissipaton are the same as those of the input beam, i. e.,

[TABLE]

Also, the asymptotic expression for the radial component of the intensity current can be evaluated to be $j_{\rho}\simeq\mbox{Im}\left\{\tilde{A}^{\star}\partial_{\rho}\tilde{A}\right\}=\left[|b_{\rm out}(\varphi)|^{2}-|b_{\rm in}(\varphi)|^{2}\right]/(2\pi\rho)$. From Eq. (10), we obtain the constant inward radial flux at large radius offsetting the total power losses as

[TABLE]

and upon integration, the relation

[TABLE]

between the amplitudes of the inward and outward Hänkel beam components.

The validity of all previous relations is supported by the numerical analysis. As an example, Figs. 5 (a) and (b) show the intensity profile of a dissipaton [the same as in Figs. 3 (a) and (e)], and the fitted intensity profile (for $\rho>5$) using Eq. (13) with $b_{{\rm in},j}$ and $b_{{\rm out},j}$ satisfying Eqs. (20) and (22) (see caption for detailed values). The radial profiles of the dissipaton and of the asymptotic form at two particular azimuthal angles are more clearly seen to match better and better with increasing radius in Fig. 5 (c). With the values of $b_{{\rm in},j}$ and $b_{{\rm out},j}$, one can evaluate $b_{\rm in}(\varphi)$ and $b_{\rm out}(\varphi)$ from Eqs. (15), and the contrast $C(\varphi)$ and average oscillation value $R(\varphi)$ of the radial oscillations at each azimuthal direction from Eqs. (18) and (19). These are plotted in Figs. 5(d) and (e) in the above example. Compared to the same quantities for the input Bessel beam superposition (dashed curves), the average value of the oscillations is significantly smaller due to the initial stage of strong absorption. It is also characteristics of dissipatons that the contrast of the oscillations is significantly diminished along the azimuthal directions of high intensity, i. e., high $R(\varphi)$, as for the angle $\varphi=3\pi/2$ in this example.

Additional confirmation of the generality of the above description is given in Sec. V for dissipatons with $n$-fold rotational symmetry. By analogy with nonlinear BVBs, each dissipaton is determined by the set of parameters $M$, $\alpha$, $s_{j}$ and $b_{L,j}$. We can conjecture, following the same analogy, that there exist dissipatons with any value of these parameters, as for nonlinear BVBs PORRASJosaB , and that they may be stable or unstable for different choices of the parameters, e. g., unstable above $\alpha_{\rm max}$ once the other parameters are fixed. This hypothesis is suggested by the fact that when the propagation does not lead to a stationary state, as in Fig. 3(d), the beam is not completely dissipated or dispersed, but its peak intensity, power losses, and the entire beam profile oscillates more or less randomly about constant values, as in the circularly symmetric case PORRASPra2 .

V Dissipatons with $2s$-fold rotational symmetry

Of particular interest are the dissipatons with $2s$-fold rotational symmetry formed from two Bessel beams

[TABLE]

of opposite topological charges $\pm s$. Without loss of generality, we can choose real and positive $b_{L,s}$ and $b_{L,-s}$, since the only effect of a phase shift between them is a rigid beam rotation. It is also convenient to fix their values as $b_{L,s}^{2}=b_{L}^{2}[1+\sqrt{1-C_{L}^{2}}]/2$ and $b_{L,-s}^{2}=b_{L}^{2}[1-\sqrt{1-C_{L}^{2}}]/2$, with $0\leq C_{L}\leq 1$, yielding the circularly symmetric case $b_{L}J_{s}(\rho)e^{is\varphi}$ for $C_{L}=0$ ($b_{L,s}=b_{L}$, $b_{L,-s}=0$). The choice $C_{L}>0$ ($b_{L,s}>b_{L-s}$) introduces an azimuthal modulation of period $\pi/s$ and depth $C_{L}$ in the intensity pattern, but it continues to have a vortex of charge $s$ in the origin. In the limit case with $C_{L}=1$ ($b_{L,s}=b_{L,-s}$), the intensity vanishes at $2s$ azimuthal directions. Two examples of these input beams with respective low and high modulation depths ($C_{L}=0.3$ and $C_{L}=0.98$) are shown in Figs. 6(a-c) and 7(a-c) for $s=2$. With increasing $C_{L}$, the azimuthal change in $2\pi s$ of the phase about the vortex approaches an step-like variation, reaching $2s$ discontinuous jumps of $\pi$ in the limit $C_{L}=1$, the intensity current $\mathbf{j}=(s/\rho)(b_{L,s}^{2}-b_{L,-s}^{2})J_{s}^{2}(\rho)\mathbf{u}_{\varphi}$ is always purely azimuthal, but increasingly weaker, and the associated angular momentum density $L=\tilde{a}^{2}\partial\Phi/\partial\varphi=\mbox{Im}\{\tilde{A}^{\star}\partial\tilde{A}/\partial\varphi\}=s(b_{L,s}^{2}-b_{L,-s}^{2})J_{s}^{2}(\rho)$ diminishes, vanishing in the limit $C_{L}=1$. We focus on these linear Bessel superpositions as input beams because the formed dissipatons share a number of their properties, but we stress that the same dissipatons are formed from other beams with equivalent asymptotic inward currents, that is, from $\tilde{A}(\rho,\varphi,0)=b_{L,s}H_{s}^{(2)}(\rho)e^{is\varphi}+b_{L,-s}H_{-s}^{(2)}(\rho)e^{-is\varphi}+f(\rho,\varphi)$, where $f(\rho,\varphi)$ does not carry any additional inward radial current at large radius, that do not have any of these properties around the beam center. There may not even be light in the central region of the input beam.

The intensity, phase and current patterns of the respective dissipatons formed after a transient propagation stage are depicted in Figs. 6 (d-f) and Figs. 7 (d-f). The corresponding dissipation patterns with four localized dissipation spots forming a square are shown in the insets of Figs. 6 (d) and 7 (d). The dissipathons always preserve the $2s$-fold rotational symmetry of the input Bessel beam superposition, the vortex of the same charge $s$ at $\rho=0$, and an azimuthal modulation of the intensity profile of similar depth as that of the input beam. We stress that symmetrization into a circularly symmetric nonlinear BVB is never observed, irrespective of the smallness of $C_{L}$ of the input beam. As seen below, this is a consequence of the conservation of the inward Hänkel amplitudes. The step-like azimuthal variation of the phase profiles also preserves approximately the same sharpness as that of the input beam, although it is slightly changing with radial distance $\rho$. In the limit of an input beam with $C_{L}=1$, the formed dissipaton continues to have $2s$ null directions, and the azimuthal variation of the phase $2s$ discontinuous jumps. In all cases, a substantial difference from the input beam is that the azimuthal variation of the phase is accompanied by a pronounced radial variation, associated with the permanently converging wave fronts and the permanent radial inward component of the intensity current sustaining the stationarity. For increasing azimuthal input contrast, the radial inward component is accompanied by decreasing azimuthal component and angular momentum, that vanish in the limit $C_{L}=1$.

In the above examples, we chose $\alpha\simeq 0$ for simplicity, but similar properties hold for $0\leq\alpha\leq\alpha_{\rm max}$, where $\alpha_{\rm max}$ depends now on $s$, $b_{L}$, $C_{L}$ and $M$. For $\pm s=\pm 1$ and $\pm s=\pm 3$, the intensity pattern of the dissipatons with $2$-fold and $6$-fold rotational symmetry formed in a medium with $M=5$ and cone angles and Kerr nonlinearities such that $\alpha\simeq 0$ and $\alpha_{\rm max}$ are depicted in Fig. 8. The insets show the corresponding two-spot and six-spot dissipation patterns. We note that even if the intensity profile is weakly localized and becomes increasingly irregular with increasing $\alpha$, these dissipation patterns continue to be extremely localized and ordered.

The above structure of these dissipatons is reflected in the structure beyond the dissipative region, where it admits an analytical description. The asymptotic form in Eq. (13) for these dissipatons reads as

[TABLE]

Equivalently, Eq. (14) also holds, where now

[TABLE]

represent harmonic azimuthal oscillations of period $2\pi/s$. Consequently, the contrast $C(\varphi)$ of the radial oscillations in the intensity profile along each azimuthal angle $\varphi$, and its average value $R(\varphi)$, repeat themselves each period $\pi/s$, as expected. For the dissipatons in Figs. 6 and 7, Figs. 9 (a,b) and 10 (a,b) depict the intensity profiles and their asymptotic forms satisfying the conservation of the inward Hänkel amplitudes, $|b_{\rm in,s}|=b_{L,s}$, $|b_{\rm in,-s}|=b_{L,-s}$, and the replenishment condition, $|b_{\rm in,s}|^{2}+|b_{\rm in,-s}|^{2}-|b_{\rm out,s}|-|b_{\rm out,-s}|^{2}=N_{\infty}$. From these values, the contrast of the radial oscillations at each azimuthal direction, $C(\varphi)$, are evaluated from Eqs. (25) and (18), and are plotted in Figs. 9 (c) and 10(c) as solid blue curves. As for nonlinear BVBs, dissipation manifests in a loss of radial oscillation contrast (from unity for the input beam), with the peculiarity that the decrease depends on the azimuthal direction and is more pronounced along those of higher intensity (e. g., at $\varphi=0$).

For dissipatons with $2s$-fold rotational symmetry, using the asymptotic forms of Hänkel functions in Eq. (12) again, but grouping now the terms with the same azimuthal dependence, Eq. (24) we can alternatively be written as

[TABLE]

where we have defined

[TABLE]

In this form, the intensity profile at large radius is also given by

[TABLE]

As a function of $\varphi$, these are harmonic oscillations of period $\pi/s$ of $\rho$-dependent contrast

[TABLE]

about a $\rho$-dependent average value $R(\rho)=|b_{s}(\rho)|^{2}+|b_{-s}(\rho)|^{2}$. The azimuthal contrast as a function of the radial distance for the dissipatons in Figs. 9 (a) and 10(a) are plotted as black curves in Figs. 9(c) and 10(c). As seen, the azimuthal contrast, which was independent of $\rho$ and equal to $C_{L}$ for the input Bessel beam superposition, oscillates with radial distance, remaining however of the same order as for the input beam, $C_{L}$, as observed from the numerical simulations. In the particular dissipatons formed from $C_{L}=1$ ($b_{L,s}=b_{-L,s}=b_{L}/\sqrt{2}$) symmetry considerations ($b_{{\rm in},s}=b_{{\rm in},-s}$, and $b_{{\rm,out},s}=b_{{\rm out},-s}$) imply $b_{-s}(\rho)=(-1)^{s}b_{s}(\rho)$, and therefore $C(\rho)=1$ too, i. e., the dissipaton preserves the $2s$ directions of zero intensity of the input Bessel beam superposition.

In addition to the inward radial component of the intensity flux, these dissipatons with $2s$-fold symmetry carry the azimuthal component $j_{\varphi}\simeq(s/2\pi\rho)[|b_{s}(\rho)|^{2}-|b_{-s}(\rho)|^{2}]/\rho$ at large radius. Interestingly, while the intensity and phase profiles exhibit pronounced azimuthal variations, the azimuthal flux is circularly symmetric. Thus, the vector expression of the intensity current is

[TABLE]

Accordingly, the angular momentum density is given by

[TABLE]

and is also circularly symmetric. These properties are also of interest in order to evaluate the azimuthal derivative of the phase as $\partial\Phi/\partial\varphi=L/|\tilde{A}|^{2}$. Upon integration with the above asymptotic expressions of $L$ and $|\tilde{A}|^{2}$, the phase is asymptotically given by

[TABLE]

with $\kappa(\rho)=[\mbox{arg}\,b_{s}(\rho)-\mbox{arg}\,b_{-s}(\rho)]/2$, and aside the term $-\zeta$ on propagation. From Eqs. (28) and (31), the above relation can also be expressed as

[TABLE]

relating the azimuthal variation of the phase with that of the intensity. The azimuthal phase variation of the dissipatons of Figs. 9 (a) and 10 (a) are plotted in the respective Figs. 9 (d) and 10 (d) at a particular radius. These explain the increasing sharpness of the azimuthal phase variation with increasing azimuthal contrast of the intensity $C(\rho)$, and that in the limit $C(\rho)=1$ (obtained from $C_{L}=1$), the phase of the azimuthons continues to display $2s$ discontinuous $\pi$ jumps of the input Bessel beam superposition.

Dissipatons with $2s$-fold rotational symmetry share with azimuthons, and their limiting cases of vortex solitons and soliton clusters, similar azimuthal modulation of the intensity patterns and similar stair-like azimuthal variation of the phase DESYATNIKOVPrl1 ; MINOVICHOe , but as pointed out in the introduction, they are supported by completely different mechanisms. Also, azimuthons are generally rotating, and the number of azimuthal maxima and vortex charge in their center are independent properties, while these dissipathons are static, and the number of maxima is $2s$ if the charge is $s$.

VI Conclusions

In conclusion, we have described a broad family of propagation-invariant and robust conical beams in media with Kerr-type nonlinearities and nonlinear absorption, with the property of continuously dissipating power along narrow and ordered channels, whose number and geometrical disposition can be easily controlled. We have described semi-analytically their spatial structure and found the laws that govern their formation from superpositions of Bessel beams.

In filamentation experiments, the new degrees of control offered by these beams with regard to the number and disposition of plasma channels may open new possibilities in short-range applications such as material laser ablation, wave guide writing, and in long-range applications such as filamentation guiding of weak waves in the atmosphere COURVOISIERApplPhys ; BHUYAN ; KRISNALaserPhysics ; ARNOLDJPhysB ; COURVOISIER2 .

M.A.P. acknowledges support from Projects of the Spanish Ministerio de Economía y Competitividad No. MTM2012-39101-C02-01, MTM2015-63914-P, and No. FIS2013-41709-P.