Small Winding-Number Expansion: Vortex Solutions at Critical Coupling

We study an axially symmetric solution of a vortex in the Abelian-Higgs model at critical coupling in detail. Here we propose a new idea for a perturbative expansion of a solution, where the winding number of a vortex is naturally extended to be a real number and the solution is expanded with respect to it around its origin. We test this idea on three typical constants contained in the solution and confirm that this expansion works well with the help of the Pad\'e approximation. For instance, we analytically reproduce the value of the scalar charge of the vortex with an error of $O(10^{-6})$. This expansion is also powerful even for large winding numbers.


Introduction
A significant feature of many gauge theories is the existence of topological solitons which may appear when the gauge and/or global symmetries are spontaneously broken. Monopoles, vortices and domain walls are by now familiar, and have found important applications in vast areas of modern physics, such as cosmology, condensed matter physics and particle physics. From another view point, the topological solitons can be seen as nontrivial solutions of nonlinear differential equations. A direct way to study topological solitons is solving such nonlinear equations exactly. For instance, a beautiful systematic method to construct exact solutions for instantons has been well-established and is widely known as the ADHM construction [1]. This is, however, a special case and for many other types of solitons numerical calculations are needed to study solutions.
The present work concerns the so-called Abrikosov-Nielsen-Olesen (ANO) vortex as the simplest topological soliton with finite energy in the (1 + 2)dimensional theory. This vortex appears as a topological defect [3] in Ginzburg-Landau theory [2] and may be viewed as a static solution to the equations describing the 1+2 dimensional Abelian Higgs model [4]. In this theory all the vortex features depend on one dimensionless parameter λ = m s /m v 1 : the ratio of the Higgs boson mass m s to the vector boson mass m v . The intervortex force is, roughly speaking, a superposition of an attractive force caused by the Higgs boson and a repulsive force caused by the vector boson as seen in a scalar potential [10] for a well-separated pair of vortices with a large distance R. Here, q s and q m stand for a vortex scalar charge and a magnetic dipole moment, respectively. Therefore the force with the longest correlation length is dominant and the vortices attract (repel) each other for λ < 1 (λ > 1) [6]. The critical coupling λ = 1 is a rather special case where net intervortex forces are exactly canceled thanks to the coincidence of the two coefficients, q s = q m ≡ 2πC 1 . From a mathematical viewpoint, the Euler-Lagrange equations reduce to the first order differential equation called the Bogomol'nyi-Prasad-Sommerfield (BPS) equations for vortices saturating Bogomol'nyi bound, whose total energy is quantized as E k = |k|πv 2 with the winding number k ∈ Z. In this critical case, the constant C 1 appears, for instance, in a potential for a pair of moving vortices [12] U λ=1 (R, u) ≃ πv 2 × C 2 with a relative velocity u, since only the magnetic field accepts a Lorentz boost and the two forces are not canceled out. Unlike the remarkable cases of instantons and monopoles, no analytic solutions for this BPS equation in flat spacetime have been found even at this critical coupling. Thus only a few quantities are exactly calculable and a detailed study of the vortices, for instance, the calculation of a value of C 1 requires numerical analysis. In this paper, to complement the numerical analysis, we propose a simple and straightforward, but new idea for analyzing vortices at critical coupling, where fields are expanded perturbatively with respect to the winding number k ∈ Z around its origin k = 0. To justify this perturbative expansion, (let us call it "small winding-number expansion"), the BPS equations must be extended so that they allow a real winding number k ∈ R. Since the BPS equations with an infinitesimal winding number |k| ≪ 1 can be exactly solved, we can systematically perform perturbation calculations without tuning any parameters and this perturbative expansion is supposed to work well as a practical tool. Here, we calculate values of three typical quantities with λ = 1 including C 1 as the most simple examples to check this idea.
The constant C 1 has often been calculated in the literature. De Vega & Schaposnik [5] gave a semi-analytical study for axially-symmetric solutions with an arbitrary winding number k ∈ Z >0 , and constructed power-series expansions around a center of a vortex and asymptotic expressions for the opposite side. These two can be determined by only one constant D k+1 k for the power-series expansion and C k (Z k in their notation) for the asymptotic expression. Comparing these parameters in a middle region, they obtained the values: C 1 = 1.7079... and D 2 1 = 0.72791... . These values now seems to be widely accepted in literature, for instance, C 1 = 1.7079 appears in Refs. [7,11,13,15] and also in a standard textbook of Vilenkin & Shellard [8]. However, we encounter a different value for C 1 : C 1 = 10.58/2π ≃ 10.57/2π ≃ 1.682 ∼ 1.684 which was obtained by Speight [10] about twenty years later than de Vega & Schaposnik [5]. Furthermore, Tong [11] gave the supergravity prediction C 1 = 8 1/4 ≃ 1.68179... which seems to agree well with Speight's C 1 . These values also seem to be accepted in literature, for instance, Ref. [12] and another standard textbook by Manton & Sutcliffe [9]. There exists a 1.5 % discrepancy between old and new results.
In Sec.2.6, we shall conclude that the correct value is the old one C 1 = 1.7079 by using two different kinds of numerical calculations with higher accuracy. In Sec.4.2.3, we reproduce this value by using the small winding-number expansion to verify its power. This paper is organized as follows. In Sec.2, we review the BPS vortex in the Abelian-Higgs theory, and define an extended vortex function which allows the winding number of non-integer, as a solution of the BPS equations. There a non-trivial integral formula including the vortex function is derived and three typical constants C ν , D ν and S ν for a vortex solution are introduced and their analytical and numerical properties are discussed. In Sec.3 we perform a small winding-number expansion of the vortex function and the three constants using Feynman-like diagrams. Results obtained there are modified in Sec.4, using the Padé approximation to overcome problems with finite convergent radii of the expansions. Summary and discussion are given in Sec.5, and some useful inequalities and details of the calculations are summarized in the Appendices.
where φ is a complex scalar field, metric is η µν = diag.(+1, −1, −1) and covariant derivative is which has a vacuum |φ| = v where the U (1) gauge symmetry is spontaneously broken. The Higgs mechanism makes the scalar and the gauge fields massive. Their masses are given by, m s = λ ev, m v = ev respectively. The spontaneously broken U (1) symmetry gives rise to a soliton which is topologically stable object supported by π 1 (U (1)), of which element is called a winding number. To require vanishing of the kinetic term |D i φ| 2 = 0 at the spatial infinity connects this winding number with the first Chern class This topological defects are called the Abrikosov-Nielsen-Olesen vortex.
In this paper, we take the critical coupling constant, λ = 1, as the simplest model, where the two masses are identical, m v = m s ≡ m. Then we can perform the Bogomol'nyi completion of an energy density H for static configurations as and a total mass (tension in higher dimension) of vortices, T , has a lower bound The inequality is saturated by BPS states which satisfy the BPS equations Without loss of generality we will consider the BPS equations with the upper sign. In order to find general solutions of the BPS equations, it is useful to solve the second equation in Eq. (2.6) at first and, it can be solved with the complex coordinate z = x 1 +ix 2 and introducing a smooth real function ψ reg = ψ reg (z,z) where an arbitrary holomorphic function P (z) can be set to be a monic polynomial without loss of generality. Here zeros {z I ≡ x 1 I + ix 2 I ∈ C} of the Higgs field φ are topological defects and identified as positions of vortices. One of the important features of BPS vortices is that they feel no interactions since the attractive and repulsive force are exactly canceled. So we can put BPS vortices anywhere as many as we like. Note that the smooth field ψ reg must behave as ψ reg ≈ log |P (z)| 2 at the spatial infinity to obtain a finite energy, it is convenient and more familiar to rewrite ψ reg in terms of a singular field ψ, so that ψ vanishes at the spatial infinity. With this singular field, then, the first equation in Eq. (2.6) can be rewritten to be, so called, Taubes' equation (2.10) Here we used that the magnetic field can be rewritten as, which coincides with Eq.(2.3) and k is the total winding number. Existence and uniqueness of a solution for Taubes' equation with a given arbitrary J have been established by [14]. With this solution, therefore, we obtain a complete solution for φ and A i . In terms of a solution of ψ and the source J, the energy density H BPS for BPS vortices can be rewritten to which gives the lower bound in Eq.(2.5). There is, however, no known exact solution for this equation, even in the simplest case with k = 1.

Extension of Taubes' equation and particle description
In a case that k I vortices coincide at x = x I for each I, the source terms are replaced with where k I indicates the winding number at x = x I . A request that the winding number k I is positive integer is to give the single-valued Higgs field φ and Profiles of ψ and the magnetic field in Eq.(2.11) and the energy density in Eq.(2.12) can be calculated without constructing φ. If we omit constructing φ, therefore, we can formally extend Taubes' equation with the generalized source terms (2.14) Here the winding number k I is renamed ν I to stress that ν I can be non-integer and the lower bound of the winding numbers will be discussed in Sec.2.4. A 'total mass' of this extended object is formally calculated as which is no longer an element of π 1 (U (1)). In the rest of this paper, we will study this extended Taubes Furthermore we assume that the solution is differentiable with respect to {ν I }.
Under this assumption, we can derive, for each I, from Taubes equation Eq.(2.9) with the source Eq.(2.14). According to Appendix A.1 the above equation show that the solution ψ is strictly increasing with respect to each ν I , ∂ψ/∂ν I > 0. In the limit of the vanishing source J = 0, furthermore we find where the modified Bessel function of the second kind K 0 (x) emerges as a twodimensional Green's function. That is, in this limit a vortex solution is exactly solved and treated as a linear combination of free massive particles and for small |ν I | ≪ 1 at least, ψ is approximated well everywhere as This is the starting point of the small winding-number expansion which will be discussed in Sec.3.
In this particle description, it will be convenient to rewrite Taubes' equation as with σ[ψ] as dimensionless self-interaction terms. Then, by applying the Green's function method to Taubes' equation, we obtain 2 an integral equation for ψ with Green's function G( x) = K 0 (m| x|), Since σ[ψ] ≥ 0 and K 0 (x) > 0 are always hold, we find that the solution of Taubes' equation must satisfy a fundamental inequality

Scaling argument and a physical size of a vortex
Let us consider the following Lagrangian in a two-dimensional Euclidean spacetime which induces Taubes' equation as an equation of motion of ψ, and an action 3 is where K ghost is introduced to cancel UV divergences of the kinetic term and the source term and we set K ghost as, for instance, After this regularization we can apply the scaling argument to this action. For simplicity, let us consider an axially symmetric case with the source J = πv 2 At the limit m 0 → 0, 1/m 0 gives a IR cut-off and K ghost can be eliminated by Kähler transformations. Actually one can confirm that the above Kähler potential gives Samols' metric [17]. 4πνδ 2 ( x). Since K is a dimensionless quantity, the dimensional argument tells us By using equations of motion for ψ and ϕ, derivatives of K with respect to masses can be calculated by where we used Eq.(2.16). Therefore, we find the following formula [18] d 2 x ψ = 2π m 2 × ν(ν + 2). (2.30) As we seen the above this exact formula does not come from topological argument, but from the scaling argument. To check numerical calculations we use this formula in this paper. Thanks to this non-trivial identity combining Eq.(2.16), the following integral is calculated as and a size of the vortex with the positive winding number ν > 0 can be naturally defined with the energy density H BPS given in Eq.(2.12) and calculated as, which turns out to be a key point in Sec.4. It is natural for the scaling argument to determine a typical size of a soliton.

Axially symmetric solution
Let us consider a single vortex sitting the origin with the winding number ν, that is, we consider a solution with the source term J = 4πνδ 2 ( x). Its configuration is axially symmetric and described by a function ψ = ψ(mr, ν) with respect to a radial coordinate r = | x| and the winding number ν. The partial differential equation (2.9), therefore, reduces to an ordinary differential equation with the following two boundary conditions Even for the non-integer number ν, a set of the differential equation and the boundary conditions defines an unique solution under assumption of its existence. Especially for small |ν| ≪ 1, ψ is approximated in the full range of r ∈ R >0 as We just assume the existence of the solution with ν > −1 in this paper. Note that we can show the following inequalities although we have no exact solution. Applying the discussion in Appendix.A.1 to Taubes' equation with the source J in Eq.(2.14), we find the solution ψ must be positive for ν > 0 and be negative for ν < 0, and Eq.(2.33) tells us that r dψ dr is strictly increasing (decreasing) with respect to r for ν > 0 (ν < 0), and therefore the boundary conditions Eq.(2.34) give lower and upper bounds as, According to Appendix A.1, the following inequality implies that ψ is a downward-convex function, Combining Eq.(2.35) with this fact, we find that ψ/ν is strictly increasing with respect to ν and furthermore we obtain where the reduced Green's function G F (r, s) takes the following form with the step function Θ(x) and the modified Bessel function of the first kind I 0 (x).

Observable parameters
To define the solution ψ of Taubes' equation even with the positive non-integer winding number ν, we have to consider a behavior of the solution around the core of the vortex seriously. Note that in the massless limit m → 0, Taubes' equation has a general solution 4 with a positive real arbitrary constant R in as, and with the finite mass m > 0, therefore, ψ can be expanded by m and we find an expansion of ψ around the origin r = 0 in an unfamiliar form, where we treated mr and Y as if they were independent of each other, and a function F n (Y ) is independent of m and turns out to be a polynomial of order n with respect to Y determined sequentially by solving Taubes' equation as, which must vanish in the limit ν → 0 for a finite radius r due to Eq.(2.17). The dimensionless constant D ν appeared in the expansion is related to R in as 5 (2.47) Therefore the expansion of ψ can be defined by a pair of parameters {ν, R in }.
The uniqueness of the solution with a given ν means, however, that to satisfy the boundary condition at the spatial infinity, the constant R in must take a certain value corresponding to each value of ν, that is, a function R in = R in (ν), otherwise a function defined by the expansion always glows up at a large r. In Appendix A.2 this feature is analytically discussed and at the present we find a pair of lower and upper bounds of R in as According to Eq.(2.40) R in and D ν /ν turn out to be strictly increasing functions with respect to ν and take values at with Euler's gamma γ. In Fig.2, we plot a profile of D ν /ν. Note that there is an another way to calculate D ν using the integral form Eq.(2.41) as, These different two definitions of D ν will be used to double-check numerical calculations of D ν . Since the axially symmetric vortex solution we consider has the only one mass parameter m, we expect that the dimensionfull parameter R in controlling a profile of the solution should be the same order of the vortex size R BPS given in Eq.(2.32). Thanks to Eq.(2.48), roughly speaking, we find actually R BPS ≈ R in for large ν. We call R in an internal size. On the other hand D ν is directly related to a value of the action K with J = 4πνδ 2 ( x) in the previous subsection. In the same way of Eq.(2.29), we can calculate a derivative of K with respective to ν, and by setting the mass of the ghost m 0 to be m 0 = 2e −γ m, we obtain the following simple relations, Let us take m large conversely, that is, consider a infrared region r ≫ R in ≈ 2 √ ν/m. There, an asymptotic behavior of ψ can be treated as a fluctuation of a free massive scalar field around the vacuum. Due to the axial symmetry, such a fluctuation is written with a certain constant C ν ∈ R >0 as There is the similarity between this asymptotic form and the form of Eq.(2.35) and the uniqueness of the solution of Taubes' equation indicates that the two constants C ν and ν are in one-to-one correspondence. Actually, to satisfy the boundary condition at the origin r = 0, the constant C ν must be a function with respect to ν and according to Eq.(2.17), Eq(2.19) and Eq(2.39) we find These property tell us that C ν /ν is strictly increasing with respect to ν and a lower bound of C ν is given as C ν > ν. A profile of this function is shown in Fig.3. According to the integral equation Eq.(2.41), C ν can be calculated by Bringing this identity back, we can remove the explicit ν-dependence from the integral equation Eq.(2.41) as A one of purposes of this paper is to confirm the true value of C 1 .

Total scalar potential S ν
Finally let us consider the following definite integral 7 which is dimensionless and proportional to a total potential energy of the Abelian-Higgs model at critical coupling, This quantity with ν > 0 satisfies Thanks to Eq.(2.40) we find that S ν is also an increasing function with respect to ν and according to the profile of S ν shown in Fig.4 an 'energy' per an unit winding number S ν /ν is also an increasing function with respect to ν, and this property gives This inequality is consistent with the well known property of type II (type I) vortices, that is, intervortex forces are repulsive (attractive) for the coupling λ > 1(λ < 1) 8 .

Numerical Data
We numerically calculate values of C ν , D ν , S ν in most of the range of ν as ν = 1/20, 1/10, · · · , 500, 1000 using mainly the shooting method. These data are listed in Table.1. We will denote these data as N sht [C ν ], N sht [D ν ] and N sht [S ν ] for C ν , D ν , S ν respectively. In Sec.4, we use these data as references to show how the winding-number expansion introduced in Sec.3 works well. The other purpose of this subsection is to settle the problem on the numerical value of C 1 . We need, therefore, numerical calculations with high accuracy. To show accuracy of our numerical data to readers, let us enter into details of the numerical calculations we performed. Note that there exist two kinds of strategies in the shooting method and we observe a big difference in usability between them. We calculate numerical solutions of ψ in a region {r |ǫ ≤ r ≤ L} where we set m = 1 and take ǫ = 10 −2n+1 and L = 2 √ ν + p log 10 with p, n = 8 ∼ 9 referring to the flux size R flux given in Eq.(2.32). The first strategy is to take r = ǫ as the initial point of the calculation and fine-tune the parameter D ν so that a profile of ψ satisfies the boundary condition at r = L and read C ν from a profile of ψ 8 It is natural to expect the following inequalities on values of total energies E k for axiallysymmetric vortex-solutions, at r = L. Since the initial conditions are given by a pair {ν, D ν }, an incorrect pair always makes a profile function blow up at large r. The second one is to take r = L as the initial point and fine-tune the parameter C ν so that ν = −(rψ ′ /2) at r = ǫ and read D ν at r = ǫ. In this strategy the profile function is controlled by the only one initial parameter C ν which is related to ν in one-to-one correspondence thanks to dC ν /dν > 1. With the sufficiently large L, therefore, a profile function with an arbitrary C ν always gives a certain solution corresponding to a certain ν, without the profile blowing up, and thus this strategy gives a function ν = f (C ν ). Thanks to this property, it is easy to create a computer program for tuning C ν automatically with a given ν and arbitrary precision. We take the second strategy in this paper although the first strategy was taken 9 in Speight's paper [10]. As we explained above, numerical data N sht [C ν ], N sht [D ν ] for C ν , D ν are directly obtained. To double-check those data, we also use the integral formulas Eq.(2.55) and Eq.(2.50) for C ν , D ν respectively, to obtain different data as errors of these data and plot them in the right panel of Fig.5. For instance, we obtain as double-checked numbers, for ν = 1 and the numerical data listed in Table.1 have been double-checked in this sense. Therefore we conclude that the numerical result C 1 = 1.7079 given by de Vega and Schaposnik is correct. Thanks to the non-trivial identity in Eq.(2.30), we can estimate accuracy of the profile functions itself by calculating the following quantity and we plotted this in the right panel of Fig.5. Note that we observe that the precision of N sht [C ν ] generally get worse than those of δ, D ν as shown in Fig.5.
The precision of calculations in Speight's paper seems to be less than six digits and we guess that his result C 1 ≈ 1.683 has an error of O(10 −2 ) ∼ O(10 −3 ) which is consistent with the other numerical results including ours. We obtain also a stable numerical value of S 1 with long digits by the shooting method. To perform double check of the values of S ν , we also use the relaxation method as the other numerical calculation. In the relaxation method, we introduce a relaxation time τ and extend ψ( x) to be dependent on τ , ψ = ψ( x, τ ), and modify the equations of motion by adding a friction term ∂ψ/∂τ with an appropriate signature. With an appropriate initial function of ψ, ψ(r, τ = 0) = 2νK 0 (r) for instance, this friction term defines the time evolution of ψ and decreases an 'energy' of this system defined in Eq.(2.26).
In principle, therefore, the true solution could be obtained with an infinite τ as ψ(r) = lim τ →∞ ψ(r, τ ). As larger τ , we will get better accuracy in many cases. In reality, beyond a certain finite τ , we observe stability of values of the observables with small noises, since those accuracy can not be better than the calculation accuracy. For instance we stopped the time evolutions at τ ≈ 4 × 10 4 . The relaxation method is convenient and powerful to solve (simultaneous) nonlinear (partial) differential equations numerically. We need no fine-tuning of any parameters there. In the simple system we are considering, however, the shooting method is more powerful to get precision. Generally speaking, numerical data N rlx [X] for X = C ν , D ν , S ν calculated by the relaxation method get worse precision as shown as Fig.5

Small Winding-Number Expansion
In the paper [5], de Vega and Schaposnik calculated C 1 and D 1 by a semianalytical study. Their strategy was essentially as follows. Let us divide the integrals in Eq.(2.50) and Eq.(2.55) as The former integral is calculated by inserting the expansion Eq.(2.44) which depends on D ν and the latter is calculated by the expansion Eq.(2.58) which depends on C ν . Then we obtain simultaneous equations for C ν and D ν , and thus, approximate the values of C ν , D ν as their solution.
In this section we will give a different expansion of the solution ψ using Eq.(2.41) and calculate them more straightforwardly and more systematically.

ν-expansion of the vortex function ψ
In the normal case, we can not define an expansion of ψ with respect to the winding number as a topological quantum number. In the previous section, we relax the winding number ν from an integer to a real number and assume smoothness at ν = 0, and thus, we can consider a Taylor expansion of the solution for ψ with respect to the winding number as, with Eq.(2.17) where expansion coefficients σ n = σ n ( x) in the interaction terms σ[ψ] are The absence of the solution for ν ≤ −1 shown in Eq.(2.36) might indicate that a radius of convergence for the ν-expansion of ψ is less than 1. In Sec.4, we will discuss that this fact is not a big problem. We can perform calculations of the ν-expansion of ψ with the familiar technic using Feynman diagrams. The ν-expansion of ψ( x) is given concretely as using conventions for Feynman diagrams, Here diagrams of the order n have n external legs coming from the point-like vortex at the origin x = 0.

E n [C ν ]
Let us approximate C ν analytically by using the ν-expansion, In principle, its coefficients c n can be obtained by taking the ν-expansions of the both sides of Eq.(2.55) and inserting ψ n obtained in Eq.(3.7) into the right hand side. Comparing Eq.(2.41) and Eq.(2.55), however, we find that the coefficient c n can be calculated by only replacing the propagator with I 0 (m| x|) as where the triangle symbol stands for I 0 (m| x|) = . (3.11) For instance, coefficients c 2 , c 3 are calculated as See Appendix B for details. Finally we obtain As shown in Fig.6, we observe that as the order n is larger, an error of E Unfortunately the accuracy of this value is worse than that of the value C 1 ≈ 1.7079 given by de Vega and Schaposnik. According to Fig.6 a radius of convergence of the infinite series, ν c , is obviously finite and smaller than ten, ν c < 10 and we can not judge whether ν c is larger than one or not. In Sec.4, we will overcome these problems.

E n [D ν ]
Next, let us consider the ν-expansion of D ν , According to Eq.(2.50), the expansion coefficient d n for n ≥ 2, is calculated by reducing diagrams in Eq.(3.7) as, We find therefore, by performing integrals numerically, Note that this quantity is known to have the lower bound 4e −1 ν and the second coefficient is near to this bound as 1.47777 > 4e −1 = 1.47152. Finite series were expected to be good approximations, but we find their slow convergence as seen in Fig.2.

The ν-expansion of the formula Eq.(2.30)
To check consistency of the ν-expansion of the formula Eq.(2.30), we need some unfamiliar formulas. There is a non-trivial identity as, (3.21) and using Eq.(B.1) we find Using the above formula, we also find with σ 1 = 4πδ 2 (x)/m 2 , and since ψ 3 is a dimensionless quantity we can confirm To check Eq.(2.30) for more higher order, similarly we must need the dimensional argument again. Checking Eq.(2.30) is, therefore, tautological in this sense.

E n [S ν ]
To calculate the ν-expansion of S ν , at first we rewrite the definition of S ν by inserting the identity in Eq.(2.31) Here we canceled a ψ 2 term to avoid complicated and redundant calculations such as those in Sec.3.4, and thus, substituting Eq.(3.7) we easily find the following expansion, 10 10 Here a diagram of order n has n + 1 external legs.
and then, we obtain by reusing the calculations of integrals in Eq.(3.18) A finite series of order n for S ν is defined as Unfortunately we find, however, that these finite series do not work as approximations even at ν = 1 as shown in Fig.4 and it is inevitable to use some technique for obtaining good approximations.

Padé approximations and Large ν behaviors 4.1 The bag model for large ν
The result of the vortex size R BPS in Eq.(2.32) implies that the total magnetic flux of a vortex is proportional to an area occupied by the flux for ν > 0, where m 2 /2 = e 2 v 2 /2 is the maximum of the magnetic field allowed by the BPS equations Eq.(2.6) for ν > 0. This fact evokes the liquid droplet model of nuclear structure, and gives an intuitive explanation in our axially symmetric case for the Bradlow bound [20], which means just that the area πR 2 BPS must be less than the total area if we considered a closed two-dimensional base space.
In a paper [23], the size R BPS was obtained by a physically intuitive way using the bag model proposed in [21,22] for the large winding number ν. In the bag model, a vortex configuration consists of an inside Coulomb phase, the outside vacuum in the Higgs phase, and a thin domain-wall at r = R interpolating their phases. In the Coulomb phase, the magnetic field takes a non-vanishing constant determined by the total magnetic flux in Eq.(2.3) with ν = k, and vanishes in the vacuum. By omitting a thickness of the domain-wall, profiles of the Higgs field and the magnetic fields are approximated by for r < R 0 for r > R , (4.2) of which the total energy is calculated as This energy is minimized just at R 2 = 4ν/e 2 v 2 = R 2 BPS . Actually, we numerically observe profiles of the magnetic field for large ν in Fig.7. A profile of the domain-wall is almost invariant with various values of ν. For large ν, therefore, a contribution to the total energy T bag form the domain-wall can be negligible.
Since a vortex configuration for large ν drastically changes around the domainwall at r ≈ R ≫ 1/m, we expect that the approximation for r ≪ R in in Eq.(2.43) is applicable for and similarly the asymptotic behavior in Eq.(2.53) is applicable for r = R + ǫ (4.5) Inserting R = R BPS = 2 √ ν/m, these estimations give large-ν behaviors of C ν and D ν as We also estimate S ν as Not that the term proportional to √ ν comes from contribution of surface of the vortex and the coefficient β must be positive due to Eq.(2.62). The above estimations for large ν will become important clues to modify the approximations using the ν-expansion.

(Global) Padé approximations
Let us assume that we know only a finite series of order n, as a part of a certain infinite series F (ν) and it behaves as almost an alternating series like F (ν) = |f 0 | − |f 1 |ν + |f 2 |ν 2 − ..., and it seems to have a small radius of convergence ν ≈ ν c . To get a good approximation for ν > ν c with such a series, it is powerful to use the Padé approximation which replace the series by some rational functions, with n = m + l, where a Padé approximant of F (ν) is given by where coefficients of these rational functions are determined so that they satisfies for k = 0, 1, · · · , m + l. Here the two sets {a n } and {b n } are determined uniquely from the finite set {f 0 , f 2 , · · · , f n+m }.
There is arbitrariness in a choice of (m, l) for the order n. The approximant P (m,l) [F (ν)] behaves for large ν as Note that if we fix p = m − l to remove that arbitrariness, then n is restricted so that n − p = 2l. In the case of p = 1 for example, we arrange the Padé approximants for all of the order n as The series expansion for S ν seems to be almost alternative series, and according to configurations for the finite series E n [S ν ] shown in the left panel of Fig.4 we guess that the radius of convergence is around |ν| ≈ 0.5 which implies, for instance, that the function S ν has a singularity at ν ≈ −0.5. The Padé approximation can avoid such a singularity and enlarge the radius of convergence. Let us take the following rational functions P n [S ν ] with respect to ν, as Padé approximants of the order n for S ν , (4,2) [S ν ] = ν 2 + 1.34536ν 3 + 0.0556454ν 4 1 + 2.90796ν + 1.86162ν 2 , P 6 [S ν ] = P (4,3) [S ν ] = ν 2 + 1.94979ν 3 + 0.69144ν 4 1 + 3.5124ν + 3.44191ν 2 + 0.814411ν 3 . (4.14) Here we have fixed arbitrariness on choice of the Padé approximants P (m,n) [S ν ] so that all coefficients of the above are positive. 11 As a result poles and zeros of these functions turn out to sit only on the negative real axis of ν as shown in Fig.8 and the rational functions P n [S ν ] have poles around ν ≈ −0.5 in common. Actually these functions give good approximations in a wider range of ν as shown in Fig.9. Note that these rational functions behave as   Fig.8 can be regarded as disturbances for large-ν behaviors. Let us consider the large-ν behavior more seriously. The large-ν behavior in Eq.(4.7) does not always mean that the function S ν has a branch cut. For an example, a function √ ν tanh( √ ν) has no branch cut anywhere although it behaves √ ν for large ν ∈ R >0 . Here, we just assume existence of a branch cut. For instance, a function has a branch point at ν = −1/2 and desirable behaviors as and consequently it works as a quite good approximation for the full range of ν as shown in Fig.4. The Padé approximation taking account of informations for large ν is called the global Padé approximation [24]. Note that an expansion of the following quantity is also alternative series due to the singularity, +18.0239ν 4 − 37.857ν 5 + 79.5748ν 6 + O(ν 7 ), (4.18) Let us apply the Padé approximation to the above series or its squared quantity. According to Eq.(4.7), the above quantity behaves as O(ν −1 ) for large ν and this property fixes the arbitrariness of Padé approximants completely. Addition to P 1 [S ν ] in the above, then, we obtain the following functions as the global Padé approximants of S ν , 1 + 0.697034ν 1 + 4.69703ν + 6.53772ν 2 + 2.31356ν 3 , P 5 [S ν ] = ν − ν 1 + 1.11774ν + 0.064997ν 2 1 + 3.11774ν + 2.17527ν 2 + 0.0904502ν 3 , P 6 [S ν ] = ν − ν 4 1 + 1.81492ν + 0.525555ν 2 1 + 5.81492ν + 11.5348ν 2 + 8.60739ν 3 + 1.63522ν 4 , (4.19) which behave for large ν as At this stage we do not know whether β n converges to a true value of β. As we see in Fig.10, the global Padé approximation works well and P 6 [S ν ] has a quite small errors less than 10 −3 in the full range of ν. Even for small ν, the global

P n [D ν ]
The ν-expansion of D ν given in Eq.(3.18) also seems to be almost an alternating series and have a finite radius of convergence as shown in Fig.2. Hence let us 12 We wish, although, to modify a slow convergence of the large-ν behavior if possible. Note that a natural and probable expansion of Sν around the infinity ν = ∞ is which have a pole ν ≈ −1 in common as seen in Fig.8. As shown in Fig let us apply the Padé approximation not to D ν it self, but to exp(2nD ν /ν) with n = 1, 2, 4, then we obtain Hence, as shown in Fig.13 and we obtained gives which reproduces the numerical result presented by de Vega and Schaposnik. This value with the similar accuracy was also obtained analytically in Ref. [25].@ Figure 14: Errors of D 1 in the left panel, and D 5 in the right panel.

P n [C ν ]
The ν-expansion of C ν in Eq.(3.13) gives a quite good approximation for C ν , at least for ν = 1 and we do not know whether the radius of convergence is larger than one or not. In this stage, therefore, it is not useful to apply the (ordinary) Padé approximation to E n [C ν ]. Once we take account of the large-ν behavior of C ν given in Eq.(4.6), however, we notice that there exists a singularity at the infinity and we have to remove this at the first stage. Let us consider the following function which has an infinite number of zeros on a negative real axis of ν and regular everywhere except for an essential singularity at the infinity. The nearest next zero to the origin is ν = −π 2 /4 ≈ −2.47. It is, therefore, natural to assume that a quantity F ν ≡ (C ν / C ν ) 4 has an infinite number of poles (and zeros) on the negative real axis of ν. Actually we find that an expansion of F ν gives an almost alternative series as, According to Eq.(4.6), F ν must behave for large ν as which means that we removed the singularity at the infinity in success. Next, let us apply the Padé approximation to the series in Eq.(4.29) or its squared quantity, satisfying the property Eq.(4.30). We obtain, 13 In Fig.15, we observe that these functions give nice approximants in the full range of ν and modify E n [C ν ] as shown in Fig.6. Resultantly, even for C 1 , we succeeded in reproducing the numerical result C 1 = 1.7079 given by de Vega and Schaposnik as (4.35) 13 There exists still arbitrariness on a choice of a functionCν . We can choose, for example,

Summary and Discussion
We considered the small winding-number expansion (the ν-expansion) of the solution of the Taubes equation by extending the winding number, which is a topological quantum number, to be a real number larger than −1. We confirmed that the ν-expansion is useful to give good approximations of axially-symmetric vortex solutions in most of the range allowed for the winding number. Finally we found that for the scalar charge C 1 the best approximate value in terms of the ν-expansion with the help of the Padé approximation is P 6 [C 1 ] = 1.7078629 . . ., which coincides with a value N sht [C 1 ] = 1.707864175, obtained numerically by the shooting method. We judged that the result given by de Vega and Schaposnik is correct, and Tong's conjecture giving C 1 = 8 1 4 ≈ 1.68 from superstring theory perspective is incorrect as a vortex solution in the Abelian-Higgs model. Their numerical similarity might suggest a certain universality.
The Abelian-Higgs model of critical coupling is just the simplest toy model to test and establish usefulness of the ν-expansion. The idea of the ν-expansion is rather simple and more straightforward than the strategy taken by de Vega & Schaposnik. As for BPS states of vortices in further complicated systems like non-Abelian gauge-Higgs models or of separated (parallel) multi-vortices, therefore, it is expected that the ν-expansion can be straightforwardly applied to their analytical approximations. Since it is difficult to apply the shooting method to such complicated systems, we guess that the role of the ν-expansion will become more important there. The ν-expansion is also expected to be powerful to analyze dependence on dimensionless parameters of solutions, like dependence on the number N and a ratio of two gauge couplings of U (N ) = [U (1) × SU (N )]/Z N for an U (N ) vortex.
We expect that the ν-expansion can be applied to systems of non-critical coupling, although it might not be a straightforward extension. Our final goal is to establish a systematic tool to study the dynamics of vortices quantitatively without taking the critical coupling limit. Since in the ν-expansion vortices are treated as singular particles (strings) in a three(four)-dimensional spacetime, it will become possible to treat vortices of arbitrary shapes and discuss their dynamics analytically and quantitatively if we can consider such an extended ν-expansion.

A.2 Sequence of sets of upper and lower bounds
Here let us modify the inequality Eq.(2.37) for ν > 0.
to obtain a stronger set of upper and lower bounds of them. By integrating Taubes equation and P = rψ ′ , we find relations between P and ψ using integrals as, with Y = (r/R in ) 2ν and setting m = 1, · · · ≥ f M n−1 ≥ f M n > f m n ≥ f m n−1 ≥ · · · ≥ f m 0 = 0, 0 = g M 0 ≥ · · · ≥ g M n−1 ≥ g M n > g m n ≥ g m n−1 ≥ · · · ≥ g m 0 = −2ν. and therefore we find the followings are required , (A. 15) and that is, R in must satisfy otherwise a function ψ can not satisfy the set of inequalities I 0 and thus blows up at large r. With R in satisfying the above set of inequalities, the next set of inequalities I 2 can be consistently obtained as for r > 2 √ ν , In principle, you can calculate I 3 , I 4 , . . . , sequentially as you like.

B Some Integrals
Since the modified Bessel function of the second kind is a two dimensional Green's function, we can find the following relations By using the integral formulas