1. Introduction

Spatial trapping of light in nonlinear media is a result of a balance between diffraction and the self-induced nonlinear index change [1, 2]. As the rate of diffraction is determined by the spatial scale of the light beam independently of its dimensionality, the stability properties of the self-trapped beams critically depend on the particular model of the nonlinear response of the medium. Commonly encountered in the context of nonlinear optics the, so called, Kerr nonlinearity with nonlinear index change being proportional to the light intensity leads to stable spatial solitons only in one transverse dimension. For two dimensional beams the same model does not permit stable beam localization. Instead, it predicts either eventual diffraction or catastrophic collapse for input power below or above a certain threshold value, respectively [3]. Hence, in order to ensure the existence of stable self-localized beams in two dimensions a different nonlinear response is necessary. It is known that the saturation of nonlinearity at high intensities is sufficient to prevent collapse and provide stabilization of finite beams (see [4] and references therein). However, saturable nonlinearity cannot stabilize complex localized structures such as vortex beams which are prone to azimuthal instability [5, 6]. It has been shown recently that stabilization of localized waves is greatly enhanced in nonlocal nonlinear media. In such media the nonlinear response in a particular spatial location is typically determined by the wave intensity in a certain neighborhood of this location. Nonlocality often results from certain transport processes such as atomic diffusion [7], ballistic transport of hot atoms or molecules [8] or heat transfer [9, 10]. It can be also a signature of a long-range interparticle interaction such as in nematic liquid crystals [11, 12, 13]. Spatially nonlocal response is also naturally present in atomic condensates where it describes a noncontact bosonic interaction [14, 15, 16]. Extensive studies of beam propagation in nonlocal nonlinear media revealed a range of interesting features. In particular, it has been shown that nonlocality may affect modulational instability of plane waves [17] and prevent catastrophic collapse of finite beams [18, 19] as well as stabilize complex one-, two- and three-dimensional beams including vortices [20, 21, 22, 23, 24, 25, 26]. Recently, it has been shown that nonlocal media can also support stable propagation of rotating solitons, the so called azimuthons [27]. Azimuthons have been originally proposed in the context of local nonlinear media [28]. They are azimuthally modulated beams with nontrivial phase structures exhibiting steady angular rotation upon propagation. They can be considered as azimuthally perturbed optical vortices, i.e. beams with singular phase structure [29, 30]. Recent theoretical studies demonstrated both, stable and unstable evolution of azimuthons [31, 32, 33]. In the latter case it has been shown that in a highly nonlocal regimes azimuthons can undergo structural transformation resulting from the energetic coexistence of solitons of different symmetries. Numerical and analytical studies revealed that the angular velocity of azimuthons is governed by two contributions. The linear component, determined solely by the spatial structure of the beam (akin to rotation of complex wave structures resulting from beating of their constituent modes [34]) as well as the nonlinear which is brought about by the nonlinearity [31].

In this work we investigate in details rotating nonlocal solitons. We consider various non-local models and show how they determine dynamical properties of azimuthons. In particular, we will show that both, rotational frequency and intensity profile of the azimuthons can be uniquely determined by analyzing eigenmodes of the linearized version of the corresponding nonlocal problem. The paper is organized as follows: In Sec. 2 we introduce the model equations, the general ansatz for rotating solitons (azimuthons) and an exact expression for their rotation frequency. In Sec. 3, we derive an approximative, variational formula for the rotation frequency of the dipole azimuthons. In Sec. 4 we confront our semi-analytical results with exact numerical solutions of the ensuing nonlinear equations. One of our main findings is that azimuthons emerge from internal modes of stationary nonlinear soliton solutions like vortex or dipole, which is detailed in Sec. 5.

2. Rotating solutions

We consider physical systems governed by the two-dimensional nonlocal nonlinear Schrödinger equation

where θ represents the spatially nonlocal nonlinear response of the medium. Its form depends on the details of a particular physical system.

In the so called Gaussian model of nonlocality, θ is given as

where r⃗=xe⃗_x+ye⃗_y denotes the transverse coordinates.

If θ is governed by the following diffusion-type equation

solving formally Eq. (3) in the Fourier space yields

where 𝔎₀ is the modified Bessel function of the second kind.

In the following, we will assume that the nonlinear response θ can be expressed in terms of the nonlocal response function R(r)

and Eqs. (2) and (4) will serve as prominent examples.

Azimuthons are a straightforward generalization of the usual ansatz for stationary solutions (solitons) [28]. They represent spatially rotating structures and hence involve an additional parameter, the angular frequency Ω (see also [35])

where U is the complex amplitude function and λ the propagation constant. For Ω=0, azimuthons become ordinary (nonrotating) solitons. The most simple family of azimuthons is the one connecting the dipole soliton with the single charged vortex soliton [27]. A single charged vortex consists of two equal-amplitude dipole-shaped structures representing real and imaginary part of U. If these two components differ in amplitudes the resulting structure forms a “rotating dipole” azimuthon. If one of the components is zero we deal with the (nonrotating) dipole soliton. In the following we will denote the ratio of these two amplitudes by α, which also determines the angular modulation depth of the resulting ring-like structure by “1-α”.

After inserting the ansatz (6) into Eqs. (1) and (5), multiplying with U* and ∂_ϕ U* resp., and integrating over the transverse coordinates we end up with

This system relates the propagation constant λ and the rotation frequency Ω of the azimuthons to integrals over their stationary amplitude profiles, namely

The first two quantities have straightforward physical meanings, namely “mass” (M) and “angular momentum” (L_z). We can formally solve for the rotation frequency and obtain (for an alternative derivation see [36])

Note that this expression is undetermined for a vortex beam (see discussion below).

While the above formula looks simple its use requires detailed knowledge of the actual solution which can only be obtained numerically. However, it turns out that it is still possible to analyze main features of rotating solitons using a simple approximate approach. To this end let us consider the “rotating dipole”. Obviously, we have Ω=0 for α=0 (dipole soliton). On the other hand, for α=1 [vortex soliton V(r)exp(iϕ+iλ ₀ z)], we can assume any value for Ω by just shifting the propagation constant λ=λ ₀+Ω accordingly (λ ₀ accounts for the propagation constant in the non-rotating laboratory frame). However, with respect to the azimuthon in the limit α→1, the value of Ω is fixed. In what follows, we denote this value by Ω|_α=1.

3. Rotation frequency - approximate analysis

In a first attempt, let us assume we know the radial shape of the vortex soliton V(r). Then, a straight forward approximative ansatz for the dipole azimuthons is [33, 31]

where A is an amplitude factor. The family of the “rotating dipole” azimuthons has two parameters, e.g. λ and Ω. Here, for technical reasons, we prefer to use the mass M and the amplitude ratio α. We take the radial shape of the vortex for a given mass M ₀ and fix the amplitude factor A in ansatz (9) to $A=\sqrt{\frac{2}{\left({\alpha }^2+1\right)}}$ . With this choice the mass of the azimuthon equals M ₀ for all values of α. Then, we can compute all integrals in Eqs. (7) as functions of α: M=M ₀, L_z=2αM ₀/(α ² +1), I=I ₀, N=[4(1+α ⁴)N _cc +8α ² N _cs]/(α ² +1)², M′=M ₀, I′=2αI ₀/(α ²+1), and N′=4α(N _cc+N _cs)/(α ²+1). The index “0” indicates the value of the respective integral for the vortex (i.e., for α=1), and

Finally, we obtain the following simple expression for the rotation frequency from Eq. (8),

Let us have a closer look at this approximation. First of all, we have N _cc≥N _cs if we assume R(r) to be monotonically decreasing in r > 0. This immediately implies that the sense of the amplitude rotation is opposite to that of the phase rotation. In general, we need to know the actual shape of the vortex soliton to compute Eq. (11). One possibility would be to use exact numerical solutions of the vortex shape V, a time consuming task because V depends in nontrivial way on λ ₀. Moreover, Eq. (11) is only an approximation to Ω and should not require the exact shapes. Luckily, it is possible to compute quite good approximative vortex solitons using the so called variational approach [37]. We employ the following ansatz V=Ar exp(-r ²/2σ ²), U=V exp(iϕ), ψ=U exp(iλ ₀ z) and look for critical points of the Lagrangian ℒ=-λ ₀ M+I+N/2 with respect to the variables A and σ. In the following, we give results for our two exemplary nonlocal nonlinearities.

 

Fig. 1. Nonlocal model and variational approach: Solid lines show (a) width σ, (b) amplitude A, and (c) frequency Ω|_α=1 for the Bessel response as functions of the propagation constant λ ₀. Dashed lines represent the same quantities for the Gaussian nonlocal response. The inset in (c) shows Ω|_α=1 for the latter response on more appropriate scales. While both, width σ and amplitude A, exhibit similar behavior for the two nonlocal models the rotation frequency Ω|_α=1 shows a completely different dependence on λ ₀ and drastic difference in magnitude for Gaussian and Bessel responses.

Download Full Size | PDF

For the Gaussian response, Eq. (2), all integrals appearing in the variational approach can be evaluated analytically. Hence, we find the width σ of our approximate solution as

its amplitude

and the rotation frequency Ω

We show these three quantities as functions of λ ₀ in Fig. 1 (dashed lines). If λ ₀ is large, σ≪1 (nonlocal limit) and we can neglect the first two terms in Eq. (14). Then we get σ _nl=(2/λ ₀)^1/4, A _nl=λ ₀(σ ² _nl+1)^3/2 and, subsequently, Ω_nl=α/α ²+1). The immediate conclusion drawn from this formula is that dipole azimuthons in the Gaussian nonlocal model do not rotate faster than Ω~0.5. Moreover, in the nonlocal regime the rotating frequency does not change much with the propagation constant λ or the mass M of the azimuthon, but is mainly determined by the structural parameter α.

For the Bessel nonlocal response, Eq. (4), we have to compute the relevant integrals numerically. The results for both, Bessel and Gaussian models are presented in Fig. 1 using solid and dashed lines, respectively. It is clear that amplitude A as well as soliton width σ, follow similar dependencies as a function of the propagation constant λ ₀ [see Fig. 1(a–b)]. Totally different dependence shows the rotation frequency Ω [see Fig. 1(c)]. In particular, not only azimuthons in the Bessel nonlocal model rotate much faster than those predicted by the Gaussian model but also the former exhibit strong sensitivity of Ω on the propagation constant λ and mass M.

 

Fig. 2. (a) Mass M and (b) quantity N _cc - N _cs for Bessel (solid lines) and Gaussian (dashed lines) nonlocal responses as functions of the propagation constant λ ₀. It is obvious that the difference in the rotating frequency Ω|_α=1 shown in Fig. 1(c) is due to N _cc -N _cs [see Eq. (11)]. (c) Spatial Fourier Transform (FT) of the Bessel (solid lines) and Gaussian (dashed lines) nonlocal response R as a function of the transverse wave-vector $k_{\perp }=\sqrt{k_x^2+k_y^2}$ . The degree of spatial averaging is obviously larger for the Gaussian than for the Lorentzian.

Download Full Size | PDF

To shed light on the origin of the difference in the rotating frequencies of the two models let us have a look at the quantities M and N _cc-N _cs, which directly determine Ω|_α=1 in the present approximation [see Eq. (11)]. Figures 2(a,b) reveal that as the masses M are comparable for both models, the important differences come from N _cc-N _cs. Hence, we can conclude that it is not the actual shape of the azimuthon which determines Ω, but the convolution integrals N _cc and N _cs with the response function. In fact, N _cc and N _cs represent the overlaps of the nonlocal response of one dipole component with either itself or its π/2-rotated counterpart. Hence, the more of the non-rotational-symmetric shape is preserved by the nonlocal response, the larger is N _cc -N _cs. Looking at both response functions considered here in Fourier domain [Fig. 2(c)] it is obvious that the degree of nonlocality and consequently the degree of spatial averaging is greater for the Gaussian response leading to smaller difference between N _cc and N _cs. In this context, it is worth mentioning previous works dealing with spiraling waves in dissipative systems, where, in contrast to our case, the angular frequency is uniquely determined by system parameters, namely gain and dissipation [38, 39].

4. Rotation frequency - numerical analysis

How good are our approximations of the rotating frequency Ω? Is the dependency in the structural parameter α really approximatively α/(α ²+1) for a given λ ₀? To answer these questions we compute the azimuthons U numerically and propagate them over some distance in z to observe their rotation and estimate Ω. One typical example for the Gaussian nonlocal model is shown in Fig. 3. The corresponding vortex soliton has a propagation constant λ ₀=11.5 and mass M=200. As its width σ≈0.7 is smaller than the width of the response function we are definitely in the nonlocal regime. We can clearly see the reasonable agreement between approximate and fully numerical results [compare with Fig. 1(c)]. We also tested azimuthons with larger masses and therefore corresponding to higher degrees of nonlocality. For instance, for λ ₀=247, M=2000 and α=0.95 we find Ω=0.54 from direct numerical simulations, in excellent agreement with Ω=0.52 from the approximate model [see also Fig. 1(c)]. Figure 4 shows another example, this time for the Bessel nonlocal model. Again, the agreement is reasonable. Here, we have to use higher λ ₀ and larger mass in order to ensure the stability of the vortex soliton. Like in the previous example, we operate in a nonlocal regime (σ=0.17). To summarize, our simple approximate model for Ω, which involves nothing but a variational solution for the vortex soliton, allows us to predict the correct range and sign of Ω for a family of dipole azimuthons with given mass. Moreover, the dependency on the structural parameter α is also correctly modeled.

 

Fig. 3. Gaussian nonlocal model, M=200, λ ₀=11.5: Solid line shows the result of Eq. (11) computed from the variational vortex approximation. Crosses depict results obtained from direct numerical solutions. Dashed lines represent results of the linear perturbation analysis of dipole and vortex soliton, respectively (see Sec. 5).

Download Full Size | PDF

 

Fig. 4. Bessel nonlocal model, M=2000, λ ₀=418: Same coding as in Fig. 3.

Download Full Size | PDF

However, there is still serious discrepancy in the values of Ω|_α=1. One could think of two possible reasons. Firstly, it could be the deviations between the actual vortex soliton and its variational approximation. However, it turns out that 2(N _cc -N _cs)/M computed from the numerically obtained vortex profile is almost constant (within a few percents). Hence it has no effect on the overall error in Ω|_α=1. However, postulating the specific ansatz Eq. (6) we implicitly assumed a certain shape of the deformation to the soliton profile while going over from the vortex to azimuthons, in the limit α→1. Therefore it has to be this shape of vortex deformation which determines the rotation frequency Ω, since a vortex formally allows for all possible rotations (see discussion on shifting λ in end of Sec. 2). In the next section we will therefore discuss the formation of azimuthons by considering it as a process of bifurcation from the stationary non-rotating soliton solutions.

 

Fig. 5. Dipole D and deformation d for the Gaussian nonlocal model, M=200, λ _D=11.8. The deformation d has a dipole shape, rotated by π/2 with respect to D. We can construct an approximate azimuthon U=D+iΩd for small rotation frequencies Ω. The amplitude ratio max d/max D gives the slope of Ω(α) at α=0.

Download Full Size | PDF

5. Internal modes and azimuthons

Let us start with the dipole soliton D(r,ϕ), because here the resulting linear problem is simpler than in the vortex case. We can directly use the rotation frequency Ω as the linearization parameter, since for the dipole we have Ω_D=0. We use the following ansatz

where D and d are real functions representing soliton and its small deformation, respectively. We then plug this ansatz into Eqs. (1), (5) and (6) and linearize with respect to Ω. The resulting equation of order 𝒪(Ω⁰) is fulfilled by the dipole soliton D. The terms of the order 𝒪(Ω¹) lead to the following linear equation for the deformation d:

In deriving the above equation we made use of the fact that the propagation constant λ of the azimuthon depends on Ω as λ=λ _D+𝒪(Ω²), where λ _D is the propagation constant of the dipole soliton D. This can be seen easily from Eq. (7a) and M=M _D+𝒪(Ω²), L_z=0+𝒪(Ω), I=I _D+𝒪(Ω²), and N=N _D+𝒪(Ω²), where the index “D” indicates the respective integral for the dipole soliton D.

The solution of Eq.(16) for the Gaussian nonlocal model and the dipole soliton with M=200, and λ_D=11.8 is shown in Figure 5. The left plot shows the dipole soliton D(x,y) while the right shows its deformation d(x,y). As we can see from the right subplot, the deformation d has a shape of a dipole as well, but rotated by π/2 with respect to D. The width of d is about 0.85 that of D. Hence, for small rotation frequency Ω the dipole azimuthon consists of two orthogonally oriented dipoles, with the the relative phase shift of π/2 and amplitude ratio α=Ωmaxd/maxD. In effect, the whole rotating structure is slightly elliptic, with the principal axis given by D. In our example, we find maxd/maxD=1.53, hence Ω=1.53α, which agrees perfectly with exact numerical results obtained from the solution of the original nonlinear problem Eq.(1) for α →0 (see Fig. 3).

 

Fig. 6. Same as Fig. 5 but for the Bessel nonlocal model, M=2000, λ_D=444.

Download Full Size | PDF

We obtain similar results for our second example, the dipole soliton in the Bessel nonlocal model with M=2000, λ_D=444 (see Fig. 6). The ratio of the widths of a dipole D and its deformation d is 0.76, and Ω=62.1α. Again, the agreement with the full numerical results of Fig. 4 is excellent.

Let us now look at the azimuthon originating (bifurcating) from the vortex soliton. For this purpose, we recall the eigenvalue problem for the internal modes of the nonlinear potential θ which is usually treated in the context of linear stability of nonlinear (soliton) solutions. We introduce a small perturbation δV to the vortex soliton V, and plug

into Eqs. (1) and (5) and linearize with respect to the perturbation. Please note that the perturbation δV(r,ϕ, z) is complex, whereas the vortex profile V(r) can be chosen real. The resulting evolution equation for the perturbation δV is given by

With the ansatz

we then derive the eigenvalue problem for the internal modes

where

 

Fig. 7. Gaussian nonlocal model, M=200, λ ₀=11.5: Solid line shows the component δV ₁, dashed line the component δV ₂, for m=2 and eigenvalue κ=-1.4 computed from Eq. (20). The movie (4 MB or 15 MB) shows the evolution of intensity and phase of the vortex perturbed with the eigenvector shown above, U|_z=0=V exp(iϕ)+δV ₁ exp(i3ϕ)+δV ₂ exp(-iϕ). The chosen amplitudes of the perturbation are such that α=0.85, the rotating frequency is close to Ω|_α=1=0.7. [Media 1][Media 2]

Download Full Size | PDF

Please not that since $\mid \overrightarrow{r}-\overrightarrow{r}\prime \mid =\sqrt{r^2+r^{\prime 2}-2\mathrm{rr}\prime \cos (\varphi -\varphi \prime )}$ , all integrals in (21) are independent of ϕ. Real eigenvalues of Eq. (20) (κ=κ*) are termed orbitally stable and the corresponding eigenvector (δV ₁,δV ₂) can be chosen as real function. If we perturb the vortex V with an orbitally stable eigenvector, the resulting wave-function ψ can be written in the form of Eq. (6) with Ω=-κ/m and λ=λ ₀-κ/m. Hence, we expect to find an orbitally stable internal mode with m=2 corresponding to the azimuthon, where the eigenvalue κ gives Ω|_α=1=-κ/2.

Indeed, when we solve problem (20) for our two test cases numerically we find eigenvalues at the expected positions, κ=-1.4 and κ=-62 respectively. The resulting estimates for Ω|_α=1 fit perfectly with values obtained from direct simulations (see Figs. 3 and 4). The radial profiles of the corresponding eigenvectors are shown in Figs. 7 and 8). The total vorticity of the component ~δV ₁ is 2+1=3, the one of ~δV ₂ is -2+1=-1. These total vorticities manifest themselves in the shapes of the radial profiles near the origin r=0.

As for the perturbation of the stationary dipole, it is possible to construct azimuthons in the vicinity of the vortex (α≈1) from δV:

 

Fig. 8. Bessel nonlocal model, M=2000, λ ₀=418: Same as in Fig. 7, m=2, κ=-62. The movie (4 MB or 15 MB) shows the evolution of the perturbed vortex with α=0.9, the rotating frequency is close to Ω|_α=1=31. [Media 3][Media 4]

Download Full Size | PDF

Used as an initial condition to the propagation equation this object is expected to rotate with Ω|_α=1=-κ/m. Unlike in the previous case of the stationary dipole, the amplitude of the perturbation is not fixed (just the ratio between the components δV ₁ and δV ₁ is prescribed), but will eventually determine the value of the structural parameter α. Generally speaking, the smaller the resulting α the greater the error in the constructed initial condition. However, the great robustness of the azimuthons allows one to use the initial condition (22) for quite large perturbation amplitudes (resulting α~0.5). Those “bad” initial conditions result in oscillations of the azimuthon upon propagation. However, the azimuthon is structurally stable and does not decay into other solutions like the single-hump ground state. Moreover, such initial conditions play a role of excellent “initial guesses” for solver routines to find numerically exact azimuthons. In the movies shown in Figs. 7 and 8 we choose moderate amplitudes of the perturbations in order to get almost perfect azimuthons. The corresponding structural parameters α≥0.85 also guarantee that the observable rotating frequencies are very close to their limiting values Ω| _α=1.

The results of the linear stability analysis of the vortex soliton offer an easy explanation for the observed systematic deviation of the estimation for Ω|_α=1 by Eq. (11) from its true value. As mentioned before, the shape of the perturbation to the vortex soliton determines Ω| _α=1. In our ansatz (9) we only account for a perturbation similar to the component ~δV ₂ and “forget” the other component ~δV ₁ with total vorticity 3, which leads to the observed deviation. The larger the “forgotten” component, the larger is the deviation (compare Figs. 7, 3 to Figs. 8, 4).

 

Fig. 9. Gaussian nonlocal model, M=200, λ ₀=11.5: Same as in Fig. 7, m=3, κ=3.8. The movie (4 MB or 15 MB) shows the evolution of perturbed vortex U|_z=0=V exp(iϕ)+δV ₁ exp(i4ϕ)+δV ₂ exp(-i2ϕ), α=0.9, the rotating frequency is close to Ω|_α=1=-1.27. [Media 5][Media 6]

Download Full Size | PDF

There is more to say about the internal modes of the vortex soliton. On one hand, solving the eigenvalue problem (20) for other vorticities m offers an exhaustive method for finding families of azimuthons which originate from a vortex soliton. For example, with our method it was straightforward to identify a rotating triple-hump azimuthon (m=3, κ=3.8), which rotates with Ω=-κ/m=-1.27. This solution bifurcates from the M=200 vortex in the Gaussian nonlocal model. The components of the generating eigenvector are shown in Fig. 9. Again, we can easily construct various azimuthons with α≈1 from this eigenvector. The movie in Fig. 9 shows an example of the propagation of the triple-hump azimuthon with α=0.9.

It should be stressed that not all orbitally stable eigenvalues can be linked to a family of azimuthons. This is obvious for eigenvalues with |κ | > |λ ₀| in the continuous part of the spectrum. Hence, we can conclude that |Ω|_α=1|<|λ ₀|/m for azimuthons bifurcating from a single-charged vortex. Please note that the parameter m determines the number of humps of the rotating structure. However, bound orbitally stable eigenvalues do not necessary indicate an emanating family of azimuthons. For example, the M=200 vortex in the Gaussian nonlocal model possesses several other bounded internal modes for m=2, e.g. κ=7. If we perturb the vortex with the corresponding eigenvector the resulting structure is a double-hump rotating with Ω=-3.5. Closer inspection reveals that this structure decays and hence does not belong to another family of dipole azimuthons. Nevertheless, since those decaying rotating structures can survive over large propagation distances they may be indistinguishable from true azimuthons in the experimental conditions. This issue is currently under detailed investigation.

6. Conclusion

In conclusion, we analyzed properties of azimuthons, i.e. localized rotating nonlinear waves in nonlocal nonlinear media. We showed that the frequency of angular rotation crucially depends on the nonlocal response function. Starting with the single charge vortex soliton and by employing a simple variational ansatz we were able to predict accurately the frequencies of the rotating dipole azimuthon. In addition, we could explain how the actual shape of the response function determines the rotation frequencies. Further on, we computed exact dipole azimuthon solutions and their rotation frequencies numerically and showed that in the limits of maximal and/or minimal azimuthal amplitude modulation, i.e., close to the dipole or vortex soliton, the rotation frequency is determined uniquely by eigenvalues of the bound modes of the linearized version of the respective stationary nonlocal solution. Moreover, the intensity profile of the resulting azimuthon can be constructed from the corresponding linear eigen-solution. This offers a straightforward and exhaustive method to identify rotating soliton solutions in a given nonlinear medium.