Vortex Structures inside Spherical Mesoscopic Superconductor Plus Magnetic Dipole

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.


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][5][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 giantvortex 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][10][11][12][13][14][15][16][17][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 2 Advances in Mathematical Physics 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][10][11][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., , 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 volume .Outside of this volume the Cooper pair density Ψ(  →  ) describing the superconducting phase, called the order parameter, is zero [5,6,[9][10][11][12][13][14][15][16][17][18].The free energy functional is given in the GL model by [10,11]  ] V, (1) where  is the Cooper pairs mass,  →  is the quantum momentum operator in the presence of magnetic field, and  → ℎ is the intensity of the magnetic field.The temperature dependent coefficient function () < 0 and the nonlinear term coupling constant  > 0 are typical Landau second-order phase transition parameters [9][10][11][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][10][11][12].
The traditional procedure for finding Ψ by minimizing the free energy functional Eq. ( 1) for constant volume  consists in expanding the order parameter Ψ in a basis of eigenfunctions of the corresponding linearized GL problem [10,11,[16][17][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 Section 2 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.
where the electromagnetic momentum  →  is given by where  →  is the vector potential,  is the charge of the Cooper pair, and  →  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) We rewrite Eq. ( 6) in the London electromagnetic gauge, ∇ ⋅  →  = 0, and we rescale the variables  →   →  →  √−2()/ℎ,  →   →  → /(√−2()) and Ψ → Ψ√−/().We follow the GL phenomenological model from [9][10][11][12][13][14][15][16][17][18] and use the temperature dependence where  0 is the critical temperature in absence of magnetic field and () is the coherence length.In the new variables the LGL Eq. ( 6) can be written in the form of a linear operator equation Since the dipole magnetic field has axial symmetry we can use any axially symmetric orthogonal coordinates in azimuth angle  in the form  →  → (,,).In this case the  dependence in Eq. ( 8) can be separated, and Ô has exp () as eigenfunctions labeled by the angular momentum .This decomposes the space of solutions Ψ(  →  ) into subspaces parameterized by with Φ ∈ L 2 ( 2 ) square integrable function.The decomposition in Eq. ( 9) reduces Eq. ( 6) to a two-dimensional PDE.On each of these  subspaces the operator Ô has a particular form depending on  denoted by Ô , and for each such restricted operator we can attach an eigenvalue problem where  is a provisional degeneracy label, which later on becomes the radial quantum number.Of course,  and Φ , carry the  dependence.
The vector potential associated with the infinitesimal (point-like) magnetic dipole can be expressed in a simple form in spherical coordinates where  →   = (−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][10][11][12][15][16][17][18]] where  →  represents the outward normal to the spherical surface Σ defined by |  →  | = .
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 R 3 \ {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.

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  = 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][26][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/ 6 , 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][31][32][33][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 where we denoted  = 2/(ℎ), 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 Equipotential contours of the effective potential (, ) at  = 0, as given by Eq. ( 15).The dipole is along -axis.Several values for dipole strength  and .For  = 0 the potential geometry is the same for any dipole strength, but for higher magnetic field the equipotential shapes become complicated and self-intersecting.
Eq. ( 9), meaning that the wave functions have good quantum numbers for the angular momentum .We can write these wave functions in the form with  being an arbitrary positive integer.Eq. ( 13) reduces to a 2-dimensional Schrödinger equation with an effective potential in the form In Figures 1 and 2 we present some equipotential contours of the effective potential (, ) at  = 0 and  = 0.5.This potential has a repulsive tail that drops to zero towards infinity like  −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 Equipotential contours of the effective potential (, ) at  = 0.5, as given by Eq. ( 15).Dipole along -axis.Far from center the isopotential surfaces are deformed spheres.
the z-axis.Next to the origin  = 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.
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   .Since this equation contains nonhomogenous terms as powers of , even if we use Laurent polynomials instead of these monomials, the coefficients of   will always contain irretrievable powers in .This situation enforces the coefficients of spherical harmonics to include all powers of .However, even if we use arbitrary (linear independent) functions   () instead of constant coefficients for each   , the differential equations resulting for each   will contain infinite many terms.This happens because the term sin 2  multiplied with spherical harmonics of different orders decomposes (through the Wigner 3-symbols) in a sum of spherical harmonics of orders smaller than .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 (, , ).These dipole coordinates are related to the spherical ones by the relations 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 Ψ(  →  ) 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 Appendix A. From there we notice that the field lines of the magnetic infinitesimal dipole, are given by the family of curves  constant.Surfaces defined by constant  follow the field lines of the magnetic dipole and are orthogonal to the surfaces of constant .Consequently, the order parameter has a stronger dependence on  than on , such that the surfaces |Ψ| =const.are very close to the dipole coordinates surfaces  =constant.Along this reasonable approximation we can neglect the  dependence of the order parameter solutions.
In the dipole coordinates the electromagnetic momentum has a simpler expression By introducing Eqs. ( 16) and (18) in Eq. ( 6), under the hypothesis Eq. ( 14), we obtain the following partial differential equation in dipole variables (, ) for Φ(, ) → Φ(, ): where  ∈N,  are free parameters.We notice that it is not possible to fully eliminate the coordinate  from Eq. ( 19), because the inverse transformation (, ) → (, ) 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 (, ) in Eq. ( 19) we investigate solutions in the approximation  ∼  ( being the radial spherical coordinate), which implies  2 / << 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  =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  as the product between  and a power series in the parameter  =  2 /, 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  ∼  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  =const.that are very close to the spherical surface  =const, Figure 11.
In the approximation  ∼  Eq. ( 19) reduces to the form and we can factorize further the -dependence from the dependence by the substitution where ,  0 are arbitrary coupling constants between the equations in  and .Consequently, we can obtain an equation for () in the same way it was obtained in Eq. 910) in [39] This equation is associated with the boundary conditions Eq. ( 12) 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  ∼  the external boundary condition can be written as  =  ∞ = .At this point we describe the mesoscopic volume taken for the present study.The equation  =  sin 2  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 ( ≈ 0) but not at the equator region ( = ).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  = 0, ∞, [42], and consequently there are no convergent solutions in terms of Frobenius series.We need to solve this equation in the interval  ∈ (0,  ∞ ], where  ∞ =  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][44][45].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() =   () for some real values of ) with some known physical significance.For example, by the substitution

Solutions of the Dipole Equation
which is a 1-dimensional Schrödinger equation with an effective potential Here the last term represents the traditional centrifugal term plus an additional constant 2. This potential is repulsive for large values of , and it becomes attractive in a short range close to the origin because of the coupling with the exp (−/) +  0 exp (/) part of the wave function, through .
Another possible approach is provided by the substitution Ã() = (), and we obtain the radial part of a Schrödinger equation in spherical coordinates 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  ℎ (+∞) = 0, and  ℎ (0) → ±∞ depending on the value of .In the simpler case  = 0, for example, if  ≥ 4 this potential has a valley centered at This valley creates a barrier of height and relative depth with respect to the top of the valley From the aspect of this potential it results that for weak  −  coupling ( ∼ 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  = 5,  = 10.
For the substitution Q =  3/2 () we obtain the radial part of the Schrödinger equation in cylindrical coordinates Eq. ( 30) is a Fuchsian differential equation with singularities at the point at infinity, and order six singularity at  = 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  is very small, we will further reduce this equation to a fourthorder singularity Fuchsian equation.
As we mentioned in Section 3, Eqs.( 21), (22), and (30), we introduce the second hypothesis (in addition to  ∽ ) by considering a weak coupling between the  and the  coordinates in the order parameter dependence.This constraint is equivalent to  ∼ 0, and it arises somehow naturally since the order parameter Ψ(, , ) has its dominant dependence on the  coordinate.This happens because the |Ψ| 2 iso-surfaces are basically laying along the dipole surface coordinates  =constant.Under this hypothesis, Eq. ( 22) becomes and we will refer further to Eq. ( 31) as the dipole 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/ 4 , 1/ 2 , , in the potential energy.The resulting equation is mapped into the Mathieu equation for () through the transformation  = (− 2 /) 1/4 exp ().It is interesting that even holding all the terms in the powers of 1/ 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 in 1/ has a deep meaning.Although not obvious, Eq. ( 31) has an intrinsic symmetry that allows mapping it into a very symmetrical form, called the symmetrical-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.

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  coordinate, that is, for  = 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  →  = /.The resulting equation is a general double confluent Heun equation (gDCHE [27]) in () of the form where the operator  = (/) and   ,   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][48][49][50][51][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 2 and  2 = , where  = / 2 .By using a transformation of the dependent variable of the form any 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 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:  = ( √ / 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 (z) where D = z(/z), and 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  ∞ 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  and .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 at  ∞ , 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  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  ∞ .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 and 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]: where   is the holomorphic local Heun function as a solution of the DCHE (35) given by the series in the range on arg(z) ∈ (−/2, 5/2), with the coefficients   =   (,  1 ,  −1 , ) uniquely given by the three-term recursion relation with  ∈N,  ± −1 = 0, ± 0 = 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  2 is generated from the first one by application of the group operator  1 on  1 , 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 , , , .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 energies  < 0 the solutions  1 (z), that is,  1 (), 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 = − and only for odd values for .For each such odd  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/ 3 .For  2 , however, the majority of solutions have  > 0, so they cannot be reduced to quasipolynomials.
In the following, we present two examples of quasitruncated solutions  1 () obtained from the abstract solutions  1 (z): with  = −9 2 /4, and with  = − 2  2 1 /4.In order to have a three-dimensional visual on the behavior of these solutions we plot in

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  → ∞, 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 for () and  from Eqs. ( 22) and ( 23) we substitute  =  ∞ , () =  ∞ ( ∞ ), where  ∞ is the upper bound of the  variable.The resulting equation is a regular Sturm-Liouville problem.The self-adjoint form of the corresponding linearized dipole equation (22) is The linear solutions |Ψ| 2 = 0.5, 0.6, 0.75, and 0.9 isoplot surfaces for  = 0.5,  ∞ = 10,  = 1, calculated with the quasipolynomial solutions Eq. ( 42) for  = −5.3.The surfaces are torilike, as expected from the behavior of  =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. with In this () representation we have  ∈ (0, 1], and () > 0, () > 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  of Eq. ( 44) are real and form an infinite sequence  1 <  2 < ⋅ ⋅ ⋅ <   < . . .with lim →∞   = ∞.According to the physical restriction  < 1, it results that the spectrum of interest of the linear dipole boundary problem is also finite.All eigenvalues   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   ̸ =   we have Normalized with respect to the weight (), the eigenfunctions form a complete orthonormal basis of piecewise continuous functions, with piecewise square integrable continuous derivative in  ∈ [0,  ∞ ],  ≥ 0.
In order to solve explicitly the eigenvalue problem for  we apply the boundary condition at  ∞ , Eq. ( 23).Since the physical range of interest consists in large values for  ∞ (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 where we denoted The boundary condition Eq. ( 47) contains two unknowns,  1,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  1,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 where we denote On the other hand, we calculate the asymptotic behavior towards infinity of the two linear independent analytic Heun solutions  1,2 () from  1,2 (z), Eqs.(39a) and (39b).It results in an asymptotic term of the form to be matched with Eq. (49).The solution of this matching condition results in where  , has the same meaning as above, Eq. ( 50).With the coefficients  1,2 obtained in Eqs.(52a) and (52b), with a choice for  ∞ , ,  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  =   ( ∞,, ), where  is the "radial" quantum number.
The corresponding eigenfunctions are  1 (;  ∞ , , ,   ( ∞ , , )).The resulting solutions of the linearized GL equation are We present some examples of such solutions in Figure 5 for  ∞ = 10 and two values of .For each situation we obtain a discrete spectrum for , as predicted at the beginning of this section.In the investigated range we obtained between two and seven positive eigenvalues for .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  increases.
The final form of this equation is For any set of values for  ∞ ,  and  the solutions of Eq. ( 54) provide the  spectrum.For every  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  < 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 Figures 6 and 7.In Figure 8 we present the eigenfunctions and spectra for small values of energy in the range 0 <  < 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 <  < 1.The maximum number of spectral levels   depends on  and it can be calculated from the condition |Ψ| 2 < 1.The next spectral domain is within the negative energies and it also depends on all the parameters.In the normal range of values ( = 1, . . ., 6,  = 0.2 − 10) this condition offers a wide range of energies (it can go down to  ∼ −75) so this is the main range for negative energies.In this domain, because the argument of the tan function in Eq. ( 54)  47), or Eq. ( 54), for small dipole strength,  = 0.6.We present two cases of volume radii,  =  ∞ = 10, 60.The coupling of energy with angular momentum is reduced, so we notice a weak dependence of .For larger volumes the spectrum is richer, note the difference in scale of left and right frames.becomes imaginary, this tangent transforms into hyperbolic tangent.Consequently, for each value of  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.
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 Figure 9 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 ( ∞ = 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.
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 with integer values of  such that  < 1.The number of eigenvalues for positive  in this approximation can be estimated and we have such bounded states only if the following criterion is fulfilled: This lower bound for  ∞ (or equivalently, an upper bound for ) gives a criterium for the critical size, or dipole strength that  can trigger the generation of multivortex structures creation.For example, if we choose  ∞ = 20 we have 5 eigenstates for  = 0, 1 and 4 eigenstates for  = 2, 3, 4.
For larger radii of the volume the energy spectrum obtained as solution of Eq. ( 54) can be approximated with with  integer ("radial" quantum number) labeling the energy spectrum for fixed .
The energy spectra presented so far depend only on two quantum numbers: , .Even if this is a full threedimensional 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  =const.describe the vortices surfaces.By this hypothesis we can neglect the  coordinate dependence in the order parameter, which eliminates the third quantum numbers, basically like taking into account just ground states with respect to  variable excitations.

Nonlinear Vortex Patterns
The goal of this paper is to obtain the nonlinear order parameter as the solution Ψ(  →  ) 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 Ψ  , 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 (, , ), 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 volume  ∞ ∼ .In the approximation of weak coupling with the  orthogonal variable (i.e., () =constant) described above we introduce the following notation for the linear basis: with  , 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  , [7,[9][10][11][12][13][14][15][16][17][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  V ∈ (0,  ∞ ) of the dipole coordinate such that the equation |Ψ  ( V , )| = 0 has a number of distinct solutions for  ∈ [0, 2).If we denote by   ,  = 1, 2, . . .,  these solutions, by analytic prolongation theorem we can always find a neighborhood of  V on which |Ψ  | is arbitrary small for each of the roots   .The surfaces described by the values of  in these neighborhoods are enveloping surfaces of the multivortex states.In this case the value of   is called the vorticity of the state.The center-lines of the vortices are the curves obtained by the intersections of the  =  V and  =   surfaces.
In literature [15,16] there are examples of such multivortex states obtained with as little as   (  ) = 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 (  = 2,  ∈ { with  = 1, 2, . . ., | 1 −  2 | and all functions  are evaluated at  V .To visualize an example of multivortex state, we choose  ∞ = 100, and  1 = 0, 1 = 0,  2 = 0, . . ., 4. We plot in Figure 10 four isoplot surfaces with |Ψ  | = 0.9 for different vorticity numbers.Nevertheless, the topology of the multivortex surfaces is controlled mainly by the linear

B. Solutions of the Dipole Equation as Expansion in Confluent Hypergeometric Functions
The Heun functions solutions represented as series are analytic in the complex plane within some sectors, which brings some constraints on the solutions.Also, being power series they are not efficiently fast convergent.An alternative to this problem is to expand these solutions in terms of confluent hypergeometric functions [27].The advantage is that the new series provide a full analytic behavior.They are convergent absolutely and locally uniformly all over the complex plane.Also, the solutions built with these hypergeometric functions are of Floquet type; that is, they represent a series expansion around any point in the complex plane, not only about the singularities at 0 and ∞.We can write these series in the form  ( ) ,, l = −2 (   ) where ,  ± and  are the reduced parameters Eq. ( 36) for the DCHE associated with the linearized dipole equation (35), and ,  ∈ {1, 1} independently.

Figure 3 :
Figure 3: The effective potential  ℎ () in the radial part of the Schrödinger equation in spherical coordinates, Eq. (26).Solid lines represent  = 0 potentials, and dotted lines represent  = 5 potentials, for the same , .
Figure 4 few surfaces of equation |Ψ(, )| = |Ψ(/sin 2 , )| =  1 (/sin 2 )  const.The three surfaces presented in Figure 4 follow the shape of the surfaces of coordinates  =const., proving the hypothesis according to which the  dependence of the order parameter can be neglected.

Figure 6 :
Figure 6: Energy spectrum  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,  =  ∞ = 10, 60.The coupling of energy with angular momentum is reduced, so we notice a weak dependence of .For larger volumes the spectrum is richer, note the difference in scale of left and right frames.

Figure 8 :Figure 9 :
Figure8: Eigenvalues (0 <  < 1) and eigenfunctions for the linearized dipole equation for  = 0.5 and  = 0 (upper frames) and  = 3 (lower frames) with boundary condition at  ∞ = 10.In the left frames we zoom in all the eigenfunctions for a given  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  ∞ , ,  values.