# Vortex Structures inside Spherical Mesoscopic Superconductor Plus Magnetic Dipole

Authors: A. Ludu
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

We investigate the existence of multivortex states in a superconducting mesoscopic sphere with a magnetic dipole placed at the center. We obtain analytic solutions for the order parameterΨ(r→) inside the sphere through the linearized Ginzburg-Landau (GL) model, coupled with mixed boundary conditions, and under regularity conditions and decoupling coordinates approximation. The solutions of the linear GL equation are obtained in terms of Heun double confluent functions, in dipole coordinates symmetry. The analyticity of the solutions and the associated eigenproblem are discussed thoroughly. We minimize the free energy for the fully nonlinear GL system by using linear combinations of linear analytic solutions, and we provide the conditions of occurring multivortex states. The results are not restricted to the particular spherical geometry, since the present formalism can be extended for large samples, up to infinite superconducting space plus magnetic dipole.


---

## Body

## 1. Introduction

The rapidly growing field of quantum computation requires nanoscale miniaturization of electronic circuits, way beyond the silicon era type of devices [1]. Therefore, the mesoscopic superconductors, having the size comparable to the coherence length or the magnetic penetration depth [2], are the prime candidates for construction of nanodevices among all other superconducting systems. Mesoscopic physics revealed a number of open fundamental problems like quantum confinement, quantum vortices and loops, spintronics, quantum dots, etc. [3] whose solutions can bring significant knowledge in fields like nanotechnology, synthesis of new materials, novel sensors, modern lithography, or molecular biology [4–6].The most important feature of a mesoscopic superconductor is that its shape and size have an important effect on the interplay of the magnetic field and superconducting condensate. The properties of mesoscopic superconductors are very different compared to those of bulk superconductors. While in bulk superconductors penetrating vortices form a lattice due to the vortex-vortex repulsion, in mesoscopic superconductors there is a competition between the vortex lattice and the boundary which tries to impose its geometry on the vortex lattice. It is observed experimentally that flux quantum configurations have the same symmetry as the symmetry of the shape of the sample in a homogeneous magnetic field [7]. Such systems are studied via the Ginzburg-Landau (GL) model. The GL equations arise from the Euler-Lagrange equations for the free Gibbs energy for a mesoscopic superconductor sample in magnetic field. These equations must be solved under specific boundary conditions: the normal component of the superconducting current is equal to zero [8]. Near and below the transition temperature, theoretical calculations show that anti-vortex and giant-vortex can appear in order to maintain the symmetric vortex configuration.The response of mesoscopic superconductor samples of different shapes (thin discs, spheres, cones, and rings) to an external magnetic field, as well as the effect of the geometry, has been theoretically [9–18] and experimentally [7] studied. In all these cases a constant external magnetic field is applied along the the revolution axis. The small volume to surface ratio of these mesoscopic structures brings new features not found in the bulk: there are two kinds of superconducting states, depending on the strength of the magnetic field, the sample geometry (external surface), and its size: giant and multivortex states. The giant vortex state has cylindrical symmetry and is the only kind stable in small size superconductors due to the confinement effect [13, 14]. If the size of the sample increases, such giant vortex states can break up into multivortices through saddle-point transitions [17, 18]. For three-dimensional objects (sphere or cone) the vortex lines need to intersect the surface perpendicularly in order to cancel the outward supercurrent component [9–12]. Consequently, the shape of the lines is strongly affected by the sample surface. For example, in the case of a mesoscopic sphere placed in uniform external magnetic field, the vortex lines are curved inside, packing denser in the equatorial plane, and spreading out towards the poles [10, 11].In this paper we consider a different situation where the magnetic field is not anymore externally generated but is generated from inside the sample,e.g., an infinitesimal magnetic dipole placed at the center of a mesoscopic sphere. Such a configuration can generate confined vortex loops. The topological transition between open and closed vortex loops is controlled by the geometry,i.e., R, and the central magnetic dipole strength.The goal of this paper is to demonstrate the occurrence of multivortex structures in the superconducting sphere plus magnetic dipole configuration, especially below and around the transition temperature. Our approach is based on the GL model of free energy for a finite volumeV. Outside of this volume the Cooper pair density Ψ(r→) describing the superconducting phase, called the order parameter, is zero [5, 6, 9–18]. The free energy functional is given in the GL model by [10, 11](1)F=1V∫VP→Ψ22m+αΨ2+β2Ψ4+h→28πdv,where m is the Cooper pairs mass, P→ is the quantum momentum operator in the presence of magnetic field, and h→ is the intensity of the magnetic field. The temperature dependent coefficient function α(T)<0 and the nonlinear term coupling constant β>0 are typical Landau second-order phase transition parameters [9–12]. For mesoscopic samples, one can neglect the term responsible of the expulsion of magnetic flux from the superconductor, that is last term in Eq. (1) [9–12].The traditional procedure for findingΨ by minimizing the free energy functional Eq. (1) for constant volume V consists in expanding the order parameter Ψ in a basis of eigenfunctions of the corresponding linearized GL problem [10, 11, 16–18] and numerically evaluating the expansion coefficients which minimize the full nonlinear functional Eq. (1).In this paper we solve the GL linear problem analytically and investigate the properties of the eigenfunctions and spectrum. The contribution of the infinitesimal magnetic dipole will be approached in the dipole system of coordinates. The linearized GL equation factorizes in two ordinary differential equations, for the two orthogonal dipole coordinates, respectively. From the physical point of view such a factorization seems natural because far enough from the sphere surface, the abstract surfaces containing vortex lines follow the magnetic field stream lines, but these are exactly the dipole magnetic field lines. Along these surfaces the order parameter has slow variation or is practically constant. This gives sufficient physical reason for neglecting higher order terms in the dipole variable going along the magnetic field lines. This approximation allows integrating the two ordinary differential equations. The linear solution consists in a product of angular momentum eigenfunctions in the azymuthal coordinate, exponential function for one of the dipole coordinates, and a double confluent Heun function in the third coordinate. The final step is to come back to the fully nonlinear GL problem, and writeΨ as a linear combination of eigenfunctions of the linear GL problem, with arbitrary coefficients, and then minimize the free energy with respect to these coefficients.We dedicate a large part of the present calculations to solve exactly the linearized GL equation and to ensure the completeness and orthogonality of the linear basis because near and below the transition temperature, even the linear GL equation is sufficient to describe multivortex states. The order parameterΨ is very small in this range, and higher order terms of Ψ can be neglected. Nevertheless, at lower temperatures, the vortex configuration does not have to match the symmetry of the system, and higher order terms of GL equation cannot be negligible [19]. It is the contribution of these nonlinear terms which generates the multivortex states at lower temperatures.The paper is organized as follows. In Section2 we formulate the nonlinear and auxiliary linear GL problem and write the partial differential equation associated with the GL problem. In Section 3 we discuss the importance of the infinitesimal central magnetic dipole from a potential aspect, introduce the dipole coordinates, and obtain general form of the dipole equation in azimuthal symmetry plus dipole coordinates. In Section 4 we reduce the general dipole equation to a double confluent Heun equation (DCHE) by the help of a geometric approximation and a decoupling of the dipole orthogonal modes. We obtain analytic solutions for the DCHE as Heun series around the point at infinity, and we present some examples. In Section 5 we show how the dipole equation plus the physical boundary conditions can be brought to a Sturm-Liouville problem, and we solve the associated GL linear eigenproblem and find the eigenvalues spectrum. Examples of spectra for different configurations are presented. In Section 6 we describe how one can use the exact solutions and spectra of the linear GL problem to build approximate solutions for full nonlinear GL problem, by minimizing the free energy. We describe the procedure to identify multivortex states, give the sufficient criteria, and provide an example of equipotential surfaces with vortex structure inside the sphere.


## 2. Ginzburg-Landau Equation for Dipole inside Sphere

In this section we obtain the governing ODEs for the GL problem given in Eq. (1) for a mesoscopic superconducting sphere in a local magnetic field h→=∇→×A→ produced by a infinitesimal dipole placed at its center r→=0. The Euler-Lagrange equations for the free energy functional Eq. (1), written in the absence of the last term in h→2 as explained above, are known as the GL equation [9–12](2)12mP→2Ψ+αTΨ+βΨ2Ψ=0,(3)∇×h→=4πJ→c,(4)J→=q2mΨ∗P→Ψ+P→Ψ∗Ψ,where the electromagnetic momentum P→ is given by(5)P→=-iℏ∇→-qcA→,where A→ is the vector potential, q is the charge of the Cooper pair, and J→ is the current density. The main approach for the nonlinear problem is to minimize the free energy, Eq. (1), with test functions chosen from a basis of eigenfunctions of the associated linear problem. The linearized GL equation (LGL) is obtained by neglecting the last term of Eq. (2)(6)12m-iℏ∇→-qcA→2Ψ=αTΨ.We rewrite Eq. (6) in the London electromagnetic gauge, ∇·A→=0, and we rescale the variables r→→r→-2mα(T)/ℏ,A→→A→q/(c-2mα(T)) and Ψ→Ψ-β/α(T). We follow the GL phenomenological model from [9–18] and use the temperature dependence(7)α=α01-TTc0,ξT=ℏ-2mαT,where Tc0 is the critical temperature in absence of magnetic field and ξ(T) is the coherence length. In the new variables the LGL Eq. (6) can be written in the form of a linear operator equation(8)O^Ψ=△Ψ-2iA→·∇Ψ-A2Ψ+Ψ=0.Since the dipole magnetic field has axial symmetry we can use any axially symmetric orthogonal coordinates in azimuth angle φ in the form r→→(φ,a,b). In this case the φ dependence in Eq. (8) can be separated, and O^ has exp(iLφ) as eigenfunctions labeled by the angular momentum L. This decomposes the space of solutions Ψ(r→) into subspaces parameterized by L(9)L2R3=⨁L∈N⁡VL,VL=ΨΨr→=Φa,beiLφ,with Φ∈L2(R2) square integrable function. The decomposition in Eq. (9) reduces Eq. (6) to a two-dimensional PDE. On each of these L subspaces the operator O^ has a particular form depending on L denoted by O^L, and for each such restricted operator we can attach an eigenvalue problem(10)O^LΦE,κ=EΦE,κ,where κ is a provisional degeneracy label, which later on becomes the radial quantum number. Of course, E and ΦE,λ carry the L dependence.The vector potential associated with the infinitesimal (point-like) magnetic dipole can be expressed in a simple form in spherical coordinates(11)A→=μsinθr2e→φ,where e→φ=(-sinφ,cosφ,0) is the unit vector in the φ direction. The boundary conditions associated with this problem request the order parameter wave function to cancel at the origin and the flow at the mesoscopic volume surface Σ to be zero (the Saint-James-de Gennes conditions) [9–12, 15–18](12)Ψr=0=0,∂Ψ∂nΣ=0,where n→ represents the outward normal to the spherical surface Σ defined by r→=R.We can make coordinate free general statements concerning the spectrum of the eigenfunction problem in Eq. (10). Following the Leinfelder-Simader theorem [20] and taking into account that the integral over R3∖{0} of the forth power of the magnitude of the dipole vector potential is finite, we find out that the operator in Eq. (6) is essentially self-adjoint if applied on analytic functions. Moreover, from the Miller-Simon theorem [21], we know the eigenvalue problem of such type of magnetic multipole potential vector operators is given by a positive, unbounded from above, continuous spectrum. However, given the two-point boundary conditions, Eqs. (12), we expect the reduction of this spectrum to a discrete and bounded spectrum for energy.


## 3. The Generalized Dipole Equation

In this section we present analytic solutions for the system Eqs. (6), (8), and (10) with boundary conditions given by Eq. (12). The presence of the magnetic dipole renders the spherical or the cylindrical coordinates unsuitable to factorize the LGL Eqs. (8) and (10) because of the terms mixing inside the square of the electromagnetic momentum. Attempts to solve similar equations in spherical coordinates in the presence of an electric dipole have been made. In [22], for example, the Schrödinger equation is separable, and the authors find exact wave functions in terms of Bessel and Mathieu functions for m=0 spherical modes in the presence of a dipolar external electric field. A similar generalized Coulomb problem for a class of general Natanzon confluent potentials is exactly solved in [23] by reducing the corresponding system to confluent hypergeometric differential equations. More recently, in [24], the authors succeeded to solve the eigenvalue wave equation for an electron in the field of a molecule with an electric dipole moment by expanding the solutions of a second-order Fuchsian differential equation with regular singularities in cosθ=1,-1,0 in a series of Jacobi polynomials with “dipole polynomial” coefficients.In all these situations separation and integration of the Schrödinger equations are possible because in spherical (cylindrical) coordinates the resulting second-order ordinary linear differential equations are of Fuchsian type, with maximum three regular singularities. Such equations can be mapped into several types of hypergeometric differential equations.However, the presence of the dipole magnetic field makes the situation more complicated because the order of singularity growths above the hypergeometric one, and consequently the dynamics is governed by differential equations of Heun, Hill, or Riemann type [25–27]. A similar situation occurs in the case of a charged particle moving on a sphere under a radial magnetic field and Coulomb force. The corresponding Schrödinger equation is transformed into a Heun equation in canonical form, and exact solutions are obtained in terms of series of hypergeometric functions. The separation of variables is still possible since the vector potential depends only on one spherical variable, θ in this case [28]. Another example of the same type of difficulty is the two-dimensional case of interaction of three particles with a perpendicular magnetic field [29]. Here the Schrödinger equations map into biconfluent Heun’s equations, too, through the higher order singularities induced by the Coulombian interaction between particles. Similar problems related to magnetic field (finite-gap potentials, Fokker-Planck, central two-point connection, generalized central potentials up to order 1/r6, Hawking radiation, etc.) were approached in the literature and usually the resulting leading differential equation for the wave function reduces to one of Heun’s differential equations [30–34]. An interesting review and study on the use of the Heun’s type of differential equations as generalizations of the hypergeometric ones, in relation to supersymmetry, are given in [35], where a two-Coulomb-center problem is solved based on a self-adjoint separation of coordinates in prolate spheroidal coordinates.Before moving to the dipole coordinates we perform a qualitative analysis of Eq. (6) in spherical coordinates. In spherical coordinates Eqs. (6) and (10) reduce to a linear Schrödinger equation in the form(13)△r+1r2△Ω+1r2-2iχr∂∂φ-χ2sin2θr2+EΨ=0,where we denoted χ=2πqμ/(ℏc), and the first two terms in the LHS are the radial and angular parts of the Laplace operator in spherical coordinates, respectively. We look for solutions of Eq. (13) in the subspaces described by Eq. (9), meaning that the wave functions have good quantum numbers for the angular momentum L. We can write these wave functions in the form(14)Ψr,θ,φ=Φr,θeiLφ,with L being an arbitrary positive integer. Eq. (13) reduces to a 2-dimensional Schrödinger equation with an effective potential in the form(15)Vr,θ=-L2sin2θ+2χLr3-χ2sin2θr4.In Figures 1 and 2 we present some equipotential contours of the effective potential V(x,z) at y=0 and y=0.5. This potential has a repulsive tail that drops to zero towards infinity like r-4. Depending on θ the potential generates a finite barrier in a neighborhood of the origin. The resulting effective potential valley produces a finite spectrum of energies for the linearized equation, hence a finite space of eigenfunctions available for the nonlinear GL equation. The parameter χ controls the barrier height. For θ~π/4 the barrier reduces to zero, and actually the potential is almost everywhere attractive. This is normal, since without magnetic dipole the spectrum is continuous, and the only possible state is spherical isotropic without any vortex. For θ~0,π the barrier height increases and allows the accumulation of more eigenstates in the linear spectrum. This denser spectrum generates a higher probability of formation of open vortices that spring towards the sphere surface around the poles of the z-axis. Next to the origin r=0 the potential has an (attractive) infinite depth valley. For small dipole strength values the potential becomes pure repulsive, and for large values of the dipole strength the potential becomes almost attractive everywhere.Figure 1
Equipotential contours of the effective potentialV(x,z) at y=0, as given by Eq. (15). The dipole is along z-axis. Several values for dipole strength χ and L. For L=0 the potential geometry is the same for any dipole strength, but for higher magnetic field the equipotential shapes become complicated and self-intersecting.Figure 2
Equipotential contours of the effective potentialV(x,z) at y=0.5, as given by Eq. (15). Dipole along z-axis. Far from center the isopotential surfaces are deformed spheres.The Schrödinger equation with the potential in Eq. (15) cannot be resolved by traditional expansion in spherical harmonics with coefficients given by the monomials rl. Since this equation contains nonhomogenous terms as powers of r, even if we use Laurent polynomials instead of these monomials, the coefficients of YlL will always contain irretrievable powers in r. This situation enforces the coefficients of spherical harmonics to include all powers of r. However, even if we use arbitrary (linear independent) functions fl(r) instead of constant coefficients for each YlL, the differential equations resulting for each fl will contain infinite many terms. This happens because the term sin2θ multiplied with spherical harmonics of different orders decomposes (through the Wigner 3j-symbols) in a sum of spherical harmonics of orders smaller than l. Consequently, spherical coordinates do not map Eq. (6) into an integrable system.Given the symmetry of the magnetic dipole we introduce a suitable type of orthogonal curvilinear coordinates, namely,dipole coordinates denoted by (a,b,φ). These dipole coordinates are related to the spherical ones by the relations(16)a=rsin2θ,b=r2cosθ,φ=φ.Dipole coordinates can be introduced in a variety of ways [36], and they are useful in a couple of physical systems controlled by magnetic dipolar terms. That is where the lines of force have strong anisotropy, like terrestrial ionosphere, solar corona, magnetostars, toroidal magnetic moments in atomic physics, etc. [37, 38]. In the following, we use the dipole coordinates in order to provide exact analytic solutions for the order parameter Ψ(r→) in the configuration magnetic dipole plus superconducting sphere. The dipole coordinates were used in [39] to solve the LGL equation for a magnetic dipole placed in a superconducting space. In that study we demonstrate that analytic solutions of the LGL equation in dipole coordinates can be obtained by transforming it in a double confluent Heun’s differential equation. Combination of such linear solutions can numerically minimize the nonlinear GL free energy functional for multivortex states as confined vortex loops.In the present paper we develop this method even further and, in addition, we apply Dirichlet type of boundary conditions on the surface of a finite superconducting sphere. This procedure maps the dipole equation plus boundary conditions into a regular Sturm-Liouville problem, for which we obtain and discuss in detail the resulting spectra and eigenfunctions. This analytic method was previously used and tested in comparison with numerical calculations and experimental results for the simpler geometry of a finite superconducting cylinder in exterior magnetic field [40, 41].The dipole coordinates, their coordinate surfaces, and other properties are described in more detail in AppendixA. From there we notice that the field lines of the magnetic infinitesimal dipole,(17)h→=μ→·r^r^-μ→r3,are given by the family of curves a constant. Surfaces defined by constant a follow the field lines of the magnetic dipole and are orthogonal to the surfaces of constant b. Consequently, the order parameter has a stronger dependence on a than on b, such that the surfaces Ψ=const. are very close to the dipole coordinates surfaces a=constant. Along this reasonable approximation we can neglect the b dependence of the order parameter solutions.In the dipole coordinates the electromagnetic momentum has a simpler expression(18)P→=ℏi∇→a,b,c-qμcae→φ.By introducing Eqs. (16) and (18) in Eq. (6), under the hypothesis Eq. (14), we obtain the following partial differential equation in dipole variables (a,b) for Φ(r,θ)→Φ(a,b):(19)ar31+3r4b2d2Φda2+41+3r4/b21adΦda+1+3r4b2b2r6ddbb2dΦdb-ar3L-χa2Φ+EΦ=0,where L∈N, E are free parameters. We notice that it is not possible to fully eliminate the coordinate r from Eq. (19), because the inverse transformation (a,b)→(r,θ) involves an algebraic equation of order 4 whose solutions would introduce prohibitive long expressions in Eq. (19); see Eq. (A.3) in the Appendix A.In order to obtain a separation of the dipole variables(a,b) in Eq. (19) we investigate solutions in the approximation a~r (r being the radial spherical coordinate), which implies r2/b<<1 from Eq. (A.3). This approximation describes the order parameter in neighborhoods of θ=0, and θ=π, that is around the poles of the bispherical surfaces a=constant; see  Figure 11. Actually, this approximation is valid within a very large part of the superconducting sphere. If we perform the inverse transformation from dipole to spherical coordinates we can express r as the product between a and a power series in the parameter ϵ=a2/b, according to Eq. (A.7) in Appendix A. For example, if we just choose ϵ=0.2 we have a relative error of only 4% or less for a~r within a range of the azimuth angle θ∈[0,80°]∪[100°,180°], which represents that 79% of the total volume of the sphere provides very good agreement with this approximation. From the geometrical perspective, this approximation limits our calculations to surfaces a=const. that are very close to the spherical surface r=const,  Figure 11.In the approximationa~r Eq. (19) reduces to the form(20)∂2Φ∂a2+4a∂Φ∂a+b2a6∂∂bb2∂Φ∂b-1a2L-χa2Φ+EΦ=0,and we can factorize further the b-dependence from the a-dependence by the substitution(21)Φa,b=Qae-c/b+C0ec/b,where c,C0 are arbitrary coupling constants between the equations in a and b. Consequently, we can obtain an equation for Q(a) in the same way it was obtained in Eq. 910) in [39](22)a6d2Qda2+4a5dQda+c2-χ2a2+2Lχa3-L2a4+Ea6Q=0,This equation is associated with the boundary conditions Eq. (12)(23)lima→0Qa=0,dQdaa∞=0,namely, the order parameter needs to cancel at the center of the mesoscopic volume, and the normal component of the superconducting current has to cancel at the external boundary of this volume. In agreement with the approximation r~a the external boundary condition can be written as a=a∞=R. At this point we describe the mesoscopic volume taken for the present study. The equation r=Rsin2θ is a special case of the so-called rose curve, firstly described by the seventeenth century Italian mathematician, Guido Grandi. Therefore the mesoscopic volume has an “apple” shape, such that at the north (θ≈0°) and south (θ≈180°) poles the volume is depleted (r≈0) but not at the equator region (r=R). This choice of volume is key to achieve separation of variables in these coordinates. However this choice does not compromise our major results, namely, to prove the existence of confined loops in any mesoscopic volume. The boundary conditions Eq. (23) are called “separated” since each involves only one of the ends of the domain of definition. Eq. (22) is an ordinary differential equation of order two with meromorphic coefficients. It has two irregular singularities of rank 2 at a=0,∞, [42], and consequently there are no convergent solutions in terms of Frobenius series. We need to solve this equation in the interval a∈0,a∞, where a∞=R is the radius of the superconducting volume.The approach used to solve the homogenous two-point boundary value problem Eqs. (22) and (23) consists in performing a series of transformations that map Eq. (22) into a symmetric canonical form of a double confluent Heun’s equation [27, 43–45].


## 4. Solutions of the Dipole Equation

### 4.1. Qualitative Analysis of the Dipole Equation

In order to perform a qualitative analysis of Eq. (22) we can transform it into other types of differential equations (e.g., Q(a)→Q~(a)=apQ(a) for some real values of p) with some known physical significance. For example, by the substitution Q~(a)=a2Q(a) we have(24)-d2Q~da2+V1daQ~+EQ~=0,which is a 1-dimensional Schrödinger equation with an effective potential(25)V1d=-c2a6+χ2a4-2Lχa3+L2+2a2.Here the last term represents the traditional centrifugal term plus an additional constant 2. This potential is repulsive for large values of a, and it becomes attractive in a short range close to the origin because of the coupling with the exp(-c/b)+C0exp(c/b) part of the wave function, through c.Another possible approach is provided by the substitutionA~(a)=aA(a), and we obtain the radial part of a Schrödinger equation in spherical coordinates(26)-d2Q~da2-2adQ~da+-c2a6+χ2a4-2Lχa3+L2+2a2-EQ~=0.The effective potential for this case (the first four terms in the parenthesis of Eq. (26)) is presented in Figure 3. It has a repulsive asymptotic behavior for Vsph(+∞)=0, and Vsph(0)→±∞ depending on the value of c. In the simpler case c=0, for example, if L≥4 this potential has a valley centered at(27)a1=χ3L-L2-162L2+2.This valley creates a barrier of height(28)Vbar=8L2+23-4+L2+LL2-16χL3L+L2-164,and relative depth with respect to the top of the valley(29)V0=LL2-163/216χ2.From the aspect of this potential it results that for weak a-b coupling (c~0) the expected spectrum, in the absence of the boundary conditions, is continuous and positive. This is the case of very large superconducting volumes. By its richness of eigenstates, this configuration increases the chances for vortex creation. Indeed, for large superconducting samples it is known that more magnetic flux is expelled in the Meissner phase; hence superconductivity is suppressed at the boundaries. Consequently, more weak points are created that facilitate the entry of many vortices [9]. At higher values of angular quantum number and dipole intensity some resonant states may occur, like in the potential valley shown in Figure 3 for L=5,χ=10.Figure 3
The effective potentialVsph(a) in the radial part of the Schrödinger equation in spherical coordinates, Eq. (26). Solid lines represent c=0 potentials, and dotted lines represent c=5 potentials, for the same L,χ.For the substitutionQ~=a3/2Q(a) we obtain the radial part of the Schrödinger equation in cylindrical coordinates(30)-d2Q~da2-1adQ~da+-c2a6+χ2a4-2Lχa3+L2+9/4a2-EQ~=0.Eq. (30) is a Fuchsian differential equation with singularities at the point at infinity, and order six singularity at a=0, presenting some similarities with the symmetries of the Bessel equation. Since the existence and uniqueness of this type of equation will be investigated elsewhere and since later on we will justify that the value of the decoupling parameter c is very small, we will further reduce this equation to a fourth-order singularity Fuchsian equation.As we mentioned in Section3, Eqs. (21), (22), and (30), we introduce the second hypothesis (in addition to r∽a) by considering a weak coupling between the a and the b coordinates in the order parameter dependence. This constraint is equivalent to c~0, and it arises somehow naturally since the order parameter Ψ(a,b,φ) has its dominant dependence on the a coordinate. This happens because the Ψ2 iso-surfaces are basically laying along the dipole surface coordinates a=constant. Under this hypothesis, Eq. (22) becomes(31)d2Qda2+4adQda+-χ2a4+2Lχa3-L2a2+EQ=0,and we will refer further to Eq. (31) as thedipole equation. The asymptotic solutions in zero and at infinity for Eq. (31) are described in [39] in terms of Bessel, Kummer, and Tricomi special functions, but here we are not interested in the asymptotical behavior of the solution at zero or infinity. An interesting limiting case is obtained if we keep only the even power terms 1/a4,1/a2, E, in the potential energy. The resulting equation is mapped into the Mathieu equation for A(z) through the transformation a=(-χ2/E)1/4exp(iz). It is interesting that even holding all the terms in the powers of 1/a we still can map the dipole equation into Hill’s differential equation [46].The existence of situations when one can map the dipole equation in traditional equations of mathematical physics (Bessel, confluent hypergeometric, Mathieu, etc.) by neglecting higher or lower powers in1/a has a deep meaning. Although not obvious, Eq. (31) has an intrinsic symmetry that allows mapping it into a very symmetrical form, called thesymmetrical-canonical double confluent Heun equation [27, 43, 44]. We will discuss in Section 4.2 multiple advantages of expressing the dipole equation in this form, especially related to the existence of simple analytic solutions. When we neglect some of the terms in the dipole equation, certain conditions of nondegeneracy fail, and mapping the dipole equation into the symmetrical-canonical DCHE becomes impossible. Hence, under such approximations, Eq. (31) degenerates into those classical equations mentioned above.


### 4.2. Analytic Solutions for the Dipole Equation

In the following we construct analytic solutions for the dipole equation, in the form of Eq. (22), or Eq. (31), in the limit of weak coupling with the orthogonal b coordinate, that is, for c=0. Also in this section we will not take into account yet the boundary conditions and just obtain the most general solution sets.The first step is to reduce Eq. (22) to a dimensionless form by the transformation a→z=a/χ. The resulting equation is a general double confluent Heun equation (gDCHE [27]) in Q(z) of the form(32)D2Q+pzDQ+qzQ=0,pz=∑i=-11pizi,qz=∑i=-22qizi,where the operator D=z(d/dz) and pi,qi are arbitrary complex coefficients. The properties and solutions of different confluent types of Heun’s equation have been studied extensively in the last decade [31, 47–52]. There is also an extensive volume of applications of Heun functions and equations in physics; see, for example, [28]. In our case for the dipole equation we have p-1=p1=q1=0,p0=3,q-2=-1,q-1=2L,q0=-L2 and q2=λ, where λ=E/χ2. By using a transformation of the dependent variable of the form(33)Qz=zξ0eξ1z-ξ-1/zyzany gDCHE (for example, Eq. (32)) can be transformed into the so-called (nondegenerate) canonical form of the DCHE [27, 43]. For the dipole equation there are two possible transformations that map it into a DCHE canonical form, namely, ξ0=-3/2,ξ1=±λ,ξ-1=±1. It is worth mentioning that one obtains a nondegenerate (a “good”) canonical DCHE form if one has the condition(34)q2-p122q-2-p-122≠0,fulfilled (see Eq. (32)). This is actually the case for the dipole equation, but if we neglect some of its terms this condition fails, and we have a degenerated canonical DCHE that reduces to the traditional differential equations mentioned in Section 4.1.The last step towards integrating the equation is to apply one of the similarity transformations in the independent variable:z=(i/λ1/4)z~, where the multiplicity comes from the square roots of the imaginary/negative parameters. After this transformation, the dipole equation has the form of the symmetric canonical DCHE for the function y(z~)(35)D~2y+αz~+1z~D~y+β1+12αz~+α22-γ+β-1-12αz~y=0,where D~=z~(d/dz~), and(36)α=±2iiλ1/4,β1=0,β-1=∓L,γ=L2+94,are the singular parameters of the symmetric canonical DCHE associated with the dipole equation. As one can note, there are several (more than two, but less than eight) representations of the dipole equation into the symmetrical form Eq. (35).Using this symmetrical canonical form has several advantages [27]. First, it generates a recursion relation for the series coefficients of the analytic solutions in only three terms. This is a major advantage over all other representations. Second, it depends on less parameters, just four instead of eight, like in the case of original gDCHE. These four parameters control the leading terms of the asymptotical expansion of the solutions. The singular parameters α,β±1 determine the leading coefficients of the asymptotical series of the solutions at ∞ and 0, respectively. However, in the case of the dipole equation we have β1=0, and only the reduced energy (λ) controls the asymptotical behavior of the solution at infinity, through α. This behavior is in agreement with the fact that the boundary condition at a∞ is the only responsible constraint for energy quantization in prescribed levels. The singular parameters α,β-1 determine the leading coefficients of the asymptotic solutions around zero, which means that in the origin the solution depends on both E and L. The γ “accessory” parameter just influences the monodromy behavior of the solutions, and it is also responsible for the eigenvalue problem related to some appropriate boundary conditions. Last advantage of the symmetric canonical form is that the central connection problem, i.e., the relations between asymptotical solutions at different singularities, becomes trivial. In that, changing the asymptotic series from the singularity 0 to ∞ just reduces to permuting the β parameters, modulo some gauge of multiplication with a constant phase [44].The symmetrical canonical form of the DCHE emerging from the dipole equation, Eq. (35), also has two irregular singularities of rank 1, namely, z~=0,∞. That means that, according to the general theory of complex ordinary differential equations with meromorphic coefficients, all solutions can be continued analytically along any path within C∗=C∖{0}. Hence,we can obtain fundamental analytic solutions as series around the singularities with a certain prescribed asymptotic behavior. This asymptotic behavior will also provide the Schmidt-Wolf uniqueness of the fundamental set. Either set of fundamental solutions can be continued analytically over C∗, though because C∗ is not simply connected the Poincaré lemma conditions are not fulfilled, and these resulting global solutions will not be single valued. This problem can be fixed by specifying the asymptotic behavior within special sectors of C∗, bounded by the corresponding Stokes rays, depending on the argument of αz~ [45].Since the physical boundary conditions ata∞, Eq. (12), are important for the eigenvalue problem for energy we will construct the solutions of Eq. (35) based on the asymptotic behavior for z~→∞. We choose the specific asymptotic behavior in such a way that will match the asymptotic behavior of the dipole equation in the large a approximation and the corresponding Bessel solutions. This asymptotic behavior secures both uniqueness of the series solution through its Laplace transform and the application of Watson’s lemma [27] and the possibility of applying the von Neumann boundary condition at a∞. Then, this solution can be continued analytically towards the solution around z~=0. The continuation is performed through the central connection relations between asymptotic solutions, such that the behavior at zero of the fundamental solution fulfills the boundary condition in zero, too, Eq. (12).The symmetric canonical DCHE form for the dipole equation is also useful because it has a simple group of transformations which maps any solution into another solution. In this way we can generate from any certain fundamental solution, constructed a series, another linear independent solution, however not necessarily linear independent. The group is infinite and it is generated by the transformations(37)T1yα,β1,β-1,γ;z=exp-αzyiα,β-1,-β1,γ;iz,T-1y=expαzyiα,-β-1,β1,γ;-iz,and(38)T0y=exp-αz+αzyα,β1,β-1,γ;1z.The two fundamental solutions for DCHE are built by using the asymptotic construction of solutions of Eq. (35), following the receipt in [27], through the procedure of expanding around the point at infinity described in detail in [39]:(39a)y1z~=Hlα,β1,γ;z~,(39b)y2z~=ie-iπβ1e-αz~Hleiπ/2α,β-1-β1,γ;eiπ/2z~,where Hl is the holomorphic local Heun function as a solution of the DCHE (35) given by the series(40)Hlα,β,γ;z~=αz~-β-1/2∑n=0∞Cnαz~n,in the range on arg(αz~)∈(-π/2,5π/2), with the coefficients Cn=Cn(α,β1,β-1,γ) uniquely given by the three-term recursion relation(41)Cn±=1nα22-γ+n+β1-122+Cn-1±-α2n+β1-β-1-1Cn-2±,with n∈N, C-1±=0,C0±=1. The signs ± represent the choice for one of the two solutions for the reduced parameters α,β-1 in Eq. (36). We mention that the second solution y2 is generated from the first one by application of the group operator T1 on y1, and it is linear independent from this one. Finally, one can perform the variable substitutions in reverse order and obtain the solutions Eqs. (39a) and (39b) in terms of the physical parameters a,χ,E,L. The Heun functions in Eqs. (39a) and (39b) represent the fundamental set of exact solutions of the linearized dipole equation (31).An interesting feature of these solutions is that in the range of negative energiesE<0 the solutions y1(z~), that is, Q1(a), are bounded and have no oscillations, yet they can be truncated to quasipolynomials. Such special truncations, very useful for analytic and numeric calculation, are possible only for the series coefficients that have β-1=-L and only for odd values for L. For each such odd L value there is one unique value of the energy which provides the truncation of the series into a quasipolynomial (the series reduces to a polynomial, but the exponential and power in front of the polynomial do not vanish). This limitation occurs because of the special structure of the dipole ODE, namely, because of the linear term 2χL/a3. For Q2, however, the majority of solutions have E>0, so they cannot be reduced to quasipolynomials.In the following, we present two examples of quasitruncated solutionsQ1(a) obtained from the abstract solutions y1(z~):(42)Q1L=1,a=χ23a2e-3a/2χ-χ/a,with E=-9χ2/4, and(43)Q1L=3,a=χ2w1a2e-w1a/2χ-χ/a1+C1-χw1a+C2-χw1a2,with E=-χ2w12/4. In order to have a three-dimensional visual on the behavior of these solutions we plot in Figure 4 few surfaces of equation Ψa,φ=Ψr/sin2θ,φ=Q1r/sin2θeiLφconst. The three surfaces presented in Figure 4 follow the shape of the surfaces of coordinates a=const., proving the hypothesis according to which the b dependence of the order parameter can be neglected.Figure 4
The linear solutionsΨ2=0.5, 0.6, 0.75, and 0.9 isoplot surfaces for χ=0.5, a∞=10, L=1, calculated with the quasi-polynomial solutions Eq. (42) for E=-5.3. The surfaces are tori-like, as expected from the behavior of a=const. surface coordinates,  Figure 11. Some of the exterior surfaces are plotted only partially, in order to make visible the torus shape of the inner surface.


## 4.1. Qualitative Analysis of the Dipole Equation

In order to perform a qualitative analysis of Eq. (22) we can transform it into other types of differential equations (e.g., Q(a)→Q~(a)=apQ(a) for some real values of p) with some known physical significance. For example, by the substitution Q~(a)=a2Q(a) we have(24)-d2Q~da2+V1daQ~+EQ~=0,which is a 1-dimensional Schrödinger equation with an effective potential(25)V1d=-c2a6+χ2a4-2Lχa3+L2+2a2.Here the last term represents the traditional centrifugal term plus an additional constant 2. This potential is repulsive for large values of a, and it becomes attractive in a short range close to the origin because of the coupling with the exp(-c/b)+C0exp(c/b) part of the wave function, through c.Another possible approach is provided by the substitutionA~(a)=aA(a), and we obtain the radial part of a Schrödinger equation in spherical coordinates(26)-d2Q~da2-2adQ~da+-c2a6+χ2a4-2Lχa3+L2+2a2-EQ~=0.The effective potential for this case (the first four terms in the parenthesis of Eq. (26)) is presented in Figure 3. It has a repulsive asymptotic behavior for Vsph(+∞)=0, and Vsph(0)→±∞ depending on the value of c. In the simpler case c=0, for example, if L≥4 this potential has a valley centered at(27)a1=χ3L-L2-162L2+2.This valley creates a barrier of height(28)Vbar=8L2+23-4+L2+LL2-16χL3L+L2-164,and relative depth with respect to the top of the valley(29)V0=LL2-163/216χ2.From the aspect of this potential it results that for weak a-b coupling (c~0) the expected spectrum, in the absence of the boundary conditions, is continuous and positive. This is the case of very large superconducting volumes. By its richness of eigenstates, this configuration increases the chances for vortex creation. Indeed, for large superconducting samples it is known that more magnetic flux is expelled in the Meissner phase; hence superconductivity is suppressed at the boundaries. Consequently, more weak points are created that facilitate the entry of many vortices [9]. At higher values of angular quantum number and dipole intensity some resonant states may occur, like in the potential valley shown in Figure 3 for L=5,χ=10.Figure 3
The effective potentialVsph(a) in the radial part of the Schrödinger equation in spherical coordinates, Eq. (26). Solid lines represent c=0 potentials, and dotted lines represent c=5 potentials, for the same L,χ.For the substitutionQ~=a3/2Q(a) we obtain the radial part of the Schrödinger equation in cylindrical coordinates(30)-d2Q~da2-1adQ~da+-c2a6+χ2a4-2Lχa3+L2+9/4a2-EQ~=0.Eq. (30) is a Fuchsian differential equation with singularities at the point at infinity, and order six singularity at a=0, presenting some similarities with the symmetries of the Bessel equation. Since the existence and uniqueness of this type of equation will be investigated elsewhere and since later on we will justify that the value of the decoupling parameter c is very small, we will further reduce this equation to a fourth-order singularity Fuchsian equation.As we mentioned in Section3, Eqs. (21), (22), and (30), we introduce the second hypothesis (in addition to r∽a) by considering a weak coupling between the a and the b coordinates in the order parameter dependence. This constraint is equivalent to c~0, and it arises somehow naturally since the order parameter Ψ(a,b,φ) has its dominant dependence on the a coordinate. This happens because the Ψ2 iso-surfaces are basically laying along the dipole surface coordinates a=constant. Under this hypothesis, Eq. (22) becomes(31)d2Qda2+4adQda+-χ2a4+2Lχa3-L2a2+EQ=0,and we will refer further to Eq. (31) as thedipole equation. The asymptotic solutions in zero and at infinity for Eq. (31) are described in [39] in terms of Bessel, Kummer, and Tricomi special functions, but here we are not interested in the asymptotical behavior of the solution at zero or infinity. An interesting limiting case is obtained if we keep only the even power terms 1/a4,1/a2, E, in the potential energy. The resulting equation is mapped into the Mathieu equation for A(z) through the transformation a=(-χ2/E)1/4exp(iz). It is interesting that even holding all the terms in the powers of 1/a we still can map the dipole equation into Hill’s differential equation [46].The existence of situations when one can map the dipole equation in traditional equations of mathematical physics (Bessel, confluent hypergeometric, Mathieu, etc.) by neglecting higher or lower powers in1/a has a deep meaning. Although not obvious, Eq. (31) has an intrinsic symmetry that allows mapping it into a very symmetrical form, called thesymmetrical-canonical double confluent Heun equation [27, 43, 44]. We will discuss in Section 4.2 multiple advantages of expressing the dipole equation in this form, especially related to the existence of simple analytic solutions. When we neglect some of the terms in the dipole equation, certain conditions of nondegeneracy fail, and mapping the dipole equation into the symmetrical-canonical DCHE becomes impossible. Hence, under such approximations, Eq. (31) degenerates into those classical equations mentioned above.


## 4.2. Analytic Solutions for the Dipole Equation

In the following we construct analytic solutions for the dipole equation, in the form of Eq. (22), or Eq. (31), in the limit of weak coupling with the orthogonal b coordinate, that is, for c=0. Also in this section we will not take into account yet the boundary conditions and just obtain the most general solution sets.The first step is to reduce Eq. (22) to a dimensionless form by the transformation a→z=a/χ. The resulting equation is a general double confluent Heun equation (gDCHE [27]) in Q(z) of the form(32)D2Q+pzDQ+qzQ=0,pz=∑i=-11pizi,qz=∑i=-22qizi,where the operator D=z(d/dz) and pi,qi are arbitrary complex coefficients. The properties and solutions of different confluent types of Heun’s equation have been studied extensively in the last decade [31, 47–52]. There is also an extensive volume of applications of Heun functions and equations in physics; see, for example, [28]. In our case for the dipole equation we have p-1=p1=q1=0,p0=3,q-2=-1,q-1=2L,q0=-L2 and q2=λ, where λ=E/χ2. By using a transformation of the dependent variable of the form(33)Qz=zξ0eξ1z-ξ-1/zyzany gDCHE (for example, Eq. (32)) can be transformed into the so-called (nondegenerate) canonical form of the DCHE [27, 43]. For the dipole equation there are two possible transformations that map it into a DCHE canonical form, namely, ξ0=-3/2,ξ1=±λ,ξ-1=±1. It is worth mentioning that one obtains a nondegenerate (a “good”) canonical DCHE form if one has the condition(34)q2-p122q-2-p-122≠0,fulfilled (see Eq. (32)). This is actually the case for the dipole equation, but if we neglect some of its terms this condition fails, and we have a degenerated canonical DCHE that reduces to the traditional differential equations mentioned in Section 4.1.The last step towards integrating the equation is to apply one of the similarity transformations in the independent variable:z=(i/λ1/4)z~, where the multiplicity comes from the square roots of the imaginary/negative parameters. After this transformation, the dipole equation has the form of the symmetric canonical DCHE for the function y(z~)(35)D~2y+αz~+1z~D~y+β1+12αz~+α22-γ+β-1-12αz~y=0,where D~=z~(d/dz~), and(36)α=±2iiλ1/4,β1=0,β-1=∓L,γ=L2+94,are the singular parameters of the symmetric canonical DCHE associated with the dipole equation. As one can note, there are several (more than two, but less than eight) representations of the dipole equation into the symmetrical form Eq. (35).Using this symmetrical canonical form has several advantages [27]. First, it generates a recursion relation for the series coefficients of the analytic solutions in only three terms. This is a major advantage over all other representations. Second, it depends on less parameters, just four instead of eight, like in the case of original gDCHE. These four parameters control the leading terms of the asymptotical expansion of the solutions. The singular parameters α,β±1 determine the leading coefficients of the asymptotical series of the solutions at ∞ and 0, respectively. However, in the case of the dipole equation we have β1=0, and only the reduced energy (λ) controls the asymptotical behavior of the solution at infinity, through α. This behavior is in agreement with the fact that the boundary condition at a∞ is the only responsible constraint for energy quantization in prescribed levels. The singular parameters α,β-1 determine the leading coefficients of the asymptotic solutions around zero, which means that in the origin the solution depends on both E and L. The γ “accessory” parameter just influences the monodromy behavior of the solutions, and it is also responsible for the eigenvalue problem related to some appropriate boundary conditions. Last advantage of the symmetric canonical form is that the central connection problem, i.e., the relations between asymptotical solutions at different singularities, becomes trivial. In that, changing the asymptotic series from the singularity 0 to ∞ just reduces to permuting the β parameters, modulo some gauge of multiplication with a constant phase [44].The symmetrical canonical form of the DCHE emerging from the dipole equation, Eq. (35), also has two irregular singularities of rank 1, namely, z~=0,∞. That means that, according to the general theory of complex ordinary differential equations with meromorphic coefficients, all solutions can be continued analytically along any path within C∗=C∖{0}. Hence,we can obtain fundamental analytic solutions as series around the singularities with a certain prescribed asymptotic behavior. This asymptotic behavior will also provide the Schmidt-Wolf uniqueness of the fundamental set. Either set of fundamental solutions can be continued analytically over C∗, though because C∗ is not simply connected the Poincaré lemma conditions are not fulfilled, and these resulting global solutions will not be single valued. This problem can be fixed by specifying the asymptotic behavior within special sectors of C∗, bounded by the corresponding Stokes rays, depending on the argument of αz~ [45].Since the physical boundary conditions ata∞, Eq. (12), are important for the eigenvalue problem for energy we will construct the solutions of Eq. (35) based on the asymptotic behavior for z~→∞. We choose the specific asymptotic behavior in such a way that will match the asymptotic behavior of the dipole equation in the large a approximation and the corresponding Bessel solutions. This asymptotic behavior secures both uniqueness of the series solution through its Laplace transform and the application of Watson’s lemma [27] and the possibility of applying the von Neumann boundary condition at a∞. Then, this solution can be continued analytically towards the solution around z~=0. The continuation is performed through the central connection relations between asymptotic solutions, such that the behavior at zero of the fundamental solution fulfills the boundary condition in zero, too, Eq. (12).The symmetric canonical DCHE form for the dipole equation is also useful because it has a simple group of transformations which maps any solution into another solution. In this way we can generate from any certain fundamental solution, constructed a series, another linear independent solution, however not necessarily linear independent. The group is infinite and it is generated by the transformations(37)T1yα,β1,β-1,γ;z=exp-αzyiα,β-1,-β1,γ;iz,T-1y=expαzyiα,-β-1,β1,γ;-iz,and(38)T0y=exp-αz+αzyα,β1,β-1,γ;1z.The two fundamental solutions for DCHE are built by using the asymptotic construction of solutions of Eq. (35), following the receipt in [27], through the procedure of expanding around the point at infinity described in detail in [39]:(39a)y1z~=Hlα,β1,γ;z~,(39b)y2z~=ie-iπβ1e-αz~Hleiπ/2α,β-1-β1,γ;eiπ/2z~,where Hl is the holomorphic local Heun function as a solution of the DCHE (35) given by the series(40)Hlα,β,γ;z~=αz~-β-1/2∑n=0∞Cnαz~n,in the range on arg(αz~)∈(-π/2,5π/2), with the coefficients Cn=Cn(α,β1,β-1,γ) uniquely given by the three-term recursion relation(41)Cn±=1nα22-γ+n+β1-122+Cn-1±-α2n+β1-β-1-1Cn-2±,with n∈N, C-1±=0,C0±=1. The signs ± represent the choice for one of the two solutions for the reduced parameters α,β-1 in Eq. (36). We mention that the second solution y2 is generated from the first one by application of the group operator T1 on y1, and it is linear independent from this one. Finally, one can perform the variable substitutions in reverse order and obtain the solutions Eqs. (39a) and (39b) in terms of the physical parameters a,χ,E,L. The Heun functions in Eqs. (39a) and (39b) represent the fundamental set of exact solutions of the linearized dipole equation (31).An interesting feature of these solutions is that in the range of negative energiesE<0 the solutions y1(z~), that is, Q1(a), are bounded and have no oscillations, yet they can be truncated to quasipolynomials. Such special truncations, very useful for analytic and numeric calculation, are possible only for the series coefficients that have β-1=-L and only for odd values for L. For each such odd L value there is one unique value of the energy which provides the truncation of the series into a quasipolynomial (the series reduces to a polynomial, but the exponential and power in front of the polynomial do not vanish). This limitation occurs because of the special structure of the dipole ODE, namely, because of the linear term 2χL/a3. For Q2, however, the majority of solutions have E>0, so they cannot be reduced to quasipolynomials.In the following, we present two examples of quasitruncated solutionsQ1(a) obtained from the abstract solutions y1(z~):(42)Q1L=1,a=χ23a2e-3a/2χ-χ/a,with E=-9χ2/4, and(43)Q1L=3,a=χ2w1a2e-w1a/2χ-χ/a1+C1-χw1a+C2-χw1a2,with E=-χ2w12/4. In order to have a three-dimensional visual on the behavior of these solutions we plot in Figure 4 few surfaces of equation Ψa,φ=Ψr/sin2θ,φ=Q1r/sin2θeiLφconst. The three surfaces presented in Figure 4 follow the shape of the surfaces of coordinates a=const., proving the hypothesis according to which the b dependence of the order parameter can be neglected.Figure 4
The linear solutionsΨ2=0.5, 0.6, 0.75, and 0.9 isoplot surfaces for χ=0.5, a∞=10, L=1, calculated with the quasi-polynomial solutions Eq. (42) for E=-5.3. The surfaces are tori-like, as expected from the behavior of a=const. surface coordinates,  Figure 11. Some of the exterior surfaces are plotted only partially, in order to make visible the torus shape of the inner surface.


## 5. Eigenvalue Problem for the Dipole Equation

Typical study of the quantization of the energy spectrum of Heun’s equation involves theγ accessory parameter. Such an analysis is performed, for example, in [45] where the authors study the eigenvalues of a double confluent Heun equation with possible applications in gravitational theory. The boundary conditions used in this paper differ from our case, since they ask the solution to be bounded in origin, and to approach either a finite value, or infinity when z→∞, but not faster than an exponential. Consequently, the [45] boundary conditions quantize the accessory parameter γ, and the recursion relation for the coefficients of the series becomes a four-term Poincaré-Perron type, which is more difficult to handle than our three-term recursion relations.In order to solve our eigenvalue problem forQ(a) and E from Eqs. (22) and (23) we substitute a=a∞x, y(x)=a∞xQ(a∞x), where a∞ is the upper bound of the a variable. The resulting equation is a regular Sturm-Liouville problem. The self-adjoint form of the corresponding linearized dipole equation (22) is(44)Ly=Erxy,with(45a)Ly=-ddxpxdydx+1+χ22a∞2x2-Lχa∞x+L22y,(45b)rx=a∞2x22,px=x22.In this y(x) representation we have x∈(0,1], and p(x)>0,r(x)>0 in this range. In this case [21], the general theory of Sturm-Liouville two-point boundary value problem provides the following results. All eigenvalues E of Eq. (44) are real and form an infinite sequence E1<E2<⋯<En<… with limn→∞⁡En=∞. According to the physical restriction E<1, it results that the spectrum of interest of the linear dipole boundary problem is also finite. All eigenvalues En are not degenerate, and to every two distinct eigenvalues there correspond two linear independent eigenfunctions, orthogonal with respect to the weighted scalar product. That is, for Ek≠En we have(46)yk,yn=a∞22∫01x2ykxynxdx=0.Normalized with respect to the weight r(x), the eigenfunctions form a complete orthonormal basis of piecewise continuous functions, with piecewise square integrable continuous derivative in a∈[0,a∞],b≥0.In order to solve explicitly the eigenvalue problem forE we apply the boundary condition at a∞, Eq. (23). Since the physical range of interest consists in large values for a∞ (at least three times, to ten times larger than the coherence length, in order to identify multivortex structures) we use the asymptotic solution for the order parameter from [39], as we mentioned previously. In terms of Bessel functions, we can write the boundary conditions Eq. (23) in the form(47)dQ∞daa=a∞=N1a∞2E1/2JL~-1w-JL~+1w-χ2JL~w3a∞-2χ+N2a∞2E1/2YL~-1w-YL~+1w-χ2YL~w3a∞-2χ=0,where we denoted(48)L~=94+L2+2iE1/2χ,w=E1/2aχ2.The boundary condition Eq. (47) contains two unknowns, N1,2. One of them will be eliminated from the condition of smooth match between the Heun analytic solutions with the Bessel asymptotic approximation. Since this procedure will provide a proportionality between N1,2 they will be actually eliminated from the homogenous boundary condition. Consequently the only thing provided by the boundary condition equation Eq. (47) will be the value of the energy.On the one hand we use the Bessel asymptotic expansion [53] which gives(49)QBessel∽2χπcosE1/2aχ2cosϕL,EN1+N2-sinϕL,EN2-N1,where we denote(50)ϕL,E=π294+L2+2iE1/2χ+12.On the other hand, we calculate the asymptotic behavior towards infinity of the two linear independent analytic Heun solutions Q1,2(a) from y1,2(z~), Eqs. (39a) and (39b). It results in an asymptotic term of the form(51)2cosE1/2aχ22iE1/2χ2-1/2,to be matched with Eq. (49). The solution of this matching condition results in(52a)C1=cosϕL,E+sinϕL,E2πχiE1/2,(52b)C2=cosϕL,E-sinϕL,E2πχiE1/2,where ϕL,E has the same meaning as above, Eq. (50).With the coefficientsN1,2 obtained in Eqs. (52a) and (52b), with a choice for a∞,L,χ as parameters, we solve the boundary condition relation Eq. (47) and find the corresponding spectrum for each such set of parameters. The solution is not unique, and we can write the resulting spectrum as E=Ek(a∞,L,χ), where k is the “radial” quantum number.The corresponding eigenfunctions areQ1(a;a∞,L,χ,Ek(a∞,L,χ)). The resulting solutions of the linearized GL equation are(53)Ψk,La,φ=Q1a;a∞,L,χ,Eka∞,L,χeiLφ.We present some examples of such solutions in Figure 5 for a∞=10 and two values of χ. For each situation we obtain a discrete spectrum for E, as predicted at the beginning of this section. In the investigated range we obtained between two and seven positive eigenvalues for E. For higher dipole strength we have less positive eigenvalues. This is because the oscillatory character of the free solution (Bessel asymptotic) is reduced by the real exponential in the solution. Also, the spectrum has more eigenvalues when L increases.Figure 5
Examples of the application of the boundary conditions ata∞=10, for the dipole equation analytic solution. The positive part of the spectrum 0<E<1, for L=0,3, and χ=0.5,1. In the main window we present the whole wave function, and in the smaller window we zoom in to the zero derivative point at a=a∞. The number of eigenstates (oscillations) increases with L, and decreases with χ, because for smaller dipole strength the solution is closer to the Bessel solutions. The eigenvalues of energy are E(L,χ)=E(0,0.5)=0.06,0.32,0.78,0.92; E(3,0.5)=0.02,0.16,0.23,0.41,0.52,0.64,0.77; E(0,1)=0.085,0.44,0.985, and E(3,1)=0.32,0.88. Higher energy states have more oscillations.The final form of this equation is(54)tanπ94+L2+2iE1/2χ2-χ23a∞+2χ2a∞2E1/2=0.For any set of values for a∞,χ and L the solutions of Eq. (54) provide the E spectrum. For every L value the spectrum is discrete and bounded, according to the general Sturm-Liouville results mentioned in the beginning of this section, and the physical constraint E<1. We also mention that the finiteness of the spectrum is in agreement with the finite depth of the potential valley noticed in Figure 3.Examples of energy spectra for different values of the parameters are given in Figures6 and 7. In Figure 8 we present the eigenfunctions and spectra for small values of energy in the range 0<E<1. We note that the spectrum becomes richer with the increase of the volume radius. Also, for large enough dipole strength, the energy of a certain level decreases with the increasing of the angular momentum. The increasing of the dipole strength also produces a dilation of the spectral levels. From Eq. (54) we note that the energy spectrum for the linear problem has three distinct domains. The positive energy domain is rich in states and is the most important source of potential multivortex structures. We define it by 0<E<1. The maximum number of spectral levels kmax depends on L and it can be calculated from the condition Ψ2<1. The next spectral domain is within the negative energies(55)-χ294+L2<E<0,and it also depends on all the parameters. In the normal range of values (L=1,…,6, χ=0.2-10) this condition offers a wide range of energies (it can go down to E~-75) so this is the main range for negative energies. In this domain, because the argument of the tan function in Eq. (54) becomes imaginary, this tangent transforms into hyperbolic tangent. Consequently, for each value of L the spectrum contains just one negative energy level. The last domain is for energy less than the limit written above, which represents unstable physical situations, and we will ignore it in the following.Figure 6
Energy spectrumE for the linear dipole problem with boundary conditions on the volume from Eq. (47), or Eq. (54), for small dipole strength, χ=0.6. We present two cases of volume radii, R=a∞=10,60. The coupling of energy with angular momentum is reduced, so we notice a weak dependence of L. For larger volumes the spectrum is richer, note the difference in scale of left and right frames.Figure 7
Energy spectrumE for the linear dipole problem with boundary conditions on the volume from Eq. (47) for two dipole strength, χ=1.5,3 and two volume radii, R=a∞=10,60.Figure 8
Eigenvalues (0<E<1) and eigenfunctions for the linearized dipole equation for χ=0.5 and L=0 (upper frames) and L=3 (lower frames) with boundary condition at a∞=10. In the left frames we zoom in all the eigenfunctions for a given L around the exterior boundary. In the right frames, we plot the wave functions for the first two lowest energy values. The little windows inside the right frames show the whole spectrum of energies for these a∞,χ,L values.For weak dipole strength (basis of frames) the spectrum is rich and represents all the asymptotic oscillations of the Bessel solutions. There is always a range of the distance to the centrum after which the spectrum becomes sparser. At this linear level, it appears that the energy spectrum is not strongly dependent on the values of the angular momentum. In Figure9 we present the action of the boundary condition on the analytic solutions, for positive energies, for example. We choose a fixed size of the sample (a∞=10) and an average value of the relative dipole strength χ=0.5. The behavior of the eigenfunctions around the boundary limit is presented in the left frames, and the general aspect of the wave function all over its range is presented in the right frames, together within the spectrum.Figure 9
The number of eigenvaluesE in the positive spectrum versus a∞ for L=0,…,5. The number increases with L. The upper band of lines is for χ=1, and the lower one for χ=5. Obviously the number of eigenvalues increases with the radius of the volume since larger samples allow more oscillations inside.For small values of the magnitude of the energy, and also for large enoughχ, we can approximate the solutions of the spectral equation Eq. (54) with the analytical expression(56)Ek1/2~χ2πa∞94+L2+k,with integer values of k such that E<1. The number of eigenvalues for positive E in this approximation can be estimated(57)kmax=NeigenvaluesL=a∞χ2π-94+L2,and we have such bounded states only if the following criterion is fulfilled:(58)a∞>χ2π94+L2.This lower bound for a∞ (or equivalently, an upper bound for L) gives a criterium for the critical size, or dipole strength that can trigger the generation of multivortex structures creation. For example, if we choose a∞=20 we have 5 eigenstates for L=0,1 and 4 eigenstates for L=2,3,4.For larger radii of the volume the energy spectrum obtained as solution of Eq. (54) can be approximated with(59)Ea∞~-χ44k2-L2-942,with k integer (“radial” quantum number) labeling the energy spectrum for fixed L.The energy spectra presented so far depend only on two quantum numbers:L,k. Even if this is a full three-dimensional problem, a third quantum number was eliminated by the supposition that ∇Ψ is mainly directed orthogonal to the magnetic dipole field lines. Consequently, the order parameter is almost constant along these lines, and the coordinate surfaces a=const. describe the vortices surfaces. By this hypothesis we can neglect the b coordinate dependence in the order parameter, which eliminates the third quantum numbers, basically like taking into account just ground states with respect to b variable excitations.


## 6. Nonlinear Vortex Patterns

The goal of this paper is to obtain the nonlinear order parameter as the solutionΨ(r→) of the full nonlinear Ghinzburg-Landau problem Eq. (2), which minimizes the free energy functional F for the mesoscopic sphere with magnetic dipole, Eq. (1). This nonlinear solution, denoted ΨNL, is constructed in the following as a linear combination of exact analytic solutions of the associated linear GL problem, Eq. (6), or presented in a more rigorous form by Eqs. (8) and (10). These linear solutions form a complete orthogonal system of eigenfunctions, Eq. (53), over the space of dipole coordinates (a,b,φ), under boundary conditions, Eq. (12). The corresponding eigenvalues spectrum is described in Section 5.The spherical surface of our mesoscopic sample is expressed in dipole coordinates as a volumea∞~R. In the approximation of weak coupling with the b orthogonal variable (i.e., B(b)=constant) described above we introduce the following notation for the linear basis:(60)ΨL,Ea,b,φ→ΨL,ka,φ=Q1L,k;aeiLφ,with the connection between the physical parameter E and the spectrum label k being provided by Eq. (54).We search nonlinear solutions of the form(61)ΨNLa,φ=∑L=0Lmax∑k=1kmaxLCL,kΨL,ka,φ,with CL,k parameters that have to be determined. In order to generate physical solutions for the full nonlinear GL problem, the expression Eq. (61) is plugged in the Gibbs free energy expression Eq. (1) and this integral is minimized in the space of parameters CL,k[7, 9–18]. In order to identify the vortex structures inside the superconducting sphere, in this model, one needs to find that it exists at least a value av∈(0,a∞) of the dipole coordinate such that the equation ΨNLav,φ=0 has a number of distinct solutions for φ∈[0,2π). If we denote by φj,j=1,2,…,n these solutions, by analytic prolongation theorem we can always find a neighborhood of av on which ΨNL is arbitrary small for each of the roots φj. The surfaces described by the values of a in these neighborhoods are enveloping surfaces of the multivortex states. In this case the value of Lmax is called the vorticity of the state. The center-lines of the vortices are the curves obtained by the intersections of the a=av and φ=φj surfaces.In literature [15, 16] there are examples of such multivortex states obtained with as little as kmax(Lmax)=2 linear wave functions in Eq. (61). In [39], for example, we obtained multivortex states from the linear combination of two linear wavefunctions. We can apply here the same procedure, except we have to replace the solutions in the full space with the present solutions bounded by the spherical surface. Thus, we obtain for the angular positions of the vortices central lines (kmax=2,L∈{L1,L2}) the expression(62)φj=arccos1+C2Q1,L2,k2/C1Q1,L1,k122C2Q1,L2,k2/C1Q1,L1,k1+argC1Q1,L1,k1-argC1Q1,L2,k2±jπ·L1-L2-1with j=1,2,…,L1-L2 and all functions Q are evaluated at av. To visualize an example of multivortex state, we choose a∞=100, and L1=0,k1=0, L2=0,…,4. We plot in Figure 10 four isoplot surfaces with ΨNL=0.9 for different vorticity numbers. Nevertheless, the topology of the multivortex surfaces is controlled mainly by the linear GL equation. This happens because close to the multivortex surface the order parameter is very small, so the nonlinear becomes negligible.Figure 10
Isoplot surfacesΨNL=0.9 for L=1,2,3,4 from left upper frame, clockwise. We choose a∞=100,χ=1,L1=0,k1=1 and the minimal energy was obtained for k2=3.Figure 11
Left: coordinate surfaces of constanta=1,2,3 (tori-like); Right: coordinate surfaces of constant b=1,2,3 (double-volumes). The dipole is directed along z axis.


## 7. Conclusions

We investigate the existence of multivortex superconducting states in a mesoscopic sphere where the magnetic field is generated by an infinitesimal magnetic dipole placed at the center. We use the Ginzburg-Landau model in the approximation of weak magnetic field. The Euler-Lagrange equations emerging from the GL free energy functional reduce to a magnetic Schrödinger equation with mixed boundary conditions at the center and the surface of the superconducting sample. The goal of this paper is to obtain the nonlinear order parameter, i.e., the solutionΨ(r→) of the full nonlinear Ghinzburg-Landau problem Eq. (2), which minimizes the free energy functional F for the mesoscopic sphere with magnetic dipole, Eq. (1). The nonlinear solution is constructed as a linear combination of exact analytic solutions of the associated linear GL problem, Eqs. (6), (8), and (10). These linear solutions form a complete orthogonal Sturm-Liouville system of eigenfunctions, Eq. (53), over the space of dipole coordinates (a,b,φ), under the boundary conditions Eq. (12). As a first step, we use cylindrical symmetric coordinates. The linear solution can be separated in Ψ=Φ(a,b)eiLφ, Eq. (9), and it maps into the complete or general dipole equation, Eq. (19). Further on, by using a justified geometric approximation, we can separate further the dipole coordinates a and b, and the solution can be written as Φ(a,b)=Q(a)(e-c/b+C0ec/b), Eqs. (20), (21), (22), and (30). We reduce further the resulting equation by the mode decoupling condition c=0, and we obtain the dipole equation (31) for Q which is a double confluent Heun equation, Eq. (35). By performing a substitution, Eq. (33), in the form z=a/χ and Q(a)=zξ0eξ1z-ξ-1/zy(z) we can finally solve the dipole equation and write its analytic solutions, Eqs. (39a) and (39b).The two approximations introduced here are possible because we confine our calculations in the neighborhoods of the poles of the sphere. Since the gradient ofΨ is mainly directed radial and orthogonal to the magnetic dipole field lines, it results that Ψ is almost constant along these lines, and the dipole coordinates surfaces can describe the multivortex isosurfaces.We build the nonlinear solution for the GL problem as linear combinations of the exact solutions of the linearized dipole equation and tune the combination’s coefficients that minimize the free energy functional. Further on, we demonstrate that multivortex states are possible in this sphere plus dipole configuration, even without the presence of an external magnetic field. Our results are in very good agreement with the numerical calculation of the same model [10, 11].



---

*Source: 8652151-2018-11-01.xml*