# Vortex Structures inside Spherical Mesoscopic Superconductor Plus Magnetic Dipole.

1. IntroductionThe 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 volume V. Outside of this volume the Cooper pair density [mathematical expression not reproducible] 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]

[mathematical expression not reproducible], (1)

where m 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 [alpha](T) < 0 and the nonlinear term coupling constant [beta] > 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 [PSI] by minimizing the free energy functional Eq. (1) for constant volume V consists in expanding the order parameter [PSI] 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 [PSI] 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 [PSI] is very small in this range, and higher order terms of [PSI] 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.

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 [mathematical expression not reproducible] produced by a infinitesimal dipole placed at its center [mathematical expression not reproducible]. The Euler-Lagrange equations for the free energy functional Eq. (1), written in the absence of the last term in [mathematical expression not reproducible] as explained above, are known as the GL equation [9-12]

[mathematical expression not reproducible] (2)

[mathematical expression not reproducible] (3)

[mathematical expression not reproducible], (4)

where the electromagnetic momentum [??] is given by

[mathematical expression not reproducible] (5)

where [??] is the vector potential, q 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)

[mathematical expression not reproducible]. (6)

We rewrite Eq. (6) in the London electromagnetic gauge, [nabla] x [??] = 0, and we rescale the variables [mathematical expression not reproducible]. We follow the GL phenomenological model from [9-18] and use the temperature dependence

[alpha] = [[alpha].sub.0] (1 - T/[T.sub.c0]),

[xi](T) = [??]/[square root of -2m[alpha](T)], (7)

where [T.sub.c0] is the critical temperature in absence of magnetic field and [xi](T) is the coherence length. In the new variables the LGL Eq. (6) can be written in the form of a linear operator equation

[mathematical expression not reproducible]. (8)

Since the dipole magnetic field has axial symmetry we can use any axially symmetric orthogonal coordinates in azimuth angle [phi] in the form [??] [right arrow] ([phi],a,b). In this case the [phi] dependence in Eq. (8) can be separated, and [??] has exp (iL[phi]) as eigenfunctions labeled by the angular momentum L. This decomposes the space of solutions [mathematical expression not reproducible] into subspaces parameterized by L

[mathematical expression not reproducible], (9)

with [PHI] [member of] [L.sup.2]([R.sup.2]) square integrable function. The decomposition in Eq. (9) reduces Eq. (6) to a two-dimensional PDE. On each of these L subspaces the operator [??] has a particular form depending on L denoted by [[??].sub.L], and for each such restricted operator we can attach an eigenvalue problem

[[??].sub.L][[PHI].sub.E,[kappa]] = E[[PHI].sub.E,[kappa]], (10)

where [kappa] is a provisional degeneracy label, which later on becomes the radial quantum number. Of course, E and [[PHI].sub.E,[lambda]] 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

[mathematical expression not reproducible], (11)

where [[??].sub.[phi]] = (- sin [phi], cos [phi], 0) is the unit vector in the [phi] 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 [summation] to be zero (the Saint-James-de Gennes conditions) [9-12,15-18]

[mathematical expression not reproducible], (12)

where [??] represents the outward normal to the spherical surface [summation] defined by [absolute value of ([??])] = 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 [R.sup.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.

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 Schrodinger 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 [theta] = 1,-1,0 in a series of Jacobi polynomials with "dipole polynomial" coefficients.

In all these situations separation and integration of the Schrodinger 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 Schrodinger 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, [theta] 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 Schrodinger 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/[r.sup.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-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 Schroodinger equation in the form

[mathematical expression not reproducible], (13)

where we denoted [chi] = 2[pi]q[mu]/([??]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

[PSI](r,[theta][phi]) = [PHI](r,[theta])[e.sup.iL[phi]], (14)

with L being an arbitrary positive integer. Eq. (13) reduces to a 2-dimensional Schrodinger equation with an effective potential in the form

[mathematical expression not reproducible]. (15)

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.sup.-4]. Depending on [theta] 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 [chi] controls the barrier height. For [theta] ~ [pi]/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 [theta] ~ 0,[pi] 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.

The Schroodinger equation with the potential in Eq. (15) cannot be resolved by traditional expansion in spherical harmonics with coefficients given by the monomials [r.sup.l]. Since this equation contains nonhomogenous terms as powers of r, even if we use Laurent polynomials instead of these monomials, the coefficients of [Y.sub.lL] 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 [f.sub.l](r) instead of constant coefficients for each [Y.sub.lL], the differential equations resulting for each [f.sub.l] will contain infinite many terms. This happens because the term [sin.sup.2] [theta] multiplied with spherical harmonics of different orders decomposes (through the Wigner 3 j-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, [phi]). These dipole coordinates are related to the spherical ones by the relations

[mathematical expression not reproducible]. (16)

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 [PSI]([??]) 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 non- linear 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,

[mathematical expression not reproducible], (17)

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 [absolute value of ([psi])] =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

[mathematical expression not reproducible]. (18)

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 [PHI](r, [theta]) [right arrow] [PHI](a, b):

[mathematical expression not reproducible], (19)

where L [member or] 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) [right arrow] (r,[theta]) 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 [r.sup.2]/b << 1 from Eq. (A.3). This approximation describes the order parameter in neighborhoods of [theta] = 0, and [theta] = [pi], 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 [epsilon] = [a.sup.2]/b, according to Eq. (A.7) in Appendix A. For example, if we just choose [epsilon] = 0.2 we have a relative error of only 4% or less for a ~ r within a range of the azimuth angle [theta] [member of] [0,80[degrees]] [union] [100[degrees], 180[degrees]], 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 approximation a ~ r Eq. (19) reduces to the form

[mathematical expression not reproducible], (20)

and we can factorize further the b-dependence from the a-dependence by the substitution

[PHI](a,b) = Q(a)([e.sup.-c/b] + [C.sub.0][e.sup.c/b]), (21)

where c, [C.sub.0] 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]

[mathematical expression not reproducible], (22)

This equation is associated with the boundary conditions Eq. (12)

[mathematical expression not reproducible], (23)

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.sub.[infinity]] = R. At this point we describe the mesoscopic volume taken for the present study. The equation r = R [sin.sup.2] [theta] 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 ([theta] [approximately equal to] 0[degrees]) and south ([theta] [approximately equal to] 180[degrees]) poles the volume is depleted (r [approximately equal to] 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, [infinity], [42], and consequently there are no convergent solutions in terms of Frobenius series. We need to solve this equation in the interval a [member of] (0, [a.sub.[infinity]]], where [a.sub.[infinity]] = 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) [right arrow] [??](a) = [a.sup.p]Q(a) for some real values of p) with some known physical significance. For example, by the substitution [??](a) = [a.sup.2]Q(a) we have

[d.sup.2][??]/d[a.sup.2] + [V.sub.1d](a)[??] + E[??] = 0, (24)

which is a 1-dimensional Schrodinger equation with an effective potential

[V.sub.1d] = -[c.sup.2]/[a.sup.6] + [[chi].sup.2]/[a.sup.4] - 2[L.sub.[chi]]/[a.sup.3] + [L.sup.2] + 2/[a.sup.2]. (25)

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) + [C.sub.0] exp (c/b) part of the wave function, through c.

Another possible approach is provided by the substitution [??](a) = aA(a), and we obtain the radial part of a Schrodinger equation in spherical coordinates

[mathematical expression not reproducible]. (26)

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 [V.sub.sph](+[infinity]) = 0, and [V.sub.sph](0) [right arrow] [+ or -] [infinity] depending on the value of c. In the simpler case c = 0, for example, if L [greater than or equal to] 4 this potential has a valley centered at

[a.sub.1] = [chi] 3L - [square root of [L.sup.2] - 16/2([L.sup.2] + 2). (27)

This valley creates a barrier of height

[mathematical expression not reproducible], (28)

and relative depth with respect to the top of the valley

[V.sub.0] = L [([L.sup.2] - 16).sup.3/2]/16[chi square]. (29)

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, [chi] = 10.

For the substitution [??] = [a.sup.3/2]Q(a) we obtain the radial part of the Schroodinger equation in cylindrical coordinates

[mathematical expression not reproducible]. (30)

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 Section 3, 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 [PSI](a, b, [phi]) has its dominant dependence on the a coordinate. This happens because the [[absolute value of ([PSI]].sup.2] iso-surfaces are basically laying along the dipole surface coordinates a =constant. Under this hypothesis, Eq. (22) becomes

[mathematical expression not reproducible], (31)

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/[a.sup.4], 1/[a.sup.2], E, in the potential energy. The resulting equation is mapped into the Mathieu equation for A(z) through the transformation a = [(-[chi square]/E).sup.1/4] exp (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 in 1/a 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.

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 [right arrow] z = a/[chi]. The resulting equation is a general double confluent Heun equation (gDCHE [27]) in Q(z) of the form

[mathematical expression not reproducible], (32)

where the operator D = z(d/dz) and [p.sub.i],[q.sub.j] 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.sub.-1] = [p.sub.1] = [q.sub.1] = 0, [p.sub.0] = 3, [q.sub.-2] = -1, [q.sub.-1] = 2L,[q.sub.0] = -[L.sup.2] and [q.sub.2] = [lambda], where [lambda] = E/[chi square]. By using a transformation of the dependent variable of the form

[mathematical expression not reproducible] (33)

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, [[xi].sub.0] = -3/2, [[xi].sub.1] = [+ or -] [square root of [lambda]], [[xi].sub.-1] = [+ or -] 1. It is worth mentioning that one obtains a nondegenerate (a "good") canonical DCHE form if one has the condition

([q.sub.2] - [p.sup.2.sub.1]/2)([q.sub.-2] - [p.sup.2].sub.-1]/2) [not equal to] 0, (34)

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 = ([square root of i]/[[lambda].sup.1/4])[??], 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([??])

[mathematical expression not reproducible], (35)

where [mathematical expression not reproducible], and

[mathematical expression not reproducible], (36)

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 [alpha], [[beta].sub. [+ or -] 1] determine the leading coefficients of the asymptotical series of the solutions at [infinity] and 0, respectively. However, in the case of the dipole equation we have [[beta].sub.1] = 0, and only the reduced energy ([lambda]) controls the asymptotical behavior of the solution at infinity, through [alpha]. This behavior is in agreement with the fact that the boundary condition at [a.sub.[infinity]] is the only responsible constraint for energy quantization in prescribed levels. The singular parameters [alpha], [[beta].sub.-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 [gamma] "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 [infinity] just reduces to permuting the [beta] 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, [??] = 0, [infinity]. 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.sup.*] =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.sup.*], though because [C.sup.*] is not simply connected the Poincare 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.sup.*], bounded by the corresponding Stokes rays, depending on the argument of a[??] [45].

Since the physical boundary conditions at a[infinity], Eq. (12), are important for the eigenvalue problem for energy we will construct the solutions of Eq. (35) based on the asymptotic behavior for [??] [right arrow] [infinity]. 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.sub.[infinity]]. Then, this solution can be continued analytically towards the solution around [??] = 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 trans- formations 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

[mathematical expression not reproducible], (37)

and

[T.sub.0][y] = exp(- [alpha]z + [alpha]/z) y ([alpha], [[beta].sub.1],[[beta].sub.-1],[gamma];1/2). (38)

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]:

[y.sub.1]([??]) = [H.sub.l]([alpha],[[beta].sub.1], [gamma];[??]), (39a)

[mathematical expression not reproducible], (39b)

where [H.sub.l] is the holomorphic local Heun function as a solution of the DCHE (35) given by the series

[mathematical expression not reproducible], (40)

in the range on arg([alpha][??]) [member of] (-[pi]/2, 5[pi]/2), with the coefficients [C.sub.n] = [C.sub.n]([alpha], [[beta].sub.1], [[beta].sub.-1],[gamma]) uniquely given by the three-term recursion relation

[mathematical expression not reproducible], (41)

with n [member of] N, [C.sup.[+ or -].sub.-1] = 0, [C.sup.[+ or -].sub.0] = 1. The signs [+ or -] represent the choice for one of the two solutions for the reduced parameters [alpha], [[beta].sub.-1] in Eq. (36). We mention that the second solution [y.sub.2] is generated from the first one by application of the group operator [T.sub.1] on [y.sub.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 a, [chi], 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 energies E < 0 the solutions [y.sub.1]([??]), that is, [Q.sub.1](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 [[beta].sub.-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[chi]L/[a.sup.3]. For [Q.sub.2], however, the majority of solutions have E > 0, so they cannot be reduced to quasipolynomials.

In the following, we present two examples of quasitruncated solutions [Q.sub.1](a) obtained from the abstract solutions [y.sub.1]([??]):

[mathematical expression not reproducible], (42)

with E = -9 [chi square]/4, and

[mathematical expression not reproducible], (43)

with E = -[chi square][w.sup.2.sub.1]/4. In order to have a three-dimensional visual on the behavior of these solutions we plot in Figure 4 few surfaces of equation [absolute value of [PSI](a,[phi])] = [absolute value of [PSI](r/[sin.sup.2][theta],[phi])] = [Q.sub.1](r/[sin.sup.2][theta])[e.sup.iL[phi]] 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.

5. Eigenvalue Problem for the Dipole Equation

Typical study of the quantization of the energy spectrum of Heun's equation involves the y 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 [right arrow] [infinity], but not faster than an exponential. Consequently, the [45] boundary conditions quantize the accessory parameter [gamma], and the recursion relation for the coefficients of the series becomes a four-term Poincare-Perron type, which is more difficult to handle than our three-term recursion relations.

In order to solve our eigenvalue problem for Q(a) and E from Eqs. (22) and (23) we substitute a = [a.sub.[infinity]]x, y(x) = [a.sub.[infinity]]xQ([a.sub.[infinity]]x), where [a.sub.[infinity]] 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

L [y] = Er (x) y, (44)

with

[mathematical expression not reproducible], (45a)

[mathematical expression not reproducible]. (45b)

In this y(x) representation we have x [member of] (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 [E.sub.1] < [E.sub.2] < ... < [E.sub.n] < ... with [lim.sub.n[right arrow][infinity]] [E.sub.n] = [infinity]. 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 [E.sub.n] 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 [E.sub.k] [not equal to] [E.sub.n] we have

[mathematical expression not reproducible]. (46)

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 [member of] [0,[a.sub.[infinity]]], b [greater than or equal to] 0.

In order to solve explicitly the eigenvalue problem for E we apply the boundary condition at [a.sub.[infinity]], Eq. (23). Since the physical range of interest consists in large values for [a.sub.[infinity]] (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

[mathematical expression not reproducible], (47)

where we denoted

[mathematical expression not reproducible]. (48)

The boundary condition Eq. (47) contains two unknowns, [N.sub.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 [N.sub.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

[mathematical expression not reproducible], (49)

where we denote

[mathematical expression not reproducible]. (50)

On the other hand, we calculate the asymptotic behavior towards infinity of the two linear independent analytic Heun solutions [Q.sub.1,2](a) from [y.sub.1,2]([??]), Eqs. (39a) and (39b). It results in an asymptotic term of the form

[mathematical expression not reproducible], (51)

to be matched with Eq. (49). The solution of this matching condition results in

[mathematical expression not reproducible], (52a)

[mathematical expression not reproducible], (52b)

where [[phi].sub.L,E] has the same meaning as above, Eq. (50).

With the coefficients [N.sub.1,2] obtained in Eqs. (52a) and (52b), with a choice for [a.sub.[infinity]], L, [chi] 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 = [E.sub.k]([a.sub.[infinity]L,[chi]]), where k is the "radial" quantum number.

The corresponding eigenfunctions are [Q.sub.1](a;[a.sub.[infinity]], [chi], [E.sub.k]([a.sub.[infinity]],L,[chi])). The resulting solutions of the linearized GL equation are

[[PSI].sub.k,L](a,[phi]) = [Q.sub.1](a;[a.sub.[infinity]],L,[chi],[E.sub.k]([a.sub.[infinity]],L, [chi]))[e.sup.iL[phi]]. (53)

We present some examples of such solutions in Figure 5 for [a.sub.[infinity]] = 10 and two values of [chi]. 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.

The final form of this equation is

[mathematical expression not reproducible]. (54)

For any set of values for [a.sub.[infinity]], [chi] 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 Figures 6 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 [k.sub.max] depends on L and it can be calculated from the condition [[absolute value of ([PSI])].sup.2] < 1. The next spectral domain is within the negative energies

-[chi]/2 (9/4 + [L.sup.2]) < E < 0, (55)

and it also depends on all the parameters. In the normal range of values (L = 1, ..., 6, [chi] = 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.

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 ([a.sub.[infinity]] = 10) and an average value of the relative dipole strength [chi] = 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 [chi], we can approximate the solutions of the spectral equation Eq. (54) with the analytical expression

[mathematical expression not reproducible], (56)

with integer values of k such that E < 1. The number of eigenvalues for positive E in this approximation can be estimated

[mathematical expression not reproducible], (57)

and we have such bounded states only if the following criterion is fulfilled:

[a.sub.[infinity]] > [chi square] [pi] [square root of 9/4 + [L.sup.2]]. (58)

This lower bound for [a.sub.[infinity]] (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.sub.[infinity]] = 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

[mathematical expression not reproducible], (59)

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 [nabla][PSI] 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 [PSI]([??]) 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 [[PSI].sup.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,[phi]), 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 [a.sub.[infinity]] ~ 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:

[[PSI].sub.L,E] (a, b, [phi]) [right arrow] [[PSI].sub.L,k] (a, [phi]) = [Q.sub.1] (L, k; a)[e.sup.iL[phi]], (60)

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

[mathematical expression not reproducible], (61)

with [C.sub.L,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 [C.sub.L,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 [a.sub.v] [member of] (0, [a.sub.[infinity]]) of the dipole coordinate such that the equation [absolute value of ([[PSI].sup.NL]([a.sub.v],[phi])] = 0 has a number of distinct solutions for [phi] [member of] [0,2[pi]). If we denote by [[phi].sub.j], j = 1,2, ..., n these solutions, by analytic prolongation theorem we can always find a neighborhood of [a.sub.v] on which [absolute value of ([[PSI].sup.NL])] is arbitrary small for each of the roots [[phi].sub.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 [L.sub.max] is called the vorticity of the state. The center-lines of the vortices are the curves obtained by the intersections of the a = [a.sub.v] and [phi] = [[phi].sub.j] surfaces.

In literature [15, 16] there are examples of such multivortex states obtained with as little as [k.sub.max]([L.sub.max]) = 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 ([k.sub.max] = 2, L [member of] {[L.sub.1], [L.sub.2]}) the expression

[mathematical expression not reproducible] (62)

with j = 1,2, ..., [absolute value of ([L.sub.1] - [L.sub.2])] and all functions Q are evaluated at [a.sub.v]. To visualize an example of multivortex state, we choose [a.sub.[infinity]] = 100, and [L.sub.1] = 0,[k.sub.1] = 0, [L.sub.2] = 0, ..., 4. We plot in Figure 10 four isoplot surfaces with [absolute value of ([[PSI].sup.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.

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 Schroodinger 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 [PSI]([??]) 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,[phi]), under the boundary conditions Eq. (12). As a first step, we use cylindrical sym- metric coordinates. The linear solution can be separated in [PSI] = [PHI](a,b)[e.sup.iL[phi]], 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 [PHI](a,b) = Q(a)([e.sup.-c/b] + [C.sub.0][e.sup.c/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 [mathematical expression not reproducible] 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 [PSI] is mainly directed radial and orthogonal to the magnetic dipole field lines, it results that [PSI] 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].

Appendix

A. Dipole Coordinate System

The dipole coordinates defined in Eq. (16) have as domain of definition a [member of] [0, [infinity]), b [member of] (-[infinity], [infinity]), b [not equal to] 0 for r [member of] (0, [infinity]), [theta] [member of] (0,[pi]), [theta] [not equal to] [pi]/2. The Jacobian of the transformation from spherical coordinates and the element of volume are given by

J = [r.sup.2] 4 + [tan.sup.2] [theta]/[sin.sup.3][theta],

dV = [a.sup.4]/[b.sup.2] dadbd[phi] (A.1)

The dipole coordinates are actually different from other apparently similar coordinates like bipolar, bispherical, or toroidal coordinates [54-56]. In the notation used above ([x.sup.i]) = (x,y,z) [right arrow] (a, b, [phi]), these traditional orthogonal curvilinear systems of coordinates have symmetrical forms for the three Euclidean coordinates of the general form

[mathematical expression not reproducible], (A.2)

where the functions f are trigonometric or hyperbolic functions.

Consequently, the two parameters a, b are in general an angle and a distance. In the case of dipole coordinates, the parameters a, b have different dimensions, but both are related to distances, [a] = m, [b] = [m.sup.2]. The symmetry mentioned above for these traditional coordinates induces a certain symmetry in their coordinates curves and surfaces. All these curves represent circles either secant or not intersecting. On the contrary, in the case of dipole coordinates the coordinate curves are not circles, but either glued deformed circles (for a = const.) or tangent closed curves with a common singular point (for b = const). The coordinate surfaces, Figure 11, are orthogonal curvilinear coordinates. In order to obtain the inverse transformation, back to spherical coordinates, we have to solve the equation

r/a + [r.sup.4]/[b.sup.2] = 1. (A.3)

This equation is a particular case of a depressed quartic equation and it can be solved by the Ferrari method, hence reducing it to a depressed cubic equation, and then use Cardano's formulas. Since all coefficients of Eq. (A.3) are positive, this equation always has two complex conjugate solutions, a real negative solution and a nonnegative solution which represents the correct geometric one, as being the value of the magnitude of the position vector r. This solution reads

[mathematical expression not reproducible], (A.4)

where we define

f ([epsilon]) = [(9[[epsilon].sup.2] + [square root of 768[[epsilon].sup.6] + 81[[epsilon].sup.4]).sup.1/3], (A.5)

and

[epsilon] = [a.sup.2]/b. (A.6)

We expand the above expression r([epsilon], b) in Taylor series with respect to e around zero, and we obtain

[mathematical expression not reproducible], (A.7)

where obviously the dominant term is r ~ [square root of b[epsilon]] = a. This approximation is valid around the North and South poles of the glued, vertically aligned, deformed circles of equation a = const., Figure 11. The solutions of a similar type of quartic equation in d,

[a.sup.2] [cos.sup.4][theta] - 2[a.sup.2] [cos.sup.2][theta] - b cos [theta] + [a.sup.2] = 0, (A.8)

can provide the inverse [theta] = [theta](a,b).

A more compact expression can be obtained if we solve the inverse transformation Eq. (A.3) in cartesian coordinates, namely, x = r sin [phi], z = r cos [phi] as functions of (a,b), where [phi] = [pi]/2 - [theta]. Beginning with relation Eq. (A.3) and the above definitions for x, z, with [r.sup.2] = [x.sup.2] + [z.sup.2] we obtain

[mathematical expression not reproducible]. (A.9)

After some elementary algebra we obtain the inverse in the form

[mathematical expression not reproducible], (A.10)

which can be compressed even more into

[mathematical expression not reproducible], (A.11)

where

W = [a.sup.1/3] [b.sup.1/2]. (A.12)

The equations above prove the fact that dipole coordinates have the dimension unbalance between coordinates, as opposed to the other bipolar, bispherical, or toroidal coordinates.

The Lamme coefficients of the dipole coordinate system are

[mathematical expression not reproducible]. (A.13)

Consequently, the Laplacian in (a,b) coordinates has the form

[mathematical expression not reproducible]. (A.14)

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 [infinity]. We can write these series in the form

[mathematical expression not reproducible], (B.1)

where [PHI](a; b; z) is the Kummer function, the indices j, l, [??] [member of] {1,2} independently (which generates 8 solutions), and [mu] is the summation label which is an integer plus an arbitrary complex parameter v called the Floquet exponent of the DCHE. This exponent (is unique modulo Z) should be conveniently chosen according to the point about which the series expands. Finally, the series coefficients [[eta].sup.j.l.sub.[mu]] are obtained from the three-term recursion relation

[mathematical expression not reproducible], (B.2)

where [alpha], [[beta].sub.[+ or -]] and y are the reduced parameters Eq. (36) for the DCHE associated with the linearized dipole equation (35), and j,l member of] {1,1} independently.

https://doi.org/10.1155/2018/8652151

Data Availability

The analytic and numerical calculations and plots data used to support the findings of this study are available from the corresponding author upon request.

Disclosure

An earlier version of this paper was presented at the 10th International Conference on Geometry, Integrability and Quantization 2010, Varna, Bulgaria; see [39].

Conflicts of Interest

The author declares that there are no conflicts of interest.

Acknowledgments

This work was inspired by several discussions with M. Doria, visiting Professor at Condensed Matter Theory Group at University of Antwerp. His ideas and comments are acknowledged as the motivation of this paper. The author also acknowledges financial support from the Condensed Matter Theory Group at University of Antwerp and he is grateful to F. M. Peeters for support through sabbatical year.

References

[1] G. M. Whitesides, "Nanoscience, nanotechnology, and chemistry," Small, vol. 1, no. 2, pp. 172-179, 2005.

[2] S. Kruchinin, H. Nagao, and S. Aono, Modern Aspects of Superconductivity, World Scientific, New Jersey, USA, 2010.

[3] M. Safronova, D. Budker, D. DeMille, D. F. Kimball, A. Derevianko, and C. W. Clark, "Search for new physics with atoms and molecules," Reviews of Modern Physics, vol. 90, no. 2, 2018.

[4] D. Bimberg and U. W. Pohl, "Quantum dots: Promises and accomplishments," Materials Today, vol. 14, no. 9, pp. 388-397, 2011.

[5] H. Kupfer, G. Linker, G. Ravikumar et al., "Correlation of the vortex order-disorder transition with the symmetry of the crystal lattice in V3Si," Physical Review B: Condensed Matter and Materials Physics, vol. 67, no. 6, 2003.

[6] E. T. Filby, A. A. Zhukov, P. A. de Groot, M. A. Ghanem, P. N. Bartlett, and V. V. Metlushko, "Shape induced anomalies in vortex pinning and dynamics of superconducting antidot arrays with spherical cavities," Applied Physics Letters, vol. 89, no. 9, p. 092503, 2006.

[7] A. K. Geim, I. V. Grigorieva, S. V. Dubonos et al., "Phase transitions in individual sub-micrometre superconductors," Letter to Nature, vol. 390, no. 6657, pp. 259-262, 1997.

[8] V. R. Misko, "Recent Advances in Superconductivity and Vortex Matter: Selected Topics," Reports in Advances of Physical Sciences, vol. 1, no. 4, p. 1750009, 2017.

[9] B. Xu, M. V. Milosevic, and F. M. Peeters, "Magnetic properties of vortex states in spherical superconductors," Physical Review B: Condensed Matter and Materials Physics, vol. 77, no. 14, 2008.

[10] M. M. Doria, A. R. Romaguera, and F. M. Peeters, "Effect of the boundary condition on the vortex patterns in mesoscopic three-dimensional superconductors: Disk and sphere," Physical Review B, vol. 75, Article ID 064505, 2007.

[11] M. M. Doria, A. R. Romaguera, M. V. Milosevic, and F. M. Peeters, "Threefold onset of vortex loops in superconductors with a magnetic core," EPL (Europhysics Letters), vol. 79, p. 47006, 2007.

[12] B. J. Baelus, D. Sun, and F. M. Peeters, "Vortex structures in mesoscopic superconducting spheres," Physical Review B: Condensed Matter and Materials Physics, vol. 75, no. 17, 2007.

[13] G. Zha, S. Zhou, B. Zhu, Y. Shi, and H. Zhao, "Superconducting phase transitions in thin mesoscopic rings with enhanced surface superconductivity," Physical Review B: Condensed Matter and Materials Physics, vol. 74, no. 2, 2006.

[14] G. Zha, S. Zhou, B. Zhu, Y. Shi, and H. Zhao, "Charge distribution in thin mesoscopic superconducting rings with enhanced surface superconductivity," Physical Review B: Condensed Matter and Materials Physics, vol. 73, 2006.

[15] V. A. Schweigert and F. M. Peeters, "Phase transitions in thin mesoscopic superconducting disks," Physical Review B: Condensed Matter and Materials Physics, vol. 57, no. 21, pp. 13817-13832, 1998.

[16] V. A. Schweigert, F. M. Peeters, and P. S. Deo, "Vortex Phase Diagram for Mesoscopic Superconducting Disks," Physical Review Letters, vol. 81, no. 13, pp. 2783-2786, 1998.

[17] J. J. Palacios, "Metastability and Paramagnetism in Superconducting Mesoscopic Disks," Physical Review Letters, vol. 84, no. 8, pp. 1796-1799, 2000.

[18] J. J. Palacios, "Vortex matter in superconducting mesoscopic disks: Structure, magnetization, and phase transitions," Physical Review B: Condensed Matter and Materials Physics, vol. 58, no. 10, pp. R5948-R5951, 1998.

[19] O. Sato, M. Kato, and J. Physics 871, "A variety of vortex state solutions of Ginzburg-Landau equation on superconducting mesoscopic plates," Journal of Physics: Conference Series (JPCS), Article ID 012029, pp. 1-4, 2017.

[20] H. L. Cycon, R. G. Froese, W. Kirsch, and B. Simon, Schrodinger Operators, Schrodinger Operators, Berlin, Heidelberg, 1987.

[21] K. Miller and B. Simon, "Quantum magnetic Hamiltonians with remarkable spectral properties," Physical Review Letters, vol. 44, no. 25, pp. 1706-1707, 1980.

[22] D. Picca, "The dipolar potential," Journal of Physics A: Mathematical and General, vol. 15, no. 9, pp. 2801-2808, 1982.

[23] G. Levai and B. W. Williams, "The generalized Coulomb problem: an exactly solvable model," Journal of Physics A: Mathematical and General, vol. 26, no. 13, pp. 3301-3306, 1993.

[24] A. D. Alhaidari and H. Bahlouli, "Publisher's Note: Electron in the Field of a Molecule with an Electric Dipole Moment," Physical Review Letters, vol. 100, no. 12, 2008.

[25] D. Zwillinger, Handbook of Differential Equations, Academic Press, San Diego, Calif, USA, 1989.

[26] P. Moon and D. E. Spencer, Field Theory Handbook, Springer, Berlin, Germany, 1971.

[27] A. Ronveaux, Heuns Differential Equations, University Press, Oxford, England, 1955.

[28] A. Ralko and T. T. Truong, "Heun functions and the energy spectrum of a charged particle on a sphere under a magnetic field and Coulomb force," Journal of Physics A: Mathematical and General, vol. 35, no. 45, pp. 9573-9583, 2002.

[29] A. Ralko and T. T. Truong, "Behaviour of three charged particles on a plane under perpendicular magnetic field," Journal of Physics A: Mathematical and General, vol. 35, no. 45, pp. 9671-9683, 2002.

[30] A. Debosscher, "Unification of one-dimensional Fokker-Planck equations beyond hypergeometrics: factorizer solution method and eigenvalue schemes," Physical Review E: Statistical, Nonlinear, and Soft Matter Physics, vol. 57, no. 1, pp. 252-275, 1998.

[31] A. Y. Kazakov, "The central two-point connection problem for the reduced confluent Heun equation," Journal of Physics A: Mathematical and General, vol. 39, no. 10, pp. 2339-2348, 2006.

[32] C. W. Bauer, Z. Ligeti, M. Luke, A. V. Manohar, and M. Trott, "Global analysis of inclusive B decays," Physical Review D: Particles, Fields, Gravitation and Cosmology, vol. 70, Article ID 094017, 2004.

[33] D. Bouaziz and M. Bawin, "Regularization of the singular inverse square potential in quantum mechanics with a minimal length," Physical Review A: Atomic, Molecular and Optical Physics, vol. 76, Article ID 032112, 2007.

[34] B. Gao, "Solutions of the Schrodinger equation for an attractive 1/[r.sup.6] potential," Physical Review A, vol. 58, p. 1728, 1998.

[35] A. A. Stahlhofen, "Susy, Gauss, Heun and physics: a magic square?" Journal of Physics A: Mathematical and General, vol. 37, no. 43, pp. 10129-10138, 2004.

[36] M. Swisdak, "Notes on the Dipole Coordinate System," https:// arxiv.org/pdf/physics/0606044.pdf.

[37] M. N. Fatkullin and Y. S. Sitnov, "Dipole coordinate system and some of its characteristics," Geomegnatism and Aeronomy, vol. 12, pp. 293-295, 1972.

[38] A. Kageyama, T. Sugiyama, and K. Watanabe, "A Note on the Dipole Coordinates," Computers & Geosciences, vol. 32, no. 3, pp. 265-269, 2006.

[39] A. Ludu, "Vortex patterns beyond hypergeometric," Journal of Geometry and Symmetry in Physics, vol. 26, pp. 85-103, 2012.

[40] A. Ludu, M. V. Milosevic, and F. M. Peeters, "Vortex states in axially symmetric superconductors in applied magnetic field," Mathematics and Computers in Simulation, vol. 82, no. 7, pp. 1258-1270, 2012.

[41] A. Ludu, J. Van Deun, M. V. Milosevic, A. Cuyt, and F. M. Peeters, "Analytic treatment of vortex states in cylindrical superconductors in applied axial magnetic field," Journal of Mathematical Physics, vol. 51, no. 8, 082903, 29 pages, 2010.

[42] E. Hille, Ordinary Differential Equations in The Complex Domain, Wiiley, NewYork, NY, USA, 1976.

[43] M. N. Hounkonnou, A. Ronveaux, and K. Sodoga, "Factorization of some confiuent Heun's differential equations," Applied Mathematics and Computation, vol. 189, pp. 816-820, 2007.

[44] W. Buhring, "The double confluent Heun equation: characteristic exponent and connection formulae," Methods and Applications of Analysis, vol. 1, no. 3, pp. 348-370, 1994.

[45] W. Lay, K. Bay, and S. Y. Slavyanov, "Asymptotic and numeric study of eigenvalues of the double confluent Heun equation," Journal of Physics A: Mathematical and General, vol. 31, no. 42, pp. 8521-8531, 1998.

[46] W. Magnus and S. Winkler, Hill's Equation, Inderscience Publishers, New York, NY, USA, 1966.

[47] K. Kuiken, "Heun's equation and the hypergeometric equation," SIAM Journal on Mathematical Analysis, vol. 10, no. 3, pp. 655-657, 1979.

[48] P. A. Clarkson and P. J. Olver, "Symmetry and the Chazy Equation," Journal of Differential Equations, vol. 124, no. 1, pp. 225-246, 1996.

[49] P. A. Becker, "Normalization integrals of orthogonal Heun functions," Journal of Mathematical Physics, vol. 38, no. 7, pp. 3692-3699, 1997.

[50] G. Valent, "Heun functions versus elliptic functions," Mathematical Physics, 2005.

[51] D. D. B. Figueiredo, "Ince's limits for confluent and double-confluent Heun equations," Mathematical Physics, 2005.

[52] P. Malits and J. Math, "On certain Heun functions, associated functions of discrete variables, and applications," International Journal of Mathematics and Mathematical Sciences, vol. 59, p. 3717, 2003.

[53] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, DC, USA, 1964, http://www.math.sfu.ca/.

[54] G. Arfken in, Mathematical Methods in Physics, Academic Press, Orlando, FL, USA, pp. 112-117, 1970.

[55] H. Bateman, "Spheroidal and bipolar coordinates," Duke Mathematical Journal, vol. 4, no. 1, pp. 39-50, 1938.

[56] E. H. Lockwood, A book of curves, Cambridge University Press, Cambridge, England, pp. 186-190, 1967.

A. Ludu (iD)

Department of Mathematics, Embry-Riddle Aeronautical University, Daytona Beach, FL 32114, USA

Correspondence should be addressed to A. Ludu; ludua@erau.edu

Received 30 July 2018; Accepted 11 October 2018; Published 1 November 2018

Academic Editor: Alexander Iomin

Caption: Figure 1: Equipotential contours of the effective potential V(x, z) at y = 0, as given by Eq. (15). The dipole is along z-axis. Several values for dipole strength [chi] 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.

Caption: Figure 2: Equipotential contours of the effective potential V(x, z) at y = 0.5, as given by Eq. (15). Dipole along z-axis. Far from center the isopotential surfaces are deformed spheres.

Caption: Figure 3: The effective potential [V.sub.sph](a) in the radial part of the Schrodinger equation in spherical coordinates, Eq. (26). Solid lines represent c = 0 potentials, and dotted lines represent c = 5 potentials, for the same [L.sub.[chi]],

Caption: Figure 4: The linear solutions |Y|2 = 0.5, 0.6, 0.75, and 0.9 isoplot surfaces for [chi] = 0.5, [a.sub.[infinity]] = 10, L = 1, calculated with the quasipolynomial 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.

Caption: Figure 5: Examples of the application of the boundary conditions at [a.sub.[infinity]] = 10, for the dipole equation analytic solution. The positive part of the spectrum 0 < E < 1, for L = 0,3, and [chi] = 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.sub.[infinity]]. The number of eigenstates (oscillations) increases with L, and decreases with [chi], because for smaller dipole strength the solution is closer to the Bessel solutions. The eigenvalues of energy are E(L, [chi]) = 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.

Caption: Figure 6: Energy spectrum E for the linear dipole problem with boundary conditions on the volume from Eq. (47), or Eq. (54), for small dipole strength, [chi] = 0.6. We present two cases of volume radii, R = [a.sub.[infinity]] = 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.

Caption: Figure 7: Energy spectrum E for the linear dipole problem with boundary conditions on the volume from Eq. (47) for two dipole strength, [chi] = 1.5,3 and two volume radii, R = [a.sub.[infinity]] = 10,60.

Caption: Figure 8: Eigenvalues (0 < E < 1) and eigenfunctions for the linearized dipole equation for [chi] = 0.5 and L = 0 (upper frames) and L = 3 (lower frames) with boundary condition at [a.sub.[infinity]] = 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.sub.[infinity]], [chi],L values.

Caption: Figure 9: The number of eigenvalues E in the positive spectrum versus [a.sub.[infinity]] for L = 0, ..., 5. The number increases with L. The upper band of lines is for [chi] = 1, and the lower one for [chi] = 5. Obviously the number of eigenvalues increases with the radius of the volume since larger samples allow more oscillations inside.

Caption: Figure 10: Isoplot surfaces [absolute value of ([[PSI].sup.NL])] = 0.9 for L = 1,2,3,4 from left upper frame, clockwise. We choose [a.sub.[infinity]] = 100, [chi] = 1, [L.sub.1] =0,[k.sub.1] = 1 and the minimal energy was obtained for [k.sub.2] = 3.

Caption: Figure 11: Left: coordinate surfaces of constant a = 1,2,3 (tori-like); Right: coordinate surfaces of constant b = 1,2,3 (double-volumes). The dipole is directed along z axis.

Printer friendly Cite/link Email Feedback | |

Title Annotation: | Research Article |
---|---|

Author: | Ludu, A. |

Publication: | Advances in Mathematical Physics |

Article Type: | Report |

Geographic Code: | 1USA |

Date: | Jan 1, 2018 |

Words: | 12707 |

Previous Article: | Autonomous Jerk Oscillator with Cosine Hyperbolic Nonlinearity: Analysis, FPGA Implementation, and Synchronization. |

Next Article: | Regularization of the Boundary-Saddle-Node Bifurcation. |

Topics: |