Spectral stability of vortices in two-dimensional Bose-Einstein condensates via the Evans function and Krein signature

We investigate spectral stability of vortex solutions of the Gross-Pitaevskii equation, a mean-field approximation for Bose-Einstein condensates (BEC) in an effectively two-dimensional axisymmetric harmonic trap. We study eigenvalues of the linearization both rigorously and through computation of the Evans function, a sensitive and robust technique whose use we justify mathematically. The absence of unstable eigenvalues is justified a posteriori through use of the Krein signature of purely imaginary eigenvalues, which also can be used to significantly reduce computational effort. In particular, we prove general basic continuation results on Krein signature for finite systems of eigenvalues in infinite-dimensional problems.


Introduction
Background. Since the experimental creation of Bose-Einstein Condensates (BEC) in alkali vapors in 1995 [5,13], BEC are one of the most active areas of modern condensed-matter physics. A general overview of the subject can be found in [12,58], and particularly in the review book [42]. In the Hartree-Fock mean-field approximation, BEC are modeled by the nonlinear Schrödinger equation (NLS) with a non-local nonlinearity. A traditional simplification, replacing the non-local interaction potential with a localized short-range interaction proportional to the delta function, leads to the Gross-Pitaevskii equation where is Planck's constant, M is the atomic mass of atoms in the condensate, θ is an azimuthal angle in cylindrical coordinates, and g is an interaction strength parameter. The total number of particles N in the condensate is given by the integral and is conserved during the evolution of the system. Equation (1) is a nonlinear Schrödinger equation with a cubic nonlinearity (focusing or defocusing depending on whether the interaction is attractive or repulsive, respectively) and with a spatially dependent trap potential V (x) stationary in a frame rotating with frequency Ω about the vertical axis. A rigorous mathematical justification of the Gross-Pitaevskii model for the BEC ground state under various conditions directly from many-body Schrödinger equations was done in a series of papers of Lieb et al. [45]- [49]. From the point of view of nonlinear waves, the interesting phenomena is that the Gross-Pitaevskii equation, similarly to some other nonlinear Schrödinger equations, supports the existence of various types of solitary wave solutions. In the two-dimensional setting we will study in particular, there are vortex solutions which have the form ψ(t, r, θ) = e −iµt e imθ w(r), where r, θ are polar coordinates, m is vortex degree, µ is the vortex rotation frequency (physically, chemical potential) and w(r) is the radial vortex profile.
Problems of stability for vortex solutions to various forms of nonlinear Schrödinger equations have drawn much attention in recent years. Related questions for models arising from nonlinear optics, micromagnetics and Bose-Einstein condensation have been considered extensively in the mathematical and physical literature. For recent work concerning spectral stability questions for various types of matter-waves, including vortices, vortex rings, multi-poles, soliton and vortex necklaces in the presence of magnetic traps and optical lattices, see for example [36,37,38,42,43]. Rigorous mathematical results on these questions are rather few, however, due to the strong nonlinearity and complexity of the system.
This work. In the present work we study the spectral stability of a single two-dimensional axisymmetric vortex trapped in an axisymmetric harmonic trap. For this simple well-studied physical setting, we develop an approach that involves a combination of analytical and numerical tools which allows us to obtain reliable results for large particle number, well into the Thomas-Fermi regime (e.g. N ∼ 10 6 atoms of 23 Na). It is important to note that in the present axisymmetric setting, rotation of the trap does not influence the dynamic stability of vortices, as the rotation term can be removed by transforming to an appropriately rotating coordinate frame.
Analytical results. Our study involves a number of analytical results that we prove in sections 5 and 6 concerning the spectrum of the operator one obtains by linearizing about a vortex solution. In section 5, we prove that due to the harmonic trapping potential the essential spectrum of this operator is emptyhence the spectrum consists entirely of isolated eigenvalues of finite multiplicity. The real part of any eigenvalue satisfies an explicit bound depending only on the (non-dimensionalized) frequence µ and degree m of the vortex, namely | Re λ| < 3(µ − m) . ( The eigenvalue problem breaks into an infinite system of coupled pairs of ordinary differential equations (ODEs) for azimuthal Fourier modes indexed by j ∈ Z. As is known from [18], only finitely many of these are relevant for possible instability, namely the ones satisfying 0 < |j| < 2m.
We construct a globally analytic Evans function [3,15] associated with each of the ODE pairs, whose zeros correspond to the eigenvalues. These results extend the approach of [54] for focusing-defocusing nonlinear Schrödinger equations to handle a harmonic trapping potential. One of the weaknesses in the numerical investigations in [54] was that one could not be confident that all unstable eigenvalues were detected, due to the absence of any bound on their imaginary part. No such bound is available in the present problem either. Nevertheless, we will show how one can indeed account for all possible instabilities through use of Krein signature (the sign of the linearized energy of associated eigenmodes).
The key property of Krein signature that makes it useful is its invariance under continuous variation of parameters such as the standing-wave frequency µ and the size of the condensate. In section 6 we extend results that are well-known for finite-dimensional systems to establish a general continuation property for any family of operators that is resolvent-continuous, a natural notion of weak continuity that one expects to hold in many applications to infinite-dimensional systems. Preservation of a definite signature is proved for any finite system of imaginary eigenvalues for such a family. By consequence, the only way eigenvalues can leave the imaginary axis is through collision with eigenvalues of opposite Krein signature.
There are three reasons why Krein signature is extremely helpful and a powerful tool in the present work. First, it allows us to explain the collisions (or avoided collisions) of eigenvalues found in our numerical computations. Moreover, in combination with numerical plots, it provides a numerically convincing a posteriori justification that there are no unstable eigenvalues outside a certain fixed box in the complex plane. Finally, it can be used to significantly reduce the amount of necessary computation, as we shall discuss in section 7.2.
Numerical methods. In order to locate zeros of the Evans functions using the argument principle, we design and implement a numerical method somewhat more robust than that used in [54]. We use a path-following technique and multiple shooting to compute the nonlinear wave profile, and a rescaled exteriorproduct representation of the Evans function to handle problems of rapid growth and numerical dependence. A numerical path-following technique for radially symmetric profiles was also used by Edwards et al. [14].
Numerical results. In agreement with previous studies in the physical literature [18,37,60,68], we find a singly-quantized vortex (m = 1) spectrally stable while the stability of multiply-quantized vortices (with m = 2 and 3) depends on the diluteness of the condensate, with alternating intervals of stability and instability as N g varies. Pu et al. [60] in their analysis use a different numerical method (finite elements) which is a priori less reliable and sensitive than the approach used here. Moreover, our results account for the appearance of all instabilities through the tracking of all eigenvalues of negative Krein signature. The presence of unstable eigenvalues (eigenvalues with positive real part) for certain parameters is also corroborated by direct simulation of time-dependent dynamics based on the splitting scheme proposed in [6].
Symmetries. In the computations one observes a special set of eigenvalues which remain constant under variation of standing-wave frequency and condensate size. These eigenvalues arise from the symmetries of the Gross-Pitaevskii equation, as we discuss in an appendix. One symmetry is particularly interesting -a breather boost, which we found first described in [59]. This is self-transformation of the Gross-Pitaevskii equation related to the Talanov lens transformation, involving a time-periodic dilation of space with an appropriate radial phase adjustment. This symmetry corresponds to eigenvalues ±2iω of the Gross-Pitaevskii equation with the harmonic potential V (x) = 1 2 ω 2 |x| 2 , linearized about a central vortex, with mode shapes corresponding to infinitesimal breathing oscillations.
Related literature. Let us now discuss some known results on stability which are relevant to our problem. In general, two different concepts of stability are distinguished in the literature: energetic and dynamic stability. A solution is energetically stable if it minimizes an associated energy functional within a certain class of functions.
The simplest energetic stability approach, where the minimization takes into account only single vortex solutions with different charges, indicates that a high enough trap rotation frequency can eventually stabilize a vortex of any degree [8]. On the other hand, without external trap rotation, the energy of a single multi-quantized vortex of charge m is larger than the energy of m singlyquantized vortices, and thus multi-quantum vortices are believed to be unstable. The total energy in this case also depends on the relative location of vortices as they tend to form regular hexagonal arrays in harmonic traps.
A mathematical framework for a rigorous variational approach was discussed by Aftalion and Du [2]. Their method for effectively 2D condensates is parallel to the Ginzburg-Landau theory of superconductors. In [66] the authors claim that sufficiently fast rotation in combination with a strong pinning potential is capable of making even multi-quantum vortices energetically stable.
A detailed rigorous analysis was conducted by Seiringer [65], who studied regimes when a vortex solution can be a global energy minimizer. He proves that for any 0 < Ω < Ω c there exists m Ω (independent of an interaction potential) such that all vortices with charge m > m Ω are energetically unstable; i.e., they are not global minimizers (ground states) of the energy functional (see also [30]). Moreover, he proves that all multi-quantized vortices, m ≥ 2, become energetically unstable for a large enough value of the chemical potential of the condensate. Finally, he proves that symmetry breaking of the axisymmetric vortex solution is inevitable for any m, even for a singly-quantized vortex for a large enough interaction strength, since no ground state is an eigenfunction of the angular momentum. The symmetry breaking of a one-dimensional ground state is also demonstrated by a dynamical systems analysis in [33] in the case of a double-well trapping potential.
Energetic stability provides a sufficient condition for dynamic stabilityground states of the energy functional are nonlinearly orbitally Lyapunov stable, i.e., if the initial data are "close" to the ground state solution then the perturbed solution remains "close" to the ground state solution for all times. (See [33] for a detailed dynamic stability study of a one-dimensional model and a sketch of the proof of well-posedness for the initial-value problem for the Gross-Pitaevskii equation.) As dynamic stability need not imply energetic stability, however, it may not be possible to draw conclusions on dynamic instability directly from considerations of the energy functional. The study of linear or spectral stability can be helpful since one may detect possible instabilities due to the presence of eigenvalues in the right half-plane.
In the physical literature, García-Ripoll and Pérez-García [18] and Pu et al. [60] have studied linear stability of single multi-quantized vortices using equations equivalent to those here. The numerical results of [60] agree substantially with those of the present paper. In [18] the search for instability is restricted only to the so-called anomalous modes -modes with a negative linearized energy. As pointed out in [31] for example, these modes are not intrinsically unstable in the sense that some dissipation mechanism must be introduced into the system for them to become relevant. The numerical techniques used in [18] and [60] rely on Galerkin-type approximations. A finite-temperature generalization is further studied in [71].
Use of the Krein signature as a tool to study stability of nonlinear waves has appeared recently in various studies, see [9,26,34,39,56] and references therein. Skryabin in [67] (also see [68]) studies a binary mixture of trapped condensates using such information. Kapitula et al. [37] study stability of various types of matter-waves including localized vortices. Their perturbation argument, combined with topological methods based on Krein signature, describes in detail transition to instability in the limit of weak atomic interactions. Finally, Kapitula and Kevrekidis [35,36] have studied Bose-Einstein condensates in the presence of a magnetic trap and optical lattice, making efficient use of information on the Krein signature of relevant eigenvalues.
Organization of this paper. First, Sections 2 and 3 introduce notation and re-call background results regarding the Gross-Pitaevskii equation and vortex solutions. Section 4 contains most of our analytical results, concerning linearization, essential spectrum, bounds on eigenvalues, and reduction to ODEs. Moreover, we establish a precise asymptotic description of eigenfunctions necessary for construction of the Evans function. The Evans function itself is constructed in Section 5. In Section 6, we discuss the Krein signature. We describe in detail our numerical procedure and discuss the numerical results in Section 7. Finally, in an Appendix we discuss symmetries and boosts of the problem and relate them to the eigenvalues which do not change as the standing wave frequency varies.

The Gross-Pitaevskii Equation
The behavior of low-temperature Bose-Einstein condensates (BEC) trapped in the harmonic potential V (x) rotating with the angular velocity Ω about the z-axis is well described by the time-dependent Gross-Pitaevskii equation. The wave function ψ(x, t) satisfies (1) (in three dimensions) with trapping potential V (x) = V (x, y, z) given by The interaction strength parameter g is where a is the s-wave scattering length [45]. The term Ω∂ θ corresponds to the angular momentum Ω · (r × ∇) of the condensate caused by a rotating frame of coordinates. The total number of particles in the condensate N is given by the integral and is conserved during the evolution of the system. For disc-shaped (pancake) traps (ω 2 z ≫ ω 2 x , ω 2 z ≫ ω 2 y ) it was justified [6] that the system is well approximated by a planar two-dimensional reduced model. The equation (1) formally does not change; one only needs to set where ω x = ω tr and ω y = ω tr λ tr . For purpose of numerical investigations in this work the same values of parameters were used as in [60]: a condensate consisting of atoms of 23 Na is considered with a = 2.75 nm, ω z = 2π × 200 Hz, ω tr = 2π × 10 Hz (ω z ≫ ω tr ), M = 10 −26 kg and the Planck constant = 6.6261 × 10 −34 Js. A similar set of parameters used in the experiment as cited in [32,66] with a 87 Rb condensate is: M = 3.81 × 10 −26 kg, a = 5.77 nm, ω = ω z = ω tr = 2π × 200 Hz. In these experiments the number of particles was approximately N = 2 × 10 5 , horizontal and vertical condensate sizes were R = 20 µm, L = 10 µm and the temperature T c = 1 µK. Both 23 Na and 87 Rb represent alkali gases with a repulsive interaction potential. As an example of an attractive interaction, 7 Li with a = −1.45 nm can serve. Note that the parameter a can be tuned via Feshbach resonance. In this work the following assumptions will be made. We assume that the magnetic trap is axisymmetric (λ tr = 1). Thus only two-dimensional wave functions of the form ψ = ψ(t, r, θ) will be considered. The symmetry of the trap also eliminates dependence of stability of axisymmetric vortices on the trap rotation. Hence, we will set Ω = 0. Moreover, although the model includes both attractive and repulsive interaction interparticle potential (sign of the nonlinear term), for simplicity only the more interesting case of repulsive potential will be considered here (some results concerning stability in transition between repulsive and attractive potential can be found in [60]).
To nondimensionalize the equation (1) we use the following scaling The Gross-Pitaevskii equation is then expressed (dropping the primes) as The energy functional is given by Note that the Thomas-Fermi regime [16,45] N a/d 0 >> 1 (here d 0 is the mean oscillator length, d 0 = /M ω 0 and ω 0 is the mean trap frequency ω 3 0 = ω x ω y ω z ) corresponds to K → ∞ since K = 2|a|N 2πM ω z / , i.e., This is the limit under which Lieb and Seiringer [45] justified the Gross-Pitaevskii energy functional to be a good approximation for the N -body quantum system. Note, that the only free parameter which stays in (6)-(7) is K, the L 2 -norm of the wave function ψ.

3 Vortex Solutions
In this section we describe the structure of vortex solutions [51] to the Gross-Pitaevskii equation (6), of the form and describe the numerical procedure which allows us to approximate them with high precision. Here m is an integer, and in this paper it suffices to always assume that m ≥ 1 due to reflection symmetry. The radial profile function w(r) of a vortex solution satisfies the equation We will require that the radial profiles be non-negative and spatially localized, i.e., satisfy the boundary conditions w(r) is bounded as r → 0 + , w(r) → 0 + as r → ∞.
We note that the boundary condition as r → 0 + implies w(r) → 0 + for m ≥ 1.
Existence. The existence of positive solutions to (11) for some µ, corresponding to any given K > 0 in (7), can be proved using a well-known variational argument [65]. One minimizes (8) among functions with the spatial dependence in (10), with m and K fixed, and a minimizer may be found that is positive. For any positive solution of (11) with finite energy, necessarily µ > m + 1, since multiplying (11) by 2πrw and integrating in r yields µK > π ∞ 0 w 2 r + m 2 r 2 + r 2 w 2 r dr ≥ (m + 1)K.
The last inequality follows since the integral is minimized at a positive solution of which is (11) linearized at zero. Analysis of this equation (see below) yieldŝ µ = m + 1, w = cw (m) 0 from (18) below, where c is constant. The following proposition describes a global bound on any vortex solution with positive profile. A proof can be found in [44] (except for the statement that µ > m + 1). The bound (13) was also proved in [65]. Proposition 1. Let w(r) be a finite-energy positive solution to (11), satisfying (12), where m > 0 is an integer. Then µ > m + 1, and w(r) is increasing on (0, R) and decreasing on (R, ∞), for some R ∈ (m/ √ 2µ, √ 2µ). Moreover, for all r > 0 we have |w(r)| 2 < µ − m.  The asymptotic behavior of a vortex profile can be determined directly from (11) (see Fig. 1). It is not difficult to see that as r → 0 + the equation has the same character as the linear Schrödinger equation (14) (one argues as in [29]), and w(r) ∼ d 0 r m for m ≥ 1, for some positive constant d 0 . As r → ∞, the nonlinear term for a localized solution becomes negligible and the linear equation (14) is again a good approximation. The proof that the positive solution w(r) to (11) approaching 0 as r → ∞ satisfies w(r) = O(r µ−1 e −r 2 /2 ) is also given in detail in [44]. The goal of this paper is to study spectral stability of solutions to (11) both analytically and numerically. Naturally, for a careful numerical stability investigation it is crucial to obtain very precise numerical solutions of (11) first. The approach used here is based on path-following along a branch bifurcating out of the trivial solution w = 0 and is similar to the one used in [14].
For µ given by (17) both solutions w (1) (r) and w (2) (r) reduce to a constant multiple of the single solution where L  The numerical algorithm designed to find solutions to (11) is based on the following observation. It is reasonable to expect that introduction of the nonlinear term leads to the existence of a solution branch (µ(s), w µ (s)) bifurcating from the trivial solution w = 0 for µ = µ 0 . To justify such a behavior one can use the Crandall-Rabinowitz theorem [11,53], to prove the following Theorem (for details see [44]).
Note that in addition to a branch of positive solutions (µ, w) bifurcating from the trivial solution at µ = µ 0 = m + 1 in a direction (0, w Numerics. Hence it is possible to numerically trace the solution curve (µ, w) from the branching point (µ 0 , 0), see Fig. 2(a). The typical radial profile for m = 2 is on Fig. 1. Note that the Crandall-Rabinowitz bifurcation theory also provides some linear stability information for solution branches given in Theorem 1 [11,53], but one must take the stability results with caution. Although the theory predicts stability for the bifurcating branch, it is only with respect to radially symmetric perturbations in (14). This is not sufficient to determine stability of vortex solutions with respect to the full dynamics in (6).
The behavior of norms relative to the norm of the parabolic Thomas-Fermi regime approximation is illustrated in Fig. 2(b). The first point on the approximate solution curve is set to be (µ 0 , εw (m) 0 ), ε = 0.1, which is an O(ε 2 ) approximation of the exact solution. Then an implementation of a predictorcorrector path-following algorithm [4] is used to get solutions for large values of the parameter µ.
Since our stability study will require evaluations of the profile at any given point within the computational domain, the precision of calculations is improved by optimizing the already calculated profile for any given µ by a multiple shooting procedure [69]. This allows one to achieve high precision in evaluating w(r) by simple integration from a nearby mesh point. Note that the calculation is almost independent of the size of the parameter µ since the size of the computational domain, and so the number of necessary nodes, grows very slowly. Therefore it is possible to reach large values of µ. Also note, that with the growing parameter µ, the L 2 -norm (9) of profiles grows (Fig. 2) and hence states far in the physically interesting Thomas-Fermi regime for a wide range of µ's are obtained for a small computational cost. On the other hand, as pointed in [14], for computation for a single value of µ this method has significant overhead. linearized equations have the same form as the so-called Bogoliubov equations [14,58] commonly used in physics literature. Note, however, that the relation between the derivation here and the physical derivation of the Bogoliubov equations is by no means straightforward.
Neglecting nonlinear terms in (6) yields The complex character of the equation (19) complicates the analysis. Therefore, we decompose the complex wave function v as The equation (19) is then equivalent to the real system where To understand the dynamic stability of the vortex solution ψ, we will study the spectrum of the operator A as an unbounded operator on L 2 (R 2 , R 2 ) with an appropriate domain D(A). As is rather well-known, spectral stability of A (meaning absence of spectrum in the right half of the complex plane), need not necessarily imply linear stability (in the sense that the zero solution of (20) is stable), nor nonlinear stability of vortices. In this paper, we avoid these subtle issues and confine ourselves to studying spectral stability. After all, the presence or absence of eigenvalues with positive real part is interesting in itself.
The precise definition of the operator A is somewhat involved and requires the concept of a quadratic form [62]. Write (only formally for now) where I is the 2 × 2 identity matrix. Then Then define a quadratic map The quadratic form q Lc is semibounded: The quadratic form q Lc is closed if it has a closed graph, i.e., if D(q) is complete under the graph norm Ψ +1 = q Lc (Ψ, Ψ) + Ψ 2 L 2 . This is true since both H 1 and L 2 (r 2 ) can be obtained from L 2 by completing the space of C ∞ 0 functions under the H 1 and L 2 (r 2 ) norms respectively. Also observe that for some C > c > 0. (The lower bound follows from the semiboudedness; the proof of the upper bound is analogous.) Then Theorem VIII.15, pp. 278 of [62] yields that q Lc is the quadratic form of a unique self-adjoint operator L c . The domain of the operator L c , denoted by D(L c ), is dense in L 2 . Clearly for Ψ ∈ D(L c ) in the sense of distributions. It is easy to see that the operator L w in (21) is bounded on L 2 (R 2 , R 2 ) Therefore the operator A = −J(L c + L w ) with the domain D(A) = D(L c ) is closed and densely defined in L 2 (R 2 , R 2 ).

Essential spectrum
We investigate the spectrum σ(A) of the operator A given by (22), regarded as an unbounded operator on complexified space D(A) ⊂ L 2 (R 2 , C 2 ). The spectrum of such an operator in general consists of two parts: isolated eigenvalues of finite multiplicity form the discrete spectrum σ disc (A), and the remaining part -the essential spectrum σ ess (A). The latter is empty as stated in the next proposition.
Proposition 2. For any m and µ the spectrum of A consists entirely of isolated eigenvalues of finite multiplicity. I.e., the essential spectrum of the operator A is empty.
Proof. The proof has three steps. First, it is easy to see that the essential spectrum of L c is empty by Theorem XIII.16 on p. 120 of Reed and Simon [63]. Then we prove the same for JL c . Finally, the generalized Weyl theorem for non-self adjoint operators yields the same for JL.
Let us prove that the essential spectrum of JL c is empty. Since L c is a nonnegative operator, 0 is not an eigenvalue and therefore L c has bounded inverse. Moreover, L c only has discrete spectrum and its eigenvalues are isolated with the only possible accumulation point ∞. Then by Theorem XIII.64 on p. 245 of [63] the operator L −1 c is compact. Now consider the following identity: If λ / ∈ σ disc (JL c ), then 1 λ / ∈ σ disc (L −1 c J −1 ) and the right-hand side of (23) has bounded inverse given by It remains to prove that the eigenvalues of JL c are isolated and of finite multiplicity, which prohibits discrete spectra to be embedded in the essential spectrum. To show that, consider the resolvent equation The operator λL −1 c J −1 − I is a (multiple of) compact perturbation of identity and therefore it is also Fredholm. Also, it is analytic everywhere except for the discrete spectra of L −1 c J −1 . By a general result of Gohberg and Krein [22], p.21. or Kato [41], p.370, the set of values for which I −T (λ) is not invertible is at most countable with their only possible accumulation point infinity. Therefore the eigenvalues of JL c are isolated. Also, the spectral projection on the eigenspace associated with a particular eigenvalue of JL c has finite dimensional range, since it is given by an integral of a compact operator. Hence the essential spectrum of JL c is empty and consist of isolated eigenvalues of finite multiplicity with only accumulation point infinity, i.e., σ ess (JL c ) = ∅ .

15
Finally, we use the generalization of Weyl's theorem to non-self-adjoint operators to prove σ ess (JL c ) = σ ess (J(L c + L w )) .
It is enough to observe that J(L c + L w ) is a relatively compact perturbation of JL c , i.e., that JL w (λI − JL c ) −1 is compact whenever λ / ∈ σ disc (JL c ). This is true since JL w is bounded and (λI − JL c ) −1 is compact.
Let us point out that although the stability of an axisymmetric vortex in an axisymmetric trap does not depend on the trap rotation frequency Ω, it has an influence on the linearized operator A. The definition of the operator and the proof of emptiness of the essential spectra can be adjusted to account for the rotation as long as |Ω| ≤ ω tr . Beyond this threshold it is not clear how to define the operator and whether its essential spectrum stays empty in this parameter regime. This threshold may be significant if the axial symmetry is broken. A further discussion on other features in this regime can be found in Chapter 12 [42]. On the other hand, it is easy to see that the eigenvalues (the discrete spectrum) of the linearized problem suffer only a shift by a purely imaginary number depending on rotation, so stability of this part of the spectrum of A is unaffected by rotation of the trap.

Eigenvalues
By the result above, the investigation of the spectrum of the operator A is reduced to study of the eigenvalue equation which can be rewritten as Similarly as in [54] it is useful to represent solutions of (25) in the basis of the eigenvectors of the matrix J: Using the information on asymptotic decay and a simple bootstrap argument one can deduce that Φ ± ∈ H k loc for each k > 0, so Φ ± ∈ C ∞ (R 2 , C) ∩ L 2 (R 2 , C).
Furthermore, decompose Φ ± (r, θ) at each fixed r into Fourier modes (with shifted indices for notational ease) After the introduction of the Fourier modes the equation (25) transforms to an infinite-dimensional system of linear equations.
The system decouples to coupled pairs for nodes y + = y By the symbol △ r we denote the radial Laplace operator △ r = ∂ 2 ∂r 2 + 1 r ∂ ∂r . The proper boundary conditions for this system must be determined only from the system itself and the property Φ ± ∈ L 2 (R 2 , C). Here asymptotic behavior of solutions to (29)-(30) described in Theorem 2 of the next Section can be used. It implies that the appropriate boundary conditions are lim r→0 + y ± (r) exists and lim r→∞ y ± (r) = 0 .
Therefore the eigenproblem for A is decomposed into countable many problems L j y = iλRy , y = (y + , y − ) T , where Similarly as before the bootstrap argument gives y ∈ C ∞ ∩ L 2 (R + , C 2 ; r). The associated inner product is given by On the other hand, one can argue (as in [54]) that any solution (λ, y, j, m, µ) to (32) defines a solution (λ, Φ, m, µ) to (24) with Φ ∈ L 2 (R 2 , C 2 ). Note that L ± j = L ∓ −j , a fact connected to the Hamiltonian symmetries of A. The next proposition analogous to [54] summarizes the results of this section. If (y + , y − ) form a solution to (32) for some pair (λ, j), then (y − , y + ) form a solution for (−λ, −j) and (ȳ − ,ȳ + ) form a solution for (λ, −j).

Bounds on unstable eigenvalues
At a first glance it may seem impossible to solve infinitely many systems of the form (32). Fortunately, similarly as in [54] it is possible to restrict the index j for which an unstable eigenvalue may occur. The proof of the following Proposition can be found in [18] but for the sake of completeness and clarity of exposition we provide it here as well.
Proof. First we prove that |j| ≤ 2m. The main idea of the proof of this part is the same as in [18] but we present it here for clarity. Assume the contrary, i.e. |j| > 2m, namely, |j − m| > m and |j + m| > m. The operator L j is analogous to the operator L c in Section 4 and thus the integrals used here are well defined. Then for any y = 0 by integration by parts Here we also used the inequality The function w(r) is a non-negative minimizer of the linearized energy (36) within the family of functions f (r) with φ(r, θ) = e imθ f (r) ∈ D(L j ). Since E w (w) = 0 it follows that E w (f ) ≥ 0 for all e imθ f ∈ D(L j ). Hence (L j y, y) > 0 .
Therefore iλ must be real if |j| > 2m. This result can be also interpreted in terms of Krein signature (see Section 6): the signature of all eigenvalues λ for |j| > 2m is positive.
Next, we prove the stronger result that |j| < 2m and j = 0. First, let us consider the case j = 0. In that case one obtains directly The equality is valid only if there is equality in (35), i.e., if y + = −y − . Therefore (L j y, y) = 0 only if y + = −y − almost everywhere and E w (y + ) = 0. Since f (r) = w(r) is a non-negative solution of the differential equation associated with (36) This result can be also interpreted in terms of Krein signature (see Section 6): the signature of all eigenvalues λ for |j| > 2m is positive.
One can also prove the following estimate which restricts the possible unstable eigenvalues to lie in a vertical strip.
Proposition 5. The real part of every eigenvalue of the operator A is bounded: Proof. First, recall the simple bootstrap argument justifying that Φ ± of (26)- . Then, split the operator A to a vortex-profiledependent part A w and independent part A c : Since J and R are constant matrices (bounded by 1 in the matrix norm) the estimate (15) implies that the norm of the profile dependent part A w is bounded above: where · L 2 denotes the operator norm in L 2 (R 2 , C 2 ) and M (w) = max Multiply (25) by the smooth complex conjugate Φ and integrate over R 2 to obtain The second term on the right hand side of (39) can be estimated using (38) 2π Finally, the real part of 2π 0 ∞ 0 A c Φ · Φ r drdθ vanishes, which can be checked by a simple but long integration by parts (omitted here). The statement of the proposition then immediately follows by Proposition 1.
Propositions 4 and 5 restrict the search for unstable eigenvalues to a finite number of equations (with indices 0 < j < 2m) and to a vertical strip. Since one can expect infinitely many stable eigenvalues in both directions on the imaginary axis it is hopeless to prove that the imaginary part of an eigenvalue is bounded. Nevertheless, the possible number of unstable eigenvalues is limited (see Section 6).

Evans Function
While the finite element and the Galerkin approximation methods provide a fast and simple way to find eigenvalues of the problem (32), the Evans function technique [3,15] method has proved to be the most reliable and robust in certain cases. This approach will be implemented here. It is parallel to [54], where the reader can find many details of the procedure.
The main idea of this approach is to identify eigenvalues of the operator L j as zeros of an analytic function E j (λ). First, write the system (29)-(30) as a 4 × 4 system of first order ordinary differential equations: where 20 and The coefficients k + and k − are given by The asymptotic behavior of solutions to (41) is described in the next theorem.
Theorem 2. For any λ ∈ C and m > 0, µ real, j integer, there exist solutions y as r → ∞ .
The full proof of the theorem which relies on the asymptotic theory of Coppel [10] can be found in [44].
The asymptotic analysis reveals that (41) has two exponentially growing solutions asymptotically equivalent to y (r) for r ≫ 1. Note that these particular solutions are not in any way unique. From now on the notation y will always refer to solutions with the given asymptotics. The two-dimensional growing and decaying subspaces non-trivially intersect only if λ is an eigenvalue. Their intersection can be detected by vanishing of the Wronskian W (r, λ) = det(y ). By Abel's formula, this determinant satisfies a differential equation W ′ (r) = Tr (B) W (r). It is convenient to remove the dependence on r by setting This Evans function is then independent of r. It is evident that E j (λ) = 0 is a necessary and sufficient condition for the existence of an eigenvalue. A different representation of the Evans function is more convenient for analytic and numerical use, however, to avoid problems such as maintaining linear independence of solutions for extreme values of parameters and at mode collisions. It is also desirable to ensure global analyticity of the Evans function so that the presence and location of eigenvalues may be studied by means of contour integrals via the argument principle, An alternative way to construct and evaluate the Evans function involves introducing the adjoint system z ′ = −zB(r, j, λ) .
The fundamental matrices Y (z) (its columns are y i ) and Z(z) (with rows z i ) of systems (41) and (44) are related by ZY = I. Therefore Theorem 2 (by a simple direct calculation of the inverse matrix) also guarantees existence of four independent solutions of (44) as z One can also easily deduce the asymptotic behavior of z (∞) i as r → ∞. Furthermore, a simple calculation [54] shows that Constructing E j (λ) in this way still involves maintaining the independence of particular solutions y , however. A quite simple idea to overcome this difficulty has been used in a number of earlier works including [54]. Instead of considering the system (41) one can construct a larger 6 × 6 system for exterior products of solutions to (41). This exterior system has a unique solution of maximum growth rate given byŷ . Corresponding statements hold for the adjoint system. Evaluation of the Evans function is then given by the simple formula It is important to realize that the Evans function given by (46) is analytic in C and it is solution-and spatially-independent. On the other hand, the values are the same as given by (43) and (45). These Evans functions have the following symmetries. The proof is the same as in [54].

Krein Signature
In this section we establish some general basic continuation results on Krein signature for finite systems of eigenvalues in infinite-dimensional problems, and identify the eigenvalues of negative signature for the linearized Gross-Pitaevskii equation (14).
A typical context [24,39] in which Krein signature arises is in study of linearized Hamiltonian systems v t = JLv (47) on a Hilbert space X with inner product (·, ·), where v ∈ X, L : X → X is symmetric, (Lu, v) = (u, Lv), and J : X → X is invertible and skew-symmetric, (Ju, v) = −(u, Jv). We are interested in cases when L is unbounded and not a positive operator and has a finite number of negative eigenvalues. In particular, we will apply the results of this section to the operators L = L j , J =Ĵ = −iR, see (32). Note that individual operators L j that we consider here do not have all the symmetries of the Hamiltonian operator L = L c + L w of (21). To define Krein signature of a discrete eigenvalue λ of JL (or more generally, a finite collection of such eigenvalues), consider the restriction of the energy form (L·, ·) to the associated generalized eigenspace. If this restriction is positive or negative definite, the Krein signature of the eigenvalue is positive or negative, respectively. In the case when the restricted energy is indefinite, the Krein signature of the eigenvalue is indefinite as well. Note that it is easy to see that the Krein signature of each eigenvalue off the imaginary axis is zero, since if JLu = λu with λ = −λ, then (Lu, u) = 0 since λ(J −1 u, u) = (Lu, u) = (u, Lu) = (u, λJ −1 u) = −λ(J −1 u, u).
Let us point out that in [20,72] (also see [23,50]) it is rigorously proved for finite-dimensional systems that the only way eigenvalues of a system depending on a parameter can leave the imaginary axis (as a quadruplet, since the symmetry of the problem forces two complex conjugate pairs of the purely imaginary eigenvalues to collide simultaneously) is via a collision of two purely imaginary eigenvalues of opposite signature. This behavior is generic -the passing of two eigenvalues of opposite Krein signature on the imaginary axis is an event of codimension one (R. Kollár & P. Miller 2010, unpublished), as a particular quantity involving eigenvectors of the associated eigenspaces must vanish.
For simplicity, in what follows we will always assume J is invertible. Then λu = JLu is equivalent to Lu = λJ −1 u, and thus For this reason it is sometimes more convenient to use the real quadratic form (iJ −1 u, u) instead of (Lu, u).

Finite systems of eigenvalues
In this subsection we formulate precise statements of some key properties of Krein signature that apply to unbounded operators of the type we consider.
A key fact that makes the Krein signature so useful for continuation problems is that for finite systems of imaginary eigenvalues, it can never be (positive or negative) semi-definite without being definite. This is a consequence of the following.
Lemma 1. Let X be a Hilbert space, and let L be a symmetric and J a skewadjoint operator on X with a bounded inverse J −1 . Let Σ ⊂ iR be a finite set of discrete purely imaginary eigenvalues of JL and let X 1 be the corresponding spectral subspace (the span of all generalized eigenvectors for eigenvalues in Σ). Then the quadratic form (iJ −1 u, u) is non-degenerate on X 1 .
Proof. This result is well-known in the finite-dimensional case (see [72, p. 180]) and the proof here is not very different, based on the spectral decomposition of the underlying Hilbert space into the finite-dimensional JL-invariant subspace corresponding to Σ and its JL-invariant spectral complement. The spectrum of JL| X1 is Σ, and the spectrum of JL| X2 contains no point of Σ. (See Theorem III-6.17 of [41].) We claim that (J −1 u, v) = 0 whenever u ∈ X 1 and v ∈ X 2 . This is true because, if u is a generalized eigenvector for some λ ∈ Σ, then (λ − JL) n u = 0 for some n, and if v ∈ X 2 , then w = (λ − JL) −n v ∈ X 2 , and one checks that Now, if the form (iJ −1 u, u) is degenerate on X 1 , it means there exists a nontrivial u ∈ X 1 such that (J −1 u, v) = 0 for all v ∈ X 1 , hence for all v ∈ X. Hence J −1 u = 0, so u = 0, a contradiction.
A simple corollary for simple eigenvalues immediately follows. Corollary 1. Let X, L, J be as above. Let λ ∈ iR be a simple isolated eigenvalue of JL with corresponding eigenvector u. Then (iJ −1 u, u) is non-zero, and if λ = 0, so is the Krein signature sgn(Lu, u).
In the application we will consider, the operator L depends continuously (in a sense we will make precise) on a real parameter µ. By the corollary above, the only possibility for a continuously varying eigenvalue λ = λ(µ) of JL(µ) to change its Krein signature is for it to collide with another eigenvalue, or to cross zero. As a simple consequence of (49), a simple eigenvalue crossing zero will flip its Krein signature to the opposite. (If the change is from negative to positive this may correspond to symmetry-breaking instability [50].) Corollary 2. Let X be a Hilbert space, and let J be a skew-symmetric operator on X with bounded inverse. Suppose s → L(s) is a family of symmetric operators on X, for s in an interval I ⊂ R. Assume that for s ∈ I, (λ(s), u(s)) is a continuous family of isolated simple purely imaginary eigenvalues and corresponding eigenvectors for the problem λ(s)u(s) = JL(s)u(s) .
The previous two results guarantee that the Krein signature of a purely imaginary eigenvalue does not change unless the eigenvalue crosses the origin or collides with other eigenvalues. Also, Theorem 2 implies that the operators JL = iRL j we will consider have only discrete eigenvalues of finite multiplicity, so only finitely many eigenvalues may collide at one point, for any value of the parameter µ.
Such colliding eigenvalues will form a finite system of eigenvalues in the sense of Kato (see section IV-3.5 of [41]), with a corresponding spectral projection X → X 1 (µ) that varies continuously with µ. If no eigenvalue of negative or indefinite signature is present in this family before collision, then the signature is positive definite and must remain so at collision and after. Then all eigenvalues must remain on the imaginary axis after collision -no pair of eigenvalues can bifurcate off the axis, because the signature of such a pair is indefinite. This is well known for finite-dimensional systems [50]. To justify these statements about continuation of Krein signature for unbounded operators, we should first define an appropriate notion of continuity. The following definition is essentially related to the concept of convergence of closed operators in the generalized sense of [41], see Theorem 2.25, Section IV-2.6. Definition 1. Let s → T (s) be a family of closed operators in X, defined for s in an interval I ⊂ R. We say the family T (s) is resolvent-continuous on I if for each s 0 ∈ I, there exists ξ ∈ C such that the resolvent map s → (ξ − T (s)) −1 is defined and norm-continuous for s in a neighborhood of s 0 .
For such a resolvent-continuous family T (s), the resolvent map is actually (locally) continuous in s for each point ξ of the resolvent set of T (s).
Suppose now that Σ 0 is a finite set of discrete eigenvalues of T (s 0 ). This set may be continued continuously to comprise a finite system of eigenvalues Σ(s) defined on a maximal interval of existence I ′ ⊂ I in a standard way: Let Γ be any smooth contour (a collection of small circles, say) that contains Σ 0 inside, and no other point of the spectrum of T (s 0 ). The spectral projection Theorem 3. Let X be a Hilbert space and let J be a skew-symmetric operator on X with bounded inverse. Suppose s → L(s) is a family of symmetric operators for s in an interval I ⊂ R such that the family JL(s) is resolventcontinuous. Suppose Σ(s) is a finite system of discrete eigenvalues of JL(s) as described above, with associated generalized eigenspace X 1 (s), defined for s in some maximal interval I ′ .
Then, if the quadratic form (iJ −1 u, u) is (positive or negative) definite on X 1 (s 0 ) for some s 0 ∈ I ′ , it has the same definiteness on X 1 (s) for all s ∈ I ′ . By consequence, all eigenvalues in Σ(s) are purely imaginary for all s ∈ I ′ .
Proof. The proof is again a rather straightforward extension of well-known results from finite dimensions [20,72]. Whenever the form (iJ −1 u, u) is definite on X 1 (s) we must have Σ(s) ⊂ iR as follows from (48). Considering the positive definite case, let ρ(s) be the infimum of (iJ −1 u, u) over the unit sphere in X 1 (s). Then ρ(s 0 ) > 0 and ρ is continuous on I ′ . If definiteness fails at some point of I ′ , then ρ(s) > 0 on some maximal strict subinterval of I ′ and ρ(s 1 ) = 0 at an endpoint s 1 ∈ I ′ . But then Σ(s 1 ) ⊂ iR by continuity of the eigenvalues and so ρ(s 1 ) > 0 by Lemma 1, a contradiction.
Note that in our application the resolvent-continuity condition on the family of operators is satisfied, as the domain Y ⊂ X of the operators L(µ) is fixed, and the operators JL(µ) vary continuously in L(Y, X).

Negative-signature eigenvalues for the GP equation
Here we identify all imaginary eigenvalues with negative Krein signature for µ = m + 1, when the linearized problem reduces to the case of (14), the Gross-Pitaevskii equation linearized about a trivial solution. The linearized system (32) for given j and m then reduces to see (32)- (33). According to Proposition 4 we may assume j, m > 0. Since, by the discussion following (16), the eigenvalues of H k are exactly |k| + 1 + 2n for n = 0, 1, 2, . . . , with eigenfunctions w The eigenvalues of the linearized Gross-Pitaevskii equation for µ = m + 1 then come in two families for j > 0: for n = 0, 1, 2, . . .. It is easy to determine the Krein signature of these eigenvalues since in the different cases we find (L(u, v) T , (u, v) T ) = ((H j+m − µI)u, u) = iλ(u, u) = (j + 2n)(u, u) , respectively. Note that according to Proposition 4 we can restrict our search for possible unstable eigenvalues to the modes with 0 < j < 2m. In particular, this means that when µ = m + 1, the only eigenvalues with negative Krein signature are

Numerical Methods and Results
In this section we first describe in detail the way we evaluate the Evans function and determine the presence and location of eigenvalues. We then show how the results on Krein signature from Section 6 are used to justify the finding of all unstable eigenvalues and interpret eigenvalue collisions. Based on our analytical and numerical results we also suggest an improved numerical method to detect unstable eigenvalues by tracking eigenvalues of negative Krein signature. Despite we have not used the method in the present work we believe that such an algorithm offers a significant reduction of computational cost in a wide variety of numerical studies of stability problems that use the Evans function. In the further subsections, numerical results are discussed separately for singly-and multi-quantized vortices since the stability diagrams (diagrams of stable and unstable eigenvalues) reveal different patterns.

Evans function evaluation
To attain high precision and stability for all computations, exterior products were used throughout for numerical evaluation of the Evans function. Note that exterior products may be avoided, however, if one uses the algorithm introduced in [27,28]. As easily seen from the asymptotic descriptions in Theorem 2, the behavior of solutions of the system (41) is significantly different for r ≪ 1 and r ≫ 1. Consequently it is useful to rescale the solution during the integration process over the interval (0, ∞). (The implementation approximates this interval with [ε, R], where ε = 10 −7 and R = R(µ) is set in such a way that the vortex solution is negligible at r = R. This R(µ) is chosen as an increasing function with R(0) = 5, R(35) = 25.) The aim is to rescale in such a way that the matrix B(r, λ, j) (and the solution) remains appropriately bounded throughout the integration. The details of the rescaling used here can be found in [44].
The presence of unstable eigenvalues is detected by contour integration using the argument principle, similarly as in [54], for j's restricted by Proposition 4. The algorithm adaptively calculates the argument of the Evans function E j (λ) along an approximately rectangular contour Γ which encloses a bounded region in R 2 . The region is pictured on Fig. 3. Note that Proposition 6 allows us to confine the integration to the right half plane (the total argument change is twice as large) and also reduces the set of j's for which the calculation is necessary to positive values. Moreover, Proposition 5 restricts the location of unstable eigenvalues to a vertical strip. Naturally, it is not possible to perform numerical calculations without imposing also some vertical bound for the region enclosed by the contour. The vertical bound in the implementation is chosen is such a way that the behavior of the stable eigenvalues becomes predictable. We justified that all the unstable eigenvalues are included a posteriori by examining the Krein signature of eigenvalues.
Stable eigenvalues on the imaginary axis are, thanks to the symmetry of the Evans function, zeros of the real-valued function E(λ). Hence one can plot that real function and determine the location and multiplicity of its zeros within a finite interval. In the actual implementation this is done automatically -one first interpolates the real function by a cubic spline and then uses the Newton method for locating zeros. In a small neighborhood of a possible double zero, a very fine mesh was used to resolve any ambiguity.
The total number of eigenvalues of E j (λ) enclosed in the region is then determined by the difference n su = n s+u − n s between the total number of eigenvalues inside Γ, and the number of (stable) eigenvalues n s on the imaginary axis, including their algebraic multiplicity. If n su is not equal to zero it must by Proposition 6 be an even positive integer and corresponds to twice the number of pairs of stable and unstable eigenvalues (λ, λ) enclosed within Γ. The precise location of unstable eigenvalues can be theoretically determined by the generalized argument principle. Higher moments give the sum of k-th powers of positions of all eigenvalues enclosed by Γ. Given the approximate location of eigenvalues on imaginary axis and number of eigenvalues off the axis, the problem reduces to solving a set of nonlinear equations. This is particularly efficient in the case n su = 2, where only s 1 is necessary. Unfortunately, even in this case, the numerical error involved can be significant, with a major contribution coming from a finite difference approximation of E ′ j (λ). Therefore the obtained values are considered only approximate. The calculated location is further used to construct a smaller contour Γ s which lies solely in the right-half plane and encloses only a single eigenvalue. The presence of a zero of E j (λ) inside a smaller contour is again justified by the argument principle. Its location is then determined by the generalized argument principle. This process can be repeated a few times until a desired precision is attained. In the implementation the threshold for a precision was set up to be 10 −3 .
Note that this method of locating eigenvalues does not allow one to calculate the eigenfunction directly; for that another method must be used.

Numerical uses of Krein signature
Here we show how the Krein signature serves to provide evidence that no unstable eigenvalues are missed in the numerical computations. Moreover we suggest a very efficient numerical algorithm based on the theory in Section 6. Finally, we discuss an approach that does not require homotopy and avoids path-following completely.
In our numerical approach we consider the linearized eigenvalue problem for a large range of the parameter µ, continuously connecting the linearized eigenvalue problem about a vortex solution to a linearized eigenvalue problem about a trivial solution. (See [14] for a similar approach). The latter problem has only purely imaginary eigenvalues, and their Krein signatures were determined in Section 6.2. We plot the location of all imaginary eigenvalues in a bounded interval on the imaginary axis which contains the continuation of all negativesignature eigenvalues for the whole range of µ considered (µ ∈ (m + 1, 35)). By calculating the Evans function on the contour Γ and on an appropriate portion of the imaginary axis, we indirectly track signature changes using the theory of Section 6, and can thus infer restrictions on possible departures of eigenvalues from the imaginary axis. This tracking process provides strong evidence that there are no unstable eigenvalues other than the ones we find.
For a significant reduction in computational effort, we propose a new numerical method that avoids the computation of contour integrals for Evans functions. To identify all unstable eigenvalues, it is only necessary to perform the following calculations: • evaluate the Evans function on the imaginary axis close to zero to detect any eigenvalues crossing zero; • follow the eigenvalue branches starting at negative Krein signature eigenvalues at the reference value of the continuation parameter (in our case start at µ = m + 1 and follow a total of 2m − 1 different branches), to detect any collision with positive Krein signature eigenvalues; • follow any branches of eigenvalues that bifurcate into the complex plane.
We emphasize that the theory presented in Section 6 justifies that no unstable eigenvalues can be left out. No contour integration is necessary at all since one can just use a path-following algorithm (a boundary-value solver). Moreover, global analyticity of the Evans function is not needed, since only zeroes of the real Evans function on the imaginary axis need to be found initially, and one may use a continuation algorithm [27,61] to track them. A negative feature of the homotopy technique is that it has a large overhead if one is interested in only a small set of values of the parameter. Therefore we examine ways of avoiding this overhead by directly evaluating the Krein signature of a given eigenvalue on the imaginary axis. This is not feasible with the approach of evaluating Evans functions using exterior products, since the eigenfunction is unavailable and it must be calculated separately by a different method. On the other hand, the restriction of the Evans function to the imaginary axis is a real function. This means that only a continuous representation of the Evans function is required and problems of rapid growth and numerical dependence are much less severe. By consequence, it is not necessary to use exterior products to evaluate the Evans function, and it can be evaluated either (i) by calculating a relevant Wronskian of solutions to the linear eigenvalue equations which makes it easy to compute the eigenfunction, or (ii) one can use a simple continuation algorithm proposed by Sandstede [27]. By either of these methods, it is possible to evaluate the Krein signature directly, which then allows one to determine the number of unstable eigenvalues.
Then for an eigenvalue problem of the form JLu = λu the following method may be used. First, determine (somehow, by using oscillation theory, for example) the number n total of negative eigenvalues of L. If the number is zero, there are no unstable eigenvalues of JL. If n total > 0, the number gives the total number n u of pairs of unstable eigenvalues of JL plus the number n s of stable eigenvalues with negative Krein signature [39,40,56] (also see [25]). If the present method (the Evans function) detects a number of pairs of eigenvalues off the imaginary axis n u = n total , there are no other unstable eigenvalues. If n u < n total , perform a calculation of the Krein signature for eigenvalues on the imaginary axis. If n u + n s = n total there are no other unstable eigenvalues, otherwise, increase the area of search and repeat the whole process.
For small values of µ, close to µ 0 = m + 1, the eigenvalues are close to the eigenvalues of the reduced uncoupled linear problem neglecting the |w| 2 dependence. As µ increases, certain eigenvalues remain constant: a double eigenvalue λ = 0 and simple eigenvalues λ = ±2i for j = 0 and simple eigenvalues λ = ±i for j = 1. These eigenvalues originate in the symmetries and boosts of the Gross-Pitaevskii equation and are present for every m (see the Appendix).

Multi-quantized vortices m ≥ 2
The eigenvalue diagrams for multi-quantized vortices with m = 2 show more complexity. The radial vortex profile for µ ≈ 35 is illustrated in Fig. 1.
In the case m = 2 the possible unstable eigenvalues may appear for |j| = 1, 2, 3. The modes j = 0 and j = 1 demonstrate the same features as in the case of the singly-quantized vortex with the same constant eigenvalues: a  double eigenvalue λ = 0 and simple eigenvalues λ = ±2i for j = 0 and simple eigenvalues λ = ±i for j = 1 (see Fig. 5 (a)).
A different behavior appears for modes j = 2 and j = 3, as shown on Fig. 5 (b) and Fig. 7. For j = 3 there are no unstable eigenvalues present and the stable eigenvalues do not collide but rather diverge from each other when they approach each other. Collisions do occur for j = 2 and cause instability. A collision of two stable purely imaginary eigenvalues produces a pair of stable and unstable complex eigenvalues symmetric with respect to the imaginary axis. After an increase of µ these two eigenvalues return to the imaginary axis and split to two purely imaginary eigenvalues as illustrated in Fig. 6. This "collision -split -collision" process ("bubbles of instability" [50]) repeats regularly for the whole range of µ studied. This behavior is caused by the presence of a single eigenvalue with negative Krein signature [50,68]. The imaginary part of this eigenvalue is decreasing with increasing µ, and on its way it encounters eigenvalues with the opposite signature. After each collision the eigenvalues split off the imaginary axis and become eigenvalue pairs with zero Krein signature symmetric relative to the imaginary axis. Reversibility of this process suggests that the eigenvalues come back to the imaginary axis and the process repeats itself (for a larger parameter µ). This indicates a surprising thing: transitions to instability for larger µ may happen at a large frequency (Im λ large), and therefore there is no hope to confine the imaginary parts of unstable eigenvalues to a finite interval independent of µ. The behavior of the eigenvalues demonstrates strong agreement with [60]. This is also consistent with the results of Seiringer [65] where he proved that for any m ≥ 1 for large enough µ the vortex becomes energetically unstable in the sense that it cannot be a global minimizer of the energy and is subject to symmetry breaking. Also note, that the maximum real part of unstable eigenvalues grows very slowly with growing parameter µ. The maximum of the real part is much smaller than the bound we were able to obtain in Proposition 5.
In the case j = 3 the eigenvalues have the same Krein signature and therefore they cannot split off the imaginary axis The eigenvalue which originates at λ = −i for µ = 3 has negative signature at first, but after it crosses zero for µ close to 5, it changes its signature to positive according to Corollary 2. Then all eigenvalues have the same signature, making splitting impossible. Instead, eigenvalues repel each other upon approach, as expected by [50].

Further results
The presence of exponential instability was also checked by direct simulations. First, an approximate eigenvector was obtained by Galerkin approximation [18] and then the Strang-splitting scheme [6] was used for time evolution. The initial perturbation of a vortex solution by an approximated eigenvector showed a good agreement with the expected exponential growth.
For the case m = 3, we performed similar computations but omit details for brevity. There are more eigenvalues of negative Krein signature, and overlapping bubbles of instability. See [60] for the analog of Fig. 6 in this case.  Finally, it would be plausible to describe the asymptotics of the eigenvalues as µ → ∞ when the condensate approaches the Thomas-Fermi regime. We were only able to study the asymptotic behavior of purely imaginary eigenvalues numerically by plotting a loglog plot of the first order differences of λ(µ). We observed a clear linear trend implying an algebraic approach to a limit. Our numerical results agree with the recent analytical results in [57] for stability of vortex solutions in the limiting Thomas-Fermi regime (see also [17] for stability of ground states).
The approximate behavior of eigenvalues with a small imaginary part for m = 1, j = 1, 2, and m = 2, j = 2, 3, is presented in Table 1. For most eigenvalues, the asymptotic behavior as µ ≫ 1 is well approximated by On the other hand, certain eigenvalues show different rate of convergence clearly distinct from (−µ −1 ), but we do not have any explanation of this phenomena.

Appendix: Symmetries and eigenvalues
The symmetries of the Gross-Pitaevskii equation and its linearization imply the presence of a special set of eigenvalues. Phase symmetries. For any m, λ = 0 (for j = 0) is a constant double eigenvalue for all µ ≥ µ 0 . Its multiplicity comes from the symmetry of the Gross-Pitaevskii equation under a phase change (this generates an eigenvector) and under a change of standing-wave frequency (this generates a generalized eigenvector). A detailed discussion on these two symmetries and their implications for spectra is given in [54]. GGV boost. Due to the presence of the harmonic potential, the other usual invariants of the nonlinear Schrödinger equation-spatial translations-do not apply. Similarly, one cannot perform the typical Galilean boost. Instead, one can apply a boost for quadratic potentials that was described by García-Ripoll et al. [19]. This particular boost transforms any solution Ψ(r, t) (r = (x, y)) to a new solution of the form Ψ R (r, t) = Ψ(r − R(t), t) exp(iθ(r, t)) , where R(t) is the path of a classical particle moving in the potential well. For the harmonic potential V (r) = 1 2 r 2 , this requires simply R tt = −R, so that each component of R(t) can be any linear combination of cos t and sin t. According to [19], θ is given up to a constant by θ(r, t) = r · R t .
As discussed in [54], any one-parameter family u τ of solutions to (6) gives rise to a solution of the corresponding equation linearized about u 0 , given bỹ u := ∂ τ u τ | τ =0 .
Therefore the GGV boost implies the existence of the eigenvalues (for any m) λ = ±i (for j = 1).
Breather boost. Finally, we describe the source of the presence of the eigenvalues λ = ±2i (for j = 0), which appear for every m ≥ 1 and µ ≥ µ 0 . These eigenvalues corresponds to a "breathing" symmetry of the Gross-Pitaevskii equation in 2 + 1 dimensions, as explained by Pitaevskii and Rosch in [59]. Here we provide a somewhat different perspective, showing that the symmetry corresponds to the Talanov lens transformation via an exact transformation to the cubic Schrödinger equation.
Rybin et al. [64] noted that the radially symmetric Gross-Pitaevskii equation in 2+1 dimensions transforms exactly to the cubic Schrödinger equation without potential, according to a transformation described originally by Niederer [52] for the linear Schrödinger equation. More generally, they noted that this transformation works in any dimension as long as the nonlinearity is critical, of the form |u| 4/n u in n space dimensions. This symmetry was used by Carles [7] to study various mathematical aspects of the Gross-Pitaevskii equation. By transforming from GP to cubic Schrödinger, using a Talanov lens transformation [70], and transforming back, one can get a self-transformation of GP (see also [21]). The symmetry that we will describe is somewhat more general, and works as follows. Suppose that v(x, t) is any solution of the equation where t ∈ R and x ∈ R n . Let where a, b, c and τ are functions of time t. Then we find that i∂ t u + 1 2 ∆u − 1 2 ν 2 |x| 2 u − γ|u| p u = 0, provided that b and c satisfy b ′ − b 2 + ω 2 c 4 = ν 2 , c ′ = bc, and a = c n/2 , τ ′ = c 2 , γ = λc 2−np/2 .
In the other direction with constant ω 2 > 0 and ν 2 = 0 we get [7]: The consequences of this transformation for the Gross-Pitaevskii equation in 2+1 dimensions are striking. From any solution of (54), one gets from (55) a solution obtained by a breather boost -a time-periodic dilation of space with an appropriate radial phase adjustment. See [59] for more discussion and a relation to representations of the group SO(2,1).
A short calculation shows that this transformation corresponds to the eigenvalues ±2i of the normalized linearized equations (29)- (30). The relation to these eigenvalues can be seen also from a simple consideration, that this selftransformation produces small oscillations at exactly twice the trap frequency for classical oscillations in the harmonic potential 1 2 ω 2 |x| 2 . The exact formulae for this breathing boost also show that breathing oscillations need not be small. Also note that a breather boost can be applied to any solution, not only standing waves.