# Vortices in Ginzburg-Landau billiards

## I Introduction

This article is devoted to the study of the existence and stability of vortex solutions of the Ginzburg-Landau equations in finite two-dimensional domains with boundaries. For the case of the infinite plane, a large amount of work has been devoted to such questions, motivated by the ubiquitous character of the Ginzburg-Landau equations either in the study of superconductors or for the Abelian Higgs model. It is known  that for the infinite system the Ginzburg-Landau energy functional has a lower bound which is saturated for a special choice of the parameters (which corresponds to the limiting case between type I and type II superconductors). The minimizing equations at this special point, also called the dual point, are first order differential equations and they admit solutions with vortices. These different vortex solutions are classified by an integer of topological origin, namely the winding number of the order parameter which, roughly speaking, counts the number of vortices.

Although the Ginzburg-Landau equations in finite domains have been investigated , there is no generalization so far of the results concerning the dual point for finite two-dimensional domains (that we shall also call hereafter ‘billiards’). In this paper, we study such a generalization; besides a theoretical interest, our motivation has been triggered by a set of new experimental results obtained on small aluminium disks  in a mesoscopic regime where the radius $`R`$ is comparable with both the coherence length $`\xi `$ and the London penetration length $`\lambda .`$ There, the magnetization, as a function of the applied magnetic field, presents a series of jumps with an overall shape reminiscent of type-II superconductors. Such a behaviour is very different from that of macroscopic systems: in the large, aluminium is a genuine type-I superconductor. These experiments were first analyzed in the framework of the linearized Ginzburg-Landau equations . Although this approach explains some of the observed features, it fails to provide a satisfactory quantitative agreement, and can not account for the large spatial variations of the magnetic field inside the sample (e.g. vortices). Numerical solutions  of the Ginzburg-Landau equations give a much better description of this phenomenon, and emphasize the importance of the nonlinear term. However, in this context, no analytical results are available for the full Ginzburg-Landau equations. Using our results on the dual point, we shall derive an analytical expression for the free energy and for the magnetization of a mesoscopic superconductor as a function of the applied magnetic field.

The plan of this paper is as follows. In the section 2, we shall describe the main features and the known results about the dual point and the lower energy bound for an infinite system. Then, in section 3, we shall study the case of a finite domain and derive an expression of the free energy which, in addition to the previous lower bound, includes a contribution of the boundary. Section 4 contains a geometrical interpretation of these results. In section 5, our results are applied to analyze experimental data on small superconducting disks. In the conclusion, we propose some extensions of our work and a scenario for the nucleation of vortices at the boundary.

## II The dual point of Ginzburg-Landau equations in an infinite system

From now on, we study a superconducting billiard within the framework of the Ginzburg-Landau equations (this assumes that both the order parameter and the vector potential have a slow spatial variation). The expression of the Ginzburg-Landau energy density $`a`$ is given by

$$a=a_0+a_2|\psi |^2+a_4|\psi |^4+a_1|(\stackrel{}{}i\frac{2e}{\mathrm{}c}\stackrel{}{A})\psi |^2+\frac{B^2}{8\pi }$$
(1)   
where $`\psi =|\psi |e^{i\chi }`$ is the complex-valued order parameter, $`B`$ is the magnetic field and the $`a_i`$’s are real parameters. The coherence length and the London penetration length are related to these parameters as follows : $`\xi ^2=\frac{a_1}{|a_2|}`$ and $`\lambda ^2=\frac{\sqrt{2}}{4\pi }(\frac{\mathrm{}c}{2e})^2\frac{a_4}{a_1|a_2|}`$, so that we obtain for the dimensionless free energy $``$

$$=_\mathrm{\Omega }\frac{1}{2}|B|^2+\kappa ^2|1|\psi |^2|^2+|(\stackrel{}{}i\stackrel{}{A})\psi |^2$$
(2)   
where $`\psi `$ is measured in units of $`\psi _0=\sqrt{\frac{|a_2|}{2a_4}}`$, $`B`$ in units of $`\frac{\varphi _0}{4\pi \lambda ^2}`$, where $`\varphi _0=\frac{hc}{2e}`$ is the quantum of flux. The lengths are measured in units of $`\lambda \sqrt{2}`$ (the numerical factor $`\sqrt{2}`$ is for further convenience). The ratio $`\kappa =\frac{\lambda }{\xi }`$ is the only free parameter in (2) and it determines, in the limit of an infinite system, whether the sample is a type-I or type-II superconductor . The integral is over the two-dimensional domain $`\mathrm{\Omega }`$ of the superconducting sample.

The Ginzburg-Landau equations for the order parameter $`\psi `$ and for the magnetic field $`\stackrel{}{B}=\stackrel{}{}\times \stackrel{}{A}`$ are obtained from a variation of $``$:

$`(\stackrel{}{}i\stackrel{}{A})^2\psi `$ $`=`$ $`2\kappa ^2\psi (1|\psi |^2)`$ (3)   
$`\stackrel{}{}\times \stackrel{}{B}`$ $`=`$ $`2\stackrel{}{ȷ}`$ (4)   
where the current density $`\stackrel{}{ȷ}`$ is given by

$$\stackrel{}{ȷ}=Im(\psi ^{}\stackrel{}{}\psi )|\psi |^2\stackrel{}{A}$$
(5)   
The Ginzburg-Landau equations are nonlinear second order differential equations and their solutions are not known except for some particular cases. However, for the special value $`\kappa =\frac{1}{\sqrt{2}}`$, the equations for $`\psi `$ and $`\stackrel{}{A}`$ can be reduced to first order differential equations. This special point was first used by Sarma  in his discussion of type-I vs type-II superconductors and then identified by Bogomol’nyi  in the more general context of stability and integrability of classical solutions of some quantum field theories. We now review some of these properties of the Ginzburg-Landau free energy at the dual point. We shall use the following identity true for two dimensional systems

$$|(\stackrel{}{}i\stackrel{}{A})\psi |^2=|𝒟\psi |^2+\stackrel{}{}\times \stackrel{}{ȷ}+B|\psi |^2$$
(6)   
where $`\stackrel{}{ȷ}`$ is the current density and the operator $`𝒟`$ is defined as $`𝒟=_x+i_yi(A_x+iA_y)`$. At the dual point, $`\kappa =\frac{1}{\sqrt{2}},`$ the expression (2) for $``$ can be rewritten using the identity (6) as follows

$$=_\mathrm{\Omega }\frac{1}{2}|B1+|\psi |^2|^2+|𝒟\psi |^2+_\mathrm{\Omega }(\stackrel{}{ȷ}+\stackrel{}{A}).\stackrel{}{d}l$$
(7)   
where the last integral over the boundary $`\mathrm{\Omega }`$ of the system results from Stokes theorem.

For an infinite system we impose, as in , that the system is superconducting at large distance, i.e. $`|\psi |1`$ and $`\stackrel{}{ȷ}0`$ at infinity. The boundary term in (7) then becomes

$$_\mathrm{\Omega }(\stackrel{}{ȷ}+\stackrel{}{A}).\stackrel{}{d}l=_\mathrm{\Omega }(\frac{\stackrel{}{ȷ}}{|\psi |^2}+\stackrel{}{A}).\stackrel{}{d}l$$
(8)   
This last integral is known as the London fluxoid. It is quantized and using (5) one shows that it is equal to $`_\mathrm{\Omega }\stackrel{}{}\chi .\stackrel{}{d}l=2\pi n`$, where $`\chi `$ is the phase of the order parameter. The integer $`n`$ is the winding number of the order parameter $`\psi `$ and as such is a topological characteristic of the system. The extremal values of $``$, namely $`=2\pi n`$, are obtained when the bulk integral in (7) vanishes identically, i.e. when the two Bogomol’nyi  equations are satisfied:

$`𝒟\psi `$ $`=`$ $`0`$ (9)   
$`B`$ $`=`$ $`1|\psi |^2`$ (10)   
These two equations can be decoupled and one obtains that $`|\psi |`$ is a solution of the second order nonlinear equation

$$^2ln|\psi |^2=2(|\psi |^21)$$
(11)   
This equation is related to the Liouville equation. The set of equations (10,11) has been obtained without any assumption on the nature of the magnetic field and appears in various other situations, e.g. Higgs , Yang-Mills  and Chern-Simons  field theories. It was proven that these equations admit families of vortex solutions . For infinite systems, it can be shown that each vortex carries one flux quantum and that the winding number $`n`$ is equal to the number of vortices in the system. However, for an infinite system there is no mechanism to select the value of $`n`$, which only plays the role of a classifying parameter. It will be precisely the role of the boundary of a finite system to introduce such a selection mechanism and to determine $`n`$, according to the applied magnetic field.

## III The finite size system

From now on, we shall study finite size systems in an external magnetic field. The question then arises to know if such systems can sustain stable vortex solutions and how they behave as a function of the applied field. A simplified version, without applied magnetic field, was extensively studied by Bethuel et al. . In this work the mechanism for vortex creation is based on Dirichlet boundary conditions of the type $`\psi =f`$ on $`\mathrm{\Omega }`$ and where $`f`$ is a complex function of degree $`n`$. In the London limit, namely $`\kappa \mathrm{}`$, $`|\psi |`$ is 1 almost everywhere but, because of the degree $`n`$ imposed on the boundary, $`\psi `$ must vanish $`n`$ times in the bulk therefore leading to vortices. Moreover, numerical simulations of the Ginzburg-Landau equations for a long and thin parallepiped in a uniform magnetic field  show the existence of stationary vortex solutions whose number depends on the applied magnetic field. These simulations then indicate that the physical picture derived for $`\kappa =\frac{1}{\sqrt{2}}`$ remains valid for quite a large range of values of $`\kappa `$, and that the corresponding change of free energy is small (see also ). Indeed, for finite systems their size $`R`$ becomes relevant as a new length scale so that $`\kappa `$ is no more the only dimensionless parameter, and its value does not alone control the physics at the mesoscopic scale as it does for infinite systems. We shall study the case $`\kappa =\frac{1}{\sqrt{2}}`$, i.e. the dual point, and extend the previous approach to a system with finite size where boundary effects are important.

In a finite system, there are in general non-zero edge currents and the order parameter is not equal to 1 on the boundary. Hence, the identification of the boundary integral in (7) with the fluxoid (8) is not possible anymore, and the free energy can not be minimized just by imposing Bogomol’nyi equations (10). However, the currents on the boundary of the system screen the external magnetic field and therefore produce a magnetic moment (a circulation) opposite to the direction of the field, whereas vortices in the bulk of the system produce a magnetic moment along the direction of the applied field. Hence currents in the bulk circulate in a direction opposite to those at the boundary. If one assumes cylindrical symmetry, the current density $`\stackrel{}{ȷ}`$ has only an azimuthal component, with opposite signs in the bulk and on the edge of the system (the radial component of $`\stackrel{}{ȷ}`$ is zero since $`\stackrel{}{ȷ}`$ is divergenceless). Thus, there exists a circle $`\mathrm{\Gamma }`$ on which $`\stackrel{}{ȷ}`$ vanishes . This allows us to separate the domain $`\mathrm{\Omega }`$ into two concentric subdomains $`\mathrm{\Omega }=\mathrm{\Omega }_1\mathrm{\Omega }_2`$ such that the boundary $`\mathrm{\Omega }_1`$ is the curve $`\mathrm{\Gamma }`$ (see Fig.1). On $`\mathrm{\Omega }_1`$, the current density $`\stackrel{}{ȷ}`$ is zero, therefore one has

$$_{\mathrm{\Omega }_1}\stackrel{}{ȷ}.\stackrel{}{d}l=_{\mathrm{\Omega }_1}\frac{\stackrel{}{ȷ}}{|\psi |^2}.\stackrel{}{d}l=0$$
(12)   
Thus one deduces as above that Bogomol’nyi and Liouville equations are valid in the finite domain $`\mathrm{\Omega }_1`$ as in the case of the infinite plane. The existence of vortices in a finite domain such as $`\mathrm{\Omega }_1`$ was checked using a numerical solution  of (11): assuming cylindrical symmetry, one shows that $`|\psi |`$ vanishes as a power law at the center of the disk, hence there is a vortex in the center (more precisely a multi-vortex whose multiplicity is determined by the exponent of the power law); moreover, $`|\psi |`$ saturates very rapidly to a constant value close to one for lengths larger than $`\lambda `$. The same conclusion can be reached by defining $`f(r)=ln|\psi |^2`$ and linearizing (11) around $`|\psi |=1`$. Then, $`f`$ satisfies $`^2f=2f`$ whose general solution is $`f(r)=aI_0(r\sqrt{2})+bK_0(r\sqrt{2})`$. From the behaviour of the Bessel functions $`I_0`$ and $`K_0`$, one obtains that for small $`r`$, $`|\psi |`$ vanishes as a power law and saturates rapidly to a constant of order one in a finite range of $`r`$’s (in units of $`\lambda \sqrt{2}`$).

The magnetic flux $`\mathrm{\Phi }(\mathrm{\Omega }_1)`$ through $`\mathrm{\Omega }_1`$ is calculated, in units of the flux quantum $`\varphi _0`$, using the fluxoid and (12) so that

$$\mathrm{\Phi }(\mathrm{\Omega }_1)=n_{\mathrm{\Omega }_1}\frac{\stackrel{}{ȷ}.\stackrel{}{d}l}{|\psi |^2}=n$$

As before $`n`$ is the winding number, i.e. $`_{\mathrm{\Omega }_1}\stackrel{}{}\chi .\stackrel{}{d}l=2\pi n`$, as well as the number of vortices. The free energy in $`\mathrm{\Omega }_1`$ is

$$(\mathrm{\Omega }_1)=2\pi n.$$
(13)   
We consider now the contribution of $`\mathrm{\Omega }_2`$ to the free energy. It is given by (2) and can be rewritten using the phase and the modulus of the order parameter $`\psi `$, as

$$(\mathrm{\Omega }_2)=_{\mathrm{\Omega }_2}(|\psi |)^2+|\psi |^2|\stackrel{}{}\chi \stackrel{}{A}|^2+\frac{B^2}{2}+\frac{(1|\psi |^2)^2}{2}$$
(14)   
We know, from the London equation, that both the magnetic field and the vector potential decrease rapidly (exponentially) away from the boundary $`\mathrm{\Omega }`$ of the system over a distance of order 1 in units of $`\lambda \sqrt{2}`$. Over the same distance, at the dual point, $`|\psi |`$ saturates to unity. One can thus estimate the integral (14) using an elementary version of the saddle-point method. We assume cylindrical symmetry, and we neglect the term $`(|\psi |)^2`$ on the boundary because of the covariant Neumann boundary conditions at the interface between a superconductor and an insulator . We obtain for the free energy the expression

$$(\mathrm{\Omega }_2)_\mathrm{\Omega }|\psi |^2|\stackrel{}{}\chi \stackrel{}{A}|^2+\frac{B^2}{2}+\frac{(1|\psi |^2)^2}{2}$$
(15)   
where the integral is now over the boundary of the system. To go further, we need to implement boundary conditions for the magnetic field $`B(R)`$ and the vector potential $`A(R)`$. The choice $`B(R)=B_e`$, where $`B_e`$ is the external imposed field, corresponds to the geometry of an infinitely long cylinder, where the flux lines are not distorted outside the system. This boundary condition is not adapted to describe a flat thin disk. A more suitable choice is provided by demanding $`\varphi =\varphi _e,`$ which means that the total flux through the disk is identical to the flux of the external field, although flux lines are distorted by the superconducting sample. One can check, using numerical simulations , that this condition is well satisfied if the disk is thin enough. The boundary condition $`\varphi =\varphi _e`$ implies that the vector potential is identified by continuity to its external applied value $`\stackrel{}{A_e}`$ which has only the azimuthal component $`\frac{\varphi _e}{2\pi R}`$. It should be noticed that the magnetic field $`\stackrel{}{B}`$ has a non monotonous variation: it is low in the bulk, larger than $`B_e`$ near the edge of the system, because of the distortion of flux lines, and eventually equal to its applied value far outside the system .

Different choices of boundary conditions will give rise to different limits for a very large system (i.e. for $`R\mathrm{}`$). The limit of an infinitely long cylinder corresponds to a superconducting bulk sample, whereas the limit of a thin disk, which is the case we consider here, corresponds to a superconducting thin sheet.

The formula (15) for the free energy is similar to the Little-Parks expression  for a quasi one-dimensional hollow cylinder in a uniform applied magnetic field. The minimization with respect to $`|\psi |`$ gives $`1|\psi |^2=|\stackrel{}{}\chi \stackrel{}{A}|^2`$, such that (15) can be written as

$$(\mathrm{\Omega }_2)=_\mathrm{\Omega }|\stackrel{}{}\chi \stackrel{}{A_e}|^2+\frac{1}{2}B^2\frac{1}{2}|\stackrel{}{}\chi \stackrel{}{A_e}|^4$$
(16)   
Performing the integral over the boundary of the system, we obtain

$$\frac{1}{2\pi }(\mathrm{\Omega }_2)=\frac{\lambda \sqrt{2}}{R}(n\varphi _e)^2\frac{1}{2}(\frac{\lambda \sqrt{2}}{R})^3(n\varphi _e)^4$$
(17)   
We have neglected the contribution of the $`B^2`$ term, which is similar to the first term in the r.h.s. of (17) but smaller by a factor of the order $`(\lambda /R)^2`$. The integer $`n`$ which appears in (17) is the same as in (13), since the order parameter $`\psi `$ is the same function in both subdomains. The circulation of its phase $`\chi `$ (the winding number) counts the number of zeroes in the domain $`\mathrm{\Omega }`$, i.e. the number of vortices. The thermodynamic Gibbs potential $`𝒢`$ of the system is obtained from $`(\mathrm{\Omega }_1)+(\mathrm{\Omega }_2)`$ by a Legendre transformation so that

$$\frac{1}{2\pi }𝒢(n,\varphi _e)=n+\frac{\lambda \sqrt{2}}{R}(n\varphi _e)^2\frac{1}{2}(\frac{\lambda \sqrt{2}}{R})^3(n\varphi _e)^4\frac{2\lambda ^2}{R^2}\varphi _{e}^{}{}_{}{}^{2}$$
(18)   
This relation consists in a set of quartic functions indexed by the integer $`n`$. The minimum of the Gibbs potential is the envelop curve defined by the equation $`\frac{𝒢}{n}|_{\varphi _e}=0`$, i.e. the system chooses its winding number $`n`$ in order to minimize $`𝒢`$. This provides a relation between the number $`n`$ of vortices in the system and the applied magnetic field $`\varphi _e`$. In the limit of a large enough $`\frac{R}{\lambda }`$, the quartic term is negligible and the Gibbs potential reduces to a set of parabolas (Fig.2). The winding number $`n`$ is then given by the integer part

$$n=[\varphi _e\frac{R}{2\sqrt{2}\lambda }+\frac{1}{2}]$$
(19)   
The magnetization $`M=\frac{𝒢}{\varphi _e},`$ of the system, is given by

$$M=\frac{2\sqrt{2}\lambda }{R}(\varphi _en)\frac{4\lambda ^2}{R^2}\varphi _e$$
(20)   
For $`\varphi _e`$ smaller that $`\frac{R}{2\sqrt{2}\lambda }`$, we have $`n=0`$ and $`(M)`$ increases linearly with the external flux. This corresponds to the London regime before the first vortex enters the system. The field $`H_1`$ at which the first vortex enters the system corresponds to $`𝒢(n=0)=𝒢(n=1)`$, i.e. to

$$H_1=\frac{\varphi _0}{2\pi \sqrt{2}R\lambda }+\frac{\varphi _0}{2\pi R^2}$$
(21)   
The subsequent vortices enter one by one for each crossing $`𝒢(n+1)=𝒢(n)`$; this happens periodically in the applied field, with a period equal to

$$\mathrm{\Delta }H=\frac{\varphi _0}{\pi R^2}$$
(22)   
This gives rise to a discontinuity of the magnetization $`\mathrm{\Delta }M=\frac{2\sqrt{2}\lambda }{R}.`$

There is a qualitative similarity between the results we derived using the Bogomol’nyi equations within the domain $`\mathrm{\Omega }_1`$ and those obtained from a linearised version of the Ginzburg-Landau functional . But the two approaches differ in their quantitative predictions due to the importance of the non linear term. An illustration of this is given in the next section.

## IV A geometrical interpretation

In an infinite system, one shows using the boundary conditions $`|\psi |1`$ and $`\stackrel{}{ȷ}0`$ at infinity that equation (8) implies the quantization of the magnetic flux:

$$_{R^2}\stackrel{}{B}.d\stackrel{}{S}=n$$
(23)   
If one interprets $`B`$ as a curvature, this relation is analogous to the Gauss-Bonnet theorem, which states that, for a compact manifold $`M`$ without boundary, the integral of the Gaussian curvature $`K`$ over the surface is equal to the Euler-Poincaré characterics $`\chi `$ of the manifold, which is a topological invariant integer:

$$_MK=\chi $$
(24)   
This result is in fact more than an analogy and at the dual point the Ginzburg-Landau functional has a useful geometrical interpretation that we shall now highlight .

The Ginzburg-Landau functional (2) for the energy of a superconductor corresponds to a U(1) gauge symmetry. For that symmetry, one can identify an abelian one-form connexion and a two-form curvature $`\mathrm{\Omega }`$ respectively given by the vector potential $`\stackrel{}{A}`$ and by the magnetic field $`\stackrel{}{B}`$. The complex order parameter $`\psi `$ is a section of the U(1) fiber and the basis manifold is the infinite plane. The above boundary conditions allow to map the plan onto a compact manifold without boundary, namely the sphere. Topological invariants of the problem are obtained from the Chern-Weil invariant polynomial $`P(\mathrm{\Omega })=det(1+\frac{i}{2\pi }\mathrm{\Omega })`$ (see e.g. ) For the U(1) bundle over a sphere there is only one Chern class, $`\mathrm{\Omega }(=B)`$. The integral of that Chern class over the basis manifold is a topological invariant integer called a Chern number and is precisely given by (23). This Chern number plays a role similar to $`\chi `$. At the dual point, the energy is $`=B=2\pi n;`$ this fact can be translated in geometrical terms by saying that the extremal free energy is identical to a topological invariant of the problem namely its Chern number.

When the basis manifold $`M`$ has a boundary $`M`$ which is not a geodesic, the integral of its curvature is neither an integer nor a topological invariant. The Gauss-Bonnet theorem is then generalized by adding a boundary term so that the Euler-Poincaré characteristics $`\chi `$ is now given by the relation

$$\chi =\frac{1}{2\pi }_MK𝑑S+\frac{1}{2\pi }_Mk_g𝑑l$$
(25)   
where $`K`$ and $`k_g`$ are respectively the Gaussian curvature of the manifold and the geodesic curvature of the boundary. A similar result holds for the U(1) problem in a bounded domain (or in an infinite domain with boundary conditions different from those chosen above). In that case the relation equivalent to (25) is given by the fluxoid relation (8):

$$n=\frac{1}{2\pi }_\mathrm{\Omega }\stackrel{}{B}.d\stackrel{}{S}+_\mathrm{\Omega }\frac{\stackrel{}{ȷ}}{|\psi |^2}.d\stackrel{}{l}$$
(26)   
As before, $`B`$ is the curvature of the connexion and the current density $`\frac{\stackrel{}{ȷ}}{|\psi |^2}`$ plays here the role of a geodesic curvature . The expression obtained in (17) for the Bogomol’nyi free energy of a system with a boundary can be rewritten as:

$$=_\mathrm{\Omega }B+_\mathrm{\Omega }\eta (\frac{\stackrel{}{ȷ}}{|\psi |^2})K+\eta (k_g)$$
(27)   
The boundary correction is a function $`\eta `$ of the geodesic curvature. The results obtained in the preceeding section show that the geodesic curvature is given by $`n\varphi _e`$ for a cylindrically symmetric system and using the appropriate approximations we determined the function $`\eta `$ as being an even fourth order polynomial in the geodesic curvature.

This geometric interpretation makes us believe that an expression such as (27) is fairly general. It could be well suited, as an Ansatz, to describe finite systems which are known to have a topological description in the infinite limit (for instance, a suitable generalisation of (27) to other symmetry group like SU(2) could describe some phases of superfluids in a bounded domain).

## V comparison with the experimental results

In order to compare our results to the experimental data given in , we shall consider the limit where the radius $`R`$ is larger than $`\lambda `$ and $`\xi `$, (typically, $`R10\lambda `$ is considered experimentally ). The thickness of the system is assumed to be small enough so that we can neglect variations of both the magnetic field and the order parameter along the cylinder axis. We shall prove that within these approximations, the expression (18) captures the main features observed experimentally i.e. the behaviour of the magnetization at low fields (before the first discontinuity), the periodicity and the linear behaviour between the successive jumps and provides a fairly good quantitative agreement.

The magnetization curve in fig.3 (plotted for the ratio $`\frac{\lambda }{R}=0.14`$, chosen arbitrarily) agrees qualitatively with both the numerical and experimental curves of fig.3 in . Besides, taking the experimental parameters of , namely $`R=1.2\mu m`$ and $`\lambda (T)=84nm`$ at $`T=0.4K`$, we compute from our expressions $`H_1=25G`$ and $`\mathrm{\Delta }H=4.6G.`$ These values agree with the results of  to within a few percent. We emphasize that $`H_1`$ scales like $`\frac{1}{R}`$, whereas $`\mathrm{\Delta }H`$ scales like $`\frac{1}{R^2}`$ in accordance with . We calculate the ratio of the magnetization jumps to the maximum value of $`M`$ to be 0.20 as compared to a measured value of 0.22. The total number of jumps scales like $`R^2`$ and the upper critical field is independent of $`R`$ in our theory. This to agrees fairly well with the experimental data . We emphasize that this quantitative agreement cannot be obtained from the linearized Ginzburg-Landau equations.

## VI Conclusion

We have investigated the problems of the existence and the stability of vortices in a two-dimensional bounded superconducting system. To that purpose, we have used the Ginzburg-Landau energy functional at the special dual point characterized by the value $`\kappa =\frac{1}{\sqrt{2}}`$ which corresponds, for an infinite system, to a superconductor between type I and type II. We have shown for that case that it is still possible to obtain vortex solutions at the thermodynamic equilibrium. In contrast to the Bogomol’nyi solution obtained for the infinite plane, where there is no definite value for the number $`n`$ of vortices, there is a selection mechanism for a finite billiard in an external magnetic field that allows to compute the number of vortices as a function of the applied field. Our reasoning is based on the existence of a contour $`\mathrm{\Gamma }`$ which allows to separate the system in two parts, the bulk containing the vortices and the edge with screening currents. As such the system might be viewed as a kind of two-dimensional Josephson junction or weak superfluid link . Vortices enter (or are expelled from) the superconductor through $`\mathrm{\Gamma }`$, if the applied magnetic field is increased (or decreased).

Although we considered the special case of the dual point, our analysis provides a satisfactory quantitative description of the behaviour observed experimentally on such small superconducting billiards at least for the regime where the density of vortices is low enough i.e. for small applied magnetic fields. For higher fields, the expression (18) does not describe properly the tail of the magnetization curve where both the periodicity $`\mathrm{\Delta }H`$ and the amplitude $`\mathrm{\Delta }M`$ of the jumps decrease and eventually vanish. In this regime, the vortices interact both between themselves and with the edge currents. This regime can be studied by a perturbative analysis around the dual point .

So far, we have studied only equilibrium states. On the basis of numerical simulations of the time-dependent Ginzburg-Landau equations , we propose a mechanism for vortex nucleation. At each discontinuity of the winding number, namely when $`\varphi _e\frac{R}{2\sqrt{2}\lambda }`$ is a half integer, there is a nodal line, joining the center of the system to its boundary, along which the order parameter $`\psi `$ vanishes and where its phase is ill defined. This might be interpreted as an opening of the ring of the screening currents which allows a flux line of the external field to enter the system. In this case, we expect the contour $`\mathrm{\Gamma }`$ to coincide with the boundary of the system. The existence of such a nodal line has been discussed in the context of Ginzburg-Landau equations and for the related Aharonov-Bohm problem of a magnetic flux line piercing either the infinite plane  or a finite domain .

Acknowledgment It is our pleasure to thank A.Ori and T.N.Tomaras for useful and enlightening discussions. This work is supported in part by a grant from the Israel Academy of Sciences and by the fund for promotion of Research at the Technion. K.M. acknowledges support by the Lady Davies foundation. E.A. ackowledges the very kind hospitality of the Laboratoire de Physique des Solides and the LPTMS at the university of Paris (Orsay).