Dissecting zero modes and bound states on BPS vortices in Ginzburg-Landau superconductors

In this paper the zero modes of fluctuation of cylindrically symmetric self-dual vortices are analyzed and described in full detail. These BPS topological defects arise at the critical point between Type II and Type I superconductors, or, equivalently, when the masses of the Higgs particle and the vector boson in the Abelian Higgs model are equal. In addition, novel bound states of Higss and vector bosons trapped by the self-dual vortices at their core are found and investigated.


Introduction
Vortex filaments carrying a single quantum of magnetic flux were discovered by Abrikosov in the realm of the Ginzburg-Landau theory of Type II superconductors in Reference [1]. The same magnetic flux tubes reappeared in the relativistic context of the Abelian Higgs model in the paper of Nielsen and Olesen [2], where their stringy nature was emphasized. Analytic formulas are available in this Reference for the 1vortex profile near the center of the core and far away from the origin, although the behaviour at infinity was refined in [3]. An important step forward in our knowledge of the mathematical properties of these extended structures was achieved by Bogomolny, who identified in the seminal paper [4] a system of firstorder PDE such that their solutions are the ANO vortices at the transition point between Type II and Type I superconductivity phases, the critical value where the quotient of the scalar and vector particle masses is one. These Bogomolny-Prasad-Sommerfield, [5], or self-dual 1 , vortices have very interesting features: (1) The magnetic flux is a topological quantity related to the first homotopy group of the circle of degenerate vacua. (2) At the BPS limit, these line defects do not interact with one another. They are thus free to move and zero modes of BPS-vortex fluctuations exist.
The primary aim of this work is to investigate the self-dual vortex zero modes of fluctuation. The BPS vortex PDE equations admit multivortex solutions. Proof of the existence of this type of solitons was given by Jaffe and Taubes in Reference [6]. The BPS multivortex moduli space with magnetic flux equal to 2π e n, with n ∈ Z integer, is the space of |n| unordered points in the R 2 plane [7]. The freedom in the locations of the centers preludes the existence of 2|n| linearly independent zero modes of fluctuation, a fact proved by E. Weinberg in [8] with a shrewd generalization of the index theorem of elliptic operators. This computation, motivated by a physical problem, paved the way to extending the Atiyah-Singer index theorem usually observed in compact spaces, with or without boundaries, to open spaces where problems with the continuous spectrum arise, see e.g. [9,10]. The self-dual vortex solutions with cylindrical symmetry aroused special interest. In this case the vortex first-order equations reduce to an ordinary differential equation system which is solvable near the origin and very far from the vortex core. Several interpolation methods have been developed, either numerically or through some functional series expansion, to obtain the full multivortex solution, which is never expressible in terms of elementary or special functions, see [11]. In a later development, the self-dual cylindrically symmetric vortex zero mode fluctuations were studied in detail starting with E. Weinberg seminal paper [8], see also Chapter 3 in the recent monograph [12]. Given the rôle of the vortex zero modes in the analysis of the low-energy vortex dynamics as geodesic motion in the moduli space of BPS vortex solutions, see e.g. [7]- [15], better ansatzes for the analytical structure of the zero mode fluctuations of BPS vortices were proposed for this purpose in References [13,14]. This task was fully achieved in the papers just mentioned for the solutions with a low number of magnetic flux quanta, e.g., n = 2, 3. In the first half of this paper we perform a complete and detailed analysis of the structure of the zero modes of fluctuation around BPS cylindrically symmetric vortices. Relying on the Ruback-Burzlaff ansatz, we describe the vortex zero mode profiles with the same level of precision as the precision attained in the knowledge of the BPS vortices themselves. After identifying analytically the zero mode radial profile near the core and close to infinity we perform the interpolation between these two regimes by means of a shooting procedure implemented numerically. The angular dependence of the zero mode wave function is fixed analytically by Fourier analysis. The regularity of the wave function near the origin and exponential decay at infinity, all together guaranteeing normalizability, impose the existence of 2|n| linearly independent vortex zero modes in concordance with the index theorem. The interest of this study is twofold: (1) It extends the work of several authors on this subject to BPS vortices with more than three quanta of magnetic flux. (2) Recently, in [16] and [17] two of us improved on the one-loop shift calculations of kink masses and domain wall surface tensions by controlling the inaccuracies induced by zero modes in the heat kernel/zeta function regularization procedure. The new method requires precise information about the zero mode wave functions such that the information gathered in this paper is necessary to improve the results obtained in [18,19,20,21,22] by diminishing the impact of zero modes in heat kernel expansions 2 .
However, vortex zero modes exist and are influential not only in critical vortices between Type I and II superconductors. Jackiw et alli, see e.g. [24], discovered that the spectrum of the Dirac operator in a vortex background includes |n| linearly independent fermionic eigenfunctions of zero eigenvalue, where n is the vortex magnetic charge. From a mathematical point of view the existence of zero modes in the vortex-fermion system obeys an index theorem on a open space, the R 2 plane. Besides these topological roots underlying their existence, the vortex-fermion zero modes add quantum states in the middle of the mass gap of the Dirac spectrum which, in turn, induce the phenomenon of fractionary charge. Thus, the context in which fermionic zero modes in a classical vortex field are considered is completely different: there are no scalar and vector particles and the vortex external field does not fluctuate as in the Abelian Higgs model. Quite recently, this problem gained importance in condensed matter physics, for instance in the mathematics and physics of graphene, see e.g. [25], or, in class A chiral superconductors, see [26].
Soon after the discovery of vortex filaments in Type II superconductors, interest aroused in the investigation of fermionic bound states trapped at the vortex core by looking at the one-particle spectrum of the Bogolyubov-de Gennes equation near a magnetic flux line background, see [27]. This pioneer paper by de Gennes et al prompted a long search aimed at unveiling the nature of this type of bound states, although without complete success from the analytic point of view. Nevertheless, interesting effects of these bound states on the vortex core have been disclosed in a superfluid phase of the 3 He isotope, see [28]. As a secondary goal, we shall study here the bound states arising when scalar and/or vector bosons are trapped at the core of a self-dual vortex in the framework of the Abelian Higgs model, mutatis mutandis in the Ginzburg-Landau phenomenological theory of superconductivity. Contrarily to the bound states mentioned above the particles trapped by the BPS vortices are bosons rather than fermions. In the context of the Abelian Higgs model bound states of mesons by vortices were discovered by Goodman and 2 Our method applies not only to conventional topological defects but also to instantons, see Reference [23]. [29] in the mid nineties. Another papers where the rôle of these bound states in the framework of topological defects in Cosmology is emphasized are [30,31,32]. Taking profit of the supersymmetric quantum mechanical structure linked to BPS topological defects we were able in the short letter [33] to offer a quite detailed description of such bizarre bound states. We shall develop in this work a more complete analysis of the meson bound states on BPS vortices and we shall discuss their properties by comparison with the well known BPS vortex zero modes. Our approach follows the pattern found in the λφ 4 kink. Fluctuations of the domain wall defects in this model are of three types: 1) translational (zero) modes where a meson travels together with the kink center of mass without disturbing the defect profile. 2) kink internal modes of fluctuation where a meson is trapped forming a meson-kink bound state that produces an oscillating in time deformation of the defect profile. The existence of this second type of fluctuations is due to the supersymmetric quantum mechanics of the kink stability problem. 3) Scattering of mesons through the wall.

Hindmarsh in Reference
We shall address types 1) and 2) of fluctuation concerning the BPS vortices in the Abelian Higgs, a much more difficult task. In fact, the search for vortex fluctuations with frequencies greater than 0 but lower than the threshold of the continuous spectrum only differs from the search for zero modes in the fact that the eigenvalue is unsettled a priori. The shooting procedure for obtaining the form factor of the bound state in the intermediate region is thus ineffective and we shall approximate the radial ODE by means of a discretization of the radial coordinate, transforming this ODE into a linear system of difference equations. The bound state eigenvalues will be identified via diagonalization of the matrix of the linear system, after which the eigenfunctions will be found numerically.
The paper is organized as follows: In Section §.2 the Abelian Higgs model is revisited with the aim of fixing our notational conventions. Section §.3 is devoted to describing in a detailed manner the critical regime between Type I and Type II superconductivity leading to the BPS system of first-order PDE governing the static solutions of finite energy density. Also, the second-order differential operator, which is usually referred to as the Hessian, determining the small fluctuations around the vortex solutions is discussed in this Section and its factorization as the product of two first-order PD operators is explained. In Sections §.4 and §.5 the general structure of the fluctuation spectrum of cylindrically symmetric BPS vortices is developed forming the main contribution of the paper. Section §.4 offers a comprehensive analysis of the BPS vortex zero modes of magnetic flux n and unveils the general pattern of zero mode fluctuations of a cylindrically symmetric BPS vortex carrying n quanta of magnetic flux. In Section §.5 a similar picture describing the features of several boson-vortex bound states also with low magnetic charge, is developed. Finally, in Section §.6 we draw some conclusions and speculate about some future prospects.

2
Topological defects carrying quantized magnetic flux in superconducting systems We start from the action of the Abelian Higgs model that describes the minimal coupling between a U (1)gauge field and a charged scalar field in a phase where the gauge symmetry is broken spontaneously. In terms of non-dimensional coordinates, couplings and fields, the action functional for this relativistic system in R 1,2 Minkowski space-time reads: The main ingredients are one complex scalar field, φ(x) = φ 1 (x) + iφ 2 (x), the vector potential A µ (x) = (A 0 (x), A 1 (x), A 2 (x)), the covariant derivative D µ φ(x) = (∂ µ − iA µ (x))φ(x) and the electromagnetic field tensor F µν (x) = ∂ µ A ν (x) − ∂ ν A µ (x). We choose the metric tensor in Minkowski space in the form g µν = diag(1, −1, −1), with µ, ν = 0, 1, 2, and use the Einstein repeated index convention. In the temporal gauge A 0 = 0, the energy of static field configurations becomes which, in a non-relativistic context, is the free energy of a superconducting material arising in the Ginzburg-Landau theory of superconductivity, see formula (17) in [1] where the order parameter φ responds to the Cooper pairs density. The search for static configurations requires us to look at the extrema of the functional: The critical points of V [φ, A] are the static fields satisfying the second-order PDE system Solutions of (2) that comply with the asymptotic boundary conditions at the circle at infinity, i.e. when have finite energy. In fact, choosing where θ = arctan x 2 x 1 and n is an integer, as representatives of (3), one checks that the configuration space of the static fields is the union of Z topologically disconnected sectors: The fields in each sector C n are asymptotically constrained by the formula (4) where it is evident that n is the winding number of the map φ| ∞ : lim r→∞ S 1 r −→ S 1 from the circle at infinity in the x 1 : x 2 -plane to the vacuum orbit determined by the phase of the scalar field at infinity. Thus, all the field configurations in the non-trivial sectors, n = 0, are endowed with a quantized magnetic flux: Rotationally symmetrical solutions of (2) with finite energy for κ 2 = 1 and a quantum of magnetic flux, n = ±1, are vortices given that the vector field (A 1 , A 2 ) is purely vorticial. Choosing, e.g., n = 1, the ansatz, see [2], together with the asymptotic conditions, f (∞) = 1 and β(∞) = 1, and the regularity conditions, f (0) = 0 and β(0) = 0, convert the PDE system (2) into a second-order ODE system and at the same time guarantees regular behavior at the origin and appropriate fall-off at infinity. There is no available analytical solution to this ODE system. The asymptotic form of these solutions, however, is known, see [3]. Linearization of equations (6), (7) around the scalar and vector fields vacuum values reveals that: where c V and c H are integration constants. The full Nielsen-Olesen vortex profiles can only be identified by numerical methods. One might also search for vortex solutions carrying n quanta of magnetic flux. Vector and scalar mesons respectively produce repulsive and attractive forces of Yukawa type between charged objects of the same type. Thus, an effective potential arises prompting vortex solutions of unit magnetic flux to either repel, if κ > 1 (Type II superconductors), or attract each other when κ < 1 (Type I superconductivity materials). In the first case, the vortices are arranged in a triangular Abrikosov lattice, whereas in Type I superconductors the magnetic flux aggregates on slices piercing the material.

Self-dual/BPS vortices and their fluctuations
At the transition point κ = 1 between Type I and II superconductors no forces exist between the vortices, which thus, become very special. In order to investigate these critical vortices it is convenient to write V [φ, A] in the form, see Reference [4]: up to a total derivative term that integrates to zero over the whole plane if the fields tend to their vacuum values at infinity. The parameter κ, determined by the φ 4 and electromagnetic couplings as κ 2 = λ e 2 , measures the quotient between the penetration lengths of the scalar and electromagnetic fields in the superconducting medium. Values such that κ 2 > 1 characterize Type II superconductors (typically alloys) whereas type I superconductors (metals) correspond to κ 2 < 1, as explained above. In the QFT context the parameter κ is the quotient between the masses of the Higgs particle, m H = √ λ v, and the vector meson, m V = e v, after the Higgs mechanism has taken place giving to the photon a finite mass.
The critical vortices are solutions of the first-order PDE's which, written in terms of the real and imaginary parts of the complex scalar field φ = φ 1 + iφ 2 , read Moreover, the BPS vortices are subjected to the asymptotic conditions (4). It is clear that at κ 2 = 1 the energy of self-dual vortices saturates the Bogomolny topological bound: Additionally, it may be checked that self-dual vortices also solve the second-order PDE system (2). Proof of the existence of vorticial solutions of the PDE system (9) has been developed in Reference [6]. Given a positive integer n, there exists a moduli space of self-dual vortices solving the PDE system (9) characterized by 2n parameters, the centers of the magnetic flux tubes located at the zeroes of the scalar field counted with multiplicity n k , i.e., n = n 0 k=1 Φ k n k , where Φ k is the quantized magnetic flux of a vortex (the multiplicity) and n 0 is the total number of flux lines, see also Reference [7]. Behind this particular structure lies the fact that there are no forces between self-dual vortices (κ = 1) of one or several quanta of magnetic flux, which thus move freely throughout the x 1 -x 2 -plane.

The first-order fluctuation operator: hidden Supersymmetric Quantum Mechanics
Knowing that the energy of self-dual vortices is a topological quantity, there is no doubt about the stability of these topological solitons. The main theme of this paper, however, is the analysis of field fluctuations around self-dual vortices. We shall concentrate on two special modes of fluctuations: 1) the 2n vortex zero modes, those belonging to the kernel of the second-order fluctuation operator (the Hessian in variational calculus terminology), which arise because of the freedom of motion of the centers.
2) vortex internal modes of fluctuation corresponding to bound state normalizable eigenfunctions of the Hessian in the discrete spectrum and 3) scattering eigenfunctions in the continuous spectrum. Let us denote the scalar field φ V and the vector potential A V corresponding to a self-dual vortex solution of vorticity n as: The self-dual vortex fluctuations (a 1 ( x), a 2 ( x)) and ϕ( x) = ϕ 1 ( x) + iϕ 2 ( x) built around the BPS vortex fields (A 1 ( x; n), A 2 ( x; n)) = (V 1 ( x; n), V 2 ( x; n)) + (a 1 ( x), a 2 ( x)) are zero modes of fluctuation if the perturbed fields (10) are still solutions of the first-order equations (9). To discard pure gauge fluctuations, we select the "background" gauge as the gauge fixing condition on the fluctuation modes. The system of four PDE equations (9)-(11) is satisfied if and only if the four-vector field which assembles all the BPS vortex field fluctuations is annihilated by the first-order PDE operator: Note that this D operator is obtained by deforming the PDE system (9) together with the background gauge. Thus, the zero mode ξ 0 ( x) BPS vortex fluctuation fields belong to the kernel of the first-order PDE operator, Dξ 0 ( x) = 0, or, in components, are L 2 solutions of the PDE system: Important information about the particle spectrum in the Abelian Higgs model not only comes from the vortex zero mode fluctuations but also from vortex fluctuations demanding positive energy because these positive perturbations determine the dynamics of mesons in the different topological sectors. Thus, we shall investigate the spectral condition where λ is a label in either the discrete or the continuous spectrum useful to enumerate the eigenfunctions and eigenvalues, and H + is the second-order vortex small fluctuation operator coming from linearizing the field equations (in the background gauge) around the BPS vortices, see [20]. The fluctuation vectors ξ( x) belong in general to a rigged Hilbert space, such that there exist square integrable eigenfunctions ξ j ( x) ∈ L 2 (R 2 ) ⊕ R 4 belonging to the discrete spectrum, for which the norm ξ( x) is bounded: together with continuous spectrum eigenfunctions ξ ν ( x) with ν ranging in a dense set. One observes that the second-order differential operator H + is one of the two SUSY partner operators obtained from D as which are isospectral in the positive part of the spectrum 4 . Explicitly, we find The Hamiltonians H ± are superpartners in a supersymmetric quantum mechanical system built from the "supercharges" 5 which is governed by the SUSY Hamiltonian: The stability of the BPS n-vortex solutions implies that the H + -spectrum consists of non-negative eigenvalues. Indeed in [8] within the framework of index theory in open spaces, see [9]- [10], E. Weinberg proved that there are 2n linearly independent normalizable BPS vortex zero modes in the topological sector of magnetic flux n, i.e., the dimension of the algebraic kernel of D, henceforth of H + is n . By inspection one sees that H − lacks zero modes (all the potential wells are non-negative) and the index theorem dictates: which means that H + has 2n zero modes, see Appendix B in [12]. Analysis of the zero mode eigenfunctions begun in [8] and was further developed in References [13] and [14]. In the last two references the motivation to describe in detail the vortex zero modes came from the study of vortex scattering at low energies within the approach of geodesic dynamics in their moduli space, see e.g. [15].
The first goal in this paper is to seek a complete description of the vortex zero modes not fully developed in the previous references because interest was focused in sectors of very low magnetic charge. The second task that we envisage is the search for the eigenfunctions in the strictly positive spectrum of H + , in particular the bound states, i.e., the internal modes of fluctuation where the self-dual vortex captures scalar and/or vector mesons.

Zero mode fluctuations of BPS cylindrically symmetric vortices
The Nielsen-Olesen ansatz (5) generalized to the topological sector C n , in cylindrical coordinates also in field space, reads: In this ansatz we assume the radial gauge A r = 0, such that the vector field is purely vorticial. Plugging (16) into the first-order PDE system (9), the following ODE system emerges: The solutions for the radial profiles f n (r) and β n (r) are the self-dual vortex solutions 6 . The asymptotic conditions (3) demand that f n (r) → 1 and β n (r) → 1 as r → ∞. In fact, it is immediate to check that the asymptotic behaviour fits formula (8) when κ = 1, i.e., the critical value where the mass of the Higgs boson is equal to (henceforth less than twice) the mass of the vector boson. The requirement of regularity at r = 0 also fixes the behaviour of the solutions near the origin to be f n (r) = d n r n and β n (r) = e 2 r 2 , where d n and e 2 are integration constants. Choice of these constants must be tailored to fit the asymptotic behaviour. A shooting procedure implemented numerically allows us to solve the system (17) by interpolating the field profiles between their shapes in the neighborhoods of the origin and of infinity. In this way we construct the cylindrically symmetric self-dual n-vortices. Investigation of zero mode fluctuations around cylindrically symmetric self-dual vortices begins with rewriting the PDE system (13) in polar coordinates: where we recall that a 1 = a r cos θ − a θ sin θ and a 2 = a r sin θ + a θ cos θ.

Analytical investigation of the algebraic kernel of the first-order operator D
In this Section we shall prove a proposition characterizing the analytical behavior of the eigenfunctions which belong to the kernel of the operator D acting on self-dual rotationally symmetric vortices. First, we state a helpful Lemma which unveils the existence of a useful symmetry in BPS vortex fluctuation space.
Regarding vortex zero modes the main result is as follows: Proposition: There exist 2n orthogonal zero mode L 2 (R 2 )⊕R 4 -fluctuations of the self-dual cylindrically symmetric n-vortex solution of the form where k = 0, 1, 2, . . . , n − 1, and the zero mode radial form factor h n (r) satisfies the second-order ODE with the contour conditions h nk (0) = 0 and lim r→∞ h nk (r) = 0.
Proof: Because the discrete symmetry explained in the previous Lemma we are only interested in the construction of the n zero modes ξ 0 ( x, n, k), their orthogonal partners ξ ⊥ 0 ( x, n, k) follow immediately. The structure of the PDE system (18)-(21) suggests the ansatz where the radial and angular dependencies of the components of ξ 0 ( x, n, k) are separated, in the search for the BPS vortex zero modes. Plugging this ansatz into the set (18)-(21) of four ODE's, only two independent but coupled first-order ODE's remain: which govern the behaviour of the g nk (r) and t nk (r) functions. Continuity in the angular part of the solution requires k to be an integer, k ∈ Z, and the L 2 integrability of the fluctuations demands that ξ 0 ( x, n, k) 2 = 2π r dr g 2 nk (r) + t 2 nk (r) < +∞ (27) Solving for the function t nk (r) in (25), we obtain and plugging this expression into (26) we end with a single second-order ODE for g nk (r) in terms of the self-dual vortex profiles f n (r) and β n (r). Now we shall investigate the behaviour of this function: • Regularity of the function g nk (r) at the origin: The origin r = 0 is a regular singular point of the second-order differential equation (29). In this situation there exists a single analytic solution at r = 0 although the second linearly independent solution of (29) has a singularity. The analytic solution admits a series expansion around r = 0 of the form: s is chosen as the minimum value that selects c (n,k) 0 = 0, i.e., h nk (r) is regular, and does not vanish, at r = 0. Plugging (30) into (29), and taking into account that f n (r) ∼ r n ∞ =0 d n+2 r 2 and β n (r) ∼ e 2 r 2 + r 2n+2 ∞ =0 e 2n+2+2 r 2 near the origin, we obtain the identity from which we extract recurrence relations that determine the series expansion coefficients c (n,k) j up to order 2n + 1. Annihilation of the term independent of r in (31) unveils the indicial equation Because c (n,k) 0 = 0, the two characteristic exponents, the two values of s compatible with (32), are: Both possibilities are equivalent: simply redefine k, k ↔ −(k + 2). Thus, we shall stick to choice A. In terms of the function h nk (r), let us recall that g nk (r) = r n−k−1 h nk (r), the norm (27) of the vortex zero modes is Near the origin, the behaviour of the first summand of the integrand of (33) is In order to avoid singularities in the integrand coming from a pole at the origin, we demand that The term proportional to r in (31) is null if (1 + 2k)c (n,k) 1 = 0 for the characteristic exponent A. Thus, the second coefficient in the expansion vanishes: c (n,k) 1 = 0. Moreover, the recurrence relations extracted from (31) for the odd terms proportional to r 2i+1 , with i = 1, 2, . . . , [n − 1 2 ], show that all the odd coefficients c 2i−1 in the equation above are non-null. The two-term recurrence relations between the even coefficients is read from the annihilation of the terms proportional to r 2i , i = 1, 2, . . . , n, in (31): The relations (34) imply c 2k+2 not being null. Thus, in a neighborhood of the origin the function h nk (r) adopts the form: with c (n,k) 0 and c (n,k) 2k+2 arbitrary constants. Therefore, near the origin the second summand of the integrand in (33) behaves as: Singularities in the integrand coming from a pole at the origin are prevented if the inequality holds. Together with the inequality k ≤ n − 1, the univaluedness of the BPS vortex zero modes ξ 0 (r, θ, n, k) = ξ 0 (r, θ + 2π, n, k), equivalent to k ∈ Z being an integer, restricts the possible values of k to being 0 ≤ k ≤ n − 1, as stated in the proposition.
• Asymptotic behaviour of the function g nk (r): The following step is to investigate the asymptotic behaviour of the function g nk (r) far from the origin. Replacing the asymptotic vortex profiles when r → +∞, f n (r) → 1 and β n (r) → 1 into the differential equation (29) the modified Bessel equation appears. The solution is a linear combination of the modified Bessel functions with very well known asymptotic properties: By adequately tuning the integration constants c (n,k) 0 and c (n,k) 2k+2 we can get the exponentially decaying particular solution. We thus keep the exponentially decaying tail (36) in the vortex zero mode eigenfunctions that does not impose new restrictions. We conclude the finiteness of the norm (33) of ξ 0 ( x; n, k) and therefore ξ 0 ( x, n, k) ∈ L 2 (R 2 )⊕R 4 . Together with their partners ξ ⊥ 0 ( x; n, k), perpendicular to them in field space, they form a system of 2n orthogonal zero modes of fluctuation around any rotationally symmetric self-dual vortex of magnetic flux n that are normalizable. The orthogonality between the zero modes ξ 0 ( x; n, k) (22) arises directly from the angular dependence of these eigenfunctions. Now let us insert g nk (r) = r n−k−1 h nk (r) into equations (29); we end with the differential equation (24) and the general form (22) of the vortex zero modes according to the above Proposition. Notice that t nk (r) = −r n−k−1 h nk (r) fn(r) .

Weakly deformed cylindrically symmetric BPS vortices: moving around magnetic flux quanta
We have just found new solutions of the self-duality equations in the linear approximation obtained by means of perturbations of the BPS cylindrically symmetric vortices induced by adding to them the ξ 0 ( x; n, k) zero modes, mutatis mutandis the ξ ⊥ 0 ( x; n, k) zero modes. Therefore, the following infinitesimal deformations of the self-dual rotationally symmetric n-vortex are still self-dual solutions of the vortex equations: or, in complex variables in field space, ψ = ψ 1 + iψ 2 and V = V 1 + iV 2 , ψ( x; n, k) = f n (r)e inθ − r n−k−1 h nk (r) f n (r) e ikθ V ( x; n, k) = i n r 2 β n (r) e iθ + i r n−k−1 h nk (r) e i(k+1−n)θ .
The main properties of the perturbed self-dual solutions are unveiled by studying the behaviour of the Higgs field near the origin: ψ( x, n, k) = d n r n e inθ − (2k + 2)c (n,k) 2k+2 d n r k e ikθ = r k e ikθ d n r n−k e i(n−k)θ − (2k + 2)c (n,k) 2k+2 d n which tells us that the origin is a zero of ψ( x, n, k) with multiplicity k and that other n − k zeroes appear characterized by the conditions where we have used the fact that h nk (r) is a decreasing function near the origin. Therefore, n − k zeroes located at the origin in the rotationally symmetric self-dual vortices migrate under this perturbation to the positions r 0 (n, k, j)e iθ(n,k,j) = (2k + 2)|c (n,k) 2k+2 | d 2 n 1 n−k e i (2j+1)π n−k , j = 0, 1, . . . , n − k − 1 infinitesimally apart from the origin. In sum, the perturbation induced by the zero mode ξ 0 ( x; n, k) on a self-dual n-vortex leaves k quanta of magnetic flux at the origin but moves n−k Higgs field zeroes away to be placed on the vertices of an infinitesimal regular (n−k)-polygon, the (n−k)-roots of the unit multiplied by an infinitesimal factor times e i π n−k . Of course, the perturbation produced by the ξ ⊥ 0 ( x; n, k) zero modes gives rise to a phase change in the new positions shifted by − n 2 π: θ ⊥ (n, k, j) = θ(n, k, j) − π 2(n−k) . Although the locations of the zeroes of the Higgs field of a self-dual n-vortex obtained as a zero mode fluctuation of a rotationally symmetric vortex capture the essential features of these stringy objects, it is also interesting to look at the associated magnetic field. The magnetic fields of a rotationally symmetric vortex perturbed along the ξ 0 ( x; n, k) and ξ ⊥ 0 ( x; n, k) zero modes are: Because the zero mode fluctuations of rotationally symmetric self-dual vortices are still solutions of the first-order equations, their energy density per unit surface area reads which is exactly the magnetic field of the deformed solution, as it should be the case. We emphasize that the information contained in formulas (37)-(39) describes the neutral equilibrium fluctuations of any BPS cylindrically symmetric n-vortex not merely as the free motion of the Higgs zeroes, the vortex centers, but also tells us how the scalar, vector and magnetic fields are deformed in these processes in a manner as precise as permitted by the limited analytical knowledge of the BPS multivortex solutions. It is clear from (22) that the construction of the 2n zero modes of fluctuations of a rotationally symmetric self-dual n-vortex is limited to the numerical computation of the k = 0, 1, . . . , n − 1 radial form factors h nk (r).
In the previous subsection, analytical expressions of these functions were found near r = 0 and very far from the origin. It remains to describe the vortex zero modes in the intermediate region by interpolating between these two regimes. Thus, we search for numerically generated solutions h nk (r) using the same procedure as that previously employed to find the self-dual vortex profiles f n (r) and β n (r). Taking into account that (24) is a homogeneous linear differential equation, we fix the value of the function at some determined point r = r 0 , i.e., h nk (r 0 ) = c 1 and allow the value of its derivative at this point h nk (r 0 ) to vary. Using a shooting numerical technique, we tune the value of h nk (r 0 ) to obtain the exponentially decaying solution, which defines h nk (r) as a normalizable function, otherwise the solution goes to infinity. This procedure allows us to construct h nk (r) in the [r 0 , ∞) range. Starting from the same integration constants and the same scheme with a negative step, we obtain a numerical profile of the function h nk (r) also in the [0, r 0 ] range. The numerical results confirm the theoretical behavior for small and large values of r derived analytically.
We shall now offer graphic representations of the self-dual vortex zero modes ξ 0 ( x; n, k) and the corresponding perturbed fields ψ and V defined in (10) or (37) with the aim of gaining some insight into the meaning of the new solutions. We warn the reader that in order to reach an appreciable visualization we shall choose a finite instead of taking some infinitesimal value. Graphics describing the 5-vortex zero modes ξ 0 ( x; n, k) for every value of the angular momentum k = 0, 1, 2, 3, 4 are collected in Figure 1, which arranges the relevant plots in a table format in order of decreasing value of k. Vector plots of the scalar and vector field zero mode fluctuations ϕ( x) and a( x) are displayed respectively in the first and third columns of Figure 1. In the second and fourth columns of Figure 1, however, we insert respectively vector plots for the scalar field ψ and vector potential V of the perturbed vortex solutions together with a density graphics of the value of its moduli. Here the darker the color is the less the value of the modulus is in such a way that the single vortex centers emerge as shadowed regions in the plots. The first row in Figure 1 illustrates the profile of the zero mode ξ 0 ( x, 5, 4) of angular momentum k = 4 and the behavior of the corresponding perturbed fields. Under this zero mode fluctuation one of the zeroes of the n = 5-vortex scalar field ψ( x) moves along the x 1 -axes, i.e., one of the five single quanta forming the BPS 5-vortex configuration segregates from the other ones, which remain located at the origin of the plane. The corresponding perturbed vector potential V ( x) is shown in the last Figure plotted in this row. The differences between the perturbed and non-perturbed vector potentials are not so prominent as for the companion scalar fields. In the second row, Figure 1, the zero mode fluctuation of angular momentum k = 3, ξ 0 ( x, 5, 3), and its influence on the BPS 5-vortex solution is displayed. It is seen in the perturbed scalar field that two of the vortex quanta, initially amalgamated in the 5-vortex center, separate following opposite senses along the x 2 direction from the remaining ones which stay at the origin. The fluctuation for angular momentum k = 2 is illustrated in the third row, Figure 1. Graphical representations of the zero mode ξ 0 ( x, 5, 2) scalar and vector fields are respectively shown in the first and third columns. Under this zero mode fluctuation of the original 5-vortex configuration only two constituent quanta remain fixed at the origin while the other three magnetic flux quanta are ejected in directions whose relative angles are 2π 3 radians. Next, fourth row, Figure 1, the zero mode ξ 0 ( x, 5, 1) of angular momentum k = 1 fluctuation fields are plotted. Four quanta of magnetic flux are ejected following the diagonal lines of every quadrant in the plane, while the remaining single quanta stays at the original location. Finally, in the last row and second column in Figure 1, the cylindrically symmetric self-dual vortex solution perturbed by the action of the k = 0 zero mode scalar field fluctuation is plotted. We observe that all the five single quanta in the original configuration are expelled in the directions determined by the vertices of a regular pentagon.

Boson-vortex bound states. Internal fluctuation modes on BPS cylindrically symmetric vortices
We shall investigate now the existence of excited fluctuation modes in the pure point spectrum of the H + -operator. With respect to the vortex zero modes an important difference is the need of solving second-order differential equations, instead of copying with first-order differential equations systems. To overcome this difficulty we shall take profit from the supersymmetric structure of the variational problem around BPS vortices in the Abelian Higgs model.
Searching for positive eigenfunctions ξ + λ ( x) of the operator H + , we shall impose the background gauge condition in order to eliminate spurious gauge fluctuations.
Proposition: Besides the 2n zero modes, the spectrum of the operator H + governing the fluctuations of a rotationally symmetric BPS n-vortex also includes excited (positive) orthogonal eigen-modes of two generic classes: 1. Class A. There is a class of fluctuations belonging to the strictly positive spectrum of H + of the form: Fluctuations of type ξ + λ ( x, n, k) are the k-term of a cosine Fourier series, while χ + λ ( x, n, k) arises in a sine Fourier series. Moreover, ξ + λ ( x, n, k) and χ + λ ( x, n, k) are linearly independent. The radial form factor v nk (r) is determined in both cases as a solution of the 1D Sturm-Liouville problem There are excited modes of fluctuations of the form (42)-(43) built from the solutions of the radial equation (44) that belong to either the continuous or the discrete spectrum of H + .
• Bound states: there exists a discrete set of eigenvalues ω 2 j confined to the open interval ω 2 j ∈ (0, 1) such that the associated eigenfunctions of the form (42) and (43) are normalizable. In order to belong to L 2 (R 2 ) ⊕ R 4 these excited eigenfunctions must comply with the boundary conditions dv nk dr (0) = 0 and lim r→∞ v nk (r) = 0. • Scattering states: all the eigenvalues ω 2 λ (p 2 ) = p 2 + 1, where p is a positive real number, p ∈ R + − {0}, are admissible. The corresponding eigenfluctuations belong to the continuous H + -spectrum if the form factor v nk (r) satisfies scattering boundary conditions. Thus, the continuous spectrum {p 2 + 1} p∈R + −{0} emerges from the threshold value 1.
• Zero modes: thoroughly described in the previous Section.
• Scattering states: h nk (r) behaves as an scattering function at infinity.
Note that in this class of fluctuation regularity of the wave functions at the origin restricts the angular momentum to be 0 ≤ k ≤ n − 1.
Proof: In order to prove the Proposition stated in the previous subsection we shall exploit the SUSY quantum mechanical structure hidden in the operators D, D † , H + , and H − . Apart from the zero modes, the operators H + and H − are isospectral. If ω 2 λ > 0 the SUSY structure implies that the eigenfunctions of H come in pairs and are related through the supercharges, i.e., and our strategy in the search for positive eigenfunctions of H + will be to solve the spectral problem of H − and apply the D † operator to ξ − λ to find the eigenfunctions of H + . The reasons are twofold: First, the kernel of H − is null. Second, H − is block-diagonal: It consists of two 1 × 1 blocks, prompting a complete decoupling of the vector field fluctuations, and one 2 × 2 block arising because of the decoupling of the scalar field fluctuations from the vector field perturbations. In sum, the spectral problem of H − is divided into three classes of positive eigenfunctions, each class associated with one of the diagonal sub-blocks. We shall now investigate these possibilities separately: • Class A: There exist H − -eigenfunctions of the form The corresponding SUSY partner H + -eigenfunction, which shares the eigenvalue ω 2 λ with ξ A− λ ( x), is: It is immediate to check that ξ A+ λ ( x) satisfies the background gauge:   (42) and (43) given in the Proposition for the excited modes of fluctuation ξ A+ λ ( x) and χ A+ λ ( x). The norm of the ξ A+ λ ( x) eigenfunctions, for example, is simplified by means of their SUSY nature: where C k = 1 if k = 0 and C k = 2 if k = 0. Normalizability of ξ A+ λ thus requires regularity at the origin, exponential decay at infinity, and a sufficiently smooth interpolation in between, to v nk (r). We shall now investigate these regimes: • Regularity of the función v nk at the origin: Near the origin, one solves for v nr (r) by applying the Frobenius method to the ODE (44) in the vicinity of the regular singular point r = 0, i.e., one plugs the series expansion v nr (r) = r s ∞ j=0 v (n,k) j r j near r = 0 starting from the condition v (n,k) 0 = 0 into the ODE (44). This procedure trades the ODE (44) for the two term recurrence relation where the coefficients v (n,k) j are grouped by powers. Annihilation of the j = 0-term induces the indicial equation −s 2 + k 2 = 0. The characteristic exponents are s = ±k but we shall choose s = k because our solutions are such that k ≥ 0. The recurrence relations between the odd coefficients arising from (51) means that these coefficients vanish at least up to the 2n + 1 coefficient. The recurrence relations (51) between even coefficients j = 2i, i = 1, 2, · · · imply that −4i 2i−2 for i = 1, 2, . . .. Therefore, near the origin we find the solution v nk (r) ≈ v (n,k) 0 . . which is regular and presents no problems in the normalizability of ξ A+ λ ( x) .
• Asymptotic behavior of the function v nk (r): The asymptotic behavior of the function v nk (r) near infinity is governed by the differential equation which is a Bessel or modified Bessel equation. Thus, if ω λ ∈ (0, 1) where I k , K k , J k and Y k belong to the Bessel function catalogue. The asymptotic behaviour of form factor v nk (r) is accordingly The genuine Bessel function J k and Y k exhibiting a non-normalizable oscillatory asymptotic behaviour are solutions of the above ODE when ω 2 λ > 1. A continuous spectrum arises in the ω λ ∈ [1, ∞) range, i.e., for energies above the scattering threshold ω 2 λ = 1. Below this threshold, in the ω 2 λ ∈ (0, 1) range, boson-vortex bound states may exist because of the exponential asymptotic behavior of the modified Bessel functions I k and K k . To find the eigenvalues in the point spectrum of both H ± it is necessary to identify a particular solution whose behaviour near the origin is regular as indicated previously and exhibit a decreasing exponential tail (C 1 = 0 in the asymptotic solutions given before). The difficulty in finding this type of bound state fluctuations is the identification of the discrete set of ω 2 j ∈ (0, 1) enabling their existence. In sum, in the discrete spectrum of H ± we expect to find SUSY pairs of positive modes of fluctuation of the generic form (46) belonging to the L 2 (R 2 ) ⊕ R 4 Hilbert space.
Note that the second choice in our strategy of studying the decoupled blocks in H − , i.e., the search for where a 2 ( x) must be a solution of the spectral ODE (−∇ 2 + |ψ| 2 ) a 2 ( x) = ω 2 n a 2 ( x), has already been discussed because of the discrete symmetry shown in the last Lemma. This symmetry connects the eigenfunctions ξ A− λ ( x) and ζ A−⊥ λ ( x) and consequently its SUSY partners. In particular, the SUSY partner H + -eigenfunction adopt the form such that the background gauge condition (47) becomes the PDE: (−∇ 2 + |ψ| 2 )a 2 ( x) = 0. The unique regular solution is a 2 ( x) = 0 and the background gauge discards this type of fluctuation mode. The perpendicular eigenfunction ζ A+⊥ λ ( x), however, automatically satisfies the background condition (47). But by comparing (46) and (52) we conclude that the result found encompasses exactly the same type of eigenfluctuations of H + as those described before because a 1 ( x) and a 2 ( x) are both solutions of the same equation (45).
• Class B. The lower 2 × 2 block-diagonal matrix inside H − acts only on the scalar field fluctuations. There exist another class of fluctuations of the form which are a priori candidates to belong to the spectrum of H − , paired with its SUSY partners in the spectrum of H + . The scalar fluctuations must be solutions of the H − -spectral PDE system: Assuming that eigenfunctions of H − are obtained by solving (54), the corresponding SUSY partner H + -eigenfunctions are: It is easy to check that the fluctuations (55) and (56) verify the background gauge (47) provided that ∇ 2 ψ 2 = ∇ 1 ψ 1 and ∇ 2 ψ 1 = −∇ 1 ψ 2 hold. But these first-order PDE are obeyed by the self-dual vortices, and hence ξ B+ λ ( x) and ξ B+⊥ λ ( x) are admissible as candidates to be H + -eigenfunctions. Now we shall demonstrate the lack of discrete eigenfunctions of the form (53) in the spectrum of H − and consequently in the spectrum of H + . We shall also analyze the connection of this type of eigenmodes with the zero modes: • Lack of discrete spectrum of H − and H + : It is easy to check that the 2 × 2 block inside H − acting on the scalar fluctuations can be factorized as in terms of the first order partial differential operator and its adjoint. The factorization (58) means that the eigenvalues ω 2 λ of the spectral problem (54) are greater than or equal to the threshold value 1. Thus, there are no positive H − -bound states of this class and only a continuous spectrum emerges. Nevertheless, the nature of eigenfunctions of eigenvalue exactly equal to the scattering threshold ω 2 1 = 1 is unclear at this point and demands a closer investigation.
The separation in the polar coordinates ansatz converts the PDE system (54) into the spectral ODE: with the radial form factor u nk (r) as unknown. This is a radial Schrödinger differential equation with an effective potential well: The asymptotics of V B eff respectively near r = 0 or infinitely far apart of the vortex core r = ∞ are: which reinforces the previous claim about the emergence of the continuous spectrum on the threshold value ω 2 λ = 1. Indeed the ω 2 λ = 1 u {1} nk (r)-eigenfunction belong to the kernel of the second-order differential operator R a relation immediately derived from the ODE (60) with ω 2 λ = 1. The operator R {1} admits a (Infeld-Hull) factorization in the form R {1} = L † L, where L and L † are the first-order differential operators It is thus clear that the eigenfunctions with ω 2 = 1 correspond to the u n,n−1 (r) tends to a constant at infinity and go to zero at the origin as r n . These modes are thus akin to the half-bound states arising in the fluctuations of sine-Gordon and φ 4 1D kinks, respectively tanh x and 3 tanh 2 x − 1. Note that the half-bound states tend respectively to ±1 or 2 at x = ±∞ and have either one node at x = 0 or two nodes at x = ±arctan 1 √ 3 , see e.g. Reference [35]. This type of modes with k = 0, 1, 2, · · · , n − 2 decays to zero at infinity like negative powers of r, at most as r −1 , and these states could be interpreted as a new type of "fractionary"bound states by some refinement of the 2D Levinson theorem.
• BPS vortex zero modes versus class B eigenfunctions: Now we shall unveil the connection between the zero modes studied in Section §.4 and the Class B eigenfunctions discussed in this Section. In order to perform the analysis just mentioned in a form as close as possible to the description of vortex zero modes elaborated in Section §.4, we express the function u nk (r) appearing in (60) in the form: u nk (r) = g nk (r) fn(r) . The ODE (60) becomes which is exactly the ODE (29) except for the term proportional to the frequency ω 2 λ . Therefore, one sees that the positive eigenmodes of clas B in the spectrum of H + belong to the same class as the BPS vortex zero modes.
The asymptotic behavior of the radial form factor g nk (r) is determined by the Bessel ODE −r 2 g nk (r) − rg nk (r) + [(1 + k − n) 2 + r 2 (1 − ω 2 λ )]g nk (r) = 0 such that: revealing the oscillatory asymptotic behavior of g nk (r), which confirms the existence of only continuous spectrum for ω 2 λ ≥ 1. Application of the Frobenius method to the differential equation (61) In this case the radial form factor h nk (r) satisfies the ODE which differs from the vortex zero mode equation only in the ω 2 λ -term. Regularity at the origin, however, of the positive eigenfluctuations ξ B+ λ ( x) and ξ B+⊥ λ ( x) requires that 0 ≤ k ≤ n − 1, as in the zero modes of fluctuation case. This proves the second, class B, proposal in the Proposition.
Finally it remains to prove the orthogonality of the eigen-modes of H + . Orthogonality between eigenfunctions belonging to different classes is guaranteed by the conservation of the scalar product in the SUSY partnership: because class A and B eigenfunctions of H − are clearly orthogonal. Orthogonality between eigenfunctions belonging to the same class with different angular dependence is established by the Fourier theory and, in general, Sturm-Liouville theory of the radial Schrödinger differential equations guarantees that eigenfunctions having different eigenvalue are orthogonal.

Portrait of bosonic bound states on BPS cylindrically symmetric vortices
In this last Section we shall try to elucidate the existence and to understand the structure of excited fluctuations of class A belonging to the discrete spectrum of H + with positive eigenfunctions lower than 1. Physically, these fluctuation modes obey to certain combinations of the scalar and vector field fluctuations trapped by the rotationally symmetric BPS vortex. Thus, fluctuations with these properties correspond to boson-vortex bound states and our goal is to identify some of these possible bound states as well as to investigate their properties. This task is more difficult than the search for vortex zero modes because we need to find numerically not only the eigenfunctions but also the discrete eigenvalues ω 2 j in the (0, 1) interval.

Numerical procedure for finding BPS vortex excited modes of fluctuation of bound state type
The search for and the analysis of some fluctuations in the discrete spectrum of H + reduces to the numerical computation of the radial form factor v nk (r) in the ODE (44). The existence of boson-vortex bound states can be physically understood analysing the properties of the effective potential wells which arise in the appropriate Schrödinger equation. These equations can also be understood as governing the quantum planar motion of a particle moving in a central potential. Accordingly, we expect to find bound states when the angular momentum is low enough. Our strategy to achieve this finding is to employ a second-order finite-difference scheme which simulates the differential equation (44) by the recurrence relations where we have confined the problem to the interval [0, r max ] for a large r max enough. We denote v (i) nk;j = v nk;j (i∆x), with ∆x = r max /N , and choose a mesh of N points with i = 0, 1, 2, · · · , N . The eigenfunction and the eigenvalue depend on the values of the angular momentum k and the quantized magnetic flux n. The index j labels the discrete eigenfunctions. The contour conditions are: nk;j = 0 A good estimation of the discrete eigenvalues ω 2 nk;j is obtained through diagonalization of the N × N matrix in the left member of the linear system (62). We show the eigenvalues of H + for the lowest values of n and k in Table 1. In particular, we have implemented the algorithm (62) in a Mathematica environment with a choice of two grids constructed respectively with N = N 0 = 4000 and N = 2N 0 = 8000 grid points. The results achieved manifest a satisfactory reliability shown in the data displayed in Table 1. Different figures in the displayed results corresponding to the two previously mentioned meshes have been emphasized by enclosing them in parentheses. We remark that we have chosen generically the value r max = 15 in the numerical scheme in such a way that the precision of the Laplacian operator in (62) is of order O(∆x 2 ) ∼ 3.4 × 10 −6 using the mesh with 2N 0 = 8000 points. Special treatment is required for states whose eigenvalues are near the scattering threshold value 1 where we choose r max = 50 such that O(∆x 2 ) ∼ 0.00004. The reason is that eigenfunctions close to the scattering threshold exhibit quite slow decaying asymptotic behavior, which demands a greater value of r max .
There is a logical structure almost so rich as the structure of vortex zero modes. In general, we observe that the number of bound states increases with the magnetic flux n. In particular, for magnetic flux n = 1, we find only one boson vector-vortex bound state. The eigenvalue of this state is (ω A 10,1 ) 2 = 0.77747, which arises for angular momentum k = 0. For magnetic flux n = 2, we find two boson vector-vortex bound states (ω A 20,1 ) 2 = 0.53859 and (ω A 21,1 ) 2 = 0.97303 whose Fourier wave numbers are k = 0 and k = 1 respectively. For magnetic flux n = 3 the situation is completely similar to the case n = 2 with eigenvalues (ω A 30,1 ) 2 = 0.402708 and (ω A 31,1 ) 2 = 0.83025. We find three bound states for the case n = 4, two of them corresponding to the angular momentum k = 0 with eigenvalues (ω A 40,1 ) 2 = 0.319288 and (ω A 40,2 ) 2 = 0.98835 while the remaining eigenvalue (ω A 41,1 ) 2 = 0.701767 corresponds to the value k = 1. This last situation is also suitable for the case n = 5 with the eigenvalues (ω A 50,1 ) 2 = 0.263679 and (ω A 50,2 ) 2 = 0.938456 and (ω A 51,1 ) 2 = 0.601272 although a fourth state with k = 2 arises with eigenvalue (ω A 52,1 ) 2 = 0.94252.
Eigenvalues of the discrete spectrum of H +  Table 1: Numerical estimation of the discrete spectrum eigenvalues within the class A eigenfunctions ξ A+ λ ( x, n, k) of the second-order n-vortex small fluctuation operator H + . These results have been obtained using the algorithm (62) with two different meshes of respectively N 0 = 4000 and 2N 0 = 8000 points.
Parentheses are used to point out the disagreeing figures arising from either the thicker or the thinner meshes. Table 2 summarizes all the spectral information of the radial eigenvalue problem (44) analyzed in this Section for vortex solutions with vorticity n = 5. In this Table we display the effective potentials V eff (r) = k 2 r 2 + f 2 n (r) (blue lines), the boson-vortex bound states energies (red lines) and the radial eigenfunction profile for several values of n (arrayed in different rows) and k (arrayed in different columns). Two types of behavior are distinguished: (a) If k = 0, we have a central potential with a minimum at the origin that tends to 1 at r → ∞. Above this threshold, ω 2 λ > 1, there is no question about the existence of scattering states in the spectrum of H + , but we also find one boson-n-vortex bound state for n = 1, 2, 3 with an eigenvalue below the scattering threshold decreasing with n. For the n = 4 and n = 5 BPS vortices, a second boson-vortex bound state appear. (b) If k > 0 the centrifugal barrier in the effective potential gives rise to a hard core but there is still a minimum of the potential away from the core if k is small enough. All together we observe a shallower well (that becomes a wall for a big enough k) that tends to infinity at r = 0 and to 1 when r → ∞. In the second, third, and fourth columns of Table 1 it is observed that when k grows we need a higher magnetic flux n for finding bound states.

BPS vortex bound state wave functions
The boson-vortex bound state eigenfunctions are found by simulating the ODE system (44) by the finite-difference equations system (62) and searching for the eigenvectors. In Table 2 Table 2: Graphical representation of the effective potential wells/walls V eff ( x, n, k) in the differential Schrödinger equation (44) Figure 2 in a Table format. Indeed we shall use this case as a paradigmatic example of the discrete bound state spectrum displayed in Table 1 because the qualitative behavior of the eigenfunctions carrying identical angular momentum k is similar for every vorticity n.
In the first row of Figure 2 the scalar ϕ A+ and vector a A+ fluctuations corresponding to the lowest eigenvalue (ω A 50;1 ) 2 = 0.263679 are plotted, as well as the perturbed fields ψ + ϕ A+ and V + a A+ corresponding to the bound state positive fluctuation ξ A+ 1 ( x; 5, 0) of a BPS 5-vortex centered at the origin. It is interesting to examine the details: (1) Both ϕ A+ ( x; 5, 0) and a A+ ( x; 5, 0) have zeroes at the origin, grow in the middle region and tend to zero again at r = ∞. The BPS 5-vortex does not vary under this k = 0 fluctuation at its center, but grows at the core's middle, and remain fixed again near infinity, see Figure 2 (first row). (2) We see this instantaneous picture as an "inflationary"process, a grow of the vortex density at middle distances from the origin. But it is in fact an oscillatory process of inflation/deflation due to the temporal dependence induced by the periodic time-dependent term: cos ω 50,1 t 7 .
The next excited vector-boson bound state in the discrete H + -spectrum is the eigenmode with Fourier wave number k = 1 whose eigenvalue is (ω A 51,1 ) 2 = 0.601272. The radial form factor v 51,1 (r) corresponding to this eigenfunction is shown in Table 2 (last row) while the ξ A+ 1 ( x; 5, 1) vector-vortex bound state is graphically depicted in the third row of Figure 2 (third row). We observe a new class of fluctuations. The scalar fluctuation ϕ A+ ( x; 5, 1) is stronger in a strip enclosing the x 1 -axis when it differs appreciably 7 It is enlightening to compare this bound state with the excited bound state of a meson trapped by λφ 4 -kink: φK (x) = tanh x and ϕ3(x) = sinh x cosh 2 x , ω 2 3 = 3. Clearly, the kink remains unchanged at x = 0, increases at x > 0, but diminishes when x < 0, and diminishes until φ 1, far away from the kink center on the right, but increases until φ −1 if x −1. All this happens at t = t0 but this frozen picture oscillates in time, with frequency ω3 = √ 3.
from zero for intermediate distances from the origin. The net effect on ψ( x, 2) is the squeezing of the 2-vortex along the x 1 -axis both from the left and from the right, but less intensely along the x 2 -axis in both senses, see Figure 2 (third row). Capturing the scalar and vector bosons in this bound state the BPS rotationally symmetric 5-vortex loses symmetry and gets thinner along the x 1 -axis. Of course, there is an oscillation of this static picture, with frequency ω A 51,1 = √ 0.601272 with a certain dipolar character. A similar process emerges for the vector self-dual 5-vortex field V ( x, 5) and the corresponding fluctuation a A+ ( x, 5, 1), see Figure 2.
In addition to the first mentioned fluctuation mode there exists other bound state in the H + spectrum with angular momentum k = 0 whose eigenvalue is (ω A 50,2 ) 2 = 0.93845. The perturbed fields behave following a similar pattern than that of the previously described k = 0-eigenmode, see Figure 2 (second row), although the radial profile function involves a new node, see the probability function ρ * 5,0 (r) in the last row of Table 2.
The highest eigenvalue in the discrete H + -spectrum takes the value (ω A 52,1 ) 2 = 0.94252. This peculiar state ξ A+ 1 ( x; 5, 2) is described in the last row of Figure 2. Focusing our attention in ϕ A+ ( x; 5, 2) we observe that the perturbation pressures the scalar field density concentrated in a disk inwards along the x 1 -axis but outwards along the x 2 -axis, inducing a quadrupolar density distribution oscillating in time.
Animations illustrating the previously mentioned behavior of the fluctuation eigenmodes together with the energy density and magnetic field associated with them can be downloaded at the web page http://campus.usal.es/∼mpg/General/Mathematicatools.

Brief summary and outlook
An ortho-normal basis of BPS cylindrically symmetric n-vortex zero modes in the kernel of the matrix second-order PED operator H + of fluctuation has been constructed and their mathematical and physical properties thoroughly described. Several positive normalizable eigenfunctions, bound states, of this Hessian operator, as well as their corresponding eigenvalues, have been identified within an unexpected level of precision. This novel type of BPS cylindrically symmetric positive vector boson-vortex bound state fluctuations exhibits very intriguing properties and deserves further investigation.
The techniques and procedures used in this paper suggest that the concepts developed and the results attained in this basic superconducting setting will work in more sophisticated superconducting media provided that a BPS or self-dual structure is available. The first system of this type that comes to mind is the Chern-Simons-Higgs model, see [34]. One expects a similar structure of the self-dual Chern-Simons vortex zero mode fluctuation, demanding only an appropriate radial form factor h nk (r) to be estimated by the same method. More dubious, however, is the prospect of finding some boson-vortex bound state in this necessarily planar model. Another possible scenario to play around with these ideas is the broad domain of U (1)-gauged massive non-linear sigma models, see [36,37]. Severely constrained by the possibility of building extended N = 2 supersymmetry on these systems, the main difficulty in investigating vortex zero modes and bound states in this framework appears to be the precise development of an index theorem in this context. Even more tantalizing, the development of a similar analysis about the zero modes of the BPS defects in the system discussed in [38] seems to be plausible, even though the high derivative terms look harmful.