Abstract
We investigate the existence of multivortex states in a superconducting mesoscopic sphere with a magnetic dipole placed at the center. We obtain analytic solutions for the order parameter 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.
1. Introduction
The rapidly growing field of quantum computation requires nanoscale miniaturization of electronic circuits, way beyond the silicon era type of devices [1]. Therefore, the mesoscopic superconductors, having the size comparable to the coherence length or the magnetic penetration depth [2], are the prime candidates for construction of nanodevices among all other superconducting systems. Mesoscopic physics revealed a number of open fundamental problems like quantum confinement, quantum vortices and loops, spintronics, quantum dots, etc. [3] whose solutions can bring significant knowledge in fields like nanotechnology, synthesis of new materials, novel sensors, modern lithography, or molecular biology [4–6].
The most important feature of a mesoscopic superconductor is that its shape and size have an important effect on the interplay of the magnetic field and superconducting condensate. The properties of mesoscopic superconductors are very different compared to those of bulk superconductors. While in bulk superconductors penetrating vortices form a lattice due to the vortex-vortex repulsion, in mesoscopic superconductors there is a competition between the vortex lattice and the boundary which tries to impose its geometry on the vortex lattice. It is observed experimentally that flux quantum configurations have the same symmetry as the symmetry of the shape of the sample in a homogeneous magnetic field [7]. Such systems are studied via the Ginzburg-Landau (GL) model. The GL equations arise from the Euler-Lagrange equations for the free Gibbs energy for a mesoscopic superconductor sample in magnetic field. These equations must be solved under specific boundary conditions: the normal component of the superconducting current is equal to zero [8]. Near and below the transition temperature, theoretical calculations show that anti-vortex and giant-vortex can appear in order to maintain the symmetric vortex configuration.
The response of mesoscopic superconductor samples of different shapes (thin discs, spheres, cones, and rings) to an external magnetic field, as well as the effect of the geometry, has been theoretically [9–18] and experimentally [7] studied. In all these cases a constant external magnetic field is applied along the the revolution axis. The small volume to surface ratio of these mesoscopic structures brings new features not found in the bulk: there are two kinds of superconducting states, depending on the strength of the magnetic field, the sample geometry (external surface), and its size: giant and multivortex states. The giant vortex state has cylindrical symmetry and is the only kind stable in small size superconductors due to the confinement effect [13, 14]. If the size of the sample increases, such giant vortex states can break up into multivortices through saddle-point transitions [17, 18]. For three-dimensional objects (sphere or cone) the vortex lines need to intersect the surface perpendicularly in order to cancel the outward supercurrent component [9–12]. Consequently, the shape of the lines is strongly affected by the sample surface. For example, in the case of a mesoscopic sphere placed in uniform external magnetic field, the vortex lines are curved inside, packing denser in the equatorial plane, and spreading out towards the poles [10, 11].
In this paper we consider a different situation where the magnetic field is not anymore externally generated but is generated from inside the sample, e.g., an infinitesimal magnetic dipole placed at the center of a mesoscopic sphere. Such a configuration can generate confined vortex loops. The topological transition between open and closed vortex loops is controlled by the geometry, i.e., , 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–18]. The free energy functional is given in the GL model by [10, 11]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 and the nonlinear term coupling constant are typical Landau second-order phase transition parameters [9–12]. For mesoscopic samples, one can neglect the term responsible of the expulsion of magnetic flux from the superconductor, that is last term in Eq. (1) [9–12].
The traditional procedure for finding by minimizing the free energy functional Eq. (1) for constant volume consists in expanding the order parameter in a basis of eigenfunctions of the corresponding linearized GL problem [10, 11, 16–18] and numerically evaluating the expansion coefficients which minimize the full nonlinear functional Eq. (1).
In this paper we solve the GL linear problem analytically and investigate the properties of the eigenfunctions and spectrum. The contribution of the infinitesimal magnetic dipole will be approached in the dipole system of coordinates. The linearized GL equation factorizes in two ordinary differential equations, for the two orthogonal dipole coordinates, respectively. From the physical point of view such a factorization seems natural because far enough from the sphere surface, the abstract surfaces containing vortex lines follow the magnetic field stream lines, but these are exactly the dipole magnetic field lines. Along these surfaces the order parameter has slow variation or is practically constant. This gives sufficient physical reason for neglecting higher order terms in the dipole variable going along the magnetic field lines. This approximation allows integrating the two ordinary differential equations. The linear solution consists in a product of angular momentum eigenfunctions in the azymuthal coordinate, exponential function for one of the dipole coordinates, and a double confluent Heun function in the third coordinate. The final step is to come back to the fully nonlinear GL problem, and write as a linear combination of eigenfunctions of the linear GL problem, with arbitrary coefficients, and then minimize the free energy with respect to these coefficients.
We dedicate a large part of the present calculations to solve exactly the linearized GL equation and to ensure the completeness and orthogonality of the linear basis because near and below the transition temperature, even the linear GL equation is sufficient to describe multivortex states. The order parameter is very small in this range, and higher order terms of can be neglected. Nevertheless, at lower temperatures, the vortex configuration does not have to match the symmetry of the system, and higher order terms of GL equation cannot be negligible [19]. It is the contribution of these nonlinear terms which generates the multivortex states at lower temperatures.
The paper is organized as follows. In 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 produced by a infinitesimal dipole placed at its center . The Euler-Lagrange equations for the free energy functional Eq. (1), written in the absence of the last term in as explained above, are known as the GL equation [9–12]where the electromagnetic momentum is given bywhere 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, , and we rescale the variables and . We follow the GL phenomenological model from [9–18] and use the temperature dependencewhere 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 equationSince 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 as eigenfunctions labeled by the angular momentum . This decomposes the space of solutions into subspaces parameterized by with 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 problemwhere 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 coordinateswhere is the unit vector in the direction. The boundary conditions associated with this problem request the order parameter wave function to cancel at the origin and the flow at the mesoscopic volume surface to be zero (the Saint-James-de Gennes conditions) [9–12, 15–18]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 of the forth power of the magnitude of the dipole vector potential is finite, we find out that the operator in Eq. (6) is essentially self-adjoint if applied on analytic functions. Moreover, from the Miller-Simon theorem [21], we know the eigenvalue problem of such type of magnetic multipole potential vector operators is given by a positive, unbounded from above, continuous spectrum. However, given the two-point boundary conditions, Eqs. (12), we expect the reduction of this spectrum to a discrete and bounded spectrum for energy.
3. The Generalized Dipole Equation
In this section we present analytic solutions for the system Eqs. (6), (8), and (10) with boundary conditions given by Eq. (12). The presence of the magnetic dipole renders the spherical or the cylindrical coordinates unsuitable to factorize the LGL Eqs. (8) and (10) because of the terms mixing inside the square of the electromagnetic momentum. Attempts to solve similar equations in spherical coordinates in the presence of an electric dipole have been made. In [22], for example, the Schrödinger equation is separable, and the authors find exact wave functions in terms of Bessel and Mathieu functions for 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 in a series of Jacobi polynomials with “dipole polynomial” coefficients.
In all these situations separation and integration of the Schrödinger equations are possible because in spherical (cylindrical) coordinates the resulting second-order ordinary linear differential equations are of Fuchsian type, with maximum three regular singularities. Such equations can be mapped into several types of hypergeometric differential equations.
However, the presence of the dipole magnetic field makes the situation more complicated because the order of singularity growths above the hypergeometric one, and consequently the dynamics is governed by differential equations of Heun, Hill, or Riemann type [25–27]. A similar situation occurs in the case of a charged particle moving on a sphere under a radial magnetic field and Coulomb force. The corresponding Schrödinger equation is transformed into a Heun equation in canonical form, and exact solutions are obtained in terms of series of hypergeometric functions. The separation of variables is still possible since the vector potential depends only on one spherical variable, in this case [28]. Another example of the same type of difficulty is the two-dimensional case of interaction of three particles with a perpendicular magnetic field [29]. Here the Schrödinger equations map into biconfluent Heun’s equations, too, through the higher order singularities induced by the Coulombian interaction between particles. Similar problems related to magnetic field (finite-gap potentials, Fokker-Planck, central two-point connection, generalized central potentials up to order , Hawking radiation, etc.) were approached in the literature and usually the resulting leading differential equation for the wave function reduces to one of Heun’s differential equations [30–34]. An interesting review and study on the use of the Heun’s type of differential equations as generalizations of the hypergeometric ones, in relation to supersymmetry, are given in [35], where a two-Coulomb-center problem is solved based on a self-adjoint separation of coordinates in prolate spheroidal coordinates.
Before moving to the dipole coordinates we perform a qualitative analysis of Eq. (6) in spherical coordinates. In spherical coordinates Eqs. (6) and (10) reduce to a linear Schrödinger equation in the formwhere we denoted , 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 . We can write these wave functions in the formwith being an arbitrary positive integer. Eq. (13) reduces to a 2-dimensional Schrödinger equation with an effective potential in the formIn Figures 1 and 2 we present some equipotential contours of the effective potential at and . This potential has a repulsive tail that drops to zero towards infinity like . 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 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 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 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 multiplied with spherical harmonics of different orders decomposes (through the Wigner -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 relationsDipole 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 expressionBy 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 from Eq. (A.3). This approximation describes the order parameter in neighborhoods of , 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 , according to Eq. (A.7) in Appendix A. For example, if we just choose we have a relative error of only or less for within a range of the azimuth angle , which represents that 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 formand we can factorize further the -dependence from the -dependence by the substitutionwhere 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 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 () and south () poles the volume is depleted () 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 , [42], and consequently there are no convergent solutions in terms of Frobenius series. We need to solve this equation in the interval , 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–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., for some real values of ) with some known physical significance. For example, by the substitution we havewhich is a 1-dimensional Schrödinger equation with an effective potentialHere the last term represents the traditional centrifugal term plus an additional constant . 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 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 coordinatesThe 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 , and depending on the value of . In the simpler case , for example, if this potential has a valley centered atThis valley creates a barrier of heightand relative depth with respect to the top of the valleyFrom the aspect of this potential it results that for weak coupling () 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 .
For the substitution we obtain the radial part of the Schrödinger equation in cylindrical coordinatesEq. (30) is a Fuchsian differential equation with singularities at the point at infinity, and order six singularity at , 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 fourth-order 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 , and it arises somehow naturally since the order parameter has its dominant dependence on the coordinate. This happens because the iso-surfaces are basically laying along the dipole surface coordinates constant. Under this hypothesis, Eq. (22) becomesand 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 , , in the potential energy. The resulting equation is mapped into the Mathieu equation for through the transformation . It is interesting that even holding all the terms in the powers of 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 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 coordinate, that is, for . 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 formwhere 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–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 and , where . By using a transformation of the dependent variable of the formany 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, . It is worth mentioning that one obtains a nondegenerate (a “good”) canonical DCHE form if one has the conditionfulfilled (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: , 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 where , andare 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 determine the leading coefficients of the asymptotical series of the solutions at and , respectively. However, in the case of the dipole equation we have , 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 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 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, . 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. 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 , though because 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 , bounded by the corresponding Stokes rays, depending on the argument of [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 . 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 . 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 transformationsandThe 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 seriesin the range on arg, with the coefficients uniquely given by the three-term recursion relationwith N, . The signs represent the choice for one of the two solutions for the reduced parameters in Eq. (36). We mention that the second solution is generated from the first one by application of the group operator on , 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 the solutions , that is, , 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 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 . For , however, the majority of solutions have , so they cannot be reduced to quasipolynomials.
In the following, we present two examples of quasitruncated solutions obtained from the abstract solutions :with , andwith . In order to have a three-dimensional visual on the behavior of these solutions we plot in Figure 4 few surfaces of equation 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.
5. Eigenvalue Problem for the Dipole Equation
Typical study of the quantization of the energy spectrum of Heun’s equation involves the accessory parameter. Such an analysis is performed, for example, in [45] where the authors study the eigenvalues of a double confluent Heun equation with possible applications in gravitational theory. The boundary conditions used in this paper differ from our case, since they ask the solution to be bounded in origin, and to approach either a finite value, or infinity when , 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) iswithIn this representation we have , and 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 with . According to the physical restriction , 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 haveNormalized with respect to the weight , the eigenfunctions form a complete orthonormal basis of piecewise continuous functions, with piecewise square integrable continuous derivative in .
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 formwhere we denotedThe boundary condition Eq. (47) contains two unknowns, . 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 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 giveswhere we denoteOn the other hand, we calculate the asymptotic behavior towards infinity of the two linear independent analytic Heun solutions from , Eqs. (39a) and (39b). It results in an asymptotic term of the formto be matched with Eq. (49). The solution of this matching condition results inwhere has the same meaning as above, Eq. (50).
With the coefficients 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 . The resulting solutions of the linearized GL equation areWe present some examples of such solutions in Figure 5 for 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 isFor 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 . 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 . 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 . The maximum number of spectral levels depends on and it can be calculated from the condition . The next spectral domain is within the negative energiesand it also depends on all the parameters. In the normal range of values (, ) this condition offers a wide range of energies (it can go down to ) so this is the main range for negative energies. In this domain, because the argument of the function in Eq. (54) 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 () and an average value of the relative dipole strength . 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 expressionwith integer values of such that . The number of eigenvalues for positive in this approximation can be estimatedand 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 we have 5 eigenstates for and 4 eigenstates for .
For larger radii of the volume the energy spectrum obtained as solution of Eq. (54) can be approximated withwith 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 three-dimensional problem, a third quantum number was eliminated by the supposition that is mainly directed orthogonal to the magnetic dipole field lines. Consequently, the order parameter is almost constant along these lines, and the coordinate surfaces 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.
6. 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 the connection between the physical parameter and the spectrum label being provided by Eq. (54).
We search nonlinear solutions of the formwith 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–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 of the dipole coordinate such that the equation has a number of distinct solutions for . If we denote by these solutions, by analytic prolongation theorem we can always find a neighborhood of 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 and surfaces.
In literature [15, 16] there are examples of such multivortex states obtained with as little as 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 () the expressionwith and all functions are evaluated at . To visualize an example of multivortex state, we choose , and , . We plot in Figure 10 four isoplot surfaces with 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 Schrödinger equation with mixed boundary conditions at the center and the surface of the superconducting sample. The goal of this paper is to obtain the nonlinear order parameter, i.e., the solution 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 , under the boundary conditions Eq. (12). As a first step, we use cylindrical symmetric coordinates. The linear solution can be separated in , 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 and , and the solution can be written as , Eqs. (20), (21), (22), and (30). We reduce further the resulting equation by the mode decoupling condition , and we obtain the dipole equation (31) for which is a double confluent Heun equation, Eq. (35). By performing a substitution, Eq. (33), in the form and we can finally solve the dipole equation and write its analytic solutions, Eqs. (39a) and (39b).
The two approximations introduced here are possible because we confine our calculations in the neighborhoods of the poles of the sphere. Since the gradient of is mainly directed radial and orthogonal to the magnetic dipole field lines, it results that is almost constant along these lines, and the dipole coordinates surfaces can describe the multivortex isosurfaces.
We build the nonlinear solution for the GL problem as linear combinations of the exact solutions of the linearized dipole equation and tune the combination’s coefficients that minimize the free energy functional. Further on, we demonstrate that multivortex states are possible in this sphere plus dipole configuration, even without the presence of an external magnetic field. Our results are in very good agreement with the numerical calculation of the same model [10, 11].
Appendix
A. Dipole Coordinate System
The dipole coordinates defined in Eq. (16) have as domain of definition for . The Jacobian of the transformation from spherical coordinates and the element of volume are given by
The dipole coordinates are actually different from other apparently similar coordinates like bipolar, bispherical, or toroidal coordinates [54–56]. In the notation used above , these traditional orthogonal curvilinear systems of coordinates have symmetrical forms for the three Euclidean coordinates of the general formwhere the functions are trigonometric or hyperbolic functions.
Consequently, the two parameters are in general an angle and a distance. In the case of dipole coordinates, the parameters have different dimensions, but both are related to distances, . 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 ) or tangent closed curves with a common singular point (for ). 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 equationThis 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 . This solution readswhere we defineandWe expand the above expression in Taylor series with respect to around zero, and we obtainwhere obviously the dominant term is . This approximation is valid around the North and South poles of the glued, vertically aligned, deformed circles of equation , Figure 11. The solutions of a similar type of quartic equation in ,can provide the inverse .
A more compact expression can be obtained if we solve the inverse transformation Eq. (A.3) in cartesian coordinates, namely, as functions of , where . Beginning with relation Eq. (A.3) and the above definitions for , with we obtainAfter some elementary algebra we obtain the inverse in the formwhich can be compressed even more intowhereThe 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 Lammè coefficients of the dipole coordinate system areConsequently, the Laplacian in coordinates has the form
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 and . We can write these series in the formwhere is the Kummer function, the indices independently (which generates 8 solutions), and is the summation label which is an integer plus an arbitrary complex parameter 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 are obtained from the three-term recursion relationwhere and are the reduced parameters Eq. (36) for the DCHE associated with the linearized dipole equation (35), and independently.
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 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.