# Stability of vortices in rotating traps: A 3D analysis

## I Introduction

Since the first experimental realization of Bose-Einstein condensation (BEC) in weakly interacting gases , there has been a huge theoretical and experimental effort to study its properties in the framework of fully quantum theories and in the so called mean field limit (Gross-Pitaevskii -GP- equations). These equations are formally Nonlinear Schrödinger Equations (NLS)  which appear in many fields of physics, e.g. in bulk superfluids and nonlinear optics to cite only a few examples.

All of these physical systems have been long known to exhibit solutions corresponding to topological defects , one of the simplest being known as vortices (in two spatial dimensions) or vortex-lines (in three spatial dimensions). Vortices are localized phase singularities with integer topological charge which analogous in the hydrodynamic interpretation to vortices those appearing in fluid dynamics . In the framework of BEC studies it has been raised the question of whether these nonuniform clouds of condensed gases may support the existence of vortices in a stable form.

There is a huge literature on vortices and vortex properties in the framework of NLS equations (including its particular cubic version, the GP equation) and its nonconservative extensions, the Ginzburg-Landau (GL) system, and vector GL models. In particular the stability of $`m`$-charged GP vortices in two dimensions was studied in . In three dimensions the GL case has been recently considered  and vortex lines geometric instabilities have been found to strongly deform the vortex lines. However the GP equation cannot be obtained as a limit of the GL studied there since dissipation and diffusion are essential ingredients of the models studied in Ref. . This fact makes the conservative case (GP) interesting by itself. Other analysis of vortices and vortex stability in the framework of Nonlinear Optics are included in Ref. .

The current setups utilized to generate Bose-Einstein condensates use a magnetic trap to confine the atomic cloud which is modeled by a parabolic trapping potential. This is a distinctive feature from the common NLS systems, in which the vortices are free and move in an homogeneous background. The dynamics of a vortex in a spatially inhomogeneous two dimensional GP problem was studied in Ref.  using the method of matched asymptotic expansions, but the authors did not consider the stability of the 2D vortex itself. In principle, the vortex motion equations  can be used to study the motion of a single 2D point vortex in spatially inhomogeneous GP problems. However, the dynamics of the many vortex case is more complicated and by no means trivial. For simple approaches to the problem which do not include the effect of vortex cores on the background field see Ref. . The dynamics of 3D vortices is yet more complicated allowing the so-called reconnection. To our knowledge there are no analytical results but only qualitative numerical observations available . Another theoretical framework where non-homogeneous dynamics of vortices has been investigated is the possibility of pinning vortices in type-II superconductors , but it is only the dynamics has been considered through analytical approximation techniques with no comparison with numerics. In all the previously discussed cases the vortex stability is given for granted.

In the framework of Bose–Einstein condensed gases studies the problem of the vortex stability has been considered in various papers that try to solve the problem of linear and global stability, either from a purely analytical point of view, such as in , or by mixing analytic and numerical techniques, .

In Ref.  the authors solve the Gross-Pitaevskii equation and find the energies of the condensate in vortex states, for a number of particles up to $`N=10^4`$. In Ref.  the authors solve the Bogoliubov equations for an unit charge vortex in a stationary trap with axial symmetry, their results being also limited to $`N<10^4`$. In Ref.  the authors address the problem of minimizing the energy functional with a reduced basis of trial states that is only valid in the limit of small $`U`$.

In this paper we unify and substantially extend what has been done in previous works regarding these two questions: global energetic stability and local stability of vortex states. First, in Sect. II we solve the GPE for an axially symmetric harmonic potential, with or without the action of a uniform magnetic field which resembles the effect of a rotating trap. We calculate the lowest stationary solutions that have a well defined value of the third component of the angular momentum, $`m=L_z`$, and we do this for small and for very large values of the nonlinearity ($`N10^7`$). We find that there are continuous intervals of the “angular velocity”, $`[\mathrm{\Omega }_m,\mathrm{\Omega }_{m+1})`$, in which the $`m`$-charged vortex state becomes energetically stable with respect to other states of well defined vorticity. In Sect. III, we study Bogoliubov’s equations from two different points of view: as a consequence of a linear stability analysis of the Gross-Pitaevskii equation (GPE), and as the first corrections to the mean field theory of the dilute condensate. The concepts of dynamical and energetic stability are defined, and it is demonstrated that any possible destabilization of the system must be either of energetic nature, or grow polynomially with respect to time. We next solve the Bogoliubov equations for $`m=1`$ and $`m=2`$ unperturbed vortex states in stationary traps. It is found that the $`m=1`$ and $`m=2`$ vortices are only energetically unstable, which means that the lifetime of both configurations is only limited by dissipation. A similar treatment reveals that rotation can only stabilize the unit charge vortex-line if the angular speed is in a suitable range, $`\mathrm{\Omega }[\mathrm{\Omega }_1,\mathrm{\Omega }_2)`$, while outside of this range, $`\mathrm{\Omega }_2<\mathrm{\Omega }<\mathrm{\Omega }_c`$, the minimum of the energy functional is not an eigenstate of the $`L_z`$ operator –i.e, it is not symmetric under rotations. These results are confirmed by numerical simulations of the evolution of perturbed vortices. In Sect. IV we summarize our work and discuss their implications.

## II Vortex solutions of the GPE

### A Stationary states of GPE in a uniform and constant magnetic field

For small temperatures and small densities, the condensate is modeled by the Gross-Pitaevskii equation (GPE) . We will always refer to an axially symmetric trap with term that accounts for rotation around the Z axis and may be generated by a weak magnetic field. The form of the equation is

$`i\mathrm{}{\displaystyle \frac{\psi }{t}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{}\psi +{\displaystyle \frac{1}{2}}m\omega ^2\left(\gamma ^2r^2+z^2\right)\psi `$ (1)   
$`+`$ $`U_0N|\psi |^2\psi +\stackrel{~}{\mathrm{\Omega }}L_z\psi .`$ (2)   
Here $`U_0=4\pi \mathrm{}^2a/m`$ characterizes the interaction and is defined in terms of the ground state scattering length $`a`$. In all cases we will take the normalization condition to be

$$|\psi |^2d^3x=1.$$
(3)   
It is convenient to express Eq. (1) in a natural set of units which for our problem is built up from two scales: the trap size (measured by the width of the linear ground state), $`a_0=\sqrt{\mathrm{}/m\omega }`$, and its period, $`\tau =1/\omega `$. With these definitions the equation simplifies to

$$i\frac{\psi }{t}=\left[\frac{1}{2}\mathrm{}+i\mathrm{\Omega }\frac{}{\theta }+\frac{1}{2}(\gamma ^2r^2+z^2)+U|\psi |^2\right]\psi ,$$
(4)   
while maintaining the normalization.

The new parameters, $`\mathrm{\Omega }=\mathrm{}\stackrel{~}{\mathrm{\Omega }}`$ and $`U=4\pi Na/a_0`$, represent the “angular speed” of the trap and the adimensionalized interaction strength, respectively. For stability reasons (see below), $`\mathrm{\Omega }`$ will be of the order of magnitude of or smaller than the strength of the trapping, $`\omega `$. The other parameter, $`U`$, will take values from $`0`$ to $`6\times 10^4`$. As of the experiments with rubidium and sodium, this implies a minimum of $`10^6`$ and a maximum of $`10^7`$ atoms which is in the range of current and projected experiments. The shape of the trap is dictated by the geometry factor, and in this work it will typically take two possible values: $`\gamma =1`$, corresponding to a spherically symmetric trap, and $`\gamma =2`$, corresponding to an axially symmetric, elongated trap.

A stationary solution of (4) will be of the form $`\psi (\stackrel{}{x},t)=e^{i\mu t}\varphi (\stackrel{}{x})`$ , where $`\mu `$ may be interpreted both as a frequency and the chemical potential

$$\mu \varphi =\left[\frac{1}{2}\mathrm{}+i\mathrm{\Omega }\frac{}{\theta }+\frac{1}{2}(\gamma ^2r^2+z^2)+U|\varphi |^2\right]\varphi .$$
(5)   
Any solution of (4) has an energy *per* particle which is given by the functional

$`E(\psi ,N)`$ $`=`$ $`{\displaystyle \left(\frac{1}{2}|\psi |^2i\mathrm{\Omega }\overline{\psi }_\theta \psi \right)}`$ (6)   
$`+`$ $`{\displaystyle \left[\frac{1}{2}\left(\gamma ^2r^2+z^2+U\left|\psi \right|^2\right)\left|\psi \right|^2\right]}`$ (7)   
For a stationary solution it becomes

$$E(\psi ,N)=\mu \frac{U}{2}\left|\varphi \right|^4.$$
(8)   
The stationary solutions of (4) may also be interpreted as the minimization of

$$_\mu =E(\psi ,N)\mu |\psi |^2$$
(9)   
subject to the constrain of Eq. (3). In that case $`\mu `$ is nothing else but the Lagrange multiplier of the norm.

Since we are interested in single vortex solutions to the GPE we will restrict our analysis to stationary states that are also eigenstates of the $`L_z`$ operator. That is, we will look for solutions of the form $`\psi (r,z,\theta ,t)=e^{i\mu t}e^{im\theta }\varphi (r,z).`$ Summarizing, our goal will be to find the unit norm functions $`\varphi _\mu ^{(m)}(r,z)`$ and real numbers $`\mu `$ which are solutions of the equation

$$\mu \varphi _\mu ^{(m)}=\left[\frac{1}{2}\mathrm{}m\mathrm{\Omega }+\frac{1}{2}(\gamma ^2r^2+z^2)+U|\varphi _\mu ^{(m)}|^2\right]\varphi _\mu ^{(m)},$$
(10)   
Our treatment on the following sections will be fully three-dimensional and no spurious conditions (e.g. periodicity) will be imposed on the boundaries. We want to obtain at least the lowest energy state for each value of the vorticity, $`m`$. Also the dependency of spectrum with the nonlinearity and the angular velocity, $`\mathrm{\Omega }`$, is interesting since it will allow us to find whether the vortex-line states may become energetically favorable.

### B Numerical method

Due to the nonlinear nature of the problem we want to solve \[Eq. (10)\] there are not many analytical tools available. The most common (and maybe easiest) approach to the problem is to discretize the spatial part and perform time evolution in imaginary time while trying to preserve the normalization, a method which is related to the steepest descent. The precision of the solution depends of the type of spatial discretization used: finite differences (used for example in Refs. ) or spectral methods (such as the one used in Ref. ). However these common methods, such as finite differences  and similar spectral methods , have reached a maximum value of the interaction of $`U=10^3`$, which should be contrasted with the value $`U=10^5`$ that can be attained which the technique to be presented later.

Properly speaking our technique is a Galerkin type method where one performs the expansion of the unknown solution in a complete basis of the Hilbert space under consideration. For convenience we have used the basis of eigenstates of the harmonic oscillator with fixed vorticity. In this basis our stationary solution is expressed as

$$\psi _\mu ^{(m)}(\stackrel{}{x},t)=e^{i\mu t}e^{im\theta }\underset{n}{}c_nP_n^{(m)}(r,z),$$
(11)   
Here the single index, $`n`$, denotes two quantum numbers, $`(n_z,n_r)`$, regarding the axial and radial degrees of freedom, and $`P_n^{(m)}`$ is a product of a Hermite polynomial, a Laguerre polynomial and a Gaussian

$`P_n^{(m)}`$ $`=`$ $`C_nH_{n_z}(z)L_{n_r}(\rho ^2)r^me^{im\theta }e^{(\rho ^2+z^2)/2},`$ (13)   
$`C`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{\sqrt{\gamma }\sqrt{\pi }2^{n_z}n_z!}}}\sqrt{{\displaystyle \frac{n_r!}{\pi \left(n_r+m\right)!}}},`$ (14)   
with $`\rho =r/\sqrt{\gamma }`$

Next, following the same convention about the indices, we have introduced this expansion into Eq. (10) to obtain

$$\left(E_i^{ho}\mathrm{\Omega }m\mu \right)c_i+U\underset{jkl}{}A_{ijkl}^{(m)}\overline{c}_jc_kc_l=0,$$
(15)   
where $`E_i^{ho}`$ is the harmonic oscillator energy of the mode $`P_{im}`$ and the tensor $`A^{ijkl}`$ has the following definition

$$A_{ijkl}^{(m)}=2\pi \overline{P}_i^{(m)}\overline{P}_j^{(m)}P_k^{(m)}P_l^{(m)}𝑑r𝑑z.$$
(16)   
Since the $`P_i^{(m)}`$ are products of known polynomials by exponentials it could be possible, in principle, to evaluate the coefficients in the tensor exactly with a Gaussian quadrature formula of the appropriate order. This approach was used in Ref.  for the three-dimensional case. However, when one wishes to use a large number of modes (which in our case is of about $`1600`$ for each value of $`m`$) to achieve large nonlinearities, the search of the quadrature points becomes more difficult than performing a stable integration by means of some other methods, of which the simplest accurate one is Simpson’s rule.

Once we fix all of the constants, $`E^{ho}`$, $`A_{ijkl}^{(m)}`$, $`\mu `$ and a guess for the solution, it is feasible to solve (15) iteratively -e.g. by Newton’s method. However, it is wiser to perform two simplifications before implementing the algorithm. The first one is that all of the eigenfunctions, $`P_n^{(m)}`$, can be made real and thus we can impose the coefficients in the expansion, $`\{c_n\}`$, to also be real.

The second optimization is that, thanks to the symmetry of the problem, the ground state of Eq. (5) has a well defined positive parity. This allows us to eliminate redundant modes , saving memory and reaching higher energies and nonlinearities which otherwise would be computationally hard to attain. On the other hand, we have always checked that this method produced the same results as the complete one for a selected and significant set of parameter values.

And finally it is important to note that the four–index tensor (16) is indeed a product of two smaller tensors, corresponding to the integration on the $`z`$ and $`r`$ variables . This decomposition is most important when working with a large number of modes, because then the size of $`A`$ becomes extremely large (i.e., $`1600^4`$ elements for $`1600`$ modes).

Concerning the evaluation of very high order polynomials as the ones involved in our computations it is necessary to say that it is not a simple task, specially for intermediate values of the spatial variables since then there is a lot of comparable terms with usually different signs and the cancellations induce numerical instabilities. The usual procedure to avoid this difficulty is to use Horner’s method  to evaluate the polynomial, which is comparable to using FFT techniques, but in our case this is not enough and the evaluation of higher order polynomials could be done only by the recursion formulas for the Hermite and Laguerre polynomials.

We remark that the election of this spectral technique was largely influenced by the need of reaching high nonlinearities which are not achievable using the other approaches. Further details on the numerical technique as well as convergence proofs will be given elsewhere 

### C Results for the stationary states and spectrum

By using the preceding technique, we have searched the lowest states, $`(n_z,n_r=0)`$, for each branch of the spectrum with a different vorticity, $`m=0,\mathrm{},6`$. This was performed for two geometries corresponding to $`\gamma =1`$ (spherically symmetric trap) and $`\gamma =2`$ (cigar shape trap), of a static trap, $`\mathrm{\Omega }=0`$, while varying the intensity of the interaction from $`0`$ to approximately $`50000`$. The results of this study are plotted in Fig. 1.

Remarkably, in the absence of rotation, and up from the lowest states, both the spectrum and the energies can be fitted to a simple formula

$$\mu _{0m}(N)=\mu _{00}(N)+\omega _{\text{eff}}(N)m.$$
(17)   
The first term is the chemical potential of the $`m=0`$ ground state and is not relevant to the dynamics. Using the Thomas-Fermi approximation one can show that it grows proportionally to $`\mu N^{2/5}`$, a behavior which is approximately reflected in the numerical results shown in Fig. 1(c).

The second term is much more relevant to the evolution of the condensate. It grows linearly, as the energy levels of a linear harmonic oscillator, with an effective frequency, $`\omega _{eff}(N)`$, that decreases with the interaction. The fact that the higher levels of the spectrum of $`\mu `$ remain equispaced even for large interactions is the reason why the condensate exhibits an exponentially divergent response to the parametric perturbation of the trap frequencies, as it is shown in Ref.  and .

Now we want to study the stationary solutions in the presence of rotation. For $`\mathrm{\Omega }0`$ the proper functions with definite vorticity remain the same, while the chemical potential and the energy suffer a shift that depends on the vorticity of the state

$$E_{nm}(U,\mathrm{\Omega })=E_{nm}(U,0)m\mathrm{\Omega }.$$
(18)   
This shift gives rise to an ample phenomenology which is pictured in Fig. 2.First, we see that the degeneracy with respect to $`m`$ is broken. The only other possible degeneracy that remains is with respect to the $`r`$ and $`z`$ variables, but this will be removed for the case without spherical symmetry, $`\gamma 1`$.

Second, the $`m=1,2,3\mathrm{}`$ branches of the spectrum become a minimum of the energy functional with respect to other branches for continuous intervals of the angular velocity, $`[\mathrm{\Omega }_m,\mathrm{\Omega }_{m+1}]`$, where

$$\mathrm{\Omega }_m=E_{0,m+1}E_{0m}.$$
(19)   
However this does not mean that in these intervals the $`m`$-th vortex state becomes a global minimum. Indeed, in Sect. III we will only be able to prove that only the $`m=1,2`$ vortex lines achieve the status of local minimum. It still remains an open question under which situations the ground state must have a well defined vorticity.

Third, even though the separation between the $`m=0`$ and $`m=1`$ states becomes very narrow for large interactions, the stabilization frequency $`\mathrm{\Omega }_1`$ only approaches zero asymptotically with $`U`$. As a consequence, $`m=1`$ states are never a global minimum of the energy in a stationary trap, a fact that can be checked by just inspecting the energy functional.

And finally, there is a critical value of $`\mathrm{\Omega }`$ for which the energy functional becomes unbounded by below \[See Fig. 2\] and which coincides with the separation between energy levels for large values of the vorticity. This critical value of the frequency, $`\mathrm{\Omega }_c`$, is such that all of the ground states for each value of the vorticity have the same energy

$$E_{0m}(U,\mathrm{\Omega }_c)=E_{0k}(U,\mathrm{\Omega }_c),k,m.$$
(20)   
Using Eqs. (20) and a fit similar to the one in Eq. (17), one finds that it is always smaller than the critical frequency of the linear case

$$\mathrm{\Omega }_c=\omega _{\text{eff}}(U).$$
(21)   
## III Stability of stationary states

### A The linear stability equations

In the preceding section we obtained stationary solutions of the mean field model for the Bose-Einstein condensate, all of which had a well defined value of the third component of the angular momentum operator. We named those states vortices. In this section we now want to study the stability of these solutions according to several criteria of a local nature: local energetic stability and linear stability.

We begin our study from the adimensionalized Gross-Pitaevskii equation (4). First we expand the condensate wavefunction around a stationary solution with a fixed vorticity.

$`\mathrm{\Psi }(r,z,\theta ,t)`$ $`=`$ $`\psi _0+ϵ\psi _1`$ (22)   
$`=`$ $`\left[f(r,z)e^{im\theta }+ϵ\alpha (r,z,\theta ,t)\right]e^{i\mu (\mathrm{\Omega })t}.`$ (23)   
We insert this expansion in Eq. (4) and truncate the equations up to $`𝒪(ϵ^1)`$ thus getting

$`i_t\alpha `$ $`=`$ $`\left[H_0+i\mathrm{\Omega }_\theta +2Uf^2\right]\alpha +Uf^2e^{2im\theta }\overline{\alpha },`$ (25)   
$`i_t\overline{\alpha }`$ $`=`$ $`\left[H_0i\mathrm{\Omega }_\theta +2Uf^2\right]\overline{\alpha }+Uf^2e^{2im\theta }\alpha ,`$ (26)   
with $`H_0=\frac{1}{2}\mathrm{}+\frac{1}{2}(\gamma ^2r^2+z^2)\mu (\mathrm{\Omega })`$. We can also write this equation in a more compact form

$$i\frac{}{t}\stackrel{}{W}=\sigma _z(\mathrm{\Omega })\stackrel{}{W}=(\mathrm{\Omega })\stackrel{}{W}.$$
(27)   
by using the definitions

$`\stackrel{}{W}`$ $`=`$ $`\left(\begin{array}{c}\alpha \\ \overline{\alpha }\end{array}\right),`$ (31)   
$`\sigma _z`$ $`=`$ $`\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),`$ (34)   
$`(\mathrm{\Omega })`$ $`=`$ $`H_0+\left(\begin{array}{cc}i\mathrm{\Omega }_\theta +2Uf^2& Uf^2e^{2im\theta }\\ Uf^2e^{2im\theta }& i\mathrm{\Omega }_\theta +2Uf^2\end{array}\right).`$ (37)   
In the rest of this work we wish to study the dynamics that is involved in Eq. (27). The simplest way to achieve this is to find a suitable basis in which Eq. (27) becomes diagonal or almost diagonal. In other words, we want a set of vectors, $`\stackrel{}{W}_k^t=(u_k(r),v_k(r))`$ such that

$$\lambda _k\stackrel{}{W}_k=\stackrel{}{W}_k.$$
(38)   
If $``$ has such a diagonal Jordan form, then the perturbation evolves simply as

$`\stackrel{}{W}`$ $`=`$ $`{\displaystyle c_ke^{i\lambda t}\stackrel{}{W}_k},`$ (39)   
$`\alpha (\stackrel{}{r},t)`$ $`=`$ $`{\displaystyle c_ke^{i\lambda t}u_k(\stackrel{}{r},t)}.`$ (40)   
On the other hand, the lack of a diagonal form, or the existence of complex eigenvalues leads to instability in a way that we will precise later.

Associated to Eq. (27) there is an energy functional,

$$E_2(\alpha )=2\overline{\alpha }H_0\alpha +\psi _0^2\overline{\alpha }^2+\overline{\psi _0}^2\alpha ^2+4|\psi _0|^2\alpha \overline{\alpha },$$
(41)   
and a constrained energy functional

$$_2(\alpha )=E_2(\alpha )\mu |\alpha |^2$$
(42)   
which are the $`𝒪(ϵ^2)`$ terms in the expansion of (8) and (9), i.e. the energy introduced in the system by the perturbation. If a diagonal Jordan form like the one of (38) is possible, then it is easy to check that the second functional becomes diagonal, too

$`_2(\alpha )`$ $`=`$ $`{\displaystyle |c_k|^2\lambda _kG(\stackrel{}{W}_k)},`$ (44)   
$`G(\stackrel{}{W}_k)`$ $`=`$ $`{\displaystyle |u_k|^2}|v_k|^2.`$ (45)   
If the stationary state, $`\psi _0`$, is a local minimum of the energy subject to the constrain of a fixed norm (3), then $`_2`$ must be positive for all perturbations, which has serious implications for the eigenvalues and eigenstates. We will refer to this later.

When studying the condensate using tools from Quantum Field Theory, one may try a similar procedure , which is known as Bogoliubov’s theory. In that framework, the $`\overline{\alpha }`$ and $`\alpha `$ are linear operators in a Fock space, and one searches an expansion of these operators in terms of others that diagonalize the energy functional (41) and the evolution equations (27). The resulting equations for the coefficients are known as Bogoliubov’s equations and correspond to the equations (38) for $`u_k`$ and $`v_k`$.

### B Operational procedure

It is now useful to perform an expansion of $`\alpha `$ and $`\overline{\alpha }`$ into states of fixed vorticity so that the modes are separated into subspaces according to their vorticities

$$\stackrel{}{W}_i^{(n)}=\left(\begin{array}{c}u_i^{(n)}(r)e^{in\theta }\\ v_i^{(2mn)}(r)e^{i(2mn)\theta }\end{array}\right).$$
(46)   
These subspaces are not mixed by the action of the operators of (27), and we can define their restriction to these subspaces

$`_n(\mathrm{\Omega })`$ $`=`$ $`\sigma _z_n(\mathrm{\Omega }),`$ (48)   
$`_n(\mathrm{\Omega })`$ $`=`$ $`_{0n}(\mathrm{\Omega })+𝒰_n,`$ (49)   
$`_0`$ $`=`$ $`\left(\begin{array}{cc}H^n(nm)\mathrm{\Omega }& \\ & H^{2mn}(mn)\mathrm{\Omega }\end{array}\right),`$ (52)   
$`H^n`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{}+{\displaystyle \frac{1}{2}}(\gamma ^2r^2+z^2)+{\displaystyle \frac{n^2}{2r}}+f^2\mu (0)`$ (53)   
$`𝒰_n`$ $`=`$ $`Uf^2\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right).`$ (56)   
With these definitions the diagonalization procedure (38) becomes

$$\lambda _k^n\stackrel{}{W}_k^{(n)}=_n(\mathrm{\Omega })\stackrel{}{W}_k^{(n)},nm.$$
(57)   
If $`G(\stackrel{}{W}_k^{(n)})>0`$, then $`(u_k^{(n)},v_k^{(2mn)})`$ is a Bogoliubov mode with energy $`ϵ=\lambda _k^{(n)}`$ and vorticities $`(n,2mn)`$, whereas if $`G(\stackrel{}{W}_k^{(n)})<0`$ then the excitation is $`(v_k^{(2mn)},u_k^{(n)})`$ with energy $`ϵ=\lambda _k^{(n)}`$. As a rule of thumb, the $`u`$ function must always be the one with the largest contribution, which is formally stated in $`G(\stackrel{}{W})>0`$. In the following we will refer to these branches of the spectrum by the pairs of quantum numbers $`(n,2mn)`$ and $`(2mn,n)`$, respectively.

One may find, in principle, two kinds of solutions. First, the Bogoliubov operator may have complex eigenvalues or even have a non diagonal Jordan form. In both cases we speak of *dynamical instability* because an arbitrarily small perturbation departures from the original state exponentially or polynomially in time. Second, the linearized operator may have only real eigenvalues which should be interpreted as the change of energy in the condensate due to excitations \[See Eq. (42)\]. If $`\lambda >0`$ the state under study $`\psi _0`$ is a local minimum of the energy functional (6) with respect to this family of perturbations, the $`\lambda =0`$ case corresponds to the existence of degeneracy in the system, and finally if $`\lambda <0`$ the system is told to be energetically unstable -i.e, excitations are energetically favorable and the state is not a local minimum of the energy.

All of the five cases exposed above have the same implications of stability for Eq. (4), which is a simple partial derivatives equation for an order parameter, and for the more complete Bogoliubov theory, where the perturbations are regarded as many-body corrections and involve more degrees of freedom. Nevertheless, it must be remarked that of the two types of instability that can be found, i.e. dynamical and energetic instabilities, the second one is less harmful because it only affects the dynamics when there is some kind of dissipation that drives the system through the unstable branch. And even then the lifetime of the system can be significant if the intensity of the destabilizing mode is small compared to the typical times of evolution.

### C Numerical procedure

We have discretized Eq. (57) in a basis which is essentially the same that we used to solve the stationary GPE. To be more precise, the expansion is as follows

$$\stackrel{}{W}_i^{(n)}=\underset{k}{}a_k\left(\begin{array}{c}P_{kn}\\ 0\end{array}\right)+b_l\left(\begin{array}{c}0\\ P_{l,2mn}\end{array}\right).$$
(58)   
Here we have used the index convention explained above.

In this basis, the operator $`_{0n}`$ is diagonal, while the operator $`𝒰`$ can be calculated, either by means of integrals of the wavefunction itself in position representation, or by using a tensor of four indices which is similar to the one in Eq. (16). In any case, the equations are always linear, and so the study of the Bogoliubov spectrum consists in building and diagonalizing a large matrix of real numbers.

Even though the procedure is quite simple, the matrices that one must build in order to resolve the case of strong interaction are very large and tend to exhaust computational resources. To be able to reach a large value of the nonlinearity we have had to work in a subspace of states with even parity with respect to the Z axis. This way we could find the excitations with lowest energy for different vorticities, at the cost of missing those with odd parity, which are more energetic anyway .

### D Analytical results

The study explained above does not have to be performed for all of the Bogoliubov operators in all of the possible situations. Here we will show several important results regarding when Eq. (57) may imply destabilizing modes.

Lack of exponential instabilities in the Bogoliubov theory.- Any eigenvalue $`\lambda `$ satisfying Eq. (57) and $`G(\stackrel{}{W})0`$ must be real. Eigenstates with $`G(\stackrel{}{W})=0`$ may involve complex eigenvalues but they are spurious and are introduced by the linearization procedure.

This first part is shown simply by projecting the left and right hands of Eq. (57) against the vector $`\stackrel{}{W}_i^{(n)}`$. Omitting the indices the result is

$`\lambda _n`$ $`{\displaystyle (|u|^2|v|^2)}={\displaystyle (\overline{u}_nH^nu+\overline{v}H^{2mn}v)}`$ (59)   
$`+`$ $`{\displaystyle Uf^2|u+v|^2}{\displaystyle (nm)\mathrm{\Omega }(|u|^2+|v|^2)}.`$ (60)   
The second part is more subtle. To prove it we must remember that solutions to Eq. (27) are stationary points of the action , $`S=L(t)𝑑t`$ corresponding to the following Lagrangian density

$$L=\frac{i}{2}(\alpha \overline{\alpha }_t\overline{\alpha }\alpha _t)+_2(\alpha ).$$
(61)   
Using (44) it is easy to prove that the $`G(\stackrel{}{W})=0`$ modes are null modes that do not appear in the Lagrangian, and thus are not affected by the dynamics.

It must be remarked that this result characterizes possible eigenvalues, but does not grant that $`_n`$ have a Bogoliubov diagonalization.

Sufficient condition for stability.- If the linearized Hamiltonian $`_n`$ is positive definite, then $`_n`$ may be diagonalized, all of its eigenvalues are positive real numbers and there are no dynamical nor energetic instabilities.

To prove this theorem one only needs to show that there’s a one-to-one correspondence between the eigenfunctions of $`_n^{1/2}\sigma _z_n^{1/2}`$ and the eigenfunctions of $`\sigma _z_n`$ so that

$$_n^{1/2}\sigma _z_n^{1/2}|n=\lambda |n,$$
(62)   
if and only if

$$\sigma _z_n\left(_n^{(1/2)}|n\right)=\lambda \left(_n^{(1/2)}|n\right).$$
(63)   
Then one uses this result to show that the eigenvalue in (57) must be positive.

Stability in stationary traps.-  In Eq. (57), if $`\mathrm{\Omega }=0`$ and $`n>3m`$ then the linearized Hamiltonian $`_n`$ is positive, the Bogoliubov operator $`_n`$ can be diagonalized and it is also positive. Furthermore, if $`n>m`$ then any real eigenvalue is positive, $`\lambda >0`$.

The demonstration has four steps. First, one takes any value of $`n`$ that satisfies that condition and proves that $`H^{2mn}>H^m`$ and $`H^n>H^m0`$. Second, this is used to prove that $`_{0n}>_{0m}`$. Third, it shown that $`𝒰_n`$ is positive which altogether implies $`_n>0`$. The last assertion may be easily checked with the help of (59).

The preceding two theorems imply that in a stationary trap any mode with negative energy must be comprised in the $`(0,2m),\mathrm{},(m,m)`$ families, and any dynamic instability must lay in $`(0,2m),\mathrm{},(3m,0)`$. Thus we need only diagonalize a finite number of operators to make sure the system is stable or unstable. This result is an extension of the one obtained in Ref. , where a sufficient condition for stability is found to be $`n^24m^2`$, without taking into account possible complex eigenvalues.

Local stability under rotation.- In Eq. (57) the operator $`_n(\mathrm{\Omega })`$ exhibits a linear dependence with respect to $`\mathrm{\Omega }`$

$$_n(\mathrm{\Omega })=_n(0)(nm)\mathrm{\Omega }.$$
(64)   
While the wave functions of the modes are the same as the ones of the stationary traps, the energies of the excitations suffer a global shift that depends on the vorticity

$$\lambda (\mathrm{\Omega })=\lambda (0)(nm)\mathrm{\Omega }.$$
(65)   
In general, the influence of these shifts has to be checked numerically. It is easy to show, however, that the shift is positive for $`n<m`$, which means that the possibly negative eigenvalues in the range $`0<n<m`$ can be suppressed if $`\mathrm{\Omega }`$ is large enough. Even more, as the shift is a real number, if one demonstrates that there are no dynamical instabilities in the stationary trap, then there will be no dynamical instabilities in the rotating trap, neither.

### E Numerical results

Summing up, from a practical point of view, the issue of stability consists in two different steps. The first one is the search for a stationary solution of the GPE with the appropriate vorticity, which we have already performed in Sect. II, and the second one is the study of the spectrum of the Bogoliubov operators for this particular state.

Stability of the $`m=1`$ vortex-line in a stationary trap.- In this case of unit charge one only has to study a single operator, $`_0`$, to know whether the system is stable. This calculation provides us with the branch of the spectrum of excitations which is characterized by the quantum numbers $`(0,2)`$ and $`(2,0)`$, as already explained. We have done this for a wide range of nonlinearities in the absence of rotation, $`\mathrm{\Omega }=0`$, and the first conclusion is that the Bogoliubov operator has a diagonal Jordan form with all eigenvalues being real.

In Fig 4 we show a selected set of the eigenvalues of the Bogoliubov operator, both for a spherically symmetric trap and an elongated trap. In those pictures one sees several things. First, there are two constant eigenvalues $`\lambda =1`$ which correspond to oscillations of the vortex line along the Z axis. Second, there is a single neutral mode $`\lambda =0`$ only for the spherically symmetric trap, which corresponds to the symmetry of rotation of the condensate around an axis on the XY plane. The symmetry and the mode disappear when $`\gamma =2`$ \[See Fig. 4\]. And finally there is at least one negative eigenvalue $`\mu <0`$ (more in the case of an elongated trap) which is responsible for the energetic destabilization of the system. The largest contribution to this destabilizing mode is a wavefunction captured in the vortex line and has zero vorticity (i.e. it is a *core* mode) \[See Fig. 5\] as it was qualitatively predicted by Rokhsar in Ref. .

We must remark that the number of unstable modes increases with the geometry factor: the more elongated the trap is, the easier is to transfer energy from the vortex to the core plus longitudinal excitations. In other words, for $`\gamma 1`$ (spherical or “pancake” traps) there’s only one negative eigenvalue which corresponds to an excitation with a different vorticity than the unperturbed function, while for $`\gamma 1`$ we still have that mode, plus some more which are excited with respect to the Z axis. As a consequence, if the experiment is subject to dissipation and these unstable modes play a significant role in the dynamics, then the more elongated the trap is the less stable the vortex will be.

In Fig 6 we also show the lowest eigenvalues of the families $`(1,3)`$,$`(1,1)`$, $`(0,2)`$, $`(2,0)`$ and $`(2,4)`$, that is, excitations where the main contribution is an eigenstate of $`L_z`$ with eigenvalues $`m=0,\pm 1,\pm 2`$. In those pictures one sees that subspaces with excitations of the same vorticity but opposite sign have also different energy, a phenomenon which is solely due to the interaction.

Stability of $`m=1`$ vortex lines in rotating traps.- It was already proved (64) that, as the effect of rotation is gradually turned on, the modes with $`n<m`$ and with $`n>m`$ are shifted up and down in the spectrum, respectively. It remained the question of whether the shift is enough to stabilize the vortex states, and the answer is yes, according to numerical experiments.

First, as it is shown in Fig. 7, the negative eigenvalue is slightly smaller than the stabilizing frequency, $`|\lambda _0|<\mathrm{\Omega }_1`$, which implies that for $`\mathrm{\Omega }>\mathrm{\Omega }_1`$ the energetically unstable branch with vorticity $`m=0`$ disappears. And second, the eigenvalues of $`_n`$ for $`n>m`$ are found to be larger than $`(nm)\mathrm{\Omega }_1`$. In consequence for at least the interval $`[\mathrm{\Omega }_1,\mathrm{\Omega }_2)`$ all of the $`_n`$ operators are positive and the vortex with unit charge is a local minimum of the energy functional.

In any case the shifts are always real, which implies that the $`_n`$ operators remain diagonalizable with real eigenvalues and without dynamical instabilities.

Stability of the $`m=2`$ vortex line.- Another interesting configuration is the $`m=2`$ multicharged vortex line. Here one suspects that a configuration with several vortices of unit charge has less energy than a single multicharged vortex, under all circumstances. In other words, they must be always energetically unstable.

This intuitive perception is confirmed by the numerics. First the diagonalization of $`_1`$ reveals that this operator has at least one negative eigenvalue, while $`_0`$ has both negative eigenvalues and a pair of complex eigenvalues that, as we saw above, do not participate in the dynamics and must be discarded. Regarding the negative eigenvalues, they do not decrease with the nonlinearity, but are always larger in absolute value than their linear limits. This implies that there are always negative eigenvalues which cannot be suppressed with any rotation below the critical value, $`\mathrm{\Omega }_c>\mathrm{\Omega }>\mathrm{\Omega }_2`$.

The immediate consequence of this linear stability analysis is that, due to the linearization of the energy (41) not being positive, the $`m=2`$ vortex-line is never a local minimum of the energy. This is true even for the parameter interval, $`[\mathrm{\Omega }_2,\mathrm{\Omega }_3)`$, in which it has less energy than the rest of stationary states of well defined symmetry. If the $`m=2`$ ground state is not a minimum, and the other symmetric states have more energy, we can conclude that the minimum of the energy functional in the rotating trap with $`\mathrm{\Omega }[\mathrm{\Omega }_2,\mathrm{\Omega }_3)`$ must be a state which is not symmetric with respect to rotations . A similar analysis can be performed for the stationary states with $`m=3,4\mathrm{}`$ which extends this result to larger rotation frequencies, all below the critical one.

### F Lyapunov stability

Speaking roughly, a solution of Eq. (4) is Lyapunov stable when every perturbed solution which is close enough to the original wave remains close throughout the evolution. The concepts of Lyapunov stability and linear stability are close, but the latter does not imply the former as it is only defined in the limit of infinitesimal perturbations.

Studying the Lyapunov stability of Eq. (4) theoretically is a difficult task that should be subject of further investigation. In the mean time we have performed an “empirical” study of the Lyapunov stability of the stationary solutions with $`m=1`$ and $`m=2`$ vorticities, by simulating numerically how they evolve for small perturbations and long times. The simulation was performed with a three-dimensional split-step pseudospectral method like the one from Ref. , using a $`80\times 80\times 80`$ points grid to study both the $`\gamma =1`$ and $`\gamma =2`$ problems.

The main result of this complementary work is that both the unit charge vortex line and the multicharged vortex line are stable to perturbations which involve the destabilizing modes as defined by (57). For example, one may try to add a small contribution ($`0.5\%`$) of a core mode to the $`m=2`$ vortex, and with the result that the vortex line is split into two unit charge vortex lines, which rotate but remain close to the origin. We must remark that, although these simulations only work for finite times which are dictated by the precision of the scheme and the computational resources, these times are typically 20 or 30 periods of the trap, which is much larger than any of the magnitudes that one may address theoretically to the destabilization process (i.e. the negative or complex eigenvalues of Eq. (57))

In the end, what this type of simulations reveal is that the $`m=1`$ and $`m=2`$ stationary states are energetically unstable, but this has no influence on the dynamic unless some other “mixing” or dissipative terms participate in the model.

## IV Conclusions

We have studied the vortex solutions of a dilute, nonuniform Bose condensed gas as modeled by the Gross-Pitaevskii equation (4), both in a stationary, axially symmetric trap, and subject to rotation (or a uniform magnetic field).

First, we have searched solutions of Eq. (5) that have the lowest energy and which are also eigenstates of the third component of the angular momentum operator, $`\psi (r,z,\theta )=f(r,z)e^{im\theta }`$, both in a stationary trap and in a rotating trap, and from small to very large nonlinearities. It has been found that a nonzero angular speed (or magnetic field) is necessary in order to turn a vortex line state into a minimum of the energy functional with respect to other states of well-defined vorticity. However it remains open the question of whether the minimum of energy must have a well defined vorticity.

Next we have studied the stability of these stationary solutions of the GPE. We have formulated a set of coupled equations that describe both the linearization of the GPE around a stationary solution and Bogoliubov’s corrections to the mean field theory that describes the condensate. It has been proved that the problem may not exhibit dynamical instabilities of exponential nature, plus several other theorems that describe the phenomenology associated to the possible instabilities.

The perturbative equations have been solved numerically for stationary states having $`m=1`$ and $`m=2`$ vorticities. In both cases it has been found that the only instability is of energetic nature, being limited to a small number of modes whose nature had already been predicted in .

For the vortex with unit charge we have found that this instability may be suppressed by rotating the trap at a suitable speed, and even when the trap is stationary, it is expected that it plays no significant role in the dynamics unless there is enough dissipation as to take the system through the unstable branch. On the other hand, the linear stability analysis for the $`m=2`$ multicharged vortex reveals that the energetic instability may never be suppressed, and that this configuration is never a minimum of the energy functional, even though its lifetime is, once more, only conditioned by possible dissipation.

The last and probably most important conclusion of this work is that in the rotating trap, and for $`\mathrm{\Omega }>\mathrm{\Omega }_2`$, the state of minimum energy is not an eigenstate of the $`L_z`$ operator, and thus it is not symmetric with respect to rotations. A similar result has been found in Ref.  by means of a minimization procedure which is only justified in the limit of very small $`U`$, while our demonstration remains valid for all nonlinearities, as far as the linearization procedure may be carried on.

From a experimental point of view, this work has several implications. First, it is clear the conclusion that vortex lines with unit charge may be produced by rotating the trap at a suitable speed and then cooling the gas. Second, once rotation is removed, these vortices will survive for a long time if dissipation is small. Third, the multi-charged vortices are not minimum of the energy functional and thus it will be difficult to produce them by mean of cooling a rotating gas. And finally, if these multi-charged vortices are produced by some other mean such as Quantum Engineering, then we can assure that their lifetime will only depend on the intensity of dissipation, whose effect is to take the system either to the $`m=0`$ ground state if $`\mathrm{\Omega }<\mathrm{\Omega }_1`$, to the unit charge vortex-line state if $`\mathrm{\Omega }<\mathrm{\Omega }_2`$, or to a symmetry-less multicharged state if $`\mathrm{\Omega }>\mathrm{\Omega }_2`$ (A phenomenon which is regarded as splitting in the literature).

The numerical results found in the paper have been possible due to the use of a powerful Galerkin spectral method optimized to allow the consideration of thousands of modes which is a step forward with respect to the previous analysis.

###### Acknowledgements.

This work has been partially supported by DGICyT under grants PB96-0534 and PB95-0389 We thank Prof. Cirac from the Institud für Theoretische Physik of Innsbruck for proposing us the problem and helping in the development of the stability theorems.