# Flux penetration and expulsion in thin superconducting disks

## Abstract

Using an expansion of the order parameter over the eigenfunctions of the linearized first Ginzburg-Landau (GL) equation, we obtain numerically the saddle points of the free energy separating the stable states with different number of vortices. In contrast to known surface and geometrical barrier models, we find that in a wide range of magnetic fields below the penetration field, the saddle point state for flux penetration into a disk does not correspond to a vortex located nearby the sample boundary, but to a region of suppressed superconductivity at the disk edge with no winding of the current, and which is a nucleus for the following vortex creation. The height of this nucleation barrier, which determines the time of flux penetration, is calculated for different disk radii and magnetic fields.

PACS number(s): 74.24.Ha, 74.60.Ec, 73.20.Dx

The study of magnetic flux penetration and expulsion in type-II superconductors has traditionally attracted much attention in view of important technological and fundamental questions concerning hysteretic behavior and phase transitions in bounded samples. The vortex creation problem is also related to phase transitions in superfluids. It is well known  that for type-II superconductors ($`\lambda /\xi >1/\sqrt{2}`$, $`\lambda `$, $`\xi `$ \- are the penetration and coherence lengths, respectively), the Meissner state becomes energetically unfavorable with increasing magnetic field at $`H=H_{c1}`$ in comparison to the Abrikosov vortex lattice. In a finite system these two states, which correspond to minima of the superconductor free energy, are separated by a barrier. Therefore, a first-order transition between the Meissner and Abrikosov states takes some time which decreases with temperature and approximately follows an Arrhenius law $`\tau exp(U/kT)`$, where $`U`$, $`T`$ are the barrier height and sample temperature, respectively. For $`T=0`$, the Meissner state survives up to the penetration field $`H_p`$ and transits suddenly to the Abrikosov state due to dynamic instability of the order parameter . With decreasing the magnetic field at zero temperature, the vortex state remains stable up to the expulsion field $`H_e<H_{c1}`$ and then goes to the Meissner state due to vortex expulsion.

The origin of barriers for flux penetration and expulsion has been discussed during the last thirty years. According to the Bean-Livingston (BL) model, the surface barrier appears due to a competition between the vortex attraction to the sample walls by its mirror image and its repulsion by screening currents. This model was further developed for: i) cylindrical samples, where the vortex shape was assumed not to be an infinite line but a semicircle , ii) thin disks  and iii) strips , where shielding, due to finite size effects, does not decay exponentially. For samples with a non-elliptical cross section, the geometrical barrier arises because of Meissner screening currents flowing on the top and bottom surfaces of a flat strip . In addition, vortex pinning by defects can play an important role in the delay of vortex expulsion or promotion of vortex penetration. It should be stressed, that the above mentioned barrier models, which are based on the London theory, do not account for the process of vortex formation and describe only the vortex motion far from the sample boundary. The Ginzburg-Landau (GL) theory has been previously applied for the study of barriers only for the 1D cases of narrow wires and rings. The approaches based on solving time-dependent GL equations  allow to treat flux penetration (expulsion) only for magnetic fields higher (lower) than the penetration (expulsion) field. In this Letter, starting from the non-linear GL theory we present an approach for finding the saddle point states in thin disks and calculate numerically the heights of the free energy barriers separating the stable states with different number of vortices.

We consider a superconducting defect free disk with radius $`R`$ and thickness $`d`$ immersed in an insulating media in the presence of a perpendicular uniform magnetic field $`H`$. For thin disks $`Rd\lambda ^2`$ we can neglect the distortion of the magnetic field, which are induced by screening and vortex currents, and write the GL functional as

$$G=G_n+𝑑\stackrel{}{r}\left(\alpha |\mathrm{\Psi }|^2+\frac{\beta }{2}|\mathrm{\Psi }|^4+\mathrm{\Psi }^{}\widehat{L}\mathrm{\Psi }\right),$$
(1)   
where $`G`$, $`G_n`$ are the free energies of the superconducting and normal states; $`\mathrm{\Psi }`$ is the complex order parameter; $`\widehat{L}=(i\mathrm{}\stackrel{}{}e^{}\stackrel{}{A}/c)^2/2m^{}`$ is the kinetic energy operator for Cooper pairs of charge $`e^{}=2e`$ and mass $`m^{}=2m`$; $`\stackrel{}{A}=\stackrel{}{e}_\varphi H\rho /2`$ is the vector potential of the uniform magnetic field written in cylindrical coordinates $`\varphi ,\rho `$; $`\alpha `$, $`\beta `$ are the GL coefficients depending on the sample temperature. Expanding the order parameter $`\mathrm{\Psi }=_i^NC_i\psi _i`$ in the orthonormal eigenfunctions of the kinetic energy operator $`\widehat{L}\psi _i=ϵ_i\psi _i`$  we go from the functional form (1) to the free energy written in terms of complex variables

$$GG_n=(\alpha +ϵ_i)C_iC_i^{}+\frac{\beta }{2}A_{kl}^{ij}C_i^{}C_j^{}C_kC_l,$$
(2)   
where the matrix elements $`A_{kl}^{ij}=𝑑\stackrel{}{r}\psi _i^{}\psi _j^{}\psi _k\psi _l`$ are calculated numerically. Note, that the sample geometry enters in the calculations only through the eigen energies $`ϵ_i`$ and eigenfunctons $`\psi _i`$, which are well known for the disk case . In thin ($`d\xi `$) disks these eigenfunctions have the form $`\psi _{i=(l,n)}=exp(il\varphi )f_n(\rho )`$, where $`l`$ is the angular momentum and the index $`n`$ counts different states with the same $`l`$. In contrast to the approaches , we do not restrict ourselves to the lowest Landau level approximation (i.e. $`n=1`$) and expand the order parameter over all eigenfunctions with $`ϵ_i<ϵ_{}`$, where the cutting parameter $`ϵ_{}`$ is choosen such that increasing it does not influence our results. The typical number of complex components used are in the range $`N=30÷50`$. Thus the superconducting state is mapped into a 2D cluster of $`N`$ classical particles $`(x,y)(Re(C),Im(C))`$, which is governed by the Hamiltonian (2).

To find a saddle point, which presents an extremum of the free energy, we use a technique similar to the eigenvector following method . We start with some set of coefficients $`C`$. In the vicinity of this point the free energy $`\delta G=G(C^n)G(C)`$ can be represented as a quadratic form for small deviations $`\delta =C^nC`$:

$$\delta G=F_m\delta _m^{}+B_{mn}\delta _n\delta _m^{}+D_{mn}\delta _n^{}\delta _m^{}+c.c.,$$
(3)   
where $`F_m=(\alpha +ϵ_i)C_m+\beta A_{kl}^{mj}C_jC_k^{}C_l`$, $`B_{mn}=(\alpha +ϵ_m)I_{mn}+2\beta A_{kl}^{mn}C_kC_l^{}`$, $`D_{mn}=\beta A_{kl}^{mn}C_kC_l`$, and $`I_{mn}`$ is the unit matrix. The quadratic form (3), which is Hermitian, can be rewritten in normal coordinates $`\delta _m=x_kQ_m^k`$ as $`\delta G=2(\gamma _kx_k+\eta _kx_k^2)`$, where $`\gamma _k=Q_m^kF_m`$, the eigenvalues $`\eta _k`$ and eigenvectors $`Q^k`$ are found by solving numerically the following equation

$$\left|\begin{array}{cc}B+Re(D)\hfill & Im(D)\hfill \\ Im(D)\hfill & BRe(D)\hfill \end{array}\right|\left|\begin{array}{c}Re(Q^k)\\ Im(Q^k)\end{array}\right|=\eta _k\left|\begin{array}{c}Re(Q^k)\\ Im(Q^k)\end{array}\right|$$

Moving in the direction with negative free energy gradient $`\gamma _k`$ we will approach a minimum of the free energy corresponding to the ground or a metastable state. In order to find a saddle point we move to a minimum of the free energy in all directions $`x_{kl}=\gamma _k/(ϵ+\eta _k)`$ except one, which has the lowest eigenvalue, and for which we go to a maximum $`x_l=\gamma _l/(ϵ+\eta _l)`$, and find $`C_m^n=C_m+x_kQ_m^k`$ for all $`k`$. The iteration parameter $`ϵ>0`$ controls the convergency, which is always reached for any initial state close enough to a saddle point. Starting from different initial states, for which the coefficients $`C`$ are choosen randomly, we find the saddle points for different magnetic fields (Fig. 1). Due to fluctuations (i.e. thermal, …) the system will be able to reach the saddle point and can then transfer to the other superconducting state. When the magnetic field approaches the expulsion or penetration field, the attraction region of a saddle point state decreases and random searching becomes inefficient. Therefore, to trace the saddle point evolution in the vicinity of the penetration (expulsion) field we start from the saddle point state and increase (decrease) the magnetic field up to the penetration (expulsion) field, when the lowest eigenvalue goes to zero (see Fig. 2, dashed curve).

The spatial distributions of the superconducting electron density $`|\mathrm{\Psi }|^2`$ and velocity $`\stackrel{}{V}=\mathrm{}\stackrel{}{}S2e\stackrel{}{A}/c`$ ($`\mathrm{\Psi }=Fexp(iS)`$) in the saddle point state corresponding to the transition from the Meissner state to the vortex state are depicted in Figs. 3 and 4 for different magnetic fields and disk radius $`R=4.8\xi `$. These figures demonstrate two different stages in the saddle point evolution. Below the penetration field, the saddle point state corresponds to a region of suppressed superconductivity (Fig. 3(d)) with a minimum of $`|\mathrm{\Psi }|^2`$, which is located at the disk boundary. While the minimum value of the order parameter remains different from zero, the vorticity $`L=𝑑\stackrel{}{l}S/\stackrel{}{l}/2\pi `$, where integration is performed along the disk boundary, equals zero and the supervelocity distributions is similar to that of the Meissner state (Fig. 4(a)). When the order parameter reaches its zero value at the nucleation field $`H_n`$ (Fig. 2, solid curve), the vorticity transits suddenly to $`L=1`$. For lower magnetic fields $`H<H_n`$, the saddle point state presents a vortex-like state with closed velocity circulation (Fig. 4(b)). Note, that this transition is not followed by any discontinuity in the free energy or the curvature of the potential curve $`\eta `$. With further decreasing the magnetic field, the saddle point corresponds to a vortex closer to the disk center (Fig. 3(d),(c)). This physical picture of flux expulsion and penetration remains valid for other transitions $`LL+1`$ with different $`L`$, independently of the disk radius and the type (giant vortex or multivortex) of superconducting state.

The free energy, measured in the condensation energy $`G_0=\alpha ^2\pi R^2d/2\beta `$, is shown in Fig. 5 for the saddle point (dotted curves) and stable (solid curves) states for the disk radius $`R=4.8\xi `$. The difference between the free energy of the saddle point state and the nearby metastable state corresponds to the transition barrier shown in the inset of Fig. 5 for transition $`01`$. As seen from Fig. 5, the penetration barrier grows more slowly deep inside the metastable region than the expulsion barrier. Therefore, we expect a larger fluctuation of the penetration field at a finite sample temperature, which agrees with recent experimental observations by Geim .

Below the nucleation field, when the saddle point state is similar to a vortex state, the penetration and expulsion barriers (see Fig. 6) can be estimated from the London theory, which leads to the following expression for the vortex free energy 

$$\frac{G_1}{4\pi G_{}}=ln(\frac{R^2\rho ^2}{r_cR})\mathrm{\Phi }(1\frac{\rho ^2}{R^2})+\frac{1}{4}(\mathrm{\Phi }^2\frac{R^2}{\xi ^2}),$$
(4)   
where $`G_{}=\alpha ^2\xi ^2d/2\beta =G_0\xi ^2/\pi R^2`$, $`\rho `$ is the radial vortex position, $`r_c\xi `$ is the vortex core radius, $`\mathrm{\Phi }=\pi HR^2/\mathrm{\Phi }_0`$ is the unitless magnetic flux, and $`\mathrm{\Phi }_0=hc/2e`$ is the flux quantum. Note that: i) the expulsion field $`H_e=\mathrm{\Phi }_0/\pi R^2`$, ii) the vortex position $`\rho _s=R\sqrt{11/\mathrm{\Phi }}`$ in the saddle point, and iii) the BL expulsion barrier $`U=G_1(\rho _s)G_1(0)=4\pi G_{}(\mathrm{\Phi }1ln\mathrm{\Phi })`$ do not depend on the vortex core energy, which is represented by the first term in Eq. (4). As seen from Fig. 6(a), the London theory prediction for the expulsion barriers are confirmed by the GL theory in the limit of large disks $`R\xi `$. We extended the BL model to arbitrary $`R/\xi `$ by taking into account the spatial nonuniformity of the modulus of the order parameter, which obeys the first GL equation

$$\frac{\mathrm{}^2}{2m^{}}\left(\mathrm{}F(\stackrel{}{}S\frac{2e}{\mathrm{}c}\stackrel{}{A})^2F\right)=\alpha F+\beta F^3,$$
(5)   
with the boundary condition $`(F/\rho )_{\rho =R}=0`$. Following the BL model we assume that the phase distribution is created by a vortex and its mirror image, which are located at the distances $`\rho _v,R^2/\rho _v`$ from the disk center , respectively. Solving Eq. (5) numerically for different vortex positions $`\rho _v`$ we find the expulsion (solid circles) and penetration (open circles) barriers shown in the inset of Fig. 5. Below the nucleation field there is an excellent quantative agreement between our GL theory and this improved BL model. Nevertheless, this model breaks down in the range $`H_n<H<H_p`$. Note, that the barrier height at $`H=H_n`$ increases with disk radius and the role of the nucleation barrier may become even more important in macroscopic systems, where possible 3D (for $`d>\xi `$) and demagnetization (for $`Rd>\lambda ^2`$) effects must also be taken into account. In the unitless variables $`2\pi H\xi R/\mathrm{\Phi }_0`$ ($`\mathrm{\Phi }/\mathrm{\Phi }_0`$), the penetration (expulsion) barriers measured in $`G_{}`$ are proportional to the disk thickness and increase slightly with the disk radius.

In conclusion, we have demonstrated that in a wide range of magnetic fields $`H_n<H<H_p`$ the saddle point state presents a vortex nucleus, which is a region of suppressed superconductivity surrounded by a background of Meissner state, which transits to a vortex state at $`H<H_n`$. We have found the penetration field and the corresponding nucleation barriers for thin disks. For lower magnetic fields $`H_e<H<H_n`$, the saddle point state can be reasonably described by the conventional London theory. We also extended the BL model to finite disk radius.

We gratefully acknowledge discussions with A.K. Geim, R. Blossey, A. MacDonald and V. Moschchalkov. This work is supported by the Flemish Science Foundation (FWO-Vl) through project 5.0277.97 and the “Interuniversity Poles of Attraction Program - Belgian State, Prime Minister’s Office - Federal Office for Scientific, Technical and Cultural Affairs”. One of us (FMP) is a research director with the FWO-Vl.