Jun 21, 2016 - that of the normal state we call it EVS. Consider the left vortex state immediately after the applied field is switched off. Vortices h...

0 downloads 10 Views 1MB Size

Departement Fysica, Universiteit Antwerpen,

arXiv:1606.06400v1 [cond-mat.supr-con] 21 Jun 2016

Groenenborgerlaan 171, B-2020 Antwerpen, Belgium 2

Instituto de F´ısica, Universidade Federal do Rio de Janeiro, 21941-972 Rio de Janeiro, Brazil∗

3

Dipartimento di Fisica, Universit`a di Camerino, I-62032 Camerino, Italy 4

Departamento de F´ısica, Universidade Federal Rural de Pernambuco, 52171-900 Recife, Pernambuco, Brazil

Abstract We consider excited vortex states, which are vortex states left inside a superconductor once the external applied magnetic field is switched off and whose energy is lower than of the normal state. We show that this state is paramagnetic and develop here a general method to obtain its Gibbs free energy through conformal mapping. The solution for any number of vortices in any cross section geometry can be read off from the Schwarz - Christoffel mapping. The method is based on the first order equations used by A. Abrikosov to discover vortices. PACS numbers: 74.20.De, 74.25.Uv, 75.20.-g

1

INTRODUCTION

Individual vortices in type II superconductors were first seen by U. Essmann and H. Trauble through the Bitter decoration technique [1, 2]. Since then several other techniques have been developed for this purpose [3], such as scanning SQUID microscopy [4], high-resolution magneto-optical imaging [5, 6], muon spin rotation (µSR) [7, 8], scanning tunneling microscopy [9, 10], and magnetic force microscopy [9]. These advancements in the visualization of individual vortices opens the gate to investigate new properties [11, 12], such as those of the excited vortex state (EVS). The EVS brings a new paradigm to the study of the transient dynamics of vortices in superconductors with boundaries, a subject of current interest due to the onset of instabilities [13].

A type II superconductor in presence of an external applied magnetic field contains vortices in its interior whose density is fixed by the external applied magnetic field. Once the applied field is switched off this state becomes unstable and vortices must leave the superconductor. However their exit can be hindered by microscopic inhomogeneities which pin them inside the superconductor. Here we make an important distinction concerning this left vortex state according on how its energy compares with that of the normal state. Although the left vortex state is always unstable, only in case its energy is lower than that of the normal state we call it EVS. Consider the left vortex state immediately after the applied field is switched off. Vortices have topologically stability but only inside the superconducting state. Thus as long as the superconducting state exists they cannot be created or destroyed and therefore can only enter or exit through the boundary. However if the superconducting state ceases to exist vortices vanish all together inside the material. Therefore the obtainment of the Gibbs free energy difference between the superconducting state, which contains the left vortex state inside, and the normal state must be done to determine whether the superconducting state still exists. In case this energy difference is positive one expects the collapse of the superconducting state and the onset of the normal state, possibly because of non-linear thermal effects that will develop and pinning will not be able to uphold this transition. In case this energy difference is negative the EVS exists and can remain inside the superconductor provided that pinning impedes the motion of vortices towards the superconductor boundary. 2

In the London limit the energy of the vortex state is clearly positive since it is the sum of all two-vortex repulsive interactions. This is equivalent to say that the energy is the sum of the vortex lines self energies plus the sum over the repulsive interaction between different vortices. Once again it becomes clear that vortices inside the superconductor can only be sustained by the external pressure exerted by the applied field and once this pressure ceases to exist they go away leaving behind the state of no vortices. However the study of the EVS cannot be done in this context of the London theory since this theory does not describe the condensate energy. Within the London theory it is not possible to compare the energy of the vortex state with the normal state. Therefore the EVS must be studied in the context of the Ginzburg-Landau (GL) theory.

The EVS is intrinsically paramagnetic [14–16], which means that its magnetization points in the same direction of the switched off applied field. A simple way to understand this property is to analyse the two types of currents in a superconductor. In the Meissner phase there are only shielding currents at the boundary, which are diamagnetic, namely, they create an equal and opposite field to the applied one in its interior resulting into a net null field. However the current around the vortex core is opposite to the shielding current, a property that explains why the increase in the number of vortices weakens the diamagnetic response set by the shielding currents. This increase reaches the critical point where the magnetization vanishes and then the upper critical field is reached. Consider, for the sake of the argument, that after the sudden switch off the applied field the shielding currents immediate disappear. This must occur in some moment otherwise there will be a left field of opposite direction to the applied one. The vortices caught in this sudden transformation either remain pinned or start their move towards the boundary. In case the Gibbs energy difference is positive the whole superconducting state will move towards its collapse, but if negative there will be a EVS left inside, which corresponds to a vortex state without the shielding currents and therefore with a purely paramagnetic response. There is no Lorentz force to push vortices at the boundary since the shielding currents are absent. Vortices have no impediments at the boundary, there is no Bean-Livingstone barrier in zero applied field [17, 18].

3

In this paper we propose a general method to calculate the Gibbs free energy difference between the left vortex and the normal states valid near the superconducting to normal transition, where the order parameter is small. This method is general since it applies to any vortex configuration in any geometry of the cross section area. The method applies to very long superconducting cylinders such that the vortex lines are parallel to each other and oriented along its major axis. As shown here, the order parameter that describes the left vortex state can be obtained from the well-known mathematical problem of conformal mapping. Remarkably the order parameter is just an analytical function with constant modulus at the boundary of a given cross section geometry. In power of the order parameter we obtain the Gibbs energy difference, the local magnetic field and currents inside the superconductor. Besides we also obtain other interesting features, such as the magnetic field at the center of the vortex core, and also the paramagnetic magnetization. It is well-known that at low vortex density the magnetic field inside the vortex core is twice the value of the lowest critical field [19]. For higher densities a vortex lattice sets in and the magnetic field inside the vortex core varies according to the parameter κ (the ratio between the penetration and the coherence lenghts) and the vortex density [20]. Here we report the surprising result that as vortices move towards the boundary the magnetic field at their cores, and also the paramagnetic magnetization, change according to their positions. The magnetization vanishes when vortices reach the boundary. The present approach stems directly from A. Abrikosov’s seminal work [21] that led to the discovery of vortices. There he found two identities that we refer here as the first order equations. They were later rediscovered by E. Bogomolny [22] in the context of string theory and found to solve exactly the second order variational equations, that stem from the free √ energy, for a particular value of the coupling (κ = 1/ 2). These first order equations are able to determine the order parameter and the local magnetic field. A. Abrikosov used them just in case of bulk superconductor, which means a superconductor with no boundary to a non superconducting region. Therefore he took the assumption of a lattice such that periodic boundary conditions apply. Obviously the bulk is an idealized system that simplifies the theory but hinders important boundary effects. Here we essentially extend this very same treatment to the case of a superconductor with a boundary. Interestingly in case of the bulk, and of no applied external field, the first order equations predict no vortex solution. The only possible solution is that of a constant order parameter. However, as shown here, the exis4

tence of a boundary changes dramatically the scenario and vortex solutions become possible even without the presence of an applied field. Thus the profound connection between the mathematical theory of conformal mapping and the vortex solutions relies on the existence of boundaries. The theory of conformal mapping was developed in the nineteenth century. The Riemann mapping theorem of 1851 states that any simply connected region in the complex plane can be conformally mapped onto any other, provided that neither is the entire plane. The mapping useful to us is the one that takes any finite geometry into the disk [23] which is essentially a variant of the Schwarz Christoffel conformal transformation [24] of the upper half-plane onto the interior of a simple polygon. The Schwarz Christoffel mapping [25] is used in potential theory and among its applications are minimal surfaces and fluid dynamics.

Although we are mainly interested here in the EVS problem, it must be stressed that the present method also applies in presence of an applied external field. The only condition imposed is to be near to the superconducting to normal transition, i.e., near to the upper critical field line. The zero field case just corresponds to particular case near to the the critical temperature, Tc . However in the general case the connection between the order parameter solution and the conformal mapping only holds in case of a circular cross section. Recently a new method to solve the linearized GL problem for mesoscopic superconductors was proposed by means of conformal mapping [26]. Thus our approach is distinct from this one since we are also able to obtain the local magnetic field from the order parameter whereas the above method is not. We find useful to summarize this paper as follows. We start describing the GL theory in section and in section introduce the key element to our approach, which is the decomposition of the kinetic energy into a sum of three terms, namely, the ground state condition plus the magnetic interaction plus the boundary term. This decomposition is directly related to the first order equations. In subsection we discuss the role played by the boundary conditions in the context of the first order equations and obtain expressions for the magnetization and the Gibbs free energy. Dimensionless units are introduced in section , which are useful in the treatment of two examples of conformal mapping considered here. There the limit of weak order parameter and small field is shown to help in solving the first order equations. Until this point the proposed method is very general and applies in case of an applied field as well. Only in section the external field is switched off and the connection between the 5

first order equations and the conformal mapping established. In order to show the power of the present method we study two examples in section , whose solutions are obtained analytically. Both have the cross section geometry of a disk, the first one with a vorticity L at the center, subsection , and the other with vorticity one in any position is discussed in subsection . We reach conclusions in section and in the appendices and details of our analytical calculations are given.

THE GIBBS FREE ENERGY OF A LONG SUPERCONDUCTOR AND THE VARIATIONAL SECOND ORDER EQUATIONS

Effects of the top and the bottom of a very long superconductor are neglected such that symmetry along the major axis is assumed. Any cross sectional cut at a given x3 plane reveals the same area Σ and the same physical properties. The external constant magnetic field ~ = Ha xˆ3 . Hence the order parameter is only expressed is oriented along the major axis, H by coordinates in this plane, ψ(x1 , x2 ), and the only one non-zero component of the local magnetic field is perpendicular to this plane, h3 (x1 , x2 ) = ∂1 A2 (x1 , x2 ) − ∂2 A1 (x1 , x2 ). The difference between the superconducting and normal Gibbs free energy densities, defined as ∆G, is given by, ∆G ≡ Gs − Gn =

Z Σ

β d2 x α(T )|ψ|2 + |ψ|4 Σ 2 (

~ 2 (h3 − Ha )2 |Dψ| + + . 2m 8π

(1)

The normal state, Gs = Gn , is reached for ψ = 0 and h3 = Ha . The well-known GL parameters have the properties that β > 0 and that,

α(T ) = α0

< 0 for T < Tc T − 1 ⇒ α(T ) is , Tc > 0 for T > Tc

~ = where α0 is a positive constant. The vector notation is two-dimensional, such as D D1 xˆ1 + D2 xˆ2 with Dj ≡ (¯ h/i) ∂j − q Aj (x1 , x2 )/c, j = 1, 2. The well-known second order ~ and δψ ∗ , equations are obtained by doing variations with respect to the fields, namely, δ A and lead to the non-linear GL equation, ~ 2ψ D + αψ + β|ψ|2 ψ = 0, 2m 6

(2)

and to Amp`ere’s law, ~ × ~h = 4π J, ~ ∇ c

(3)

q ∗~ ψ Dψ + c.c. . J~ = 2m

(4)

where the current density is given by,

Boundary conditions must be added to find the physical solutions. They correspond to no current flowing out of the superconductor, as vacuum is assumed outside, and the local field inside must be equal to the applied field outside. Let Υ be the perimeter of area Σ, thus at the boundary it must hold that, ~ n ˆ · J| ~ x at Υ = 0, and,

(5)

h3 |~x at Υ = Ha ,

(6)

where n ˆ is perpendicular to the boundary. Since the current is given by Eq.(4) the following condition on the derivative of the order parameter, ~ n ˆ · Dψ| ~ x at Υ = 0,

(7)

is enough to guarantee the condition of Eq.(5).

DUAL VIEW OF THE KINETIC ENERGY AND THE FIRST ORDER EQUATIONS

The kinetic energy density admits a dual formulation due to the mathematical identity proven in the appendix, given by Eq.(90), ~ 2 |D+ ψ|2 h ¯q h ¯ |Dψ| = + h3 |ψ|2 + (∂1 J2 − ∂2 J1 ) , 2m 2m 2mc 2q

(8)

where D+ ≡ D1 + i D2 and the current is given by Eq.(4). This decomposition of the kinetic energy as a sum of three terms is exact, and its derivation [27] is given in the appendix . Therefore the kinetic energy density, Fk =

Z Σ

~ 2 d2 x |Dψ| , Σ 2m 7

(9)

is also given by, Fk =

Z Σ

d2 x Σ

|D+ ψ|2 h ¯q h ¯ 1I ~ ~ + h3 |ψ|2 + dl · J. 2m 2mc 2q Σ Υ !

(10) In this dual view there is a superficial (perimetric) contribution due to the current. For a bulk superconductor, where Σ → ∞, the superficial current vanishes either because the currents are localized within some region in the bulk away from the boundary or because of periodic boundary conditions, the latter case being that one considered by Abrikosov [21]. However in case of a finite area, such as for a mesoscopic superconductor, the current at the boundary must be considered and does not vanish. The most interesting property of this dual view of the kinetic energy is that the current also acquires a new formulation. The current can be simply obtained by the variation of the kinetic energy with respect to the vector potential, 1 Z d2 x ~ ~ J · δ A, δFk = − c Σ Σ

(11)

and use of the dual formulation leads to q h ¯q [(D+ ψ)∗ ψ + ψ ∗ (D+ ψ)] − ∂2 |ψ|2 , 2m 2m q h ¯q J2 = i [(D+ ψ)∗ ψ − ψ ∗ (D+ ψ)] + ∂1 |ψ|2 . 2m 2m J1 =

(12) (13)

The superficial current term does not contribute since at the edge the variation of the vector ~ = 0, for ~x at Υ. potential vanishes, namely, δ A In this paper we seek the minimum of the free energy through solutions of the first order equations, given below, D+ ψ = 0, and, h ¯q 2 h3 = H 0 − 2π |ψ| , mc

(14) (15)

instead of solutions of the second order equations, given by Eqs.(2) and (3). On Eq.(15), H 0 is a constant to be determined but whose interpretation is very clear. It is the local field at the vortex core since there ψ = 0, and so, h3 = H 0 . In the absence of vortices the order parameter is constant and H 0 = 2π(¯ hq/mc)|ψ|2 since it must hold that h3 = 0 everywhere. We show that the above equations provide an exact solution to Amp`ere’s law 8

and an approximate solution to the non-linear GL equation, and so, provided an easy and efficient method to search for the minimum of the free energy. Amp`ere’s law is exactly solved and to see it, just take the condition of Eq.(14) into the current, as given by Eqs.(12) and (13). The Amp`ere’s law, given by ∂2 h3 = 4πJ1 /c and ∂1 h3 = −4πJ2 /c, becomes ∂2 h3 = −(4π¯hq/2mc)∂2 |ψ|2 and ∂1 h3 = −(4π¯hq/2mc)∂1 |ψ|2 , respectively, since Eq.(15) holds. Then one obtains that, h ¯q ~ 2. xˆ3 × ∂|ψ| J~ = 2m

(16)

The surface term, contained in the dual formulation of the kinetic energy, is given by Eq.(8). h ¯ h ¯2 ~ 2 2 (∂1 J2 − ∂2 J1 ) = + ∂ |ψ| . q 2m

(17)

The mathematical identity of Eq.(8), becomes, ~ 2 hq 2 h h ¯2 ~ 2 2 |Dψ| ¯q = H0 − |ψ| |ψ|2 + ∂ |ψ| , 2m mc 2mc 4m !

(18)

once assumed that the first order equations are satisfied. The non-linear GL equation is approximately solved in the sense that its integrated version is exactly solved. This integrated version is obtained by multiplying the non-linear GL equation with ψ ∗ , and next integrating over the entire area Σ of the superconductor:

Z Σ

~ 2ψ d2 x ∗ D ψ + α(T )|ψ|2 + β|ψ|4 = 0. Σ 2m

(19)

Next we transform this equation into an algebraic one whose usefulness is to fix the scale of the order parameter which has remained undefined when solving the scale invariant Eq.(14). The first term of the integrated equation summed with its complex conjugate and divided by 2, can be expanded as follows,

∗

~ 2ψ ~ 2ψ ~ 2 1 ∗D D h ¯ 2 ~ 2 2 |Dψ| ψ = − ψ + ∂ |ψ| + . 2 2m 2m 4m 2m (20) Inserting Eq.(20) into the integrated equation, upon summing with its complex conjugate and dividing by 2, one obtains that,

Z Σ

~ 2 d2 x h ¯ 2 ~ 2 2 |Dψ| ∂ |ψ| + − + α(T )|ψ|2 + β|ψ|4 = 0. Σ 4m 2m

(21) 9

The use of Eq.(18) turns the integrated equation into the following one.

Z Σ

¯q d2 x H 0 h + α |ψ|2 − π Σ 2 mc !

h ¯q mc

!2

− β |ψ|4

= 0.

(22) In summary Eq.(14), together with Eq.(22), fully determines the order parameter, ψ, and from Eq.(15) one obtains the local magnetic field, h3 .

The boundary conditions, the magnetization and the Gibbs free energy

There should be no current flowing out of the superconductor and the magnetic field must be continuous at the boundary. Here we address the question on how to satisfy these boundary conditions in the context of the first order equations. The boundary conditions themselves are first order equations, as seen in Eqs.(5) and (6), and so, their fulfilment is easily understood in the context of the second order equations. For instance the derivative of the order parameter normal to the surface must vanish, according to Eq.(7), but this condition cannot be imposed on ψ, obtained through Eq.(14) because this is a first order equation itself and there is not enough free parameters in this solution. However it is possible to satisfy Eq.(5) simply by the requirement that the density |ψ|2 be constant at the boundary, which introduces a new parameter, c0 , fixed by the non-linear integrated GL equation, Eq.(22). |ψ|2 = c20 for ~x at Υ, and, h ¯q 2 H 0 = Ha + 2π c. mc 0

(23) (24)

Therefore the constant H 0 is automatically fixed by c0 according to Eq.(15). The important point is that Eq.(23) is enough to guarantee that there is no current flowing outside the superconductor. This is just a direct consequence of Eq.(16). As |ψ|2 is constant along the ~ 2 is border there is no gradient tangential to it and only perpendicular to it, namely, ∂|ψ| normal to the surface, rendering J~ always tangent to the surface. The magnetization M3 is also directly obtained from the present formalism and easily shown to be paramagnetic in the absence of an applied external field. According to Eq.(15), B3 ≡

Z Σ

Z d2 x d2 x 2 h3 = H 0 − 4πµB |ψ| , Σ Σ Σ

10

(25)

given in units of the Bohr’s magneton, µB = h ¯ q/2mc. Just take the thermodynamical relation B3 = Ha + 4πM3 and Eq.(24) to obtain that, M3 = µB

Z Σ

d2 x 2 c0 − |ψ|2 , Σ

(26)

which means that the integral has the dimension of inverse volume. In case of no applied external field, the order parameter is maximum at the boundary, namely (|ψ|2 ≤ c20 ), and so the magnetization is paramagnetic, M3 > 0. Under the condition that the first order equations are satisfied the Gibbs free energy of Eq.(1) becomes, ∆G =

Z Σ

d2 x Σ

(

1 + β − π 2 +

!

¯q Ha h α(T ) + |ψ|2 2 mc h ¯q mc

!2 |ψ|4

(H 0 − Ha )2 + 8π

h ¯ 2 I dl ~ 2 n ˆ · ∂|ψ| . 4m Υ Σ

(27)

Two direct contributions of the boundary to the free energy, which are not present in A. Abrikosov’s treatment of the GL theory [21], are the field energy due to H 0 6= Ha , according to Eq.(24), and the perimetrical contribution of the normal gradient of |ψ|2 . Near to the transition to the normal state the order parameter is weak, fact that allows for the expansion of the free energy in powers of ψ. From the present point of view this weakness also leads to the proposal of an iterative method to solve the first order equations. Firstly one seeks a solution for ψ in Eq.(14) under a known external field H0 sufficiently near to the upper critical field line which sets the order parameter in the vicinity of the normal state transition. Any solution of Eq.(14), multiplied by a constant is also a solution and this constant is fixed by the integrated equation (Eq.(22)). Eq. (14) can be rewritten as, "

∂ ∂ −i + ∂x1 ∂x2

!

#

2π − (A1 + iA2 ) ψ = 0. Φ0

(28)

Let us introduce the complex notation, z ≡ x1 + ix2 , z¯ ≡ x1 − ix2 into Eq.(28), and consider a constant external field, Ha , such that A1 = −Ha x2 /2 and A2 = Ha x1 /2. In this case the above equation becomes, ∂ψ(z, z¯) πHa =− z ψ(z, z¯), ∂ z¯ 2Φ0

11

(29)

where ∂/∂z = 1/2(∂/∂x1 − i∂/∂x2 ) and its complex conjugate is ∂/∂ z¯ = 1/2(∂/∂x1 + i∂/∂x2 ). The solution of Eq.(29) is promptly found to be, −

ψ(z, z¯) = f (z)e

πHa 2Φ0

z z¯

,

(30)

where f (z) is any function of z. The local field is equal to Ha plus a correction proportional to |ψ|2 according to Eq.(15). Very near to the normal state one expect that this correction is small such that it suffices to solve Eq.(14) and feed the solution in Eq.(15) without any ~ associated to the local further recurrence, namely, a returns to Eq.(14) with a corrected A field h3 . Thus the quest for a solution ψ(z, z¯) in a given geometry with cross sectional area Σ and boundary Υ is reduced to the search of an analytical function f (z) that will render |ψ(z, z¯)|2 constant at the boundary Υ. Assuming that in the boundary |ψ|2 = c20 this means that the analytical function f (z) must satisfy the condition |f (z)| = c0 exp

πHa |z|2 2Φ0

, where

z belongs to the boundary. For a circular disk, |z| is automatically constant at the boundary, but in case of an arbitrary geometry this characteristic is no longer satisfied. In this paper we will only consider the case that there is no external applied field, Ha = 0.

THE DIMENSIONLESS UNITS

At this point we find useful to switch to dimensionless units and rewrite all previous expressions in this fashion. As well known the GL theory has only one coupling constant, 1/2

κ = λ/ξ, the ratio between the London penetration length, λ = (mc2 /4πq 2 ψ02 )

1/2

coherence length, ξ = h ¯ 2 /2m|α|

, respectively, where ψ0 = (|α|/β)

1/2

, and the

. This renders

this ratio temperature independent, κ = (β/2π)1/2 mc/¯ hq and only dependent on material parameters. Let us refer just in this paragraph to the dimensionless units by the prime notation. For instance distance is measured in terms of the coherence length such that ~x = ξ~x0 . The magnetic field is expressed in units of the upper critical field, ~h = Hc2~h0 , √ √ Hc2 = 2κHc = Φ0 /2πξ 2 , Φ0 = hc/q, where Hc = Φ0 /2π 2λξ is the thermodynamical ~ = Hc2 ξ A ~ 0 , and the covariant field. Then the dimensionless vector potential is given by A derivative becomes Dj = Dj0 /ξ, Dj0 ≡ 1i ∂j0 − A0j since ∂j = ∂j0 /ξ. Lastly the dimensionless order parameter is obtained from ψ = ψ0 ψ 0 . The order parameter value at the boundary is also defined dimensionless, c0 = ψ0 c00 , and finally the dimensionless magnetization is M3 = M30 Hc2 /8πκ2 since ψ02 µB = Hc2 /8πκ2 . 12

Hereafter we drop the prime notation in all quantities, which means that they are all expressed in dimensionless units. The first order equations become, D+ ψ = 0, and, 1 h3 = H 0 − 2 |ψ|2 . 2κ

(31) (32)

Once the order parameter and the local magnetic field satisfy the first order equations, the kinetic energy of Eq.(8) becomes, ~ 2 = H 0 − 1 |ψ|2 |ψ|2 + |Dψ| κ2

1~ 2 2 ∂ |ψ| . 2

(33)

The integrated equation, Eq.(22), in reduced units becomes, d2 x 1 (H 0 − 1) |ψ|2 + 1 − 2 |ψ|4 = 0. Σ 2κ

Z Σ

(34) The Gibbs free energy difference of Eq.(1) is in units of Hc2 /4π = |α|2 /β , given by, ∆G ≡

d2 x 1 ~ 2 −|ψ|2 + |ψ|4 + |Dψ| 2 Σ Σ o κ2 (h3 − Ha )2 .

G−Gn = |α|2 /β

+

Z

(35)

Similarly, Eq.(27), becomes 1 1 d2 x ∆G = − (1 − Ha ) |ψ|2 + 1 − 2 |ψ|4 2 2κ Σ Σ 1~ 2 2 + ∂ |ψ| + κ2 (H 0 − Ha )2 . 2

Z

(36)

THE VORTEX EXCITED STATE AND CONFORMAL MAPPING

In the absence of an applied field outside the superconductor, Ha = 0 , and with the help of Eq.(34), the Gibbs free energy density of Eq.(36) becomes, ∆G =

Z Σ

1 d2 x − (1 + H 0 ) |ψ|2 + Σ 2

1~ 2 2 2 ∂ |ψ| + κ2 H 0 . 2

(37) Then H 0 follows directly from Eq.(24). Near to the normal state the order parameter is small enough that iteration of the first order equations is not necessary. This means that the local 13

magnetic field is small enough that it can be dropped in Eq.(38) and directly obtained from Eq.(39). In this case Eqs.(31) and (32) are reduced to, ∂+ ψ = 0, and, c2 − |ψ|2 h3 = 0 2 . 2κ

(38) (39)

The connection between vortex states and conformal mapping stems from Eq.(29), which simply becomes, ∂ψ(z, z¯) = 0. ∂ z¯

(40)

Expressing the order parameter as ψ ≡ ψR +iψI in Eq.(38), one obtains that ∂1 ψR +∂2 ψI = 0 and ∂2 ψR −∂1 ψI = 0, which are the well-known Cauchy-Riemann conditions for analyticity. We recall the so-called Maximum Modulus Theorem in Mathematics which says that for an analytical function ψ(z) in a given region Σ, the maximum of |ψ(z)| necessarily falls in its boundary Υ. This theorem is useful because from it we know that the order parameter, constant at the boundary Υ, indeed reaches its maximum value there, and therefore the magnetization is necessarily paramagnetic. The presence of vortices inside will only lead to the vanishing of the order parameter at points in its interior, ψ = 0, and still the maximum of |ψ|2 is at the boundary. This feature helps to establish fundamental differences between the finite boundary and the bulk superconductor, the latter understood as the case that only periodic solutions are sought. According to Liouville’s theorem [28] any periodic analytical function ψ with zeros must also diverge, from which one draws the conclusion that the only possible physical solution is the constant one. In other words it is not possible to find vortex solutions without the presence of an applied field in a unit cell with periodic boundary conditions. However they do exist for the long superconductor with a finite cross section. Thus we conclude from the above discussion that the present method is general and applies for any number of vortices in any cross sectional geometry. As examples of our general method, we detail in this paper two particular cases of a disk, namely, a vortex with vorticity L in its center and a vortex with vorticity one in any position inside the disk. Let us define the integrals, 1 Z d2 x 2 |ψ| , c20 Σ Σ 1 Z d2 x 4 I4 ≡ 4 |ψ| , and, c0 Σ Σ I2 ≡

14

(41) (42)

1 Z d2 x ~ 2 2 ∂ |ψ| . IΥ ≡ 2 c0 Σ Σ

(43)

From the integrated equation, Eq.(34) and the condition that the local magnetic field vanishes at the boundary, one obtains that, 1 − H 0 I2 , and , 1 − 1/2 κ2 I4 c2 H 0 = 02 . 2κ c20 =

(44) (45)

Solving these equations for c0 and H 0 , we obtain that, 2κ2 I2 /I4 , and , 2κ2 − 1 + I2 /I4 I2 /I4 H0 = 2 . 2κ − 1 + I2 /I4 c20 =

(46) (47)

Notice that I2 > I4 is indeed satisfied since the order parameter divided by c0 is always smaller than one inside the disk, and so, at each point the second power, Eq.(41), is larger √ than the fourth power, Eq.(42). The fact that I2 > I4 and κ > 1/ 2, always render a solution for c0 . The difference of the free energy density given at Eq.(37) can be expressed in terms of c0 , 1 1 ∆G = − (I2 − IΥ ) c20 − 2 (I2 − 1)c40 , 2 4κ

(48)

∆G = −κ2 (I2 − IΥ ) H 0 − κ2 (I2 − 1)H 02 .

(49)

or in terms of H 0 ,

The magnetization given by Eq.(26) becomes M3 = 2κ2 H 0 − I2 c20 , and so, equal to, M3 =

2κ2 (1 − I2 ) I2 /I4 . 2κ2 − 1 + I2 /I4

(50)

This interesting expression shows that the magnetization depends on the position of the vortices inside the superconductor, since I2 and I4 vary accordingly, as shown in the next example. In case of no vortices, I2 = I4 = 1, the magnetization vanishes. We also include the general expression for the magnetic field at the center of the vortex, as h3 (v) = H 0 : h3 (v) =

I2 /I4 . 2κ2 − 1 + I2 /I4

(51)

We consider two particular examples of a thin long cylinder with radius R, namely, a vortex L at the center and a vortex L = 1 at arbitrary position. Through them we see the general 15

aspects of the present theory such as the importance of the boundary term which makes the Gibbs energy explicitly dependent on R. Remarkably the cooper pair density at the boundary determines the local magnetic field at the vortex core. There ψ(v) = 0 and so h3 (v) = c20 /2κ2 according to Eq.(32) where v refers to the center of the vortex.

LONG CYLINDER WITH A CIRCULAR CROSS SECTION

We find useful to apply our theory to Niobium [29–31], one of the favored materials to study the characteristics of vortex matter in superconductors and also used to construct nano-engineered superconductors [32] All figures are expressed in reduced units and to retrieve the predicted values for Niobium we take the parameter values reported in Ref. 31, namely, κ = 2.1, λ(0) = 42 nm, ξ(0) = 20 nm. In particular we choose the temperature of T = 7.7 K which is near to the critical temperature, Tc = 9.3 K, such that an order parameter approach is valid. Thus for this temperature the local magnetic field, expressed in units of the upper critical field, must be multiplied by Hc2 (T = 7.7K)= 780 Oe. Similarly the magnetization must be multiplied by Hc2 (T = 7.7K)/8πκ2 = 7.0 Oe.

Vorticity L at the center

In this simple example we show that the EVS exists only in a special range of the radius and the vorticity. The search for the order parameter in a disk cross section of radius R that satisfies Eq.(38), is reduced to find an analytical function that is constant at the perimeter of the disk. This is simply given by,

ψ = c0

z R

L

= c0

r R

L

eiLθ ,

(52)

where c0 is the value of the order parameter at the boundary and L is an integer since the order parameter is assumed to be single-valued, ψ(θ + 2π) = ψ(θ). In polar coordinates √ the order parameter is expressed in terms of r ≡ z z¯ and tan θ ≡ x2 /x1 . The fulfillment of |ψ| = c0 at the circle guarantees the confinement of the current to the superconductors boundaries and ψ = 0 at the center means that the solution has non-zero vorticity for L 6= 0. For L = 0 there is no zero and so, describes the homogeneous ground state. Then one can 16

determine the integrals of Eqs.(41), (42) and (43), respectively. 1 , L+1 1 I4 = , and, 2L + 1 4L IΥ = 2 . R

(53)

I2 =

1.5

(54) (55)

L 0 1 2 3 4

|ψ(r)|2

1.0

0.5

0.0 0.0

0.2

0.4

r/R

0.6

0.8

1.0

FIG. 1. The superconducting density as a function of the relative distance to the cylinder’s center for vorticity ranging from L = 0 to L = 4 in case of a vortex fixed at the center (κ = 2.1).

From Eqs.(46) and (47), one obtains that, 2κ2 (2L + 1) , and , (56) 2κ2 (L + 1) + L (2L + 1) H0 = 2 . (57) 2κ (L + 1) + L Inserting these expressions into the Gibbs free energy density, either Eq.(48) or (49), it gives c20 =

that, "

1 1 4L L (2L + 1) ∆G = − − 2− 2 2 L+1 R L + 1 2κ (L + 1) + L 17

#

×

2κ2 (2L + 1) . 2κ2 (L + 1) + L

(58)

The condition for the existence of the EVS, ∆G ≤ 0, can only be achieved for a limited range of parameters R and κ. It means that the EVS exists for a given κ superconductor in case the radius is larger than a critical value, given by, Rc =

v u u 2tL(L + 1) ·

2κ2 (L + 1) + L . 2κ2 (L + 1) − 2L2

(59)

q

In the limit κ → ∞, the critical radius becomes Rc = 2 L(L + 1), and the Gibbs free energy for R > Rc becomes negative, ∆G → [4L(L + 1) − R2 ](2L + 1)/[2R2 (L + 1)2 ]. The order parameter and the local magnetic field are given by, ψ(r, θ) = h3 (r) =

v u u t

2κ2 (2L + 1) r L iLθ e , and, 2κ2 (L + 1) + L R

h r 2L i (2L + 1) 1 − . 2κ2 (L + 1) + L R

(60) (61)

The paramagnetic magnetization and the field at the center of the vortex, are given by, L 2κ2 (2L + 1) , and, L + 1 2κ2 (L + 1) + L (2L + 1) . h3 (v) = 2 2κ (L + 1) + L

M3 =

(62) (63)

respectively. Fig. 1 shows the superconducting density according to the distance to the center for vorticity L = 0 (the homogeneous state) to L = 4. The density |ψ|2 obtained from Eq. (52) and c20 given by Eq. (56). The superconducting density is maximum at the boundary, which implies in a paramagnetic effect, as previously shown. The superconducting density √ reaches c20 at the boundary and is a slow growing function of L for κ ≥ 1/ 2. Fig. 2 depicts the local magnetic field versus the distance to the center, for vorticity ranging from L = 0 to L = 4. This plot shows h3 (r) obtained from Eq. (39), with c20 and |ψ|2 as described in Fig. 1. Notice that the field at the core h3 (v), where v corresponds to r = 0, follows from Eq.(51). Thus it depends on L through the integrals I2 and I4 , defined by Eqs. (53) (54)). Interestingly the magnetic field at the center of the vortex and the superconducting density at the boundary are directly related to each other, as previously pointed out. With the exception of the homogeneous state, which has no internal magnetic field, for all other states L ≥ 1 the local magnetic field reaches its maximum at the center and vanishes at 18

L 0 1 2 3 4

0.25

h3(r)

0.20 0.15

0.10 0.05 0.00 0.0

0.2

0.4

r/R

0.6

0.8

1.0

FIG. 2. The local magnetic field as a function of the relative distance to the cylinder’s center for vorticity ranging from L = 0 to L = 4 in case of a vortex fixed at the center (κ = 2.1). To recover the Niobium values at the temperature of T = 7.7 K the vertical axis must be multiplied by Hc2 (T = 7.7K)= 780 Oe

the boundary, as expected. For L ≥ 2 the local field varies slowly from the center to the middle of the cylinder to then abruptly show a strong decay. The Gibbs free energy density as a function of the cylinder’s radius is shown for two cases, R = 1 and R = 10, in Fig. 3. This plot shows that the homogenous state L = 0 is the absolute ground state with the minimum energy, ∆G = −0.5. Recall that our search is for the EVS, namely for ∆G < 0 otherwise ∆G > 0 and the superconducting state can somehow decay into the normal state. For instance, the extreme value of this plot, R = 10, shows that the L = 0, 1, 2 states are EVS whereas the L = 3, 4 ones are not. The Gibbs free energy has the following values, ∆G = −0.27, −0.11, 0.003 and 0.09 for L = 1, 2, 3, 4, respectively. The EVS exists for R > Rc , Rc = 3.1, 6.0, 10.2, 17.3 for L equal to 1, 2, 3, 4, respectively. 19

1.0

L 0 1 2 3 4

0.8 0.6 0.4

∆G

0.2

0.0 0.2 0.4 0.6 0

2

4

R

6

8

10

FIG. 3. The Gibbs free energy difference between the superconducting and the normal states as a function of the cylinder’s radius for vorticity ranging from L = 0 to L = 4 in case of a vortex fixed at the center (κ = 2.1). The excited vortex state only exists in the negative range of this difference.

Vorticity L = 1 at any position

The analytical function of z that describes the order parameter of a single vortex at any position 0 ≤ a ≤ R inside a disk is given by, z R

− Ra ψ = c0 , 1 − Rz Ra

(64)

where c0 is the order parameter value at the boundary. The coordinate z and the position a can be expresses in polar form, z = r eiθ and a = a0 eiα , where we can consider α = 0 without any loss of generality. ψ = c0 R

reiθ − a0 . R2 − a0 r eiθ 20

(65)

The density of superconducting electrons is given by, r2 + a20 − 2 a0 r cos(θ) |ψ| = |c0 | · R 4 . R + a20 r2 − 2 a0 r R2 cos(θ) 2

2

2

(66)

The integrals of Eqs.(41),(42) and (43) can be exactly obtained and are given by, !2

R2 a20 R2 − 1 ln 1 − 2 , I2 = 2 − 2 − a0 a20 R 4 2 R R I4 = −4 +6 −1 a0 a0 #2 ! 2 " 2 R R a20 −4 − 1 ln 1 − 2 , and, a0 a0 R 4 IΥ = 2 . R !

(67)

(68) (69)

Next we take two special limits of these expressions. In the first one the vortex is just slightly displaced from the center, and the expressions of the previous section must be retrieved. This is the limit a0 → 0, but to obtained it, the following approximate expansion for the logarithm function must be introduced, ln(1 − x2 ) ≈ −x2 − x4 /2 − x6 /3 − x8 /4 − x10 /5 − x12 /6. Then one obtains the needed approximated expressions for the integrals I2 and I4 , respectively I2 ≈ + I4 ≈ −

1 1 a0 2 1 a0 4 1 a0 6 + + + 2 3 R 12 R 30 R 1 a0 8 2 a0 10 1 a0 12 − + , and, 60 R 15 R 6 R 2 a0 4 1 a0 6 1 1 a0 2 + + + 3 3 R 15 R 15 R 8 a0 8 2 a0 10 + . 15 R 3 R

(70)

(71)

From these approximated expressions we can easily verify that I2 = 1/2 , I4 = 1/3 and IΥ = 4/R2 , which are exactly the Eqs. (53), (54) and (55) for L = 1. The other interesting limit is the vortex very near to the boundary of the cylinder. To treat it we change to the coordinate that express the distance of the vortex to the boundary, defined by y ≡ R − a0 . Inserting this new parameter into the order parameter we obtain the superconductor’s density given by, 2 y

|ψ|2 = c20

R

1+

+ 1−

2 r R

y 2 − R 2 y

1−

R

2 1−

y R

−2 1−

y R

r R

cos(θ)

r R

.

cos(θ) (72)

21

The local magnetic field (Eq. (32)) takes into account the density given by Eq.(72). The expressions for c0 and H 0 , given by Eqs.(46) and (47), respectively, are functions of I2 and I4 , which in terms of the coordinate y are given by,

I2 = 2 − 1 − "

y R

−2

−2

#2

"

#

y 2 − − 1 ln 1 − 1 − , and, R y −4 y −2 I4 = −4 1 − +6 1− −1 R R #2 " " # y 2 y −2 y −2 − 1 ln 1 − 1 − . −4 1− 1− R R R y 1− R

(73)

(74) Similarly the Gibbs free energy density ( Eq.(48) or (49)) is a function of Eqs.(73) and (74), but also of Eq.(69), which renders it explicitly R dependent. Very near to the boundary of the disk R >> y, which means that the limit y ≈ 0 must be taken. Consider only the first order term in the expansion in y to obtain that I2 ≈ 1 − 2y/R and I4 ≈ 1 − 4y/R, and so, I2 /I4 ≈ 1 + 2y/R. Then we obtain that, 1 1 1 H ≈ 2 + 2 1− 2 2κ κ 2κ 1 y c20 ≈ 1 + 2 1 − 2 . 2κ R

0

y , and, R

(75) (76)

The paramagnetic magnetization becomes, y M3 ≈ 2 . R

(77)

The Gibbs free energy density is given by, 4 1 4 1 1 ∆G ≈ − 1 − 2 + 2 + 2 1 − 2 2 R κ R 2κ

y . R

(78)

For R → ∞, c20 = 1, H 0 = 1/2κ2 , M3 = 0 and ∆G = −1/2, which means that the homogenous state is recovered. Fig. 4 shows the Gibbs free energy density with respect to the ratio a0 /R. The plotted curves are obtained from Eq.(48) with c20 given by Eq. (46) and the values of the integrals I2 , I4 and IΥ given by Eqs. (67), (68) and (69), respectively. The two cases shown, R = 3 and R = 10 demonstrate that the existence of EVS depends on the position of the vortex. In both cases the Gibbs free energy decreases monotonically from the center to the boundary of 22

the cylinder. For R = 10 the vortex in any position is in a EVS, since the Gibbs free energy is always negative, but this is not so for the R = 3 cylinder. There the Gibbs free energy is positive for a0 < 0.24R such that only for a0 > 0.24R there is an EVS. The Gibbs free energy is null for a0 = 0.24R and decreases to reach the value ∆G = −0.27 at the boundary. Fig. 5 illustrates how the local magnetic field at the vortex’s core varies according to the vortex position. This figure shows that this field h3 (v) depends on the ratio a0 /R according to Eq. (51). Thus it only depends on the ratio I2 /I4 , defined by Eqs. (67) and (68), respectively. The local field is maximum at the vortex’s center. From its side h3 (v) reaches a maximum when the vortex is at the center and slowly decreases when the vortex is positioned at the boundary. The Fig. 5 is in agreement with Fig. 2 since h3 (0) for

R =3 R =10

0.0

∆G

0.1 0.2 0.3 0.4 0.50.0

0.2

0.4

a0/R

0.6

0.8

1.0

FIG. 4. The Gibbs free energy difference between the superconducting and the normal states as as a function of the vortex’s position inside the cylinder in case of vorticity one. Two radius of cylinders are considered (κ = 2.1). The excited vortex state only exists in the negative range of this difference.

23

a0 = 0 corresponds to h3 (r = 0) for L = 1. Fig. 6 shows the paramagnetic magnetization as a function of the ratio a0 /R as obtained from Eq. (50). Interestingly the magnetization depends on the position of the vortex, according to the integrals I2 and I4 , given by Eqs.(67) and (68), respectively. The magnetization is stronger for the vortex near to the cylinder’s center and weaker near to the cylinder’s border. For a vortex at the boundary, a0 = R, the magnetization vanishes because I2 = 1. This is consistent with the description of the exit of the vortex, although h3 (v) is not zero as can be seen in Fig. 5. The vortex next to the boundary means that the superconducting density is almost homogenous on the entire region of the cylinder and the integral I2 is almost equal to the unity what makes the magnetization goes to zero as can be seen from Eq.(50). Finally the three rows of Fig. 7 depict the vortex at three different positions, a0 /R = 0.1, 0.5 and 0.9, respectively. The columns correspond

0.16 0.15

h3(v)

0.14 0.13 0.12 0.110.0

0.2

0.4

a0/R

0.6

0.8

1.0

FIG. 5. The magnetic field at the center of the vortex as a function of the position of the vortex inside the cylinder (κ = 2.1). To recover the value for Niobium at the temperature of T = 7.7 K the vertical axis must be multiplied by Hc2 (T = 7.7K)= 780 Oe.

24

to the density, local magnetic field, current and phase, respectively. The density is obtained from Eqs.(66) and (46), with the integrals of Eqs.(73) and (74). Therefore the density ranges from zero at the center of the vortex to the maximum value c20 at the boundary. The density is shown in Figs. 7(a), (e) and (i) in color (on line) scheme ranging from low density (cyan) to high density (magenta). The local magnetic field is obtained from Eq.(32) and ranges from the value given by Eq.(63) at the center of the vortex to zero at the boundary. Thus the (color on line) scheme of Figs. 7(b), (f) and (j) ranges from the maximum at the center of the vortex (magenta) to zero at the boundary(cyan). The vortex current is obtained from Eq.(16) and depicted in Figs. 7(c), (g) and (k). Notice that there is no current flowing out of the cylinder as the vortex moves towards the boundary, which is a general property guaranteed by the present formalism. Finally Figs. 7(d), (h) and (l) show the phase of the

0.8 0.7 0.6

M3

0.5 0.4 0.3 0.2 0.1 0.0 0.0

0.2

0.4

a0/R

0.6

0.8

1.0

FIG. 6. The magnetization as a function of the position of the vortex inside the cylinder (κ = 2.1). To recover the value for Niobium for the temperature of T = 7.7 K the vertical axis must be multiplied by Hc2 (T = 7.7K)/8πκ2 = 7.0 Oe.

25

order parameter, defined as tan−1 (ψI /ψR ). The discontinuity, represented as a straight line that abruptly separates white to black ,confirms the presence of a vortex in the cylinder.

CONCLUSION

In this paper we have included the treatment of boundaries in the first order equations method used by A. Abrikosov to predict the vortex lattice. We find a direct connection between this method and the theory of conformal mapping. Using this method we have considered vortex states left inside the superconductor after the external applied field is switched off. This vortex state is unstable which means that vortices must leave the superconductor since their presence is thermodynamically forbidden. However vortices are topological stable as long as the superconducting state exists. We define here the excited vortex state, whose free energy is lower than that of the normal state. Thus the excited vortex state cannot collapse into the normal state. We calculate the paramagnetic magnetization of the excited vortex state for some particular cases in a thin long wires of Niobium. Acknowledgments: We acknowledge J. C. Lopez Ortiz for helpful discussions. M.M.D acknowledges CNPq (Brazil) support from funding (23079.014992/2015-39). R. R. G. acknowledges CNPq (Brazil) support from funding (207007/2014-4). A. R. de C. R. acknowledges FACEPE (Brazil) support from funding (APQ-0198-1.05/14).

The kinetic energy decomposition

Consider the term |D+ ψ|2 , which can be casted as, |D+ ψ|2 = [(D1 ψ)∗ − i(D2 ψ)∗ ][(D1 ψ) + i(D2 ψ)] = (D1 ψ)∗ (D1 ψ) + (D2 ψ)∗ (D2 ψ) + i[(D1 ψ)∗ (D2 ψ) − (D2 ψ)∗ (D1 ψ)].

(79)

Expanding only the derivative D1 , q h ¯ (D1 ψ)∗ (D2 ψ) = − ∂1 ψ ∗ A1 ψ ∗ (D2 ψ) i c h i h ¯ ∗ = ∂1 − ψ (D2 ψ) i h ¯ ∗ q − ( ψ )∂1 D2 ψ − A1 ψ ∗ D2 ψ. i c

26

(80)

Rearranging the terms, one obtains that, h ¯ (D1 ψ)∗ (D2 ψ) = − ∂1 (ψ ∗ D2 ψ) + ψ ∗ D1 D2 ψ. i

(81)

In the same way, we obtain that, h ¯ (D2 ψ)∗ (D1 ψ) = − ∂2 (ψ ∗ D1 ψ) + ψ ∗ D2 D1 ψ. i

(82)

The complex conjugate of these relations give us that, h ¯ ∂1 [(D2 ψ)∗ ψ] + (D1 D2 ψ)∗ ψ, i h ¯ ∗ (D1 ψ) (D2 ψ) = ∂2 [(D1 ψ)∗ ψ + (D2 D1 ψ)∗ ψ. i (D2 ψ)∗ (D1 ψ) =

(83) (84)

At this point, we have two formulations for the same identity, where one relation is the complex conjugate of the second one. The expression ~ 2 |D+ ψ|2 = |Dψ| h i h ¯ + i − ∂1 (ψ ∗ D2 ψ) + ψ ∗ D1 D2 ψ i i h h ¯ − − ∂2 (ψ ∗ D1 ψ) + ψ ∗ D2 D1 ψ , i becomes, ~ 2 + iψ ∗ [D1 , D2 ]ψ |D+ ψ|2 = |Dψ| −h ¯ ∂1 (ψ ∗ D2 ψ) + h ¯ ∂2 (ψ ∗ D1 ψ).

(85)

The commutator, h ¯ q h ¯ q [D1 , D2 ] = [ ∂1 − A1 , ∂2 − A2 ] i c i c h ¯q = − (∂1 A2 − ∂2 A1 ), ic

(86)

becomes, [D1 , D2 ] = −

h ¯q h3 . ic

(87)

Then one obtains that, ~ 2 − i([D1 , D2 ]ψ ∗ )ψ |D+ ψ|2 = |Dψ| −h ¯ ∂1 [(D2 ψ)∗ ψ] + h ¯ ∂2 (D1 ψ)∗ ψ], 27

(88)

~ 2, which gives expressions for |Dψ| h ¯q h3 |ψ|2 c +h ¯ [∂1 (ψ ∗ D2 ψ) − ∂2 (ψ ∗ D1 ψ)],

(89)

h ¯q h3 |ψ|2 c +h ¯ {∂1 [(D2 ψ)∗ ψ] − ∂2 [(D1 ψ)∗ ψ]}.

(90)

~ 2 = |D+ ψ|2 + |Dψ|

and ~ 2 = |D+ ψ|2 + |Dψ|

We sum the two expressions and divide by 2 to obtain that, h ¯q h3 |ψ|2 c ψ ∗ D2 ψ + (D2 ψ)∗ ψ ψ ∗ D1 ψ + (D1 ψ)∗ ψ +h ¯ {∂1 − ∂2 }. 2 2

~ 2 = |D+ ψ|2 + |Dψ|

(91) Introducing the definition of the current we obtain the desired dual formulation of the kinetic energy given by Eq.(8)

Calculation of the integrals

In this appendix we calculate the following integrals. R2 Z R dr r f2 (r), Σ 0 Z 2π r2 + a20 − 2 a0 r cos(θ) f2 (r) = dθ 4 , R + a20 r2 − 2 a0 r R2 cos(θ) 0 R4 Z R I4 = dr r f4 (r), Σ 0 2 Z 2π (r2 + a20 − 2 a0 r cos(θ)) f4 (r) = dθ 2, 0 (R4 + a20 r2 − 2 a0 r R2 cos(θ)) 1 Z 2 ~2 1I ~ Υ (r, θ), IΥ = d x ∂ fΥ (r, θ) = dl n ˆ · ∂f Σ Σ ∂Σ r2 + a20 − 2 a0 r cos(θ) fΥ (r, θ) = R2 4 . R + a20 r2 − 2 a0 r R2 cos(θ) I2 =

(92) (93) (94) (95) (96) (97)

The integral I2 is associated with the density of Cooper pairs. As shown in Eq.(92) the integral of |ψ|2 is obtained in two parts: first we integrate the angular part, that results in a 28

function of the radius f2 (r) and after we integrate the function f2 (r) in the radial coordinate. The angular integral has a form with a well-known result [33]. Z

dx

B Ab−aB Z 1 A + B cos(x) = x+ dx , a + b cos(x) b b a + bcos(x)

1 2 =√ 2 tan−1 2 a + bcos(x) a −b 2 2 , for a > b .

Z

dx

√

(98) a2 − b2 tan(x/2) a+b

!

(99)

To calculate this integral we consider separately the integration in the four quadrants,as shown below. A + Bcos(x) Z π/2 A + Bcos(x) dx = a + bcos(x) a + bcos(x) 0 0 Z π Z 3π/2 A + Bcos(x) A + Bcos(x) dx dx + + a + bcos(x) a + bcos(x) π/2 π Z 2π A + Bcos(x) . dx + a + bcos(x) 3π/2 Z 2π

dx

(100)

We make the following change of variables. First quadrant: y = x Z π/2

Z π/2 A + Bcos(x) A + Bcos(y) dx dy → . a + bcos(x) a + bcos(y) 0

0

(101) Second quadrant: y = x − π Z 0 A + Bcos(x) A − Bcos(y) dx → dy . a + bcos(x) a − bcos(y) π/2 −π/2

Z π

(102) Third quadrant: y = x − π Z 3π/2

dx

π

Z π/2 A + Bcos(x) A − Bcos(y) → dy . a + bcos(x) a − bcos(y) 0

(103) Fourth quadrant: y = x − 2π Z 2π 3π/2

dx

Z 0 A + Bcos(x) A + Bcos(y) → dy . a + bcos(x) a + bcos(y) −π/2

(104) 29

Inserting the variable changes in the original integral, we obtain a new interval of integration Z 2π

A + Bcos(x) Z π/2 A + Bcos(y) dx = dy a + bcos(x) a + bcos(y) 0 −π/2 Z π/2 A − Bcos(y) . + dy a − bcos(y) −π/2

(105)

We calculate the integrals in the interval −π/2 < x < +π/2 Z +π/2

A + Bcos(x) B = π a + bcos(x) b −π/2 √ ! 2 4(Ab − aB) a − b2 π −1 + √ 2 tan tan( ) , a+b 4 b a − b2 dx

A − Bcos(x) B = π a − bcos(x) b −π/2 √ ! 2 4(A b − a B) π a − b2 −1 + √ 2 tan tan( ) . a−b 4 b a − b2 Z +π/2

(106)

dx

(107)

We obtain the expressions for the integral in the interval 0 < x < 2π A + Bcos(x) B 4(A b − a B) = 2π + √ 2 a + bcos(x) b 0 b a − b2 √ √ " ! !# a2 − b 2 a2 − b 2 −1 −1 × tan + tan , a+b a−b Z 2π

dx

(108)

and verify that, √

q

1/2 (a + b)(a − b) a−b − = = , a+b (a + b)2 a+b q √ !1/2 (a + b)(a − b) a2 − b 2 a+b = = . a−b (a − b)2 a−b

a2

b2

!

(109) Using the identity tan−1 (x) + tan−1 (1/x) = π/2 if x > 0, we simplify the expression to obtain that, Z 2π 0

dx

A + Bcos(x) B (A b − a B) = 2π + 2π √ 2 . a + bcos(x) b b a − b2 (110)

At this point we have a general result for the integral (98). The next step is to apply this result for integrand of f2 (r). We consider the following identifications A = a20 + r2 , 30

B = −2a0 r, a = R4 + a20 r2 and b = −2a0 rR2 . Inserting these expressions in the Eq. (110) we obtain f2 (r) =

2π 2π (R2 − a20 )(R2 − r2 ) − . R2 R2 R4 − a20 r2

(111)

With the result obtained in Eq.(111), the integral I2 becomes, R2 Z R rf2 (r)dr Σ 0 Z R rdr 2π Z R 2π 2 2 2 = rdr − (R − a0 )R 2 4 Σ 0 Σ 0 R − a0 r 2 Z R r3 dr 2π 2 + (R − a20 ) . 2 Σ 0 R 4 − a0 r 2

I2 =

(112)

Next we calculate each one of the three integrals in Eq.(112) Z R 0

Z R 0

Z R 0

R2 , 2 rdr = 4 R − a20 r2 r3 dr = R4 − a20 r2

rdr =

1 R4 , ln 2a20 R4 − a20 R2 ! R4 R2 R4 − ln . 2a40 R4 − a20 R2 R2 − a20 !

Inserting these results in Eq.(112) and considering that Σ = πR2 we have !2

R2 R2 I2 = 2 − 2 + −1 a0 a20

a2 ln 1 − 02 . R !

(113)

Next we calculate the surface integral defined in Eq.(96). The normal vector is the radial ~ Υ = ∂r fΥ . To calculate the integral with the derivative, one, n ˆ = rˆ, such that, fΥ (r, θ): rˆ · ∂f ∂fΥ (r, θ) ∂r = 2R2 (R2 − a20 )

(R2 + a20 )r − a0 (R2 + r2 ) cos(θ) , (R4 + a20 r2 − 2a0 rR2 cos(θ))2 (114)

we use the formula [33] Z

"

A + B cos(x) 1 (a B − A b) sin(x) dx = 2 2 2 (a + b cos(x)) a −b a + b cos(x) # Z Aa−bB + dx . a + b cos(x)

(115)

To obtain the integral IΥ , we divide the interval 0 < x < 2π into four quadrants and use the following changes of variable for each of them. 31

First quadrant: y = x Z π/2 0

Z π/2 A + B cos(x) A + B cos(y) dx → dy . 2 (a + b cos(x)) (a + b cos(y))2 0

(116) Second quadrant: y = x − π Z π

dx

π/2

Z 0 A + B cos(x) A − B cos(y) → dy . (a + b cos(x))2 (a − b cos(y))2 −π/2

(117) Third quadrant: y = x − π Z 3π/2

dx

π

Z π/2 A + B cos(x) A − B cos(y) → dy . 2 (a + b cos(x)) (a − b cos(y))2 0

(118) Fourth quadrant: y = x − 2π Z 2π

Z 0 A + B cos(x) A + B cos(y) dx dy → . 2 (a + b cos(x)) (a + b cos(y))2 3π/2 −π/2

(119) We collect these results together and we obtain a new integration interval −π/2 < x < π/2, and the original integral is rewritten in two parts, with positive and negative signs in the integrand, respectively. Z 2π

Z π/2 A + B cos(x) A + B cos(y) dx dy = 2 (a + b cos(x)) (a + b cos(y))2 0 −π/2 Z π/2 A − B cos(y) dy + . (a − b cos(y))2 −π/2

(120)

For the integral with the positive sign we find that, A + B cos(y) 2(a B − A b) = 2 (a + b cos(y)) a(a2 − b2 ) −π/2 √ ! 4(A a − b B) a2 − b 2 −1 + tan , (a2 − b2 )3/2 a+b Z π/2

dy

(121)

and for the negative sign, A − B cos(y) 2(a B − A b) =− 2 (a − b cos(y)) a(a2 − b2 ) −π/2 √ ! 4(A a − b B) a2 − b 2 −1 + tan . (a2 − b2 )3/2 a−b Z π/2

dy

32

(122)

Summing the two integrals and using the same trigonometric property for tan−1 (x) shown before, we obtain that, Z 2π

dx

0

A + B cos(x) Aa−bB . = 2π 2 2 (a + b cos(x)) (a − b2 )3/2

(123)

The surface integral can be written as, 1I ∂fΥ dl Σ ∂Σ ∂r ! ∂fΥ 1 Z 2π Rdθ . = Σ 0 ∂r r=R

IΥ =

(124)

Next we insert the derivative of the function fΥ in the previous integral, and define the following parameters A = (R2 + a20 )R, B = −2a0 R2 , a = R4 + a20 R2 and b = −2a0 R3 . We consider the area Σ = πR2 to obtain that, IΥ =

4 . R2

(125)

The last integral to be calculated, I4 , defined in Eq.(94), is the square of the superconductor density over the entire area of the cross section of the cylinder. Similarly to I2 , first we calculate the angular integral and define the results as a function of the radius, as shown in Eq.(95). To calculate the latter integral we expand the numerator of the argument and divide the original integral in two parts, defined as follows: f4 (r) = f4I (r) + f4II (r), where the first part is given by, f4I (r) ≡

Z 2π

dθ

0

(r2 + a20 )2 − 4 a0 r(r2 + a20 ) cos(θ) 2 , (R4 + a20 r2 − 2 a0 r R2 cos(θ)) (126)

and the second part is given by, f4II (r) ≡

Z 2π 0

dθ

4 a20 r2 cos2 (θ) 2. (R4 + a20 r2 − 2 a0 r R2 cos(θ)) (127)

The integral (126) is of the same type of the calculated integral (123). We just need to do A = (r2 + a20 )2 , B = −4a0 r(r2 + a20 ), a = R4 + a20 r2 and b = −2a0 rR2 . The second integral we write as a derivative of an integral with a known result: Z 2π 0

dx

∂ Z 2π −B cos(x) B cos2 (x) = dx . 2 (a + b cos(x)) ∂b 0 a + b cos(x) (128) 33

This integral is found in the right hand side of the previous equation and the result is given by, Z 2π

dx

0

B aB −B cos(x) = 2π + 2π 2 , a + b cos(x) b b(a − b2 )

(129)

where A = 0. Calculating the derivative, one finds the expression, f4II = 2π

B a B(a2 − 2 b2 ) − 2π , b2 b2 (a2 − b2 )3/2

(130)

where the constants shown in the final result are related to the original expression for f4II by B = 4a20 r2 , a = R4 + a20 r2 and b = −2a0 r R2 . The final result for Eq.(95) is a function of the radius and is given by 2π 2π 1 + R4 R4 (R4 + a20 r2 )3 h R6 R8 4 4 R4 6 6 4 × − 1 − 4 a0 r + 7R − 8 2 + 4 a0 r a0 a0 a0

f4 (r) =

i

+ 7R8 − 8a20 R6 + a40 R4 )a20 r2 − R8 (R4 − a40 ) .

(131)

To obtain I4 we insert the function f4 (r) and we consider the following change of variable x = a20 r2 /R4 . With this change of variable and some simplifications one obtains that I4 = (R2 /a20 )(a20 /R2 +

P3

j=0

dj Kj ), where the dj coefficients are defined as follows: a40 , R4 a4 a2 d1 ≡ 7 − 8 02 + 04 , R R 4 2 R R d2 ≡ 4 − 8 2 + 7, a0 a0 4 R d3 ≡ 4 − 1, a0 d0 ≡ −1 +

(132) (133) (134) (135)

and the integrals in the new variable x are defined as Kj ≡

Z a2 /R2 0 0

xj dx . (1 − x)3

(136)

The calculated Kj integrals give the results 1 a2 K0 = 1 − 02 2 R

!−2

1 a2 K1 = 1 − 02 2 R

!−2

34

1 − , 2 a2 − 1 − 02 R

(137) !−1

1 + , 2

(138)

!−2

1 a2 a2 K2 = 1 − 02 − 2 1 − 02 2 R R ! 2 a − ln 1 − 02 , R !−2

1 a2 a2 K3 = 1 − 02 − 3 1 − 02 2 R R ! 2 2 a a − 02 − 3 ln 1 − 02 . R R

!−1

+

3 2 (139)

!−1

+

5 2 (140)

Then the products cj Kj are obtained and summed over in j. After some simplifications Eq.(68) is finally obtained.

∗

[email protected]

[1] U. Essmann and H. Tr¨ auble, Phys. Lett. A 24, 526 (1967). [2] U. Essmann and H. Tr¨ auble, Phys. Lett. A 27, 156 (1968). [3] S. J. Bending, Adv. Phys. 48, 449 (1999). [4] L. N. Vu, M. S. Wistrom, and D. J. Van Harlingen, Appl. Phys. Lett. 63, 1693 (1993). [5] P. E. Goa, H. Hauglin, M. Baziljevich, E. Il’yashenko, P. L. Gammel, and T. H. Johansen, Supercond. Sci. Tech. 14, 729 (2001). [6] ˚ A. Olsen, H. Hauglin, T. Johansen, P. Goa, and D. Shantsev, Physica C 408–410, 537 (2004). [7] J. E. Sonier, J. H. Brewer, R. F. Kiefl, D. A. Bonn, S. R. Dunsiger, W. N. Hardy, R. Liang, W. A. MacFarlane, R. I. Miller, T. M. Riseman, D. R. Noakes, C. E. Stronach, and M. F. White, Phys. Rev. Lett. 79, 2875 (1997). [8] C. Niedermayer, E. M. Forgan, H. Gl¨ uckler, A. Hofer, E. Morenzoni, M. Pleines, T. Prokscha, T. M. Riseman, M. Birke, T. J. Jackson, J. Litterst, M. W. Long, H. Luetkens, A. Schatz, and G. Schatz, Phys. Rev. Lett. 83, 3932 (1999). [9] E. W. J. Straver, J. E. Hoffman, O. M. Auslaender, D. Rugar, and K. A. Moler, Appl. Phys. Lett. 93, 172514 (2008), 10.1063/1.3000963. [10] H. Suderow, I. Guillam´ on, J. G. Rodrigo, and S. Vieira, Supercond. Sci. Tech. 27, 063001 (2014). [11] T. Cren, D. Fokin, F. m. c. Debontridder, V. Dubost, and D. Roditchev, Phys. Rev. Lett. 102, 127005 (2009).

35

[12] T. Cren, L. Serrier-Garcia, F. Debontridder, and D. Roditchev, Phys. Rev. Lett. 107, 097202 (2011). [13] I. Lukyanchuk, V. M. Vinokur, A. Rydh, R. Xie, M. V. Milosevic, U. Welp, M. Zach, Z. L. Xiao, G. W. Crabtree, S. J. Bending, F. M. Peeters, and W. K. Kwok, Nat. Phys. 11, 21 (2015). [14] D. J. Thompson, M. S. M. Minhaj, L. E. Wenger, and J. T. Chen, Phys. Rev. Lett. 75, 529 (1995). [15] F. T. Dias, P. Pureur, P. Rodrigues, and X. Obradors, Phys. Rev. B 70, 224519 (2004). [16] G. Li, R. R. Urbano, P. Goswami, C. Tarantini, B. Lv, P. Kuhns, A. P. Reyes, C. W. Chu, and L. Balicas, Phys. Rev. B 87, 024512 (2013). [17] P. Singha Deo, V. A. Schweigert, and F. M. Peeters, Phys. Rev. B 59, 6039 (1999). [18] A. D. Hern´ andez and D. Dom´ınguez, Phys. Rev. B 65, 144529 (2002). [19] E. H. Brandt, Rep. Prog. Phys. 58, 1465 (1995). [20] E. H. Brandt, Phys. Rev. Lett. 78, 2208 (1997). [21] A. A. Abrikosov, Sov. Phys. JETP 5, 1174 (1957). [22] E. B. Bogomolny, Sov. J. Nucl. Phys. 24, 449 (1976). [23] E. Boutsianis, S. Gupta, K. Boomsma, and D. Poulikakos, Ann. Biomed. Eng. 36, 2068 (2008). [24] T. A. Driscoll and L. N. Trefethen, Schwarz-Christoffel Mapping (Cambridge University Press, 2002) cambridge Books Online. [25] L. N. Trefethen, SIAM J. Sci. Sstat. Ccomp. 1, 82 (1980). [26] P. J. Pereira, V. V. Moshchalkov, and L. F. Chibotaru, J. Phys: Conf. Series 410, 012162 (2013). [27] A. A. Vargas-Paredes, M. M. Doria, and J. A. H. Neto, J. Math. Phys. 54, 013101 (2013), 10.1063/1.4773286. [28] B. Fine and G. Rosenberger, The Fundamental Theorem of Algebra, edited by S. Science and B. Media. (Springer, 1997). [29] D. K. Finnemore, T. F. Stromberg, and C. A. Swenson, Phys. Rev. 149, 231 (1966). [30] U. Essmann, Phys. Lett. A 41, 477 (1972). [31] P. Das, C. V. Tomy, S. S. Banerjee, H. Takeya, S. Ramakrishnan, and A. K. Grover, Phys. Rev. B 78, 214504 (2008).

36

[32] D. Bothner, C. Clauss, E. Koroknay, M. Kemmler, T. Gaber, M. Jetter, M. Scheffler, P. Michler, M. Dressel, D. Koelle, and R. Kleiner, Supercond. Sci. Tech. 25, 065020 (2012). [33] I. Gradshteyn and I. Ryzhik, Table of Integrals, Series, and Products, edited by D. Zwillinger and V. Moll (Academic Press, Elsevier, 2014).

37

a)

b)

c)

d)

e)

f)

g)

h)

i)

j)

k)

l)

FIG. 7. The superconducting density, the local magnetic field, the electric current density and the order parameter phase are depicted at three different positions of the vortex inside the cylinder, a0 /R = 0.1 ((a) to (d)), 0.5 ((e) to (h)), and 0.9 ((i) to (l)) of the vortex inside the cylinder. Figures (a), (e) and (i) display the density in (color on line) scheme ranging from low density (cyan) to high density (magenta). The (color on line) scheme is distinct for each of these three figures, as the density varies from zero at the center of the vortex to c20 at the boundary, whose value varies according to the position of the vortex inside the cylinder (see Eqs.(66), (46), (73) and (74)). The local magnetic field is shown in figures (b), (f) and (j), and the color scheme is distinct for each of these three figures. It ranges from the maximum at the center of the vortex (magenta), according to Eqs.(32) and (63), to zero at the boundary (cyan). Figures (c), (g) and (k) show the vortex current around the vortex obtained from Eq.(16). Figures (d), (h) and (l) show the phase of the order parameter. The discontinuity line splitting the white to the black region shows the existence of a single vortex inside the cylinder. The color range from white to black describes 0 to 2π phase.

38