Universal and shape dependent features of surface superconductivity

TL;DR

This paper investigates surface superconductivity in type II wires using Ginzburg-Landau theory, revealing that boundary shape minimally affects electron density, but curvature influences electron attraction.

Contribution

It demonstrates that boundary curvature affects surface superconductivity, with shape having negligible impact at first order, and curvature playing a significant role at second order.

Findings

Boundary shape has little effect on superconducting electron density.

Mean curvature of the boundary influences electron attraction.

Surface superconductivity persists near the boundary between critical fields.

Abstract

We analyze the response of a type II superconducting wire to an external magnetic field parallel to it in the framework of Ginzburg-Landau theory. We focus on the surface superconductivity regime of applied field between the second and third critical values, where the superconducting state survives only close to the sample's boundary. Our first finding is that, in first approximation, the shape of the boundary plays no role in determining the density of superconducting electrons. A second order term is however isolated, directly proportional to the mean curvature of the boundary. This demonstrates that points of higher boundary curvature (counted inwards) attract superconducting electrons.

Equations

$$\mathcal{G}_{\kappa,h_{\mathrm{ex}}}^{\mathrm{GL}}[\Psi,\mathbf{A}]=\int_{\Omega}\mathrm{d}\mathbf{r}\>\bigg{\{}\left|\left(\nabla+ih_{\rm ex}\mathbf{A}\right)\Psi\right|^{2}-\kappa^{2}|\Psi|^{2}\\ +\mbox{$\frac{1}{2}$}\kappa^{2}|\Psi|^{4}+\left(h_{\rm ex}\right)^{2}\left|\mbox{curl}\,\mathbf{A}-1\right|^{2}\bigg{\}}.$$

$$h_{\mathrm{ex}}=b{\kappa }^2={\epsilon }^{-2}$$

$$\mathcal{G}_{\varepsilon}^{\mathrm{GL}}[\Psi,\mathbf{A}]=\int_{\Omega}\bigg{\{}\bigg{|}\bigg{(}-i\nabla+\frac{\mathbf{A}}{\varepsilon^{2}}\bigg{)}\Psi\bigg{|}^{2}-\frac{1}{b\varepsilon^{2}}|\Psi|^{2}\\ +\frac{1}{2b\varepsilon^{2}}|\Psi|^{4}+\frac{1}{\varepsilon^{4}}\left|\mbox{curl}\,\mathbf{A}-1\right|^{2}\bigg{\}}.$$

$$E_{\epsilon }^{\mathrm{GL}}:=\underset{(\Psi ,\mathbf{A})}{\min }{\mathcal{G}}_{\epsilon }^{\mathrm{GL}}[\Psi ,\mathbf{A}],$$

$$1<b<{\Theta }_0^{-1}\approx 1.69$$

$${\Theta }_0:=\underset{\alpha \in \mathbb{R}}{\min }\underset{\int |u{|}^2=1}{\min }{\int }_0^{+\infty }\mathrm{d}t\left\{|{\partial }_{t}u{|}^2+(t+\alpha {)}^2|u{|}^2\right\}$$

$${\mathcal{E}}_{0,\alpha }^{1D}[f]:={\int }_0^{+\infty }\mathrm{d}t\left\{{\left|{\partial }_{t}f\right|}^2+(t+\alpha {)}^2{f}^2-\frac{1}{2b}\left(2f^2-f^4\right)\right\}$$

$$E_{\varepsilon}^{\mathrm{GL}}=-\frac{1}{2b}\bigg{(}C_{1}(b)\frac{|\partial\Omega|}{\varepsilon}\\ +C_{2}(b)\int_{\partial\Omega}k(s)\mathrm{d}s\bigg{)}+o(1),$$

$$C_1(b)$$

$$C_2(b)$$

$$\mathcal{G}[\psi]:=\int_{0}^{|\partial\Omega|}\mathrm{d}s\int_{0}^{c_{0}|\log\varepsilon|}\mathrm{d}t\left(1-\varepsilon k(s)t\right)\\ \left\{\left|\partial_{t}\psi\right|^{2}+\frac{1}{(1-\varepsilon k(s)t)^{2}}\left|\left(\varepsilon\partial_{s}-it+\mbox{$\frac{1}{2}$}i\varepsilon k(s)t^{2}\right)\psi\right|^{2}\right.\\ \left.-\frac{1}{2b}\left[2|\psi|^{2}-|\psi|^{4}\right]\right\},$$

$$\mathcal{E}^{\mathrm{1D}}_{k,\alpha}[f]:=\int_{0}^{c_{0}|\log\varepsilon|}\mathrm{d}t(1-\varepsilon kt)\big{\{}\left|\partial_{t}f\right|^{2}+V_{k,\alpha}(t)f^{2}\\ -\textstyle\frac{1}{2b}\left(2f^{2}-f^{4}\right)\big{\}},$$

$$V_{k,\alpha }(t):=\frac{(t+\alpha -\frac{1}{2}\epsilon k{t}^2{)}^2}{(1-\epsilon kt{)}^2}.$$

$$\psi (s,t)=f_0(t)e^{-i\frac{{\alpha }_{0}s}{\epsilon }}v(s,t)$$

$${\mathcal{E}}_0[v]\ge {\int }_0^{|\partial \Omega |}\mathrm{d}s{\int }_0^{+\infty }\mathrm{d}t\left(f_0^2(t)+F_0(t)\right){\left|\nabla v\right|}^2,$$

$$F_0(t)=2{\int }_0^t\mathrm{d}\eta (\eta +{\alpha }_0)f_0^2(\eta ).$$

$$\begin{gathered} {\Psi }^{\mathrm{GL}}(\mathbf{r})={\Psi }^{\mathrm{GL}}(s,\tau )\approx f_{k(s)}\left(\frac{\tau }{\epsilon }\right)\exp \left(-i\alpha (k(s))\frac{s}{\epsilon }\right) \\ \exp \left(i{\varphi }_{\epsilon }(s,\tau )\right) \end{gathered}$$

$$\mathcal{E}^{\rm{corr}}_{\alpha_{0}}[f_{0}]:=\int_{0}^{c_{0}|\log\varepsilon|}\mathrm{d}t\>t\bigg{\{}\left|\partial_{t}f_{0}\right|^{2}\\ +f_{0}^{2}\left(-\alpha_{0}(t+\alpha_{0})-\frac{1}{b}+\frac{1}{2b}f_{0}^{2}\right)\bigg{\}}.$$

$$E_{\epsilon }^{\mathrm{GL}}=\frac{|\partial \Omega |E_0^{1D}}{\epsilon }-\frac{\pi }{b}{C}_2(b)+o(1),$$

$$\begin{gathered} {\int }_D\mathrm{d}\mathbf{r}|{\Psi }^{\mathrm{GL}}{|}^4=\epsilon C_1(b)|\partial \Omega \cap \partial D| \\ +{\epsilon }^2{C}_2(b){\int }_{\partial D\cap \partial \Omega }k(s)\mathrm{d}s+o({\epsilon }^2), \end{gathered}$$

$$e_{\epsilon }^{\mathrm{GL}}\approx -\frac{1}{2b{\epsilon }^2}|{\Psi }^{\mathrm{GL}}{|}^4$$

$$\int_{-\infty}^{+\infty}\mathrm{d}t\big{\{}\left|\partial_{t}f\right|^{2}+\min\left((t+\alpha)^{2},(t-\alpha)^{2}\right)f^{2}\\ -\textstyle\frac{1}{2b}\left(2f^{2}-f^{4}\right)\big{\}}.$$

$${\int }_0^{\infty }(t+\alpha )f_{0,\alpha }^2(t)\mathrm{d}t=0$$

$$f_0^{\prime }(0)=0,f_0^{\prime \prime }(0)=f_0(0)\left(b^{-1}-{\alpha }_0^2\right)\ge 0$$

$$-f^{\prime \prime }+(t+\alpha {)}^{2}f=\lambda f$$

$$f^{\prime \prime }(t)=tf(t)\left(t-2\sqrt{{\Theta }_0}\right).$$

Full text

The response of a superconductor to an applied magnetic field is well-known to be both very rich physically and very important practically  deGennes-66 ; Tinkham-75 ; TilTil-90 ; Leggett-06 . In between the Meissner state, where the superconductor totally expels a weak magnetic field, and the normal state at very strong field, where superconductivity is destroyed, different types of mixed states may occur.

For a type-II superconductor in an external magnetic field, there are three critical values of the field marking phase transitions in the state of the material. In increasing order, one should distinguish the first critical field $H_{\mathrm{c}1}$, where bulk superconductivity starts being challenged by the emergence of quantized vortices, from the second and third critical fields, $H_{\mathrm{c}2}$ and $H_{\mathrm{c}3}$. Between these latter values, bulk superconductivity is lost altogether, but Cooper pairs of superconducting electrons can still survive close to the boundary. It is a natural question to wonder how this surface superconductivity phenomenon depends on the shape of the sample.

For a type-I superconductor however the transition from the superconducting to the normal state is more abrupt (first order transition) and only one critical value of the applied field is expected to be relevant. In fact, the behavior of a type-I superconductor can be much richer and show some connections with that of type-II superconductors Metal . In spite of some recent progress Letal , the understanding of superconductivity exactly at threshold between type-I and type-II materials, i.e., for $\kappa=1/\sqrt{2}$ is yet to be understood.

In this note we consider an infinitely long superconducting wire of smooth cross-section $\Omega\subset\mathbb{R}^{2}$, modeled by a 2D Ginzburg-Landau (GL) theory set in $\Omega$, with a perpendicular external magnetic field (thus parallel to the axis of the wire). We consider values of the magnetic field between $H_{\mathrm{c}2}$ and $H_{\mathrm{c}3}$, so that all the physics happens close to the boundary (denoted by $\partial\Omega$) of $\Omega$. We present results in two complementary directions:

$\bullet$ Universality: to leading order in the limit of large GL parameter, the physics does not depend on any detail of the surface’s shape. The order parameter’s amplitude is roughly constant in the direction tangential to $\partial\Omega$.

$\bullet$ Shape-dependence: the next-to-leading order of the energy density is directly proportional to the boundary’s local curvature. Regions of larger boundary curvature (counted inwards) attract more superconducting electrons.

Universality-type results date as far back as Saint James and de Gennes’ seminal papers  JamGen-63 ; deGennes-66 . They argued that, at the onset of surface superconductivity, the problem can be mapped to a one-dimensional one in the direction normal to the boundary. Essentially, we extend this view to the full surface superconductivity regime (where linearization of the model is not valid). Such theoretical predictions (as well as more elaborated models of surface superconductivity) have been verified experimentally several times in the past, for different superconducting materials CKKSSS05 ; MGSJ95 ; Ning-09 ; Retal03 ; RS67 ; SPSKC64 ; TGLS10 ; TLGFMA06 .

As for shape-dependent results, the influence of boundary curvature on the value of the third critical field has been known for some time: when increasing the magnetic field, superconductivity survives longer where the boundary’s curvature is maximal. This has been demonstrated SchPee-98 ; SchPee-99 ; JadRubSte-99 ; FomMisDevMos-98 ; FomDevMos-98 first in the case of domains with corners (i.e. points where boundary curvature jumps) that we do not discuss in detail here111It can be shown CorGia-17a that, as for regular samples, the energy as well as the the density of Cooper pairs is universal and unaffected by the presence of corners to leading order. The corrections due to corners are expected to show up in the next-to-leading order terms CorGia-17b .. The mathematical analysis of the third critical field for domains with corners may be found in BonDau-06 ; BonFou-07 (see (FouHel-10, , Chapter 15) for review and more references). As for smooth domains, the influence of boundary curvature was derived in BerSte-98 ; HelMor-04 (see also (FouHel-10, , Chapters 8 and 13)).

In this note, we summarize rigorous mathematical results obtained in CorRou-14 ; CorRou-16 ; CorRou-16b , discuss their physical interpretation, and complement them with numerical estimates of the shape-dependent contributions to surface superconductivity. The universality-type results we obtained had precursors in the mathematics literature Almog-04 ; AlmHel-06 ; FouHelPer-11 ; FouHel-05 ; FouHel-06 ; FouKac-11 ; LuPan-99 ; Pan-02 (see (FouHel-10, , Chapter 14) for a review). The shape-dependent results however seem to be the first of their kind to be valid in the full regime of surface superconductivity. In view of experimental results Ning-09 on the imaging of the surface superconductivity layer, it would be particularly interesting to measure the influence of the boundary’s curvature on the concentration of Cooper pairs that we predict here.

I Setting and model

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded 2D domain with smooth boundary. Our starting point is the GL energy functional for an infinite wire of cross-section $\Omega$,

[TABLE]

Here $\kappa,h_{\mathrm{ex}}$ are the GL parameter and the strength of the applied magnetic field (parallel to the wire), respectively. The order parameter is denoted $\Psi$ and $h_{\mathrm{ex}}\mathbf{A}$ is the vector potential of the induced magnetic field.

It is well-known that surface superconductivity occurs in the regime $h_{\mathrm{ex}}\propto\kappa^{2}$. To study this regime it will be convenient to change units, setting

[TABLE]

for two new parameters $b,\varepsilon$. In these units we get

[TABLE]

We consider the associated ground state problem

[TABLE]

and study the limit $\varepsilon\to 0$ with

[TABLE]

where de Gennes’ constant

[TABLE]

is the minimal ground state energy of the shifted harmonic oscillator on the half-line. In terms of the GL parameter $\kappa$ the limit $\varepsilon\to 0$ corresponds to the extreme type II case $\kappa\to\infty$. In view of (I.2) the condition (I.5) corresponds to asking that the external magnetic field lies strictly between $H_{\mathrm{c}2}$ and $H_{\mathrm{c}3}$. The numerical range for $b$ in (I.5) is set by the lowest eigenvalues for Schrödinger operators with constant magnetic field equal to $1$, in the plane and the half-plane respectively. Note that $\varepsilon$ is the only small parameter at our disposal (stressing this point is the main reason for changing units). Only when $b\to\Theta_{0}^{-1}$ (i.e. at the upper limit of the surface superconductivity regime) does the problem become linear.

In all the sequel we take the above model as our starting point. Clearly this demands that the full BCS theory that would be appropriate reduces to the simpler GL model. The main condition is that one be sufficiently close to the critical temperature (see Gorkov-59 ; deGennes-66 for standard references and FraHaiSeiSol-12 ; FraHaiSeiSol-16 for a more recent and mathematical treatment). Reduction from a 3D to a 2D model is valid in two situations:

$\bullet$ a very long wire with magnetic field parallel to the wire;

$\bullet$ a very thin film with perpendicular magnetic field.

In the latter case, one should directly impose that the induced magnetic field coincides with the external one, $\mbox{curl}\,\mathbf{A}=1$ (see e.g. AlaBroGal-10 and references therein), which only simplifies the following arguments.

II Shape independence of the Ginzburg-Landau energy

The first result we discuss is a two-term expansion for the GL ground state energy. All its parameters can be calculated by minimizing the reduced 1D energy functional

[TABLE]

both with respect to the normal density profile $f:\mathbb{R}^{+}\mapsto\mathbb{R}$ and the total phase circulation $\alpha\in\mathbb{R}$. This is a generalization of the linear problem (I.6), introduced in Pan-02 . We denote the minimal value of the 1D effective energy by $E^{\rm 1D}_{0}$ and by $f_{0},\alpha_{0}$ a minimizing pair. The above problem relates to the original one via scaling lengths by a factor of $\varepsilon$ (typical thickness of the superconducting layer) in the direction normal to the boundary and assuming a particular ansatz for the GL minimizer, discussed below.

Concerning the ground state energy, our results in CorRou-14 ; CorRou-16 ; CorRou-16b can be summarized as follows:

Theorem II.1** (GL energy expansion).**

For any $1<b<\Theta_{0}^{-1}$, in the limit $\varepsilon\to 0$,

[TABLE]

where

[TABLE]

Sketch of proof.

The methods of (FouHel-10, , Chapters 11 and 12) allow to prove that the order parameter and energy density are very small at distances larger than $\mathcal{O}(\varepsilon)$ from the boundary. Also, the induced magnetic field does not differ significantly from the applied one. Changing to scaled boundary coordinates $(s,t)=$(tangential, $\varepsilon^{-1}\times$ normal), one is then lead to consider a scaled energy functional

[TABLE]

where $c_{0}$ is a fixed, essentially arbitrary, constant. The local curvature (counted inwards) $s\mapsto k(s)$ of the boundary arises from the change of variables. The upshot is that, in first approximation (setting $\varepsilon=0$ in (II.5)), the curvature of the sample’s boundary can be neglected. The whole sample is then mapped to a half-plane. To proceed further and get more precise results, one may instead approximate the boundary at each point by the auscultating circle of radius $k(s)^{-1}$.

The particular case $|k|=R^{-1}=\mathrm{const.}$ for the above functional corresponds to the surface energy in a disk of radius $R$ (or the exterior thereof if $k<0$). Since all terms involving the $s$-coordinate in (II.5) come with $\varepsilon$ pre-factors, it makes sense to guess an ansatz of the form $f(t)e^{-i\frac{\alpha}{\varepsilon}s}$ with a constant $\alpha$ to solve the disk problem. This leads to considering the reduced functional

[TABLE]

where

[TABLE]

This is a generalization of (II.1) to the case of non-zero curvature. The main step of the proof is to vindicate the exactness of the ansatz $\psi(s,t)=f(t)e^{-i\frac{\alpha}{\varepsilon}s}$ for the disk problem ($k=\mathrm{const.}$) in (II.5).

For simplicity we explain this in the case $k=0$, corresponding to the degenerate case of a half-plane sample. The method is to write a tentative minimizer in the form

[TABLE]

and bound from below the excess energy due to a possibly non constant $v$ in the manner

[TABLE]

with a potential function

[TABLE]

The “cost function” $f_{0}^{2}(t)+F_{0}(t)$ can be interpreted as a lower bound to the kinetic energy density generated by a vortex located at distance $t$ from the boundary. In CorRou-14 we prove that $f_{0}^{2}(t)+F_{0}(t)\geq 0$ for any $t$, provided $1\leq b\leq\Theta_{0}^{-1}$, thus concluding that a vortex, or more generally any non constant $v$, would increase the energy. The ground state energy of the half-plane problem is thus simply given by $|\partial\Omega|E^{\rm 1D}_{0}$. The disk case follows similar considerations, with significant but technical additional difficulties.

Relying on the previous observations, we get a kind of “adiabatic” decoupling where the order parameter is for any $s$ in the ground state of the problem in the direction perpendicular to the boundary. This reduces matters to looking for the energy $E^{\mathrm{1D}}_{\star}\left(k\right)$ obtained by minimizing (II.6) with respect to both $\alpha\in\mathbb{R}$ and $f:\mathbb{R}^{+}\mapsto\mathbb{R}$. The rationale is that a true GL minimizer can, in a suitable gauge, be approximated in the manner

[TABLE]

with $f_{k(s)},\alpha(k(s))$ a minimizing pair for the energy functional (II.6) at curvature $k=k(s)$ and $(s,\tau)=$ (tangential coordinate, unscaled normal coordinate). Here $\phi_{\varepsilon}$ is a phase factor accounting for a uniform vortex density proportional to $\varepsilon^{-2}$ in the bulk of the sample, which is required to compensate the huge magnetic field inside the sample. Such a large circulation gets corrected by a relatively small amount $-\alpha(k(s))\varepsilon^{-1}$ in order to optimize the boundary energy. This form includes subleading order corrections, and its main feature is the decoupling of scales/coordinates.

The final expression (II.2) follows from first order perturbation theory applied to the 1D functional (II.6). Since we work in the regime $\varepsilon\to 0$, one can expand $E^{\mathrm{1D}}_{\star}\left(k(s)\right)$ in powers of $\varepsilon$. The leading order is given by the $k=0$ functional and the first correction by $-\varepsilon k$ times

[TABLE]

That $C_{1}(b)$ is equal to the expression in the right-hand side of (II.3) is a simple consequence of the variational equation satisfied by $f_{0}$, while showing that (II.9) coincides with (II.4) divided by $2b$ requires a few non trivial computations that we skip for brevity. ∎

For $b\to\Theta_{0}^{-1}$, the half-plane problem can be linearized, and the proof is much simpler. Indeed, for the linear problem, the ansatz $f_{0}(t)e^{-i\frac{\alpha_{0}}{\varepsilon}s}$ is clearly optimal (as used first in JamGen-63 and further justified in the mathematics literature (FouHel-10, , Chapter 13)). This is because the Hamiltonian of the linear problem commutes with translations in the tangential coordinate, so that one can perform a Fourier series decomposition in this direction. This nice structure is broken by the nonlinear term. In this respect, one may see Theorem II.1 as an extension to the nonlinear regime $b<\Theta_{0}^{-1}$ of the linear analysis pioneered in JamGen-63 .

The energy asymptotics (II.8) also implies estimates of the amplitude of the order parameter and the phase circulation around the boundary. Indeed, $|\Psi^{\mathrm{GL}}|^{2}$ can be shown to be close to $|f_{0}\left(\textstyle\frac{\tau}{\varepsilon}\right)|^{2}$ pointwise in the surface superconductivity layer. This indicates that no defect can be present therein, in particular no vortex.

As we have anticipated, the curvature effects emerge in the next-to-leading order corrections to the density of superconductivity. However, in the energy expansion (II.2), the second term is universal and independent of the sample’s shape: by the Gauss-Bonnet theorem, we have that

[TABLE]

for the integral of the curvature in (II.2) is just $2\pi$ times the Euler characteristic of the domain $\Omega$, equal to $1$. In the next section we turn to describing the most explicit form of curvature dependence we are able to derive.

III Shape dependence in the distribution of superconductivity

The main result we obtained in CorRou-16b is stated in terms of the distribution of $|\Psi^{\mathrm{GL}}|^{4}$, the square of the Cooper pairs’ (normalized) density. For obvious physical reasons it would be desirable to obtain instead an estimate of $|\Psi^{\mathrm{GL}}|^{2}$, but this does not follow straightforwardly from our methods.

Theorem III.1** (Curvature dependence of $\Psi^{\mathrm{GL}}$).**

Let $\Psi^{\mathrm{GL}}$ be a GL minimizer and $D\subset\Omega$ be a measurable set intersecting $\partial\Omega$ with right angles. For any $1<b<\Theta_{0}^{-1}$, in the limit $\varepsilon\to 0$,

[TABLE]

using the notation of Theorem II.1.

Sketch of proof.

Theorem II.1 is actually proved in a local fashion. Denoting $e_{\varepsilon}^{\mathrm{GL}}(\mathbf{r})$ the Ginzburg-Landau energy density of a ground state, we in fact obtain, for any domain $D\subset\Omega$, a version of (II.2) for the energy density integrated over $D$. The main idea is then that, neglecting the magnetic kinetic energy (last term in (I.3)), we essentially have

[TABLE]

and thus the result. The above (approximate) identity follows by multiplying the GL equation by $\Psi^{\mathrm{GL}}$ and integrating the result over a domain $D\subset\Omega$. Boundary terms however arise due to integrating the kinetic energy by parts, involving the derivative of $\Psi^{\mathrm{GL}}$ in the direction normal to $\partial D$. They are negligible for domains intersecting $\partial\Omega$ with $\pi/2$ angles because, for such domains, this direction is tangential to the boundary $\partial\Omega$ of the sample. The variations of $|\Psi^{\mathrm{GL}}|$ in this direction are very slow, as can be guessed from (II.8). ∎

Theorem III.1 is almost as far as we can get with an analytic approach. It is however not quite obvious to determine the sign of $C_{2}(b)$ from the expression we obtained, i.e. to conclude whether larger curvature favors superconductivity, or the other way around. The expression (II.4) can however be evaluated by numerically solving a rather simple variational problem. We next report on the conclusions of such an investigation:

Claim III.1** (Curvature favors superconductivity).**

For any $1<b<\Theta_{0}^{-1}$, we find that $C_{2}(b)>0$. Consequently, the distribution of surface superconductivity along the boundary is an increasing function of the local curvature (counted inwards).

Argument. Everything boils down to solving the 1D minimization problem (II.1) and inserting the results in (II.4). We proceed as follows:

$\bullet$ Starting from an initial value for $\alpha$, we minimize $\mathcal{E}^{\mathrm{1D}}_{0,\alpha}[f]$ as a function of the normal profile $f$. We find it simpler to consider a problem on the full real line, which is equivalent and avoids boundary conditions issues. We thus reflect the harmonic potential and minimize

[TABLE]

Up to the fact that $f$ should not be normalized, this can be seen as the computation of the ground state profile of a Bose-Einstein condensate in the double-well potential $\min\left((t+\alpha)^{2},(t-\alpha)^{2}\right)-b^{-1}$. We rely on the freely available Matlab toolbox GPELab GPELab ; AntDub-14 ; AntDub-15 , removing the normalization step imposed therein.

$\bullet$ Having computed the minimizing profile $f_{0,\alpha}$ at fixed $\alpha$, the minimization in $\alpha$ at fixed $f$ is very easy: the Feynman-Hellmann principle yields

[TABLE]

and we use this to update the value of $\alpha$.

$\bullet$ We again minimize the energy with respect to $f$, using the updated value of $\alpha$, and iteratively continue the procedure until convergence is reached.

$\bullet$ A few benchmarks can be used to ascertain the validity of the calculation. The second identity in (II.3) should always hold for the converged quantities. Next, the two expressions (II.4) and $C_{2}(b)=2b\mathcal{E}^{\rm{corr}}_{\alpha_{0}}[f_{0}]$ should always coincide. Finally, we checked that the cost function $f_{0}^{2}(t)+F_{0}(t)$ is always positive for $1<b<\Theta_{0}^{-1}$.

For $b$ very close to $\Theta_{0}^{-1}$, the numerics become rather unstable, for then $f_{0}$ takes very small values. This regime would be better analyzed by solving the lowest eigenvalue problem for the linearized equation. We have not pursued this beyond the available studies (e.g. Bonnaillie-08 ), for the result we aim at in this regime is already included in (CorRou-16b, , Lemma 2.3). Here, perturbation theory in $|b-\Theta_{0}^{-1}|$ proves that there exists $\delta>0$ such that $C_{2}(b)>0$ for any $b$ satisfying $\Theta_{0}^{-1}-\delta\leq b\leq\Theta_{0}^{-1}$. This does not settle the question analytically for the full parameter regime of interest but allows to bypass the numerical instability issue for $b\to\Theta_{0}^{-1}.$

The values of the energy corrections $C_{2}(b)$ we found are plotted in Figure 1. The result is clearly positive for $1<b<\Theta_{0}^{-1}$, as claimed in the statement. The numerical problem can also be tackled for $b<1$, but this is less relevant physically, for then bulk superconductivity starts challenging surface superconductivity, which is not taken into account in the reduced model we solve. ∎

IV Remarks on the density profile

As an illustration, we provide plots of the optimal density profile at $b=1$ and $b=1.5$ in Figure 2. According to the discussion in Section II, they give a very good approximation for the dependence of the (amplitude of the) order parameter in the direction normal to the boundary. In this respect, it is interesting to note that the density is not monotonously decreasing as a function of the distance to the boundary. This can be seen analytically because

[TABLE]

where the first equation is the Neumann boundary condition and the second uses the variational equation for $f_{0}$ (Euler-Lagrange equation associated to (II.1)) together with the known fact (CorRou-14, , Proof of Lemma 3.3) that $f_{0}^{2}(0)=2-2b\alpha_{0}^{2}\geq 0$. Thus, $f_{0}$ increases close to the origin, whereas it must ultimately decay, so that it has to reach at least one global maximum outside the origin. Since superconductivity is due to boundary effects in this regime, it is rather counter-intuitive that the maximum density of Cooper pairs is not attained exactly at the boundary.

In our simulations, this drop of density at the boundary is hardly visible for values of $b$ close to the third critical field (see the curve for $b=1.5$ in Figure 2) but becomes more pronounced for $b$ close to $1$, i.e. close to the second critical field (as per (I.2)). It might thus be a precursor of the emergence of bulk superconductivity in the sample, when the second critical field is crossed from above.

Note that the maximum density not occuring at the boundary is a nonlinear effect: the linearized model valid close to $H_{\mathrm{c}3}$ leads to a monotonously decreasing density. Indeed, close to $H_{\mathrm{c}3}$ the density is found by minimizing (II.1) without nonlinear term, under a mass unit mass constraint JamGen-63 ; FouHel-10 . This leads to the variational equation

[TABLE]

and it is known (FouHel-10, , Chapter 3) that the value of $\alpha$ giving the minimal energy is such that $\alpha=-\sqrt{\lambda}=-\sqrt{\Theta_{0}}$, so that one gets

[TABLE]

Since the minimizing $f$ must be positive, the above changes sign only once on the positive half-axis, and this is incompatible with a global maximum occurring away from the origin.

Since for $b\to\Theta_{0}^{-1}$ the minimizer of (II.1) gets closer and closer to a multiple of the linear solution, it is natural that the density maximum away from the origin is less visible in that limit.

V Conclusions

We studied the surface superconducting regime of the usual two-dimensional GL theory for a very elongated cylinder in a magnetic field perpendicular to its cross-section (also valid for a thin film with perpendicular field). Our results hold in the limit of a large GL parameter, in the full regime of applied field between the second and third critical values $H_{\mathrm{c}2}$ and $H_{\mathrm{c}3}$. In particular, we showed that:

1. To leading order, the distribution of superconducting electrons along the sample’s boundary is independent of its shape. The boundary can in first approximation be treated as a single straight line, and the density of Cooper pairs has significant variations only in the direction perpendicular to it.

2. The subleading correction to the distribution is directly proportional to the local value of the boundary’s curvature. A larger curvature (counted inwards) corresponds to a larger density of Cooper pairs. This follows by locally approximating the boundary by an auscultating circle.

These results are in accord with what was previously obtained from a linearized analysis close to the upper critical field $H_{\mathrm{c}3}$ (see FouHel-10 for review). They however constitute, as far as we know, the first complete results of their kind valid for all the relevant values of the applied field, where a fully nonlinear analysis is required.

We remark that our approach requires the boundary curvature to be smooth, i.e. that its variations happen on a scale much larger than the coherence length (proportional to $\kappa^{-1}\propto\varepsilon$). If this assumption is violated at some points of the boundary, it would be appropriate to treat these as corners (where the curvature makes sharp jumps). Such an analysis will be provided elsewhere CorGia-17a ; CorGia-17b .

As directions for future investigations we mention the case where the magnetic field is not strictly perpendicular to the cross-section (or the film), and a more complete 3D treatment in case the wire is bent for example. In the latter case, it is known that superconductivity appears first, when decreasing the external field below $H_{\mathrm{c}3}$, in regions where the magnetic field is tangent to the boundary FouHel-10 . In the regime of our interest, superconductivity should however be present along the whole surface. The first question to tackle then is certainly the influence of the angle between the magnetic field and the cross-section. Results in this direction may be found in FouKacPer-13 . In a full 3D treatment, whether the local (Gaussian) curvature of the boundary rather than the 2D cross-sectional curvature influences the Cooper pairs’ distribution remains an open question.

Acknowledgments: M.C. acknowledges the support of MIUR through the FIR grant 2013 “Condensed Matter in Mathematical Physics (Cond-Math)” (code RBFR13WAET). N.R. acknowledges the support of the ANR project Mathostaq (ANR-13-JS01-0005-01). B.D. was supported by the CNRS Inphynity project MaBoLo. We thank Romain Duboscq for helpful advice regarding the GPELab toolbox.