Existence, Stability, and Scattering of Bright Vortices in the Cubic-Quintic Nonlinear Schrödinger Equation

# Existence, Stability, and Scattering of Bright Vortices in the Cubic-Quintic Nonlinear Schrödinger Equation

## Abstract

We revisit the topic of the existence and azimuthal modulational stability of solitary vortices (alias vortex solitons) in the two-dimensional (2D) cubic-quintic nonlinear Schrödinger equation. We develop a semi-analytical approach, assuming that the vortex soliton is relatively narrow, and thus splitting the full 2D equation into radial and azimuthal 1D equations. A variational approach is used to predict the radial shape of the vortex soliton, using the radial equation, yielding results very close to those obtained from numerical solutions. Previously known existence bounds for the solitary vortices are recovered by means of this approach. The 1D azimuthal equation of motion is used to analyze the modulational instability of the vortex solitons. The semi-analytical predictions – in particular, that for the critical intrinsic frequency of the vortex soliton at the instability border – are compared to systematic 2D simulations. We also compare our findings to those reported in earlier works, which featured some discrepancies. We then perform a detailed computational study of collisions between stable vortices with different topological charges. Borders between elastic and destructive collisions are identified.

## I Introduction

The cubic-quintic nonlinear Schrödinger (CQNLS) equation is used to model a variety of physical settings. In scaled units, the CQNLS equation takes the following form:

 i∂Ψ∂t+∇2Ψ+|Ψ|2Ψ−|Ψ|4Ψ=0, (1)

where is the complex wave function, is the two-dimensional (2D) Laplacian, and the last two terms represent, respectively, the focusing cubic and defocusing quintic nonlinearities. The CQNLS equation emerges in models of light propagation in diverse optical media, such as non-Kerr crystals (20), chalcogenide glasses (34); (35), organic materials (38), colloids (2); (4), dye solutions (15), and ferroelectrics (17). It has also been predicted that this complex nonlinearity can be synthesized by means of a cascading mechanism (12). It should be noticed that, in the optics models, evolution variable is not time, but rather the propagation distance.

The competition of the focusing (cubic) and defocusing (quintic) nonlinear terms is the key feature of the CQNLS model, which allows for the existence of stable multidimensional structures which would be unstable in the focusing cubic nonlinear Schrödinger (NLS) equation. In particular, the CQNLS equation supports solitary vortices (alias vortex solitons) in 2D and 3D geometries. These are ring-shaped structures which carry a rotational angular phase. The respective integer winding number in the vortex is referred to as its topological charge (alias vorticity), . An example of a stable solitary vortex is depicted in Fig. 1.

It is well known that the cubic NLS equation supports “dark” (delocalized) and “bright” (localized) vortices, with the self-defocusing and focusing signs of the nonlinearity, respectively. Dark vortices of unitary charge are stable, while higher-order ones are unstable, splitting into unitary eddies (13). On the other hand, the “bright” modes (vortex solitons) are always azimuthally unstable in the framework of the cubic equation or its counterpart with the saturable nonlinearity (14), which leads to the breaking of the vortex soliton into a number of fragments depending on the charge (8); (24). The addition of a periodic potential to the cubic NLS equation may stabilize vortex solitons of a different type, which are built (in the simplest case) as a chain of four peaks with the phase circulation corresponding to integer vorticity (30); (5); (37); (18); (25). An example of the azimuthal breakup of an unstable vortex is displayed in Fig. 2.

In contrast to the cubic NLS, the CQNLS can support stable solitary vortices. For the first time, this remarkable fact was discovered “empirically” in simulations of the 2D CQNLS equation (31), and later was investigated in detail in a more rigorous form in Refs. (36); (23). Moreover, three-dimensional solitons with embedded vorticity also have their stability region in the 3D CQNLS equation(26) (see also reviews (10); (24)). Stable vortex solitons may find their potential applications in the design of all-optical data-processing schemes. In that respect, knowing the growth rates of unstable modes is also important, because, if the rates are small enough, the vortices may be considered as practically stable ones, as they will not exhibit an observable instability over relevant propagation distances.

One of objectives of this paper is to revisit the topic of the azimuthal modulational stability of solitary-vortex solutions in the 2D CQNLS equation. As mentioned above, numerous studies of this problem have already been published (31); (36); (10); (23); (27). The aim of the analysis was to find a critical value () of the intrinsic frequency of the vortex [for its exact definition see Eqs. (2) and (3) below], above which the vortices with a given value of the topological charge are stable. In Ref. (31), 2D azimuthally stable vortices with were shown to exist. It was found that the slope of the vortex’ profile at the pivotal point peaked at a specific value of , which was considered as the critical frequency, , while the vortices, with all values of , exist at (at , the radius of the vortex diverges, i.e., the “bright” vortex goes over into a “dark” one). The same value, is simultaneously the largest one up to which exact soliton solutions exist in the 1D version of the CQNLS equation (28). In work (31), a variational approach (VA) was developed for 2D vortex solitons, yielding results similar to those obtained in the numerical form. Full 2D simulations of the CQNLS equation reported in Ref. (31) had confirmed that vortices with were indeed stable, while those with were not.

Works (36); (10); (23) employed a more rigorous approach. They introduced small perturbations around the 2D vortex soliton and solved the resulting eigenvalue problem numerically. In this way, the critical frequencies were found as and . Also, in Ref. (29), using the Gagliardo-Nirenberg and Hölder inequalities together with Pohozaev identities, is was shown that the eigenvalues generated by the CQNLS equation possess an upper cutoff value.

In Ref. (27), 2D perturbations were considered too. Through an extensive analysis, the problem was transformed into finding, by means of numerical methods, zeros of respective Evans functions for a set of ordinary differential equations (ODEs). In this way, the existence of azimuthally stable vortices of all integer values of was predicted (in Refs. (36) and (31), the stability regions for were not identified, as they are extremely small). The predictions for the critical frequencies reported in Ref. (27), which are somewhat different from those in Refs. (36) and (31), are shown in Table 2.

The present work aims to undertake an additional study of the stability of vortices in the CQNLS equation for two reasons. Firstly, the previous studies demonstrated some (relatively small, but not negligible) disagreement in the predictions of the critical frequency. Therefore, since we will use a different approach, our results may help in comparing and verifying the previous findings. Secondly, although in Ref. (31) 2D simulations were performed for the vortices, this was only done for , and longer simulations later showed that some vortices which were originally concluded to be stable turned out to be eventually unstable (36). Therefore, stability results from additional 2D numerical simulations are needed for comparison with various predictions of the stability. In fact, our numerical results will illustrate that the predictions of Ref. (27) are the most accurate ones. In parallel to extensive direct simulations, we elaborate a semi-analytical methodology, following the lines of previous studies of the azimuthal modulational instability of vortices in the cubic NLS equation (8). The approximation amounts to splitting the full 2D equation into effectively one-dimensional radial and azimuthal equations. The former equation is used to compute the radial shape of the soliton by means of a VA to find an analytic approximation to the profile, which is then used as an initial condition to a numerical optimization routine to find the numerically ‘exact’ profile. The closeness between the numerically computed profiles and the VA approximation testify to the good accuracy of the VA. The analysis of the azimuthal equation makes it possible to predict the threshold of the onset of the splitting instability of the vortex soliton. Our results predicting the growth rates of individual azimuthal modes of unstable vortices match full simulations very well, but the predictions for the critical frequency are less accurate, when compared to the simulations. This semi-analytical technique may be quite relevant for applications in other 2D models.

Another objective of the work is to study, by means of systematic simulations, collisions between stable vortices. Except for few examples reported in Ref. (31), this problem was not studied before.

The paper is organized as follows. In Sec. 2 we derive an approximate analytical description of the vortex profile by first finding an analytical asymptotic approximation to it, which is then employed as the ansatz on which the VA is based. In Sec. 3 we use the variational ansatz as the initial guess, to generate numerically exact vortex profiles by dint of a nonlinear optimization routine. We then compare the numerical profiles with the variational ansatz to show its accuracy. In Sec. 4 we derive an approximate 1D equation for the dynamics along the azimuthal direction. A linear stability analysis is subsequently performed to find stability criteria and growth rates of unstable modes within the framework of the azimuthal equation. In Sec. 5 we present full 2D simulations of the vortices, and the respective results for their stability. These results are compared to our predictions, as well as to the predictions produced by previous works. In Sec. 6, we use direct simulations to explore collisions between stable vortices in detail. In particular, these studies allow us to identify a border between quasi-elastic and destructive collisions. Finally, in Sec. 7, we summarize our findings and formulate concluding remarks.

## Ii Approximate Analytical Profiles of Steady-State Vortices

A steady-state vortex solution to Eq. (1) is looked for as

 Ψ(r,θ,t)=f(r)A(θ,t), (2)

where real function is a stationary radial profile with azimuthal dependence given by

 A(θ,t)=ei(mθ+Ωt), (3)

with topological charge and frequency . Analytical solutions being not available for the profile , we begin the study by finding an approximate analytic expression for and identifying its existence bounds. As shown below in Sec. III, the predicted profile is very close to the true solution, therefore it can be used to predict the existence and stability regions of the vortex solitons. The developed analytical method is quite general and may be applied to other physically relevant equations.

### ii.1 The Asymptotic Vortex Profile

Inserting expression (2) into the underlying equation (1) yields the following ODE for the radial profile of the vortex with charge :

 1rddr(rdfdr)−(Ω+m2r2)f(r)+f3(r)−f5(r)=0. (4)

If we assume that the vortex takes the form of a relatively narrow ring with a large radius, then, in the region of interest, variable may be approximately replaced by a constant, , which we take to be the value of at a point where the radial profile attains its maximum. Since we assume to be large, in the lowest approximation we neglect the term in the Laplacian. An obvious consequence of dropping this term is that the approximate profiles which we are going to derive will be less accurate for smaller values of , see Sec. III.

Under the large-radius assumption, Eq. (4) becomes

 d2fdr2−Ω∗f(r)+f3(r)−f5(r)=0, (5)
 Ω∗≡Ω+m2/r2c. (6)

For a relatively narrow ring, one can treat Eq. (5) as if , in which case the equation has a well-known 1D analytical solution for the soliton (28); (9):

 f2\scriptsize asy(r)=4Ω∗1+√1−(16/3)Ω∗cosh(2√Ω∗(r−rc)). (7)

As mentioned above, this solution exists for . At point , the solution degenerates into a constant, . The family of solutions (7) parameterized by is depicted in Fig. 3.

To obtain a closed-form analytical approximation, it is now necessary to find an expression for the radius of the profile, , for given frequency . To this end, we will use the VA, and will then compare the results with numerically found profiles.

### ii.2 The Variational Approach

To apply the VA, we use the Lagrangian density corresponding to the CQNLS equation:

 L=i2(ΨΨ∗t−Ψ∗Ψt)+|Ψr|2+r−2|Ψθ|2−12|Ψ|4+13|Ψ|6, (8)

where the subscripts denote partial derivatives. Inserting here expression (2), to be used as a factorized ansatz, yields

 L=(Ω+m2r2)f2(r)+(dfdr)2−12f4(r)+13f6(r). (9)

The integration of density (9) gives rise to the full Lagrangian:

 L≡2π∫∞0L(r)dr=2π(ΩC1+C2+m2C3−12C4+13C5), (10)

where

 C1 ≡∫∞0f2(r)rdr, C2 ≡∫∞0(dfdr)2rdr, C3 ≡∫∞0f2(r)drr, C4 ≡∫∞0f4(r)rdr, (11) C5 ≡∫∞0f6(r)rdr.

Next, we use the asymptotic approximate profile (7) as an ansatz for radial profile, treating and in Eq. (7) as variational parameters. Inserting expression (7) into the integral terms in the Lagrangian, and again making use of the large-radius approximation yields

 C1 ≈2√3Trc,C2≈rc[38√Ω∗−√3T(316−Ω∗)], C3 ≈2√3Trc,C4≈−3rc[√Ω∗−√32T], (12) C5 ≈−rc[278√Ω∗−3√3T(916−Ω∗)],

where

 T≡arctanh[√316Ω∗−√316Ω∗−1]. (13)

The respective static Euler-Lagrange equations, and , yield, respectively, equation

 rc=m(Ω−316+√3Ω∗8T)−1/2, (14)

and the other one which is equivalent to Eq. (6). Eliminating from these equations, one concludes that is, for a chosen value of , a solution to the following equation,

 Ω=Ω∗2+332−√3Ω∗16T≡G(Ω∗). (15)

This yields the final form of the vortex radial profile, as predicted by the VA:

 f2\scriptsize va(r)=4Ω∗1+√1−(16/3)Ω∗cosh(2√Ω∗(r−r\scriptsize vac)), (16)
 r\scriptsize vac≡m√Ω∗−G(Ω∗). (17)

Finally, the transcendental equation (15), together with Eq. (13), was solved numerically by means of a simple bisection method with a tolerance of .

In Fig. 4 we show versus and for , along with relationship (15). We see that, with the increase of , the radius of the vortex starts from very large values (infinite at ), decreases until it reaches a minimum value, and then increases rapidly (it becomes infinite once again at ). We also see that the relationship between and does not depend on charge , and that as (the apparent gap near the right edge of the plot is due to the sensitivity of the relationship, which we discuss below).

### ii.3 Existence Bounds for the Two-Dimensional Vortex Profiles

From numerical investigations performed in Refs. (31) and (11), the existence border for the two-dimensional CQNLS equation vortex solution was found to be , while, as said above, the analytical limit, which is identical for 1D and 2D equations, is The discrepancy may be explained by the sensitivity of relationship (15) between and , as produced by the variational approximation. From Eqs. (15) and (13) we see that, as , , and so , in accordance with the exact result. However, as seen in the right panel of Fig. 4, the relationship between and seems to be extremely sensitive near this limit (see the inset in the right panel of Fig. 4). This sensitivity can also be observed in Table 1, where some values of and the corresponding are given. It is clearly seen that one quickly approaches the limit of machine double precision for as .

This effect can also be understood from the logarithmic divergence of close to the limit point, see Eq. (13). Therefore, it is not surprising then that numerical estimates of , obtained by means of a shooting method to look for profiles at different vales of , gave lower estimates of the existence bound. In Sec. III, we confirm, using accurate numerical methods to find vortex profiles for , that the actual existence bound is greater than the numerical estimates of Refs. (31) and (11).

The analytical approximation (16) of the vortex profile is used in the following sections as an input to a numerical optimization routine which solves the ODE (4) to find numerically exact radial profiles.

## Iii Numerically Exact Steady-State Vortex Profiles and Comparison to the Variational approximation

To find numerical solutions to Eq. (4), we used a modified Gauss-Newton optimization routine, with a tolerance (8). In Fig. 5, some numerically found radial profiles are displayed for charges , , and and , including as close to as . To our knowledge, the profiles above have not been shown in any previous study. It is seen that the profiles flatten out as grows, pushing the profile farther from . At the same time, the ring of the vortex becomes wider without a change in its height. It is worth mentioning that, due to the high sensitivity of the numerics mentioned in the previous sections, computing profiles past is a daunting task requiring high-precision arithmetic. Implementing high-precision arithmetics in a Gauss-Newton optimization routine (or any fixed-point iteration method) would be quite involved and very time consuming. Nonetheless, we have found that, using high-precision arithmetics (up to decimal places for frequencies as close to as ) in the calculation of in the variational equation (15) yields approximate analytical solutions that the Gauss-Newton subroutine is able to process into numerically accurate solutions for values of very close to , as shown in Fig. 5. The results clearly demonstrate the usefulness of the VA: without using this approximation for generating the initial profiles to be fed into the numerical solver, a prohibitively complex high-precision fixed-point algorithm would be needed to obtain meaningful profiles past .

The accuracy of the variational ansatz per se can be tested by comparing it to the numerical solutions computed by means of the Gauss-Newton routine. In Fig. 6 we compare the VA prediction for the vortex-ring’s radius, , to the numerically found values for . We observe that the radii provided by both methods are virtually identical, allowing one to use the VA-predicted radius in applications, that may be useful in experimental situations.

To get an even better idea of how close the variational approach is to the Gauss-Newton profile overall, in Fig. 7 we compare the relative sum of squared deviations of the variational profiles from their the numerically found counterparts, for different values of and different values of . We also plot a series of profiles produced by the VA and by the Gauss-Newton method for . We observe that, as increases, the mismatch between the variational and numerical profiles decreases. This is understandable, as we have neglected the term in the Laplacian, while formulating the variational ansatz. Therefore, since vortices with smaller have smaller inner radii, the discrepancy between the variational ansatz and the numerically exact solution are expected for lower values of . The total discrepancy is observed to be quite low overall, showing that the variational profile provides for a very accurate rendition of the true solution, especially for larger values of . Since our VA provides an accurate radial profile of the vortex solitons, we can use it to derive fully analytic azimuthal modulational stability predictions, which is done below.

## Iv The Azimuthal Modulational Stability: Analytical Results

With the profiles of the steady-state vortices available, we proceed to the study of their azimuthal modulational stability. To this end, we apply the methodology utilized in previous works, dealing with circular gap solitons in the model of the circular Bragg grating (7), and in the cubic NLS equation (8): we will derive an azimuthal equation of motion by assuming a frozen radial profile, and then perform a perturbation analysis to analyze the stability of azimuthal perturbation modes.

In Ref. (8), numerically computed eigenmodes of perturbations around the solitary vortices featured a slight coupling in the radial and azimuthal directions, hence assuming the compete separability of the radial and azimuthal directions for the evolution of the perturbations may lead to discrepancies between the predictions for the stability and numerical results. Nonetheless, the coupling between the radial and azimuthal directions is attenuated with the increase of , hence, the predictions made under the assumption of the separability should be more accurate for higher values of . The trend to the improvement in the accuracy of the analytical approximation with the increase of was observed in the numerical results presented in Sec. III.

### iv.1 The Azimuthal Equation of Motion

To derive the azimuthal equation of motion, we first insert the separable ansatz (2) into the Lagrangian density (8) to obtain

 L = i2f2(r)(AA∗t−A∗At)+(dfdr)2|A|2+ (18) 1r2f2(r)|Aθ|2−12f4(r)|A|4+13f6(r)|A|6.

Using our steady-state radial profiles, we can perform the integration of the Lagrangian density over , thus arriving at an effectively one-dimensional (in the direction of ) Lagrangian that may be used to derive the equation of motion for :

 L\scriptsize 1D=∫2π0L\scriptsize 1Ddθ, (19)
 L\scriptsize 1D ≡ i2C1(AA∗t−A∗At)+C2|A|2 (20) +C3|Aθ|2−s12C4|A|4−s23C5|A|6,

and , , are the radial integrals defined as per Eq. (11). Evaluating the variational derivative of the action functional (22), and applying a linear transformation,

 A→Aexp(−iC2t/C1), dt→C3t/C1, (21)

yields the following evolution equation for :

 iAt=−Aθθ−(C4/C3)|A|2A+(C5/C3)|A|4A. (22)

This equation provides a description of how the radially-frozen, azimuthally time-dependent solution will evolve in the CQNLS model.

### iv.2 The Stability Analysis

To study the azimuthal modulational stability of vortex-soliton solutions to the CQNLS equation, we performed a perturbation analysis in the framework of the azimuthal equation of motion (22), with the objective to compute the growth rates of small perturbations. We begin with an azimuthal “plane-wave” solution perturbed by a complex time-dependent perturbation:

 A(θ,t)=[1+u(θ,t)+iv(θ,t)]ei(mθ+Ω′t), (23)

with . Inserting this into Eq. (22), performing the linearization with respect to the perturbations, and separating the result into real and imaginary parts, we obtain a system of coupled equations (a time shift was made here, ):

 ut=−vθθ,vt=uθθ+(2C4−4C5C3)u. (24)

In order to study the azimuthal modulational stability of the solitary vortices, we expand and into Fourier series over azimuthal harmonics with integer wavenumbers :

 ^u(K,t) = ∫2π0u(θ,t)eiKθdθ, ^v(K,t) = ∫2π0v(θ,t)eiKθdθ. (25)

Applying these transforms to Eq. (24) yields two coupled equations for the amplitudes of and of each mode. In a matrix form, the equations are

 ddt[^u^v]=⎡⎢⎣0K2(2C4−4C5C3−K2)0⎤⎥⎦[^u^v]. (26)

The growth rates for each azimuthal wavenumber are simply eigenvalues of Eq. (26). Taking into account the underlying rescaling (21), they are:

 λ1,2=±C3C1√K2(K2\scriptsize crit% −K2), (27)

where is the critical value of , above which all the modes are stable:

 K\scriptsize crit≡√2(C4−2C5)C3. (28)

For the modulationally unstable vortices, it is useful to know the maximum growth rate and the mode that exhibits it. This is because, in an experiment, even an unstable vortex may be practically “stable enough” if the maximum perturbation growth rate is small enough. This information can be extracted from expression (27) by equating its derivative to zero and solving for , which reveals the fastest growing perturbation mode (and subsequently, the prediction of the number of fragments that the unstable vortices will break up into):

 K\scriptsize max=√C4−2C5C3=1√2K\scriptsize crit. (29)

We can then insert this value into Eq. (28), which yields the maximum growth rate,

 λ\scriptsize max=C4−2C5C1=C32C1K2\scriptsize crit. (30)

In fact, because must be integer, the actual fastest growing eigenmode may correspond to the integer closest to value (29). Accordingly, the actual largest growth rate may be somewhat smaller than the one given by Eq. (30).

For the vortex to be azimuthally stable against all modes, one needs either (then, there is no integer value ) or being purely imaginary. According to Eq. (28), these two stability criteria amount to the following inequalities:

 imaginary  K\scriptsize crit:  C4−2C5<0; (31)
 ∣∣K\scriptsize crit∣∣<1:  C4−2C5−C3/2<0. (32)

Since is positive, the former condition is stricter than the latter one. Although the latter condition is sufficient for the azimuthal stability, we also keep the former one (according to the above derivation, it is relevant to the infinite line) because it leads to an estimation for the critical frequency which is independent of charge (see below).

For our predictions, we compute the constants using numerical integrals of the numerically exact radial profiles. However, since the variational analytic profiles are so close to the numerically exact profiles, one can use expressions (12) to find approximate analytic values of Eqs. (31) and (32), which then yield critical values of , above which all the vortex solitons are azimuthally stable for charge . These values can then be compared to those computed by the numerically exact profiles.

### iv.3 Analytical Stability Predictions Using the Variational Profile

As we showed in Sec. III, the variational ansatz yields a useful approximation to the true radial profiles. Therefore, the VA can be employed to calculate the constants in Eq. (12), and thereby derive analytical predictions for the azimuthal modulational stability. In Sec. IV.2 it was demonstrated that, for studying the stability of the vortices, only the values of and are required,

 K\scriptsize va\scriptsize crit=m ⎷6Ω∗−(15/8)+5√3Ω∗/(4T)Ω∗−G(Ω∗), (33)
 (C3C1)\scriptsize va=Ω∗−G(Ω∗)m2, (34)

since the other coefficients can be expressed in terms of these two (here, as before, ).

The VA can predict the critical value of above which the vortices are modulationally stable. Inserting the constants defined in Eq. (12) into the stability criteria (31) and (32), we arrive at the following expressions which determine the critical frequency:

 K\scriptsize crit=0:  6Ω∗−158+5√3Ω∗4T=0; (35)
 K\scriptsize crit=1:   ⎷6Ω∗−15/8+5√3Ω∗/(4T)Ω∗−G(Ω∗)−1m=0, (36)

where is defined as per Eq. (13). We see that, at , these two criteria tend to coincide. At finite , Eqs. (35) and (36) give different critical values of , and hence of too —one which depends on charge , and the other one, independent of , which represents an upper bound on for all values of . We denote the charge-dependent critical value as , and the charge-independent one as .

Solving Eq. (35) for by means of a root finder and inserting the result into Eq. (15) yields

 Ω\scriptsize va\scriptsize st=0.144320424. (37)

Since this value is smaller than , this result predicts that azimuthally stable vortices exist for all values of , at . We also solved Eq. (36) for at various values of . The results are displayed in Fig. 8. It is seen that the critical frequency increases with the increase of and eventually converges to . Thus, the stability window in is larger for lower charges. In Ref. (27), it was also concluded that azimuthally stable vortices exist for all , and that lower charges indeed have a larger stability window. However, according to Ref. (27), , i.e., the stability window shrinks towards as . The estimate obtained in that work shows that the window shrinks as . This conclusion precludes experimental creation of higher-order stable vortices, as the respective stability interval would be too small, and the radius of the vortex too large.

According to the VA predictions, does increase with (and, as shown in Fig. 8, the difference decays proportionally to which resembles the prediction of Ref. (27)). On the other hand, since the variationally predicted value of is smaller than , there remains a stability window at all , which does not vanish at .

As shown in Table 2, the variational predictions are very close to those using the numerically-exact profiles. We then note that the numerical results reported in Ref. (27) are more accurate than our predictions, the stability window indeed shrinking to zero at high values of .

## V Azimuthal Modulational Stability: Numerical Results

In this section we report the numerical predictions and full simulation results for the azimuthal modulational (in)stability of the vortices. For the 2D simulations, we used a finite-difference scheme with a central difference in space and fourth-order Runge-Kutta in time (16). We used both polar-coordinate and Cartesian grids. The polar grid makes computing the growth of individual perturbation modes easier, therefore we used it to test our predictions for azimuthally unstable vortices. However, the polar grid forces one to use smaller time steps in the finite-difference scheme than is needed for the equivalent Cartesian grid. Therefore, for testing the critical frequency , which requires running multiple long-time simulations, we use the Cartesian grid. We used second-order differencing for the polar-coordinate simulations, and (due to the long simulation time needed) a fourth-order differencing for the Cartesian simulations.

### v.1 Unstable Vortices

In Fig. 9 we show the results for the (in)stability of vortices with charges and . In general, we see a good agreement between the numerically measured growth rates and the predictions for ; however, for values of which get closer to , the accuracy of the predictions becomes low in each case, implying that our predictions for are not precise. To further test this, we ran long-time simulations of randomly perturbed vortices.

For , we were easily able to identify the transition from unstable to stable vortices, but for we found it very difficult to detect a stable solution, because of a snake-like instability which breaks the vortex into asymmetric irregular fragments. An example of this is shown in Fig. 10. Obviously, this snake-like instability is not captured by our study of the azimuthal modulational instability. It is very plausible that it accounts for the discrepancy between the analytic and numerical growth rates observed in Fig. 9.

### v.2 Stable Vortices

To identify the values of corresponding to the stability border, we simulated the evolution of vortices with taken near the predicted values of , perturbed by a small uniformly distributed random noise in the azimuthal direction. Our aim was to find a value of that results in an unstable vortex soliton (and to observe the actual breakup), and then show that the vortex solution for is stable, by simulating its evolution for an extremely long time, compared to the time necessary to reveal the full breakup in the former case. We did this for various charges . In Fig. 11, we show the initial and final states produced by this analysis for .

Since the snaking effect hinders our ability to simulate the evolution of the vortices for extended time periods at , we were not able to make precise stability predictions in this case. However, as this effect is dynamically distinct from the azimuthal breakup, we can still give an approximate estimate of the critical frequency for higher charges.

In Table 2 we display the predictions of and their numerically found counterparts for . The column labeled “semi-numerical” are the predictions made using the numerically exact radial profiles, while the column labeled “VA” gives the predictions made using the variationally derived profiles. For comparison, in the same Table we also display findings reported in earlier works. We see that most of the predictions for are close to the numerical result. For , the results reported in Ref. (27) are most accurate ones.

It is observed that the semi-numerical and variational results are very close, the latter ones being, in fact, slightly more accurate in comparison with the numerical findings. This implies that the discrepancies in the variationally derived predictions are not due to any inherent deficiency in the VA. Rather, the discrepancies of both sets of predictions are likely due to the approximation of the vortex soliton as a separable entity composed of radial and azimuthal parts (as discussed in Sec. IV).

## Vi Collisions and Scattering of Stable Vortices

In this section we study collisions and scattering of stable vortices in the CQNLS. In Ref. (31), such collisions of vortices of unit charge were considered in a brief form. It was shown that both elastic and destructive collisions could be observed depending on the phase difference, , between the colliding vortices, as well as on their relative velocity. Here we expand this study in three directions: i) we include vortices of charge ; ii) we consider the critical velocity for elastic collisions as a function of ; iii) we study the scattering of vortices colliding at various values of impact parameters, by measuring the scattering angles and speeds.

Summarizing results of numerous simulations, we have found that, at sufficiently small collisional velocities2, the elastic or destructive character of the collision is solely determined by the phase difference at the point of contact. This fact seems to be independent of the charge of the vortices involved in the collision. In Fig. 12 we display four different cases that illustrate the main features of head-on collisions between vortices with different charges. In the first two columns of Fig. 12 we show the collision of vortices carrying charges and , with and , respectively. It can be seen that, since the vortex with has opposite phase on the collision side, with respect to the vortex with , the interaction is repulsive since it locally emulates the interaction of two out-of-phase fundamental bright solitons (21); (1). Therefore, if the velocity is small enough, the mutual repulsion determines the result of the collision. On the other hand, when the vortex with charge is phase-shifted by (or similarly, if the charge vortex were to be placed on the opposite side of the charge vortex), the adjacent sides of the colliding vortices are in phase. Therefore, the interaction is similar to that between two in-phase bright solitons, which is attractive (21); (1). This results in the merger of the two vortices and their eventual breakup into fragments.

Clearly, even when the repulsive force acts between vortices, but the collision velocity is sufficiently large, the vortices have enough momentum for causing them to merge and eventually break up as well. This case is depicted in the third column of Fig. 12.

Slowly moving vortices with opposite phases at the point of contact are quite robust against the collisions, as illustrated in the fourth column of Fig. 12. In this panel we depict a “billiard”-type example, in which one vortex collides with another, which in turn collides with a third one. All collisions in this case are elastic because the relative phases have been chosen such that the colliding sides are always out-of-phase, providing for the necessary repulsion. It is worthy to note that, although these collisions at low velocities seem to be elastic, the stable vortices clearly have internal breathing modes, that are excited during the collisions. Therefore, a small fraction of the collisional energy is transferred to these internal modes, preventing the collisions from being completely elastic. A more in-depth analysis of the excitation of these internal modes and of the degree of the inelasticity of the collisions fall outside of the scope of the present work. Nonetheless, in order to quantify the degree of the elasticity of the collisions, we have measured the critical velocity for the vortices to elastically bounce off each other: at the vortices bounce elastically, while at they merge and, most often, get destroyed. The crucial parameter that controls , for a particular combination of the charges of colliding vortices, is their phase difference . Different values of  allow to have different relative phases at the colliding sides (as discussed above), and thus (partially) repel or attract. This effect is illustrated by Fig. 13, where we depict versus for vortex pairs with charges (a) , (b) , and (c) . From panel (a) it is clear that, for charges , the most “stable” (in terms of observing elastic collisions at larger collisional velocities) collision occurs around . This corresponds exactly to the most stable case mentioned above, when the collisional sides are out of phase. In panel (b) we depict for . It is seen in this panel that is -periodic, as a function of , due to the angular symmetry of the vortex. For this vortex-charge combination, the most stable collision occurs around (or, due to the symmetry, ), which again corresponds to the adjacent out-of-phase sides. Finally, in panel (c) we present the results for two vortices with . The behavior is very similar to that in panel (b): the same periodicity, shape and value of the optimal phase difference. However, since the vortices with vortices are less stable than those with , the critical velocity is lower for the charge combination than for the case of .

In order to fully characterize the collisions between vortices with different charges, we studied the scattering of a vortex of charge impinging at velocity upon another vortex of charge at rest (), with impact parameter . The impact parameter is defined as the perpendicular distance between the trajectory of the incoming vortex and the initially quiescent one. In Fig. 14 we show several orbits for the two vortices with different impact parameters, vortex-charge combinations, and incoming velocities. It is seen, in all the panels, that all the trajectories are well defined for the chosen parameters, and clearly feature elastic scattering. The key to achieve elastic scattering is to choose, for each charge combination, the optimal phase difference, using Fig. 13 [i.e.,  for and for , , or ], and the initial velocity lower than the corresponding . For example, we chose in panel (a) of Fig. 14 and . In fact, when one uses a velocity closer to the critical one the vortices start experiencing strong deformations of the quadrupole and octupole types when they collide nearly head-on. This effect is clearly visible for the trajectory with (also and ) in panel (b) of Fig. 14, where is close to the critical value, , and the trajectory features undulations due to the internal deformation of the vortex. Even larger velocities, i.e., above , lead to elastic scattering for larger values of the impact parameter , and to annihilation of the colliding vortices at smaller (not shown here). It is interesting to note that, as for hard spheres with different masses, the scattering of vortices can also produce very different scattering angles, depending on the mass combination of the two colliding objects. This effect is made evident by the comparison of the scattering between vortices with and in panel (c), and the scattering between and in panel (d) (all the other parameters coincide in both panels). When a “lighter” vortex with collides with the “heavier” one, with , for , the former vortex bounces back after the collisions. In the reverse situation, with the “heavier” vortex impinging upon the “lighter” one with , the “heavier” vortex only loses a portion of its momentum and continues to move in a straighter trajectory, while the “lighter” vortex is pushed away at a larger velocity.

To further characterize the scattering between vortices we analyzed the scattering angle as a function of the impact parameter. In our scattering experiments, we initially define horizontally moving and quiescent vortices. After the collision, the trajectory of the moving vortex is deviated by scattering angle (in our case negative since we consider the moving vortex placed below the quiescent one), and the initially quiescent vortex is scattered at angle (in our case, positive). In Fig. 15 we depict the scattering angles corresponding to the scattering experiments displayed in Fig. 14. We depict by (blue) circles the negative of the scattering angle, , of the incoming vortex, and by (green) squares the scattering angle . In all the cases, we also depict the difference between the scattering angles, , by (red) triangles. In the case of the elastic collision between two hard spheres of equal masses it is well known that . Figure 15 reveals that the collisions between vortices with equal charges [ in panels (a) and (b), and in panel (e)] produce very close to for all values of the impact parameter. In contrast, as it is the case for hard spheres with different masses, our simulations demonstrate that for [ in panel (c)], and that for [ in panel (d)].

Another measure of the scattering can be given in terms of the initial and final velocities. In Fig. 16, we present the final velocities corresponding to the scattering events shown in Fig. 14. The figure shows, by (blue) circles and (green) squares, the post-collision velocities, and , respectively, of the initially moving and quiescent vortices. We also depict, using (red) triangles, the sum of the final velocities, , and the thin horizontal black line represents the initial velocity of the moving vortex, . It can be concluded that for relatively small values of the impact parameter, , and for larger values of .

The above results suggest that the interactions between vortices may be either elastic or destructive, depending on the parameters. In order to simulate the interplay of these interactions in a multi-vortex gas, we display in Fig. 17 a time series depicting the time evolution of a random vortex gas. The first two rows depict the local intensity plots at the indicated moments of time, while the bottom two rows depict the corresponding field phases. We use a square domain of size with periodic boundary conditions, seeding 15 vortices of charge at random positions (such that initially no two vortices were closer than twice the size of a single vortex), and with random phases. As seen from the figure, due to the random phase differences and collision angles, some vortex pairs experience elastic collisions, while others tend to merge and destroy each other. Since the destruction of vortices is an irreversible process, all the vortices get eventually destroyed, leading to a seemingly disordered pattern of interacting intensity patches with an approximately constant phase inside each one (see the correspondence between field patches in the top rows and the phase patches in the bottom rows, at the late stage of the evolution). The dynamics resembles the grain coarsening in typical chaotic spatiotemporal systems, see e.g., review (3).

## Vii Conclusions

In this work, we have revisited the issues of the existence, stability, and interactions of solitary vortices in the two-dimensional CQNLS (cubic-quintic nonlinear Schrödinger) equation, using both numerical and analytical methods. The latter one is based on the variational approach, that has been shown to be increasingly more accurate as the vortex’ topological charge increases. We have also developed the analysis of the azimuthal modulational (in)stability of vortices, using the approximation originally proposed in other contexts in Refs. (7) and (8), which postulates the separation between the frozen radial profile of the vortex soliton and free evolution of azimuthal perturbations. This approach leads to an effectively one-dimensional equation for the azimuthal dynamics. Examining the stability of the perturbed annulus within the framework of the latter equation, we were able to predict the stability of the vortices against the breakup. We then ran full 2D simulations to verify the semi-analytical predictions.

For azimuthally unstable vortices, our predictions of the largest growth rate are fairly accurate over a wide range of the vortex’ intrinsic frequency (especially for the higher-order vortices, with topological charge ). The so developed semi-analytical technique may be useful for other applications. The analytical approach is less accurate in predicting the stability border (intrinsic critical frequency of the vortex) for . This discrepancy is due to deviation from the assumption admitting the separation of the perturbed evolution into radial and azimuthal parts, while the VA proper, which was used to predict the radial shape of the stationary vortex solitons turns out to yield accurate results. Comparing our numerical results for the critical frequency with previously published ones, we have concluded that the most accurate findings were reported in Ref. (27). We have also found that higher-order vortices may be subject to a snaking radial instability.

The semi-analytical approach developed in this work may be quite useful for other 2D models —at least, for the description of the most stable vortex solitons with lowest values of the topological charge.

We have also investigated, in detail, collisions between stable vortex solitons in the two-dimensional CQNLS equation. The outcome of the collision crucially depends on the phase difference and relative velocity of the vortices. We have produced the results showing the critical velocity for elastic collisions as a function of the relative phase for all combinations of the vortices with and . We also studied the dependence on the scattering angles and post-collision velocities on the impact parameter of the collision.

One relevant direction for the extension of the present analysis is to find accurate conditions for the stabilization of the nearly isotropic vortex solitons by an external spatially periodic (lattice) potential, see Ref. (33) and references therein. A challenging problem is to extend the semi-analytical considerations to three-dimensional solitons with embedded vorticity, in the model with the same cubic-quintic nonlinearity. In fact, it is not known if such 3D vortex solitons with higher vorticity () may be stable (26); (24).

## Viii Acknowledgments

The authors thank Juan Belmonte and Jesús Cuevas for providing useful insights. RCG gratefully acknowledges support from NSF-DMS-0806762. PGK gratefully acknowledges support from NSF-DMS-0349023 (CAREER) and NSF-DMS-0806762, as well as from the Alexander von Humboldt Foundation.

### Footnotes

1. URL: http://nlds.sdsu.edu
2. A horizontal velocity is imprinted to each vortex by multiplying its corresponding steady state by with . This factor of two between and can be formally obtained by transforming CQNLS solutions into a moving reference frame with velocity (6).

### References

1. S. K. Adhikari, Mean-field model of interaction between bright vortex solitons in Bose-Einstein condensates, New J. Phys. 5 (2003) 137.1–137.13.
2. G. S. Agarwal, S. Dutta Gupta, T-matrix approach to the nonlinear susceptibilities of heterogeneous media, Phys. Rev. A 38 (1988) 5678–5687.
3. I. S. Aranson, L. S. Tsimring, Patterns and collective behavior in granular media: Theoretical concepts, Rev. Mod. Phys. 78 (2006) 641–692.
4. E. L. Falcão-Filho, C. B. de Araújo, J. J. Rodrigues, Jr, High-order nonlinearities of aqueous colloids containing silver nanoparticles, J. Opt. Soc. Am. B 24 (2007) 2948–2956.
5. B. B. Baizakov, B.A. Malomed, M. Salerno, Multidimensional solitons in periodic potentials, Europhys. Lett. 63 (2003) 642–648.
6. Juan Belmonte, personal communication.
7. J. Scheuer and B. Malomed, Annular gap solitons in Kerr media with circular gratings, Phys. Rev. A 75 (2007) 063805.
8. R. M. Caplan, Q. E. Hoq, R. Carretero-González, P. G. Kevrekidis, Azimuthal modulational instability of vortices in the nonlinear Schrödinger equation, Opt. Commun. 282 (2009) 1399–1405.
9. S. Cowan, R. H. Enns, S. S. Rangnekar, S. S. Sanghera, Quasi-soliton and other behavior of the nonlinear cubic-quintic Schrödinger equation, Can. J. Phys. 64 (1986) 311–315.
10. L. Crasovan, B.A. Malomed, D. Mihalache, Spinning solitons in cubic-quintic nonlinear media, Pramana J. Phys. 57 (2001) 1041–1059.
11. A. Desyatnikov, A. Maimistov, B.A. Malomed, Three-dimensional spinning solitons in dispersive media with the cubic-quintic nonlinearity, Phys. Rev. E 61, 3 (2000) 3107–3113.
12. K. Dolgaleva, H. Shin, R. W. Boyd, Observation of a microscopic cascaded contribution to the fifth-order nonlinear susceptibility, Phys. Rev. Lett. 103 (2009) 113902.1–113902.4.
13. A. L. Fetter, Rev. Mod. Phys. 81, 647 (2009).
14. W. J. Firth, D. V. Skryabin, Optical solitons carrying orbital angular momentum, Phys. Rev. Lett. 79 (1997) 2450–2453.
15. R. A. Ganeev, M. Baba, M. Morita, A. I. Ryasnyansky, M. Suzuki, M. Turu, H. Kuroda, Fifth-order optical nonlinearity of pseudoisocyanine solution at 529 nm, J. Opt. A: Pure Appl. Opt. 6 (2004) 282–287.
16. G. H. Golub, J. M. Ortega, Scientific computing and differential equations an introduction to numerical methods, Academic Press, San Diego, California, 2 edition, 1992.
17. B. Gu, Y. Wang, W. Ji, J. Wang, Observation of a fifth-order optical nonlinearity in BiLaFeMgO ferroelectric thin films, Appl. Phys. Lett. 95 (2009) 041114.1–041114.3.
18. G. Herring, L. D. Carr, R. Carretero-González, P. G. Kevrekidis, D. J. Frantzeskakis. Radially symmetric nonlinear states of harmonically trapped Bose-Einstein condensates, Phys. Rev. A 77 (2008) 023625.1–023625.11.
19. Yu. S. Kivshar, J. Christou, V. Tikhonenko, B. Luther-Davies, L. M. Pismen, Dynamics of optical vortex solitons, Opt. Commun. 152 (1998) 198–206.
20. B. L. Lawrence, G. I. Stegeman, Two-dimensional bright spatial solitons stable over limited intensities and ring formation in polydiacetylene para-toluene sulfonate, Opt. Lett. 23 (1998) 591–593.
21. B.A. Malomed, Potential of interaction between two- and three-dimensional solitons, Phys. Rev. E 58 (1998) 7928–7933.
22. B.A. Malomed, Variational methods in nonlinear fiber-optics and related fields, Progr. Optics 43 (2002) 71–193.
23. B.A. Malomed, L.-C. Crasovan, D. Mihalache, Stability of vortex solitons in the cubic-quintic model, Physica D 161 (2002) 187–201.
24. B.A. Malomed, D. Mihalache, F. Wise, L. Torner, Spatiotemporal optical solitons, J. Optics B: Quant. Semicl. Opt. 7 (2005) R53–R72.
25. T. Mayteevarunyoo, B.A. Malomed, B. B. Baizakov, M. Salerno, Matter-wave vortices and solitons in anisotropic optical lattices, Physica D 238 (2009) 1439–1448.
26. D. Mihalache, D. Mazilu, L.-C. Crasovan, I. Towers, A. V. Buryak, B.A. Malomed, L. Torner, J.P. Torres, F. Lederer, Stable spinning optical solitons in three dimensions, Phys. Rev. Lett. 88 (2002) 073902.1–073902.4.
27. R. L. Pego, H. A. Warchall, Spectrally stable encapsulated vortices for nonlinear Schrödinger equations, Nonlinear Science 12 (2002) 347–394.
28. Kh. I. Pushkarov, D. I. Pushkarov, I. V. Tomov, Self-action of light beams in nonlinear media: Soliton solutions, Opt. Quant. Electr. 11 (1979) 471–478.
29. V. Prytula, V. Vekslerchik, V. M. Pérez-García, Eigenvalue cutoff in the cubic-quintic nonlinear Schrödinger equation, Phys. Rev. E 78 (2008) 027601.1–027601.4.
30. H. Pu, C. K. Law, J. H. Eberly, N. P. Bigelow, Coherent disintegration and stability of vortices in trapped Bose condensates, Phys. Rev. A 59 (1999) 1533–1537.
31. M. Quiroga-Teixeiro, H. Michinel, Stable azimuthal stationary state in quintic nonlinear optical media, J. Opt. Soc. Am. B 14 (1997) 2004–2009.
32. B. Y. Rubinstein, L. M. Pismen, Vortex motion in the spatially inhomogeneous conservative Ginzburg-Landau model, Physica D 78 (1994) 1–10.
33. H. Sakaguchi and B.A. Malomed, Solitary vortices and gap solitons in rotating optical lattices, Phys. Rev. A 79 (2009) 043606.
34. F. Smektala, C. Quémard, V. Couderc, A. Barthélémy, Non linear optical properties of chalcogenide glasses measured by z-scan, J. Non-Cryst. Solids 274 (2000) 232–237.
35. G. Boudebs, S. Cherukulappurath, H. Leblond, J. Troles, F. Smektala, F. Sanchez, Experimental and theoretical study of higher order nonlinearities in chalcogenide glasses, Opt. Commun. 219 (2003) 427–433.
36. I. Towers, A. V. Buryak, R. A. Sammut, B.A. Malomed, L.-C. Crasovan, D. Mihalache, Stability of spinning ring solitons of the cubic-quintic nonlinear Schrödinger equation, Phys. Lett. A 288 (2001) 292–298.
37. J. Yang, Z. H. Musslimani, Fundamental and vortex solitons in a two-dimensional optical lattice, Opt. Lett. 28 (2003) 2094–2096.
38. C. Zhan, D. Zhang, D. Zhu, D. Wang, Y. Li, D. Li, Z. Lu, L. Zhao, Y. Nie, Third- and fifth-order optical nonlinearities in a new stilbazolium derivative, J. Opt. Soc. Am. B 19 (2002) 369–377.
121698