A Theory of Giant Vortices

We elaborate a theory of giant vortices [1] based on an asymptotic expansion in inverse powers of their winding number $n$. The theory is applied to the analysis of vortex solutions in the abelian Higgs (Ginzburg-Landau) model. Specific properties of the giant vortices for charged and neutral scalar fields as well as different integrable limits of the scalar self-coupling are discussed. Asymptotic results and the finite-$n$ corrections to the vortex solutions are derived in analytic form and the convergence region of the expansion is determined.


Introduction
Vortices, string-like topologically nontrivial solutions of field equations, naturally appear in the theory of superconductivity [2], superfluidity [3], and QCD confinement [4]. They play a crucial role in many physical concepts from cosmic strings [5] and vacuum selection [6] to mirror symmetry and dualities of supersymmetric models [7]. Giant vortices carrying large topological charge are of particular interest and are observed experimentally in a variety of quantum condensed matter systems [8][9][10]. Corresponding winding numbers range from n = 4 in mesoscopic superconductors [10] through n = 60 in Bose-Einstein condensate of cold atoms [9] and up to n = 365 in superfluid 4 He [8]. From a theoretical perspective it is quite appealing to identify characteristic features and universal properties of vortices in the limit of large n. A solution to this problem is too subtle for a straightforward numerical analysis and requires some form of analytic approach.
Though the system of vortex equations is relatively simple, its exact analytic solution is not available even for critical coupling when hidden supersymmetry reduces the order of the equations [11] and even for the lowest winding number n = 1 [12] in contrast to the apparently more complex case of magnetic monopoles [13]. In general, finding analytic solutions of higher topological charge is a very challenging problem and only a few such solutions are known in gauge models (see e.g. [14,15]). By contrast the structure of vortices vastly simplifies in the limit of large winding number n → ∞ [16,17]. In a recent letter [1] we have introduced a framework which enables a systematic expansion in inverse powers of n to obtain the asymptotic form of the axially symmetric giant vortex solution. In this framework a topological quantum number n is associated with a ratio of dynamical scales and a systematic expansion in inverse powers of n is then derived in the spirit of effective field theory. Below we present a detailed account of the method and its application to the Abrikosov-Nielsen-Olesen (charged field) and Ginzburg-Pitaevskii (neutral field) vortices.

Abrikosov-Nielsen-Olesen vortices
We consider the standard Lagrangian for the abelian Higgs (Ginzburg-Landau) model of a scalar field with abelian charge e, quartic self-coupling λ, and vacuum expectation value η in two dimensions where D µ = ∂ µ + ieA µ . Vortices are topologically non-trivial solutions of the corresponding Euclidean equations of motion. We study the axially symmetric solutions of winding number n, which in polar coordinates can be written as follows φ(r, θ) = f (r)e inθ , A θ = −na(r)/e, A r = 0. For a given winding number the solution carries n quanta of magnetic flux Φ = − F 12 d 2 r = 2πn. When the winding number grows, the characteristic size of the vortex has to grow as well to accommodate the increasing magnetic flux. Assuming a roughly uniform average distribution of the flux inside the vortex for λ, η = O(1) we can estimate its radius to be of order √ n/e. At the same time a characteristic distance of the nonlinear interaction is 1/e. Thus for large n we get a scale hierarchy and the expansion in the corresponding scale ratio is a standard tool of the effective field theory approach. Since we deal with the spatially extended classical solutions it is more convenient to perform this expansion in coordinate space at the level of the equations of motion. The analysis becomes particulary simple for critical coupling λ = e 2 , which we discuss first.

Critical vortices
In this case the vortex dynamics is governed by the first-order Bogomolny equations [11] (D 1 + iD 2 ) φ = 0 , − F 12 + e |φ| 2 − η 2 = 0 , (2.2) and the vortex energy (string tension) is proportional to the topological charge T = − Ld 2 r = 2πnη 2 . It is convenient to introduce the rescaled dimensionless quantities eηr → r, f /η → f , λ/e 2 → λ, so that in the new variables e = η = 1 and critical coupling corresponds to λ = 1. Then the Bogomolny equations in terms of the functions a(r) and f (r) take the following form with the boundary conditions f (0) = a(0) = 0 and f (∞) = a(∞) = 1. For large n the field dynamics is essentially different in three regions: the core, the boundary layer, and the tail of the vortex. Below we discuss the specifics of the dynamics and its description in each region.
The vortex core. For small r the solution of the field equations gives f (r) ∝ r n . This function is exponentially suppressed at large n for all r smaller than a critical value, which can be associated with the core boundary. For such r the contribution of f (r) can be neglected in the equation for a(r) and we get a(r) ∼ r 2 /r 2 n with r n = √ 2n, which in turn can be used in the equation for f (r). Thus in the core the dynamics is described by linearized equations in the background field Their solutions read where the integration constant F in the first line is determined by matching conditions explained below. For r n − r = O(1) we have n(1 − a(r))/r = O(1) and the equation for f (r) becomes independent of n. Hence the approximation Eq. (2.4) is not applicable anymore, the nonlinear effects become crucial, and we enter the boundary layer. Note that the magnetic flux and energy density for Eq. (2.5) are approximately 1 and η 2 , respectively, so that the core accommodates essentially all the vortex flux and energy and we can identify r n with the vortex radius. The boundary layer. In this region the field dynamics is ultimately nonlinear. However, it crucially simplifies for large n. To see this we introduce a new radial coordinate x = r−r n so that in the boundary layer x = O(1) and the expansion in x/r n converts into an expansion in 1/ √ n. In the leading order in x Eq. (2.3) reduces to a system of n-independent field equations with constant coefficients where w(x) = ln f (r n + x), γ(x) = n (a(r n + x) − 1) /r n , and prime stands for a derivative in x. The system can be resolved for w, which results in a second-order equation This equation has a first integral I = w 2 − e 2w + 2w with I = −1 corresponding to the boundary condition w(∞) = 0. Thus Eq. (2.6) can be solved in quadratures with the result which determines a unique asymptotic solution in the boundary layer. It has a Taylor expansion w(x) = ∞ m=0 w m x m where w 1 = (e 2w 0 − 2w 0 − 1) 1/2 and the higher order coefficients can be obtained recursively The asymptotic behavior of the function at x → ∞ reads where the constant is computed in Appendix A. By matching Eq. (2.11) to the core solution, Eq. (2.5), in the region 1 r n − r r n where both approximations are valid we get F = 1/ √ e. The vortex tail. For (r −r n )/r n = O(1) the boundary layer approximation breaks down and the coordinate dependence of the field equation coefficients should be restored. However, the deviation of the fields from the vacuum configuration is now exponentially small so the field equations linearize. The solution of the linearized theory is well known and reads where K m (z) is the mth modified Bessel function. It describes the field of a point-like source of scalar charge ν and magnetic dipole moment µ with ν = µ for critical coupling. Eqs. (2.8) and (2.13) should coincide in the second matching region 1 r − r n r n , which yields Next-to-leading order solution. For a finite winding number the above asymptotic formulae have a finite accuracy which may deteriorate as n decreases. To get control over the accuracy and convergence of the large-n expansion we compute the O(1/ √ n) correction to the asymptotic result. Let us consider first the boundary layer. Writing the corrections to the asymptotic solutions w and γ as δw/ √ 2n and δγ/ √ 2n, respectively, and expanding Eq. (2.3) through O(x/r n ) we get the following field equations i.e. δγ is completely determined by the solution for δw. The latter obeys the boundary conditions δw(x) ∼ x 3 /6 + O(1) at x → −∞, δw(∞) = 0, and can be found in a straightforward though not completely obvious way. By equating the variation of the first integral I to zero and changing the variable from dx to dw = w dx we get the following first-order equation for the homogeneous part of the solution δw h One solution of this equation is given by the translational zero mode w (x) and we take w (x) x −∞ dy/w 2 (y) as the second solution so that the pair has unit Wronskian. This results in the following solution of Eq. (2.15) which satisfies the boundary condition at where x 0 is an integration constant to be determined. A variation of x 0 changes Eq. (2.17) by a term proportional to w (x) and therefore does not affect its behavior at x → ∞. For x 0 → ∞ the two integrals in Eq. (2.17) can be combined into The last integral up to a term proportional to w (x) is equal to Thus the next-to-leading contribution can be written as follows where the integration constant is adjusted to satisfy the boundary condition at x → −∞ (see Appendix A). Eq. (2.20) determines the next-to-leading approximation in the boundary layer. In the core and tail regions the neglected nonlinear terms are exponentially suppressed and the corrections to the asymptotic result enter only through the matching to the boundary layer solution. To perform the matching we need the asymptotic behavior of Eq.  (9) ln (10) ln (11) ln (12) ln (13) ln (14) ln (15) ln(|ν|) n leading next-to-leading where and The first line of Eq. (2.22) in the matching region 0 r n −r r n determines the correction to the core solution integration constant Analogously the matching of the second line of Eq. (2.22) to the tail solution in the region 1 r −r n r n determines the correction to the scalar charge ν, which we write as δν/ √ n. Expanding Eq. (2.13) at large r and keeping the subleading term we get which completes the next-to-leading approximation of the critical vortex solution.
We can now scrutinise the accuracy and convergence of the large-n expansion by comparing the leading and next-to-leading approximations to numerical solutions of the exact field equations for different winding numbers. The results of the analysis for the critical vortices are presented in Figs. 1-4. In Fig. 1 the leading and next-to-leading boundary layer solutions for the scalar field f (r) are plotted against the exact (numerical) solutions with n = 1, 4, 10. Similar plots for the gauge field a(r) are given in Fig. 2. In Figs. 3 and 4 the exact numerical values of f (r n ) and ν, the natural characteristics of the vortex solution, are plotted against the leading and next-to-leading results for n = 1, . . . , 100. The expansion reveals an impressive convergence and the next-to-leading approximation works reasonably well even for n = 1.

Noncritical vortices
For noncritical scalar self-coupling λ = 1 the order of the field equations cannot be reduced and they read (2.28) Nevertheless, the general structure of the solution is quite similar to the critical case. Inside the core the contribution of the scalar potential to Eq. (2.28) is suppressed by r 2 /n 2 . Hence the core dynamics is not sensitive to λ and the core solution is given by Eq. (2.5) up to the value of the integration constants which do depend on λ through the matching to the non-linear boundary layer solution. In particular the vortex size r n is determined by the region where the two terms in the square brackets of Eq. (2.28) become comparable and the core approximation breaks down, which gives the leading order result r n = √ 2n/λ 1/4 . Note that the approximately constant energy density in the core is now λη 2 so that the total vortex energy in the large-n limit is T = 2π √ λnη 2 . In the tail solution, Eq. (2.13), the argument of K 0 gets an additional factor of √ λ to account for the variation of the scalar field mass, while the scalar charge and the magnetic dipole moment are not equal anymore and have different leading behavior at n → ∞ |ν| ∼ e 2 More accurately these parameters as well as the normalization of the scalar field in the core solution are determined by matching to the boundary layer solution. In the boundary layer by expanding r = r n + x in x/r n we get a system of n-independent equations with constant coefficients Here the first integral in the brackets is the magnetic flux 2πn while the second integral is saturated in the boundary layer with an exponential accuracy. Thus in the latter we can approximate the volume element as follows The correction to the critical vortex energy then reads where the first term originates from the expansion of the asymptotic result 2π √ λnη 2 in δλ and the second term represents the correction to the boundary layer energy. Thus the Figure 6. Same as Fig. 5 but for the normalized gauge field γ(x)/ √ λ.
energy of a near-critical giant vortex to the first order in δλ can be written as follows Large scalar self-coupling. In the limit λ → ∞, after a field rescaling γ λ (x) = γ(x)/ √ λ we may neglect the derivative term in the first line of Eq. (2.30), which becomes algebraic and can be solved for the function f (r). The field equations then become , where x 0 is the remaining integration constant formally determined by the boundary conditions at x → −∞. However, for large scalar self-coupling these boundary conditions are rather peculiar since the matching region shrinks to a point x = 0. Indeed, in the matching region the scalar field is exponentially suppressed f (r) = O(e − √ λx 2 /2 ) and at λ → ∞ it vanishes for all x < 0 and at x = 0 its derivative is infinite. At the same time for negative x the gauge field is given by the core solution and has a continuous first derivative at x = 0. These conditions are satisfied for Thus the boundary layer solution for x ≥ 0 reads while f (r) = 0, a(r) = r 2 /r 2 n inside the core for r ≤ r n + δr n , where δr n = −1. Note that now r = r n + δr n + x, i.e. the core radius which defines the origin of the x coordinate gets a O(λ 1/4 / √ n) correction. This correction makes the winding number of the above gauge field configuration an integer. We can now compute the finite-n correction to the vortex energy. As in the case of near-critical coupling it comes from two sources. Due to the variation of the radius the core energy changes by −1 in the units of 2πλ 3/4 √ 2nη 2 . The contribution of the boundary layer in the same units is 2 . In total this gives Eq. (2.37) is consistent with the known negative surface tension on the phase interface in extreme type-II superconductors [18]. The boundary layer solution, Eq. (2.36), is subject to enhanced O(λ 1/4 / √ n) corrections which can be obtained in analytic form by the method of Sec. 2.1. It is convenient to write the corrections to f (r n + δr n + x) and γ λ (x) as λ 1/4 √ 2n δf (x) and λ 1/4 √ 2n δγ λ (y), where y = x − x 0 . Then for the gauge field we get (see Appendix B) Inside the boundary layer the scalar self-coupling can be neglected in the field equations. However, the boundary conditions for the corresponding solution do depend on λ. As before matching to the core requires γ(x) ∼ √ λx at x → −∞ while matching to the half-kink solution Eq. (2.40) gives f (r n + x) ∼ λ/2x at 1 x x λ . This dependence can be eliminated by a rescaling (2.41) In the new variables the boundary layer equations read with the boundary conditions γ λ (z) ∼ z, f λ (z) = O(e −z 2 /2 ) at z → −∞ and f λ (z) ∼ z, γ λ (z) = O(e −z 2 /2 ) at z → ∞. The above system is invariant with respect to the discrete field transformations f λ → −f λ , γ λ → −γ λ , and f λ ↔ γ λ , as well as the coordinate reflection z → −z. Taking into account the boundary conditions we find that the gauge and the scalar fields have "mirror" boundary layer solutions Despite the high symmery and the existence of a first integral γ λ 2 +f λ 2 −γ 2 λ f 2 λ = 1 Eq. (2.42) cannot be integrated in quadratures, at least not in a straightforward way. However, we can conclude that inside the boundary layer of width 1/λ 1/4 the gauge and scalar fields scale as λ 1/4 , and hence the boundary layer contributes a plain O(1/ √ n) correction to the vortex energy. Thus the total correction to the energy is dominated by the scalar cloud outside the boundary layer, Eq. (2.40), which gives Eq. (2.44) agrees with a classical result for the positive surface tension in extreme type-I superconductors [19]. The scalar cloud solution itself gets the enhanced finite-n correction. It can be written as δf (y)/( √ 2nλ 1/4 ), where y = x/x λ and δf (y) = 9 + 12y − 8e −2y − e −4y 6 √ 2(2 + e 2y + e −2y ) . (2.45) The details of the derivation of this result are presented in Appendix B. We should emphasize that the accuracy and convergence of the large-n expansion for noncritical vortices crucially depend on the scale of λ, cf. Eqs. (2.37, 2.44). For large scalar self-coupling the expansion clearly breaks down at n = O( √ λ) where the ratio of scales characterizing the boundary layer and the core dynamics δr n /r n is not small anymore. At this value the dependence on the winding number qualitatively changes. Indeed, in the limit λ → ∞ with fixed n we should recover the classical result [2,19] where the scalar field core radius vanishes as n/ √ λ, a characteristic size of the gauge field core is O(1), and the energy is quadratic in the magnetic flux and logarithmic in the scalar coupling constant For small λ the convergence problem appears at n = O(1/ √ λ). For such n the expansion in the ratio x λ /r n of the scales characterizing the dynamics of the scalar cloud and the vortex core breaks down, and the structure of the solution is no longer described by the asymptotic formulae. In particular, in the limit λ → 0 at fixed n the vortex energy is dominated by the scalar cloud contribution, which can depend on the scalar self-coupling and the winding number at most logarithmically.

Ginzburg-Pitaevskii vortices
The Ginzburg-Pitaevskii vortex equation [3] for a neutral complex scalar field φ which describes a Bose-Einstein condensate of weakly interacting Bose gas can be obtained from the analysis of the previous section by setting the gauge field to zero and rescaling the scalar self-coupling to λ = 1. The corresponding radial equation reads The structure of its solution in the limit of large n is quite distinct from the Abrikosov-Nielsen-Olesen case discussed above. 1 Now the field dynamics is essentially different in the vortex core r r n , its tail r r n , and in the matching region around r = r n , where r n = n gives the leading approximation for the size of the region with unbroken symmetry phase. Let us consider the field dynamics and the vortex solution in each region. The vortex core. For r r n the solution f (r) ∝ r n is exponentially suppressed with n so that the nonlinear cubic term can be neglected. As r approaches r n the O(1) linear term in the equation becomes relevant while the exponentially suppressed nonlinear term can still be omitted and we get a linear differential equation with the regular solution proportional to the nth Bessel function J n (r). Since J where Ai(z) is the Airy function, x = r − r n , and z = (2/n) 1/3 x. The Airy function is exponentially suppressed for large positive values of the argument, hence the core approximation is valid for r n − r ∼ > n 1/3 . The vortex tail. In the opposite limit r r n the derivative term in the equation is suppressed and can be treated as a perturbation. As a result the field equation becomes algebraic. For example, in the leading order we get with the solution The higher order terms can easily be obtained by iterations to any desired order in 1/r. The first three terms read 1 After a rescaling of the radial coordinate r → r/ √ λ it describes the scalar field core of the Abrikosov-Nielsen-Olesen vortex for λ n 2 1.
The expansion Eq. (3.6) formally breaks down for x = O(n 1/3 ). Indeed, let us consider the behavior of the above series at x = κ 1−α n α /2 for some power α. Here the factor κ 1−α /2 with a constant κ is introduced for convenience. Then the series Eq.(3.6) can be written as an expansion where ρ = κ 3α−3 /n 3α−1 . Each f m can in turn be expanded in τ = (κ/n) 1−α . For the first three terms we get . . , For α ≤ 1/3 the expansion parameter ρ is not suppressed at large n and the series Eq. (3.6) is not defined (though for the boundary case α = 1/3 we have ρ = 1/κ 2 and the series may converge for sufficiently large κ). At the same time the expansion parameter τ is not suppressed for α ≥ 1. Thus the allowed interval for the expansion in negative powers of n is 1/3 < α < 1. For α = 1/2 both expansions in ρ and τ convert into a series in 1/ √ n. Note that by choosing a smaller or larger value of α we boost the convergence of the expansion in τ or ρ, respectively, and for a given α the parameter κ can be tuned to further improve the convergence of the series. The matching region. The problem now is to derive a large-n expansion in the matching region −n 1/3 ∼ < x ∼ < n α . First we note that in this interval the solution is suppressed at least as 1/n , which means that they can be polynomially matched over the corresponding interval. A linear matching of the leading core and tail solutions gives a simple analytic formula [1] f (r) =        C J n (r) , r ≤ r n , (δ/2n) 1/2 (1 + (r − r n )/δ) , r n < r < r n + δ , 1 − n 2 /r 2 1/2 , r n + δ ≤ r , (3.9) where C = 3 1/4 π 1/2 2 1/2 , δ/n 1/3 = 2 2/3 π 3 5/6 Γ 2 (2/3) , and Γ(z) is the Euler gamma-function. In the limit n → ∞ the function defined by Eq. (3.9) becomes continuous and smooth. As has been discussed above, the higher orders of the perturbative expansion of the core solution in the nonlinear term at x = O(1) as well as the perturbative expansion of the tail solution in the derivative term at x = O(n 1/3 ) are not suppressed by extra powers of n (though both expansion may de facto converge). Thus Eq. (3.9) gives the asymptotic solution up to O(1/n 1/3 ) corrections. It can be used to obtain the analytic result for the vortex energy at n → ∞. The total energy consists of the logarithmically divergent rotational part T rot and the nonrotational finite part T [3]. The rotational contribution determined by the angular The relevant solution must exponentially decay at z → ∞ and match the expansion of the next-to-leading tail solution f (r n + x) at x = √ κn/2 through O(1/n). It can only be found numerically. Note that since the expansion parameter x/n varies from 1/n to 1/ √ n through the matching region, the convergence of the series there is not uniform. In general the expansion of the tail solution in ρ matches the asymptotic expansion of χ 0 (z) at large negative z while the series in τ for each f m is connected with the corrections to the equation for χ(z).
The results of numerical analysis for the neutral scalar field vortices with n = 4, 10, 16 are presented in Fig. 7. For the leading order approximation the linear matching is used and in the next-to-leading order the numerical solution of Eq. (3.14) is matched to the first two terms of Eq. (3.6) at r n + n/2, which corresponds to κ = 2.

Summary
We have elaborated a method of expansion in inverse powers of a topological quantum number. The method is quite general and can be applied to the study of topological solitons in a theory where the corresponding quantum number is associated with a ratio of dynamical scales. When applied to axially symmetric vortices with large winding number n the expansion is in negative noninteger powers of n and in general is not uniform. In the large-n limit the complex nonlinear vortex dynamics unravels. In particular, for the Abrikosov-Nielsen-Olesen vortices at critical coupling and in the limit of large or small scalar self-coupling the field equations become integrable. At the same time for the Ginzburg-Pitaevskii vortices in a weakly interacting Bose-Einstein condensate the field equation reduces to an algebraic one. The method yields simple asymptotic formulae for the shape and parameters capturing the main features of the giant vortices. The accuracy of the asymptotic result can be systematically improved by calculating the finite-n corrections. For near-critical vortices the approximation works remarkably well all the way down to very low n after including the leading O(1/ √ n) corrections. In the case of extremely small or large scalar self-coupling λ the large-n expansion converges and gives a reliable approximation for all winding numbers satisfying the condition n √ λ, 1/ √ λ.

Acknowledgement
A.P. is grateful to Valery Rubakov for a discussion of noncritical vortices. The work of A.P. was supported in part by NSERC and the Perimeter Institute for Theoretical Physics. The work of Q.W. was supported through the NSERC USRA program.

A Calculation of the critical vortex parameters
Calculation of w ∞ . The integral in Eq. (2.8) is logarithmically divergent at w(x) → 0 and can be decomposed as follows which defines an inhomogeneous differential equation on the function δγ λ (y), where y = x − x 0 . The boundary conditions for this equation are δγ λ (∞) = 0 and δγ λ (−x 0 ) = 1/2, where the latter is determined by matching to the core solution for the gauge field at r n +δr n . Note that δf and δγ λ are not continuous at this point due to the singular character of the matching region in the limit λ → ∞ discussed previously. Due to the existence of the first integral for Eq. (2.35) the homogeneous part of the solution δγ h λ satisfies the first-order equation One of its solutions is given by the translational zero mode γ λ (y) = 2 sech( √ 2y) tanh( √ 2y). The second solution reads