Vortices on a superconducting nanoshell: phase diagram and dynamics
Abstract
In superconductors, the search for special vortex states such as giant vortices focuses on laterally confined or nanopatterned thin superconducting films, disks, rings, or polygons. We examine the possibility to realize giant vortex states and states with nonuniform vorticity on a superconducting spherical nanoshell, due to the interplay of the topology and the applied magnetic field. We derive the phase diagram and identify where, as a function of the applied magnetic field, the shell thickness and the shell radius, these different vortex phases occur. Moreover, the curved geometry allows these states (or a vortex lattice) to coexist with a Meissner state, on the same curved film. We have examined the dynamics of the decay of giant vortices or states with nonuniform vorticity into a vortex lattice, when the magnetic field is adapted so that a phase boundary is crossed.
pacs:
PACSI Introduction
Quantized vortices are a quintessential property of superfluids and superconductors. The energetically favored state when multiple quanta of vorticity are present, is a lattice of singly quantized vortices. In ultracold Fermi gases, the recent observation of such a vortex lattice formed the ‘smoking gun’ proof for superfluidity ZwierleinNAT435 (). In nanoscopic superconducting samples, controlling the vortex behavior is essential for the development of new devices based on fluxon dynamics Moshchalkov93 (). The confinement of Cooper pairs to the length scales comparable to the correlation length also offers the prospect to probe fundamentally new phase topologies predicted by the theory, such as giant MoshchalkovGiant (); MiskoPRL90 () and ringlike vortices StenuitPC332 (). This has led to renewed experimental efforts to observe giant vortex states, both in superconductors Giantvortex (); MiskoPRB64 () and in superfluid atomic gasesGiantvortexBEC ().
In this contribution, we argue that superconducting spherical nanoshells form a promising candidate for realizing giant vortex states, and for engineering phase transitions between those states and a vortex lattice. Moreover, we show that nanoshells allow the coexistence of a Meissner state and a vortex state in equilibrium on one and the same superconducting film. Nanoshells are hybrid nanostructures consisting of a dielectric core (usually a silicon oxide nanograin), coated with a thin layer of metalNanoshell (). When the metal in its bulk form is a superconductor, the nanoshell below the critical temperature will also exhibit superconductivity in the thin shell around the isolating core.
The superconducting order parameter in the nanoshell is well described by a macroscopic wave function that obeys the coupled timedependent GinzburgLandau (TDGL) equations. Vortices are characterized as topological defects in the phase (requiring a vanishing gap ). For thin shells, the description is simplified in two important ways. Firstly, when the shell thickness is much smaller than the London penetration depth, the magnetic field will be only weakly perturbed by the nanoshell. Secondly, when the shell is thinner than the coherence length, the order parameter will not vary substantially in the radial direction in the shell; that is, will only depend on the spherical angles . In the radial direction, will be constant in the shell, and zero outside it. Note that confining to the shell leads to an effective GinzburgLandau parameter that differs from its bulk value. In Section II we present the formalism, and in Sec. III the results, for thin shells. When the shell thickness is increased and becomes nonnegligible with respect to the penetration depth, the magnetic field will be more strongly perturbed, and the field gradients affect the energetics. This case and the effect on the phase diagram are discussed in Sec. IV. Finally, we summarize the results for vortices in nanoshells in Sec. V.
Ii GinzburgLandau formalism on thin shells
We assume that the shell is sufficiently thin for neglecting variations of the order parameter across the shell. In other words, the order parameter will only depend on the spherical angles . We use the spherical coordinates , , with the origin at the center of the sphere. The angle is counted from the axis parallel to the external homogeneous magnetic field. Like in Ref. Zhao, , we will make the used variables dimensionless by expressing lengths in units of , magnetic fields in units of , and the vector potential in units of , where is the penetration depth, is the magnetic flux quantum, is the Planck constant, and is the elementary charge. Thus, the dimensionless parameters , , and are linked to the radius of the nanoshell , its thickness , and the applied magnetic field by the expressions , , and , respectively.
In our numerical treatment of superconducting states on spherical shells we exploit the timedependent GinzburgLandau equation, which is known to be a powerful tool for studying both the dynamic and static properties of superconductors. For a thin shell under consideration, the behavior of the order parameter in a fixed (or slowly varying) magnetic field can be described by the TDGL equation (cp. hu72 (); kato91 ())
(1) 
where is the GinzburgLandau parameter, is the (dimensionless) vector potential, and . The dimensionless variable is linked to the time by the relation , with , the normalstate diffusion constant.
The vector potential can be represented as a sum of the contribution , related to supercurrents in the shell, and the contribution , which corresponds to the external magnetic field . The vector potential is chosen in the form
(2) 
In the case of a constant applied magnetic field , with increasing the function , given by Eq. (1), approaches one of the (meta)stable states of the system (). The thermodynamically stable state is to be found by comparing the Gibbs free energy for different solutions. The difference in the Gibbs free energy between a superconducting state and the normal state at the same magnetic field is given by the equation
(3) 
where corresponds to the superconducting state with no vortices at , i.e. the Meissner state present on the complete surface. The dimensionless density of supercurrents is denoted by and expressed in units of .
In this section, the shell is assumed to be sufficiently thin in order to make negligible the magnetic fields, induced by supercurrents. Correspondingly, we can neglect as compared to . Then, as seen from Eqs. (1) and (2), two independent parameters, which govern the solution of Eq. (1), remain:

the dimensionless size of the nanoshell , determined by the ratio of the shell radius to the GinzburgLandau coherence length , and

the parameter , equal to the number of flux quanta of the applied field that pass through the equatorial plane of the sphere.
When the magnetic field is increased beyond a critical value (computed below), a first vortex appears for nanospheres with radius large enough to sustain the vortex core, as depicted in Fig. 1, panel (a). Upon further increasing the magnetic field, more quanta of flux can penetrates the spherical surface. This can be accomodated in a variety of ways: for example as a giant vortex carrying more than one quantum , shown in panel (b) of Fig. 1. In this case, the angular momentum is uniform over the spherical surface. It is also possible to envisage states with nonuniform distributions of angular momentum: a value near the poles and a value in a band near the equator. Such states are characterized by a ringlike vortex separating the regions with different angular momentum, as illustrated in panel (c) of Fig. 1. Also states that do not have axial symmetry should be investigated: we will show that these are in many case the stablest state and that they then consist of an array of singly quantized vortices, illustrated in panel (d) of Fig. 1. Such states will be denoted by multivortex states, to emphasize that every vortex carries a single quantum of flux. In the next subsections we start by investigating the axially symmetric states: giant vortices and ringlike vortices. Then, in the next section we investigate the condition under which those states decay into multivortex states and the dynamics of this decay.
ii.1 Giant vortex states
First, let us consider superconducting states, which keep the axial symmetry of the system, so that the order parameter can be written in the form . where has the sense of the winding number (vorticity). Then for a stationary distribution the GinzburgLandau equation (1) reduces to the onedimensional equation
(4) 
with boundary conditions, determined by the requirement that the component of the current density must be zero at the axis: . Solid lines in Fig. 2 illustrate typical behavior of the freeenergy difference as a function of for cylindrically symmetric states with different vorticity . In the case of a thin spherical shell with , as illustrated in panel (a) of Fig. 2, the value of in the lowest cylindrically symmetric state increases with from 0 at to 9 at . The modulus and phase are illustated in panel (b) of Fig. 1, using hue and saturation of the color scale respectively. An increase of the applied magnetic field is seen to result also in a significant increase of the Gibbs free energy of the lowest state.
ii.2 Ringlike vortices
In Ref. Zhao (), when analyzing supercondicting states in hollow cylinders, it was suggested that – under certain conditions – cylindrically symmetrical states with changing winding number can be more energetically favorable than the states with uniform . Our calculations show that a similar situation occurs also in thin spherical shells with the dimensionless size larger than (i.e., for ), but as we will show in the next section, such states decay into a lattice of singlyquantized vortices breaking the cylindrical symmetry.
We have compared the Gibbs free energies for axially symmetric states with uniform winding number and those for states where , the winding number at and , differs from , the winding number at . The order parameter of the latter states on the sphere are illustrated in panel (c) of Fig. 1 – we will refer to such states as “ringlike vortex” states. The continuity of the order parameter as a function of requires vanishing at the boundaries between regions with different winding numbers, i.e., at and as can be seen in panel (c) of Fig. 1. At sufficiently strong magnetic fields the Gibbs free energy for states with ringlike vortices can become lower than that for states characterized by a unique winding number over the whole range. This is illustrated by Fig. 2, where the dotted lines show the calculated freeenergy difference for states with and different on a shell with in the case of . It is worth mentioning that the values , which minimize for the lowest ringlike vortex state at a given , are rather insensitive to the dimensionless nanoshell size (at least, for ). At the same time, the parameter is an increasing function of . Thus, our calculations show that this parameter changes from to when increasing from 7 to 15. However, moderate variations of around only slightly affect for the lowest state with ringlike vortices. That is why in Fig. 2 we restricted ourselves to the case of a fixed value , which coincides with at .
In Fig. 2, the curve labeled with ‘1,1’ corresponds to the state with a ringlike vortex and uniform vorticity, which is qualitativley similar to the states analyzed in Ref. StenuitPC332, . The free energy of this state is always significantly higher than for the lowest giant vortex states with pointlike core. As can be further seen from Fig. 2, at the ringlike vortex states () appear the most energetically favorable among the axially symmetric states. Thus, at , the difference in between the state with and the state with , is larger than 0.025. At even larger values of , states with three different regions of vorticity (characterized by two ringlike vortices) can become stable. However, our calculations show that in shells with , where such ringlike vortices allow for decreasing the free energy of giant vortex states as compared to the case of a giant vortex, even lower values of can be achieved by breaking up the ring like vortex (or vortices) into an array of singlyquantized vortices. A natural question arises of how stable are the aforedescribed giant vortex states with respect to decay into multiple singly quantized vortices. In order to answer this question, one has to return to the TDGL equation (1).
ii.3 Numerical treatment
The finitedifference scheme, applied here to solve Eq. (1), is similar to that of Ref. kato91, , with necessary adaptations to the case of a spherical 2Dsystem. Twodimensional grids, used in our calculations, typically have equally spaced nodes in the interval from 0 to and equally spaced nodes in the interval from 0 to . Cyclic boundary conditions for are applied at and . The boundary conditions at and are determined by the requirement . The step of the time variable is automatically adapted in the course of calculation. This adaptation is aimed to minimize the number of steps in , necessary for approaching a steady solution of Eq. (1), and – at the same time – to keep the solving procedure convergent. On average, the step in is to depending on the used grid as well as on and . When starting at from a random distribution of (with ), a (meta)stable solution of Eq. (1) is achieved typically at . When analyzing (meta)stability of states in a spherical shell, one has to keep in mind that a transition between states with different vorticity, in general, requires symmetry breaking. This means that simulations, which assume a perfectly symmetric spherical nanoshell, would tend to overestimate stability of a state with respect to a possible transition to another state with lower free energy. In order to model the effect of imperfections, inevitably present in realistic nanoshells, we consider spherical shells with small angular variations of the parameter . Importantly, for relative magnitudes ranging roughly from to , the results of simulations practically do not depend on a specific choice of the magnitude and distribution of these nonhomogeneities. An appreciable effect of those imperfections on stable distributions of the order parameter appears only for .
Iii Results and discussion for thin shells
iii.1 Decay of giant and ringlike vortices
In order to examine the stability of giant vortex states with respect to decay into multivortex states, we apply the computation scheme described in the previous section, starting at from a distribution of that corresponds to a giant or ringlike vortex state. Typical examples of the evolution of the order parameter distributions are shown in Figs. 3 and 4 for the cases when the initial state is a giant vortex and a ringlike vortex, respectively. In the case of and , the thermodynamically stable state corresponds to 7 pairs of vortices with a single quantum of flux each, and it has a relative free energy , approximately 0.07 lower than the value of for the lowest giant vortex state (see Fig. 2). Such lattices of vortices with each a singly flux quantum will be denoted as “multivortex states”.
As illustrated by Figs. 3(a) to 3(c), within a interval the initial giant vortex state with transforms into a chain of 7 singly quantized vortices, which surround each pole of the sphere. In the course of the further rearrangement of the vortex pattern, one of the vortices moves to the pole, while the remaining 6 vortices tend to form a symmetric chain around the pole (see Figs. 3(d) to 3(f)). A freeenergy gain due to this rearrangement is by one order of magnitude smaller than that due to the decay of the initial giant vortex into single vortices. Correspondingly, the interval, necessary for this rearrangement, appears to be relatively long: only at the solution reaches the equilibrium symmetric configuration of vortices (not shown in Fig. 3), similar to that found in Ref. Du04, , where multivortex states on a thin hollow sphere were studied in detail. As seen from Fig. 4, the transition from a ringlike vortex state ( for and ; for ) to a multivortex state is even faster. The ringlike vortex core, which is present in the initial state [see Fig. 4(a)], decays into a chain of 6 single vortices very quickly: clear signatures of this decay can be found already at [see Fig. 4(b)]. The equilibrium state with a vortex at the pole and 6 vortices, symmetrically surrounding the pole, is formed already at .
The results of our calculations clearly indicate that giant and ringlike vortex states are rather unstable in spherical shells with relatively large . This does not mean, however, that giant vortex states on a spherical shell are never stable. A decrease of the shell radius and/or an increase of the applied magnetic field enhance the role of the Lorenz forces, which act on the supercurrents and tend to drive vortices towards the poles of the shell. As a result, for sufficiently small and sufficiently large , the distance between vortex cores in a multivortex state becomes so small, that physically a multivortex state appears undistinguishable from the corresponding giant vortex state. A similar continuous transition from a multivortex state to a giant vortex states with increasing magnetic field was recently found when solving the linearized GinzburgLandau equation for superconducting spherical grains Baelus (). Of course, in the case of such a continuous transition, the boundary between thermodynamically stable multivortex states and giant vortex states can be drawn only approximately. As a criterion of a transition from a multivortex state to a giant vortex state, here we have chosen the condition that the angular distance of vortex cores from the pole becomes smaller than .
iii.2 Phase diagram for thin shells
Our results, related to thermodynamically stable states on thin spherical shells, are summarized in Fig. 5, where the solid lines indicate boundaries of stability regions for the normal state, the superconducting Meissner states, single vortex states, giant vortex and multivortex states. The dashed line indicates the boundary between the regions, where giant vortex states (to the left from this line) or multivortex states (to the right from this line) are thermodynamically stable. As seen from Fig. 5, formation of vortices can be energetically advantageous only on sufficiently large shells: at for states with , at for giant vortex states, and at for multivortex states.
While in Fig. 5 the phase diagram for thin spherical superconducting shells is shown in the plane, it seems interesting to analyze the phase boundaries also in more common form: in terms of the applied magnetic field and the temperature . We assume that the temperature dependence of the penetration depth is described by the empirical relation , while the (less important) temperature dependence of the GinzburgLandau parameter is roughly given by the expression (see, e.g., Ref. Tinkham, ). In Fig. 6 we plot the phase boundaries for thermodynamically stable normal states, Meissner states, single vortex, giant vortex, and multivortex states on spherical superconducting shells with different radius , measured in units of the zerotemperature GinzburgLandau coherence length . As illustrated in Fig. 6(a), in the case, where the radius is much larger than , giant vortex states are thermodynamically stable (at moderate applied magnetic fields) only in the close vicinity of the critical temperature . With decreasing , the stability range of giant vortex states gradually extends towards lower temperatures and lower values of the magnetic flux through the shell. Correspondingly, in sufficiently small shells the stability range of multivortex states is restricted to a relatively narrow interval of and to , significantly lower than [see Fig. 6(b)]. On even smaller shells (with close to ) no stable multivortex states are possible [see Fig. 6(c)]. For those small shells, the superconducting phase persists only for relatively weak (few ) magnetic fluxes through the shell. At the same time, as implied by a comparison between panels (b) and (c) in Fig. 6, the values of the applied magnetic field , which correspond to transitions with an increase of vorticity by 1, become significantly higher when decreasing the shell size.
Iv Vortex states on thick shells
iv.1 Magnetization effects in thick shells
Now, let us extend our analysis to the case of relatively thick spherical shells, where magnetic fields, induced by supercurrents, which flow in a shell, are of nonnegligible. At the same time, we assume that the thickness of a shell is still sufficiently small for neglecting variations of the order parameter and of the vector potential across the layer. For such a shell, also currents across the layer can be neglected. Expressing the vector potential through the density of current as
(5) 
the nonnegligible components of the product , which enters Eq. (1), can be written down in the following form:
(6)  
(7)  
where is the dimensionless thickness of the shell. On the other hand, in the case of constant or slowly varying magnetic fields, using the relation
(8) 
the products and , which enter Eqs. (6) and (7), can be expressed through , , and as
(9) 
(10) 
In order to find the order parameter and the corresponding vector potential, we solve selfconsistently the set of equations (1), (6), and (7), using relations (9), and (10). From Eqs. (1), (6), (7), (9), and (10), one can see that for relatively thick shells under consideration a set of independent parameters, which govern the solution, can be chosen as , , and , where the introduced additional parameter is linearly proportional to the thickness of the nanoshell and to its radius.
Figure 7 gives few examples of magneticfield distributions, which correspond to thermodynamically stable states in spherical shells with , and , . The patterns of magneticfield lines, displayed in Fig. 7, are plotted for the particular case of the GinzburgLandau parameter , the (mean) dimensionless radius of the shell , and the dimensionless thickness and . In general, none of the three mutually orthogonal components of the magnetic field is zero, so that the field lines are “threedimensional”. In Fig. 7, however, we restrict ourselves to fieldline patterns within symmetry planes, where the field lines are “flat”. As seen from Fig. 7(a), even in the case of a relatively thick shell () the magnetic fields, induced by the supercurrents in the Meissner state, are not sufficient for complete screening of the applied magnetic field inside the shell. Nevertheless, not only in the case of but also for a significantly thinner shell with [Fig. 7(b)], the net field inside the shell is much weaker than . In the case of the state with , the magnetic flux, captured by a vortex pair, is seen as an increased density of field lines at the poles of the sphere [Fig. 7(c)]. At the same time, in the depth of the sphere the magnetic field is relatively homogeneous, only slightly increasing towards the axis. Also for states with higher vorticity, a considerable local increase of the magneticflux density takes place only at the vortex cores within the superconducting shell, while in the depth of the sphere the density of magneticfield lines is considerably more homogeneous [see panels (d) to (f) in Fig. 7].
iv.2 Phase diagram for thick shells
As seen from Fig. 7, the magnetic fields, induced by supercurrents, can be considerably large even for shells with quite moderate thickness (). These fields strongly affect the stability range for superconducting states with different vorticity in a spherical shell. In Fig. 8, we present the calculated phase diagram for relatively thick spherical shells with . As follows from a comparison of Fig. 8 to Fig. 5, an increase of the thickness of a spherical shell results in a wellpronounced shift of the boundaries between states with different vorticity towards higher magnetic fields . In particular, for , the range of , where Meissner states are thermodynamically stable, is more than two times wider in the case of as compared to the case of . One can also see that for a relatively thick spherical shell () the boundary between giant vortex and multivortex states is shifted towards significantly larger values of as compared to those in the case of . The increased stability of giant vortex states agree with the results recently obtained by Baelus et al. Baelus () for the limit of a full sphere, in the framework of linearized GinzburgLandau equations.
In Fig. 9 the phase boundaries for thermodynamically stable normal states, Meissner states, single vortex, giant vortex, and multivortex state states are plotted in the plane. The results are shown for shells with different radius and thickness . In order to keep the plots more universal, it is convenient to express the thickness in units of . For thick nanoshells, there is a much more pronounced increase of the transition fields, which correspond to a change of vorticity, with lowering temperature [cp. Figs. 9(a) and 9(c) to Figs. 6(a) and 9(c)]. When comparing Fig. 9(a) to Fig. 6(a) one can also see that with increasing the nanoshell thickness the temperature range, where giant vortices are thermodynamically stable, extends towards lower temperatures. With decreasing the shell radius, this effect becomes quite pronounced even for relatively small values of [cp. Fig. 9(b) to Fig. 6(b)]. In sufficiently thick nanoshells, the temperature range, where multivortex states are thermodynamically stable, reduces to zero [see Fig. 9(c)], although in thin nanoshells of the same radius this range is relatively wide [see Fig. 6(b)]. Of course, when the temperature approaches , the phase boundaries become almost insensitive to the value of . Indeed, at the parameter always goes to zero [due to an increase of the penetration depth ], so that any nanoshell appears effectively thin.
V Conclusions
Curving a superconducting film into a spherical shell changes its vortexrelated properties drastically due to topological constraints. The hairysphere theorem Brouwer () is a straightforward example of such a constraint: it states that, in contrast to the situation on a flat film, there exists no nonvanishing continuous tangent vector field on the sphere. So, every nonvanishing supercurrent velocity field requires discontinuities, such as vortices. The interplay between the Lorentz force due to an applied field and the vortex superflow will force these vortices away from the equator (leaving an equatorial “Meissner band”) and towards the poles. This results in a ‘polar trapping potential’, which is nearly quadratic near the poles. When vortices conglomerate at the poles, they may coalesce to form giant or ringlike vortices, and these dynamics and phases are the topic of the present paper.
Three contributions to the energy should be kept in mind to interpret the phase diagrams obtained in our calculations. First, to create a vortex, the kinetic energy of the associated supercurrent (on the 2D spherical surface) should be taken into account. This contribution increases when two vortices with parallel vorticity are placed near each other, so it acts as a repulsion between the vortices. Thus, it tends to favor splitting of the giant vortices. Second, to create a vortex, the order parameter needs to be suppressed over a region typically of the size of the coherence length. The energy cost associated with this turns out to favor a multiply quantized (giant) vortex over the corresponding multivortex state. The energy cost is relatively larger for a smaller sphere, since proportionally a larger fraction of the total order parameter needs to be suppressed. The balance between these two energy contributions can be used to qualitatively understand the phase diagrams that we calculate for thin shells. Indeed, for magnetic fields corresponding to multiple quanta of vorticity, the smaller spheres will favor giant vortices, whereas the larger spheres favor the multivortex state. Note that this contribution to the energy strongly disfavors ringlike vortex states.
The third contribution to the energy is related to the gradients in the magnetic field. When the shell is much thinner than the penetration depth, the currents on the shell will not substantially perturb the applied field, and this contribution plays no role. However, for thicker shells, this contribution does become important – as can be seen from Fig. 7, the magnetic field is substantially perturbed. When a multivortex lattice is present, the magnetic field flux is concentrated near each vortex core, and shielded in between, leading to a larger magnetic contribution to the energy than for a giant vortex. Thus, for a thick shell, this contribution will favor the giant vortex state. This agrees with our phase diagram showing that the region, where the giant vortex is stable, is growing for thicker shells.
The temperature dependence of the phase diagrams was studied straightforwardly by taking temperature into account through the GinzburgLandau parameters. When multiple quanta of vorticity are present, we find that the giant vortex phase forms the preferred hightemperature phase. This offers the prospect of probing a temperaturedriven transition between a giant vortex and a multivortex state, alongside with a magneticfield driven transition. Moreover, the vortex dynamics are shown to be not sensitive to moderate imperfections in the shell; the energy contributions discussed here can overcome the pinning potential due to for example thickness inhomogeneities – such pinning potentials have in past experimental work hampered the detection of the giant vortex state. This robustness, together with the tunability of the phase diagram through a limited set of controllable parameters, makes superconducting nanoshells uniquely suited for the study of novel vortex states.
Acknowledgements.
This work was supported by the Fund for Scientific Research  Flanders projects Nos. G.0356.06, G.0115.06, G.0435.03, G.0306.00, the W.O.G. project WO.025.99N, the GOA BOF UA 2000 UA …References
 (1) M. W. Zwierlein, J. R. AboShaeer, A. Schirotzek, C. H. Schunck and W. Ketterle, Nature 435, 10471051 (23 June 2005).
 (2) V.V. Moshchalkov, L. Gielen, M. Dhallé, C. Van Haesendonck, Y. Bruynseraede, Nature 361, 617 (1993).
 (3) V. V. Moshchalkov, X. G. Qiu, and V. Bruyndoncx, Phys. Rev. B 55, 11793 (1997).
 (4) V.R. Misko, V.M. Fomin, J.T. Devreese, V.V. Moshchalkov, Phys. Rev. Lett. 90, 147003 (2003).
 (5) G. Stenuit, J. Govaerts, D. Bertrand and O. van der Aa, Physica C 332, 277 (2000); ibid. Phys. Lett. A 267, 56 (2000).
 (6) V. Bruyndoncx, J. G. Rodrigo, T. Puig, L. Van Look, and V. V. Moshchalkov, Phys. Rev. B 60, 4285  4292 (1999); R. Jonckheere et al., Phys. Rev. Lett. 85, 1528 (2000); ibid. Phys. Rev. Lett. 86, 1663 (2001); A. Kanda et al., Phys. Rev. Lett. 93, 257002 (2004).
 (7) V.R. Misko, V.M. Fomin and J.T. Devreese, Phys. Rev. B 64 (2001).
 (8) U.R. Fischer and G. Baym, Phys. Rev. Lett. 90, 140402 (2003); P. Engels, I. Coddington, P. C. Haljan, V. Schweikhard, and E. A. Cornell, Phys. Rev. Lett. 90, 170405 (2003).
 (9) R.D. Averitt, D. Sarkar, and N.J. Halas, Phys. Rev. Lett. 78, 4217 (1997); S.J. Oldenburg, R.D. Averitt, S.L. Westcott, and N.J. Halas, Chem. Phys. Lett. 288, 243 (1998).
 (10) C. R. Hu, R. S. Thomson, Phys. Rev. B 6, 110 (1972).
 (11) R. Kato, Y. Enomoto, S. Maekawa, Phys. Rev. B 44, 6916 (1991).
 (12) Hu Zhao, V. M. Fomin, J. T. Devreese, V. V. Moshchalkov, Solid State Commun. 125, 59 (2003).
 (13) Q. Du, L. Ju, J. Comp. Phys. 201, 511 (2004); ibid. Math. Comp. 74, 1257 (2004).
 (14) B. J. Baelus, D. Sun, F. M. Peeters, Phys. Rev. B 75, 174523 (2007).
 (15) M. Tinkham, Introduction to Superconductivity (2nd ed., McGrawHill, New York, 1996).
 (16) L. E. J. Brouwer, Mathematische Annalen 71, 97 (1912).