Rotating soliton solutions in nonlocal nonlinear media

We discuss generic properties of rotating nonlinear wave solutions, the so called azimuthons, in nonlocal media. Variational methods allow us to derive approximative values for the rotating frequency, which is shown to depend crucially on the nonlocal response function. Further on, we link families of azimuthons to internal modes of classical non-rotating stationary solutions, namely vortex and multipole solitons. This offers an exhaustive method to identify azimuthons in a given nonlocal medium.


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 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 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]. It has been shown recently that stabilization of localized waves is greatly enhanced in the nonlocal nonlinear media. In such media 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 [6], ballistic transport of hot atoms or molecules [7] or heat transfer [8,9]. It can be also a signature of a long-range interparticle interaction such as in nematic liquid crystals [10,11,12]. Spatially nonlocal response is also naturally present in atomic condensates where it describes a noncontact bosonic interaction [13,14,15]. In this work we investigate in details rotating nonlocal solitons. We consider various nonlocal 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.

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 = x e x + y e y denotes the transverse coordinates. If θ is governed by the following diffusion-type equation solving formally Eq. (3) in the Fourier space yields where K 0 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) [26]. They represent spatially rotating structures and hence involve an additional parameter, the angular frequency Ω 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 [25]. 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 [31]) 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 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λ 0 z)], we can assume any value for Ω by just shifting the propagation constant λ = λ 0 + Ω accordingly (λ 0 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 .

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 [28, 29] 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 0 and fix the amplitude factor A in ansatz (9) to A = 2/(α 2 + 1). With this choice the mass of the azimuthon equals M 0 for all values of α. Then, we can compute all integrals in Eqs. (7) as functions of α: 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 λ 0 . 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 [32]. We employ the following ansatz V = Ar exp(−r 2 /2σ 2 ), U = V exp(iφ ), ψ = U exp(iλ 0 z) and look for critical points of the Lagrangian L = −λ 0 M + I + N/2 with respect to the variables A and σ . In the following, we give results for our two exemplary nonlocal nonlinearities.
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 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 λ 0 and drastic difference in magnitude for Gaussian and Bessel responses.

its amplitude
and the rotation frequency Ω Ω = 2α α 2 + 1 We show these three quantities as functions of λ 0 in Fig. 1 (dashed lines). If λ 0 is large, σ ≪ 1 (nonlocal limit) and we can neglect the first two terms in Eq. (14). Then we get σ nl = (2/λ 0 ) 1/4 , A nl = λ 0 (σ 2 nl + 1) 3/2 and, subsequently, Ω nl = α/(α 2 + 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 λ 0 [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.
To shed light on the origin of the difference in the rotating frequencies of the two models let The degree of spatial averaging is obviously larger for the Gaussian than for the Lorentzian.
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 .

Rotation Frequency -Numerical Analysis
How good are our approximations of the rotating frequency Ω? Is the dependency in the structural parameter α really approximatively α/(α 2 + 1) for a given λ 0 ? 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 λ 0 = 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 λ 0 = 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 λ 0 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.
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.

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   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 α = Ω max d/ max D. In effect, the whole rotating structure is slightly elliptic, with the principal axis given by D. In our example, we find max d/ max D = 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).
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 δV = δV 1 (r)e imφ +iκz + δV * 2 (r)e −imφ −iκ * z we then derive the Eigenvalue problem for the internal modeŝ Real Eigenvalues of Eq. (20) (κ = κ * ) are termed orbitally stable and the corresponding Eigenvector (δV 1 , δV 2 ) 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 λ = λ 0 − κ/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 Please not that since | r − r ′ | = r 2 + r ′2 − 2rr ′ cos(φ − φ ′ ), all integrals in (21) are independent of φ . As for the perturbation of the stationary dipole, it is possible to construct azimuthons in the vicinity of the vortex (α ≈ 1) from δV : 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 2 , 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 2 and "forget" the other component ∼ δV 1 with total vorticity 3, which leads to the observed deviation. The larger the "forgotten" component, the larger is this deviation (compare Figs. 7, 3 and Figs. 8,4).
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. 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 |κ| > |λ 0 | in the continuous part of the spectrum. Hence, we can conclude that |Ω| α=1 | < |λ 0 |/m for azimuthons bifurcating from a singlecharged vortex 3 . 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.

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.

Acknowledgment
This research was supported by the Australian Research Council. Numerical simulations were performed on the SGI Altix 3700 Bx2 cluster of the Australian Partnership for Advanced Computing (APAC).