# Stable three-dimensional spatially modulated vortex solitons in Bose-Einstein condensates

## Abstract

We present exact numerical solutions in the form of spatially localized three-dimensional (3D) nonrotating and rotating (azimuthon) multipole solitons in the Bose-Einstein condensate (BEC) confined by a parabolic trap. We numerically show that the 3D azimuthon solutions exist as a continuous family parametrized by the angular velocity (or, equivalently, the modulational depth). By a linear stability analysis we show that 3D azimuthons with a sufficiently large phase modulational depth can be stable. The results are confirmed by direct numerical simulations of the Gross-Pitaevskii equation.

###### pacs:

03.75.Lm, 05.30.Jp, 05.45.YvRecently, a novel class of two-dimensional (2D) spatially localized vortices with a spatially modulated phase, the so called azimuthons, was introduced in Ref. Kivshar2 (). Azimuthons represent intermediate states between the radially symmetric vortices and nonrotating multipole solitons. In contrast to the linear vortex phase, the phase of the azimuthon is a staircaselike nonlinear function of the polar angle. Various kinds of azimuthons have been shown to be stable in media with a nonlocal nonlinear response Lopez1 (); Lopez2 (); Skupin (); We2 (); We3 (). Note, that azimuthon solutions were found in Refs. Kivshar2 (); Lopez1 (); Lopez2 () by using an approximate (though rather accurate) variational approach and separation of variables. The first example of exact numerically found azimuthon solutions was presented in Ref. We4 ().

The aim of this Letter is to present exact numerical solutions in the form of spatially localized nonrotating and rotating (azimuthon) multipole solitons in the three-dimensional (3D) BEC confined by a parabolic trap. We numerically show that the azimuthon solutions exist as a continuous family parametrized by the angular velocity or, equivalently, by the parameter which determines the phase modulational depth. To our knowledge, this is the first time that solutions with such a nontrivial (3D azimuthon) topology can be found with very high (up to machine precision) accuracy. Moreover, by means of a linear stability analysis, we investigate the stability of these structures and show that rotating 3D dipole solitons (azimuthons with two intensity peaks) are stable provided that the number of atoms is small enough and the phase modulational depth is large enough. Thus, we present the first example of stable rotating dipole solitons in media with local nonlinearity. The presence of the trapping potential plays a crucial role in stabilizing the solitons. The results were confirmed by direct numerical simulations of the 3D Gross-Pitaevskii equation (GPE).

We consider a condensate at zero temperature confined in an axisymmetric harmonic trap. The dynamics of the condensate is described by the normalized GPE for the wave function

(1) |

where appropriate dimensionless units are used, , where the sign corresponds to attractive (repulsive) contact interaction.

Equation (1) conserves the 3D norm (the normalized number of particles) , -component of the angular momentum , and energy

(2) |

We will seek solutions of Eq. (1) which are stationary in the frame rotating with the angular velocity . In cylindrical coordinates , such solutions of the form

(3) |

where is the chemical potential in the rotating frame, satisfy the equation

(4) |

where , and resolve the variational problem for the functional . The chemical potential in the laboratory frame is related to the chemical potential in the rotating frame by

(5) |

In what follows, we restrict ourselves to the case of a spherically symmetric trap () and attractive interaction (). Two-dimensional case corresponding (”pan-cake configuration”) was considered in Ref. We4 (). Equation (4) was rewritten in Cartesian coordinates

(6) |

and solved numerically with zero boundary conditions. The numerical method that we use to find solutions of Eq. (6) was first introduced by Petviashvili Petviashvili1 () (see also Ref. Petviashvili2 ()) and more recently in Refs. Pelinovsky1 (); Ablowitz (); Fibich (); Pelinovsky2 (); Lakoba1 (). Petviashvili employed a stabilizing factor to suppress divergence of the iteration procedure (or convergence to zero solution). Petviashvili and the authors of Refs. Pelinovsky1 (); Ablowitz (); Fibich (); Pelinovsky2 () considered periodic boundary conditions and worked in the spectral Fourier representation. We work in physical space and use zero boundary conditions so that the method should be modified. We describe it in the following way. Equation (6) can be written as , where stands for the linear part of the equation and is the nonlinear term. At each th stage of the iteration we solve the corresponding linear problem and the iteration procedure is

(7) |

where the stabilizing factor

(8) |

where , and for the cubic nonlinearity . Note, that one can also take (sometimes it yields faster convergence)

(9) |

with in Eq. (7). It can be seen that the right-hand side of Eq. (7) has homogeneity zero with respect to , and this, indeed, prevents the aforementioned divergence. The convergence can be monitored using the value , which should approach zero.

As an initial guess we use

(10) |

on a Cartesian grid, where is an integer. The parameter determines the modulation depth of the soliton intensity and without loss of generality we can assume . The ansatz Eq. (10) describes topology of the 3D azimuthon solutions. Note that the case corresponds to nonrotating multipole 3D solitons (e. g. to a dipole, to a quadrupole etc.), while the opposite case corresponds to the radially symmetric 3D vortices with the topological charge . A detailed study, including a linear stability analysis, of 3D vortices in BEC with attractive interactions and a parabolic trap was performed in Refs. Malomed1 (); Malomed2 ().The intermediate case corresponds to the rotating 3D azimuthons.

Note, that Eq. (6) with has localized solutions provided that Malomed1 (), where is the chemical potential in the laboratory frame. The fundamental solitons of Eq. (6) with and exist only if . Moreover, the radially symmetric 3D vortices with are stable if (for ) Malomed1 (); Malomed2 (), so that in the following we consider only the region .

The rate of convergence does not depend on the choice of , , and . In all runs we used , , and between and . During the iterations, the parameter (modulational depth), which is similar to the parameter in Eq. (10), can be introduced (after eliminating the constant phase factor in ) in the following way:

(11) |

We start with the case . The iteration procedure monotonically converges to a solution that consists of two dipole-shaped structures in the real and imaginary parts of (dipole azimuthon, or, azimuthon with two intensity peaks). If the amplitudes of these parts are equal (i. e. ) we get the radially symmetric 3D vortex with the topological charge . If these amplitudes are different (), we have a rotating azimuthon; if one of the amplitudes is zero (), we have the nonrotating 3D dipole soliton. The final value of depends only (for the fixed ) on the rotational velocity and does not depend on the initial guess of (for example, the initial condition with and converges to the nonrotating dipole with ). Typically we used a spatial grid with a resolution .

The progressive iterations in the relaxation method were terminated when the norm of the residual fell below (with the corresponding ). Under this, the relative differences between the successive amplitudes , where , and modulational depths were less than and respectively. Typically, from (for solutions with small ) to (for solutions with close to ) iterations were required for convergence. Increasing the number of iterations and grid points, we were able to find the solutions with accuracy up to machine precision.

By changing one can find the azimuthon solution with arbitrary and, thus, there exists a continuous family of azimuthons parametrized by the parameter . The dependence of the angular velocity of the azimuthon with two intensity peaks on the modulational depth for several different and (attractive interaction) is presented in Fig. 1. In Figs. 2(a) and (b) we demonstrate two numerically found examples of the azimuthons (nonrotating dipole and rotating azimuthon with two intensity peaks) for the model described by Eq. (1). Similar solutions can be found for (repulsive interaction) and .

The vortex solitons have the maximum angular momentum , while the dipole solitons have zero angular momentum. Structures with the intermediate values of carry out the nonzero angular momentum, and this leads to soliton rotation.

Changing topology of the initial guess one can find high-order azimuthon solutions. We took Eqs. (10) with . The iteration procedure began to converge monotonically to a solution that consists of two quadrupole-shaped structures in the real and imaginary parts of . If the amplitudes of these parts are equal () we get the radially symmetric vortex with the topological charge . In the intermediate case (), we have the azimuthon with four intensity peaks. The convergence was controlled by stopping the iteration when the value began to increase. This indicates that the iteration procedure jumps off the solution corresponding high-order azimuthon (and begins to converge to the azimuthon with two intensity peaks). Nevertheless, the high-order solutions can be found with a rather high accuracy: typically, one can reach . Then, the obtained solution can be used as an initial condition in the Yang-Lakoba iterative procedure Lakoba2 () which belongs to a family of universally-convergent iterative methods and can converge to any nonfundamental (i. e. high-order) solution of a given equation provided that the initial condition is sufficiently close to that solution. Thus, we were able to find high-order azimuthons of Eq. (1) with very high accuracy. The nonrotating quadrupole and rotating azimuthon with four intensity peaks are presented in Figs. 2(c) and (d).

The physical relevance of any stationary solution, of course, depends on whether it is stable. To study the stability of the stationary solutions, we represent the wave function in the form

(12) |

where the stationary solution is perturbed by a small perturbation . In the absence of the radial symmetry, the corresponding eigenvalue problem (i. e. after linearizing and taking ) on a spatial grid implies a complex nonsymmetric matrix and, for reasonable (say, ), represents extremely difficult task. Instead, after inserting Eq. (12) into Eq. (1), we solved the Cauchy problem for the linearized equation

(13) |

with some initial small perturbations . The final results are not sensitive to the specific form of . The chemical potential in the laboratory frame is determined from Eq. (5). If the dynamics is unstable, the corresponding solutions , undergoing, generally speaking, oscillations, grow exponentially in time and an estimate for the growth rate can be written as

(14) |

where , and are assumed to be large enough.

In Fig. 3 we plot the growth rates as functions of the modulational depth for the azimuthons with two intensity peaks and several different values of the chemical potential in the rotating frame. The linear stability analysis shows that for solutions with , the growth rate of perturbations for all so that there is no stability region. In particular, all nonrotating dipoles are unstable. The growth rate decreases as increases and can be very small if is close to . The situation, however, changes when the modulational depth exceeds the critical value . In this case, the growth rate of perturbations falls to zero, the stability window appears and the azimuthons with are stable. The critical value varies very slightly with changing the chemical potential . All high-order azimuthon solutions turn out to be unstable (for the radially symmetric vortices with and the topological charge it was shown in Ref. Malomed1 (); Malomed2 ()).

To verify the results of the linear analysis, we solved numerically the dynamical equation (1) initialized with our computed solutions with added Gaussian noise. The initial condition was taken in the form , where is the numerically calculated solution, is the white gaussian noise with variance and the parameter of perturbation . The unstable dynamics of the nonrotating dipole () with is illustrated in the top row of Fig.4. In the middle row we present unstable evolution of the rotating azimuthon with two intensity peaks and , , . A slight stochastic perturbation with was applied at . The azimuthon lost its shape after one rotational period. Stable evolution of the azimuthon with , and (i. e. in the region of stability) is shown in the bottom row of Fig.4. The initial state of the azimuthon is perturbed by a rather strong noise (). The period of rotation of the azimuthon is . The azimuthon cleans up itself from the noise and survives over many dozens of the rotational periods.

In the present paper we restricted ourselves to the case of a spherically symmetric trap (). A linear stabilty of the radially symmetric 3D vortices (i. e. azimuthons with in our case) was performed for different in Ref. Malomed2 (). By analogy with the results of Ref. Malomed2 (), for the case (”cigar-like configuration”) one could expect that the quasi-1D azimuthon solitons with are tightly confined in the corresponding cigar-like trap, which suppresses their destabilization. On the other hand, the case (”pan-cake configuration”) implies nearly 2D azimuthon solitons, with the largest energy, which are most prone to instabilty.

In conclusion, we have presented exact numerical solutions in the form of spatially localized 3D nonrotating and rotating (azimuthon) multipole solitons in the three-dimensional BEC confined by a parabolic trap. We have shown that the 3D azimuthon solutions exist as a continuous family parametrized by the angular velocity (or, equivalently, the phase modulational depth). By a linear stability analysis we have shown that 3D azimuthons with a sufficiently large modulational depth can be stable. The results are confirmed by direct numerical simulations of the Gross-Pitaevskii equation.

The author thanks A. I. Yakimenko and Yu. A. Zaliznyak for discussions.

### References

- A. S. Desyatnikov, A. A. Sukhorukov, and Yu. S. Kivshar, Phys. Rev. Lett. 95, 203904 (2005) .
- S. Lopez-Aguayo, A. S. Desyatnikov, Yu. S. Kivshar, S. Skupin, W. Krolikowski, and O. Bang, Opt. Lett. 31, 1100 (2006).
- S. Lopez-Aguayo, A. S. Desyatnikov, and Yu. S. Kivshar, Opt. Express 14, 7903 (2006).
- S. Skupin, O. Bang, D. Edmundson, and W. Krolikowski, Phys. Rev. E 73, 066603 (2006).
- V. M. Lashkin, Phys. Rev. A 75, 043607 (2007).
- V. M. Lashkin, A. I. Yakimenko, and O. O. Prikhodko, Phys. Lett. A 366, 422 (2007).
- V. M. Lashkin, Phys. Rev. A 77, 025602 (2008).
- V. I. Petviashvili, Sov. J. Plasma Phys. 2, 257 (1976).
- V. I. Petviashvili and V. V. Yan’kov, in Reviews of Plasma Physics, edited by B. B. Kadomtsev, (Consultants Bureau, New York, 1989), Vol. 14.
- D. E. Pelinovsky and Yu. A. Stepanyants, SIAM J. Numer. Anal. 42, 1110 (2004).
- M. J. Ablowitz, Z .H. Musslimani, Opt. Lett. 30, 1 (2005).
- G. Fibich, Y. Sivan, and M.I. Weinstein, Physica D 217, 31 (2006).
- M. Chugunova and D. Pelinovsky, Discrete Contin. Dyn. Syst., Ser. B 8, 773 (2007).
- T. I. Lakoba and J. Yang, J. Comput. Phys., 226, 1693 (2007).
- J. Yang and T. Lakoba, Stud. Appl. Math., 118, 153 (2007).
- D. Mihalache, D. Mazilu, B. A. Malomed, and F. Lederer, Phys. Rev. A 73, 043615 (2006).
- B. A. Malomed, F. Lederer, D. Mazilu, and D. Mihalache, Phys. Lett. A 361, 336 (2007).