Vortex nucleation barrier in superconductors beyond the Bean-Livingston approximation: A numerical approach for the sphaleron problem in a gauge theory

The knowledge of vortex nucleation barriers is crucial for applications of superconductors, such as single-photon detectors and superconductor-based qubits. Contrarily to the problem of finding energy minima and critical fields, there are no controllable methods to explore the energy landscape, identify saddle points, and compute associated barriers. Similar problems exist in high-energy physics where the saddle-point configurations are called sphalerons. Here, we present a generalization of the string method to gauge field theories, which allows the calculation of energy barriers in superconductors. We solve the problem of vortex nucleation, assessing the effects of the nonlinearity of the model, complicated geometry, surface roughness, and pinning.

In type-II superconductors, the Meissner state, characterized by the total magnetic field expulsion in the bulk, is stable up to the lower critical field H c1 . Above it, quantum vortices appear. However, the Meissner state can survive as a metastable state, causing the phenomenon known as magnetic superheating [1]. The presence of fluctuations can trigger the spontaneous decay of the metastable Meissner state through vortex formation. This phenomenon is the effect of a surface barrier, which hinders vortex nucleation from the sample boundaries [1; 2]. The barrier disappears when the applied magnetic field exceeds a critical value, the nucleation field. In contrast to the potential barrier, the nucleation field is not only amenable to analytical treatment [3][4][5][6][7][8][9][10][11], but can be calculated with a controlled numerical accuracy from iterative simulations [12][13][14].
The full problem of the vortex entry barrier is not solved even for the simplest case of a semi-infinite superconductor with an ideal surface. The Bean-Livingston estimate [2] relies on the London model, but the result depends on the choice of cutoff, as emphasized by the authors themselves. The Bean-Livingston approach has been extended to other geometries [15; 16].The effect of surface roughness is considerably more complicated. Only some approximate analytical approaches within the London model were advocated in Ref. [17]. However, roughness alters the core structure and therefore a controllable method to solve the full nonlinear model is needed.
Fundamentally, the problem of a vortex entry barrier consists of finding the sphaleron, i.e., the saddle point which separates two stable states in a gauge theory [18][19][20][21][22]. In the case of vortex nucleation, it is an energy maximum of the minimum energy path between two states with different phase windings. * alben@kth.se † andrea.maiani@nbi.ku.dk The knowledge of potential barriers for vortex nucleation is crucial for superconductor applications. The applications include current transmission, where the dissipationless state is lost if free vortices form, and superconductive rf cavities for particle accelerators where superconductive magnets often operate in a superheated state [11]. Recently, the problems of vortex entry barriers appeared in quantum technologies. In superconducting single-photon detectors, it is believed that the principle of operation consists in the creation of a currentcarrying state with a small potential barrier for vortex entry so that a single photon creates a vortex, hence a detectable signal [23; 24]. However, small barriers yield spontaneous vortex nucleation caused by fluctuations, resulting in dark counts. The ability to calculate potential barriers would allow designing devices with significantly improved performances. Likewise, the knowledge of vortex entry barriers is crucial to design superconducting topological qubits [25][26][27][28][29], where the geometry of the device plays a central role. These devices often operate at temperatures near absolute zero. However, while a microscopically derived Ginzburg-Landau model applies only close to the critical temperature, many aspects of low-temperature vortex physics may under certain conditions be fittable by effective Ginzburg-Landau-type models [30]. In this Rapid Communication, we generalize to gauge theories the simplified string method [31]. This allows us to perform surface barrier calculations in type-II superconductors for vortex nucleation (∆F n ) and escape (∆F e ), by computing the minimum energy path of the transition in the Ginzburg-Landau theory. The results differ from the Bean-Livingston theory [2], which neglects the vortex-core and nonlinear effects. More importantly, the method allows taking into account the effect of surface roughness and the presence of pinning in the computation of the vortex nucleation barrier and the superheating field. The Ginzburg-Landau free-energy functional describing the superconductor in dimensionless units is (1) The complex field ψ = |ψ|e iθ describes the state of the superconductor. The vector potential a is related to the magnetic field by b = ∇×a. The coefficient κ = λ ξ is the Ginzburg-Landau parameter, i.e., the ratio of the magnetic field penetration depth λ and the coherence length ξ (we use the definition of coherence length with a factor √ 2 absorbed as in Ref. [32]). The external magnetic field H is expressed in units of the thermodynamic critical field H c , i.e., h = H/H c . The free energy f = F/F 0 is expressed in units of F 0 , which in SI units is F 0 = µ 0 H 2 c λ 2 d, where d represents the thickness of the sample and µ 0 is the vacuum permeability. Quantities in capital letters are intended in SI units, while variables in lower case are dimensionless.
In a finite system, the Meissner state can survive in a meta-stable way for higher fields than H c1 , i.e., up to the spontaneous nucleation field H n . In the absence of fluctuations, only when this field is exceeded, vortices nucleate from the boundaries. To introduce a vortex in the system, when H < H n , we need to overcome an energy barrier due to the surface. The calculation of this barrier is a nonlinear problem that is not amenable to analytical treatment. In this work, we generalize a numerical approach, which was originally developed within the molecular dynamics field to study transitions between two metastable states through the identification of a minimum energy path [33]. The path is considered in the configuration space of the system, and it is parametrized by the transition coordinate s ∈ [0, 1]. If one denotes by q(s) = [a(s, r), ψ(s, r)], the state of the system, then q(0) and q(1) are two equilibrium solutions corresponding to the minima of the Hamiltonian in Eq. (1). By varying s from 0 to 1, we observe the transition of the system from the initial state to the final state. One can assume that a potential force g ≡ −∇f = − δf δa , δf δψ * acts on each point of the curve q(s). The minimum energy path is a trajectory in which the force g acting on each point is uniquely directed along the tangent vector ∂q/∂s , i.e., We emphasize that the optimal path defined by Eq. (2) does not correspond to the real time dynamics. It describes the most energetically favorable transformation undertaken by the system for the transition between the initial and the final state.
To identify the minimum energy path, we started from the well-known simplified string method [31] and developed a variant applicable to gauge field theories [34]. This algorithm evolves an initial guessed path in the configuration space, towards the minimal energy one. Consider a situation where we start from a Meissner state and end in the one-vortex state. To construct the initial guess, we used an ansatz for the single winded vortex state. In the initial state, for s = 0, the vortex is outside of the domain and in the final state, for s = 1, it lies in the origin. To each value of the transition coordinate s, there corresponds a particular configuration of the system. Hence, s is not equivalent to the position of the center of the vortex because static vortex-core deformations correspond in general to different values of s. This method allows us to solve the full nonlinear problem in contrast to previous approaches, and obtain exact and quantitatively valid results. A previous study of the vortex entry barrier, based on the London model, was carried out by Bean-Livingston [2]. However, this introduced the uncontrollable approximation of considering the vortex core as a rigid cylinder of radius ξ, which neglects the physics of the core and the nonlinear effects.
The energy dependence F (s) is a function with one or more maxima. At the summit of the potential barrier, the configuration of the system is a saddle point of the free-energy functional of Eq .(1), which in the context of gauge theories is called a sphaleron.
Once the minimum energy path is computed, we can define the nucleation barriers ∆F n = F sphaleron −F Meissner for a given magnitude of the external magnetic field H. Consequently, the nucleation field H n is the external field needed to nullify the nucleation barrier, i.e., ∆F n (H n ) = 0 [35]. Analogous definitions can be used for the escape barrier, i.e., the energy needed to expel one vortex from the superconductor. A quantitative description of nucleation and escape barriers is given in the Appendix. The minimum energy path contains more information than the height of the energy barrier as those paths are most likely to be followed in the nucleation process. This path helps us understand in detail how transitions between metastable states occur. In our framework, we can calculate the energy barrier free from any approximations, except fully controlled numerical errors. We begin by considering a vortex entry in a 2D superconductor with flat surfaces. We find that the vortex always enters, following the minimum energy path, from the sides of the superconductor and never from the corners. Figure 1 shows the free-energy profile of the process. The barrier presents a single maximum corresponding to the sphaleron. The substantial core deformation confirms that this kind of problem is not in general treatable with controlled accuracy in the London limit.
For complex materials, the surface can have different roughness, doping and oxidation, which strongly affects the vortex entry process. The knowledge of how impurities influence vortex physics is crucial for applications [36].
Let us consider randomly distributed inhomogeneities, The energy dependence and the saddle-point (sphaleron) configuration for the process of vortex entry from a straight edge in a two-dimensional superconductor. At the peak of the energy barrier, we have a substantial vortex-core deformation which cannot be described in the London limit, where ψ is assumed to be constant. This shows how this kind of process, in general, cannot be treated in the London model, i.e., within the Bean-Livingston approach. The variable s parametrizes the minimum energy path. In this example H = 0.6Hc and κ = 2. The inset shows |ψ| in the region when the nucleation takes place.
with a decreasing density as we enter the sample, as shown in Fig. 2. To model inhomogeneities in a superconductor, we follow a procedure similar to the one outlined in Refs. [37][38][39], where we modify the quadratic term in Eq. (1) accordingly, where σ(r) is the inhomogeneity distribution. We tested different models for the inhomogeneity distribution with similar results. In the free-energy profile, depicted in Fig. 2, there is a global minimum in between the saddle points. This has the effect of increasing the nucleation barrier for a second vortex. In fact, in a clean sample, the entry barriers for the first and the second vortex are nearly the same. However, in the presence of inhomogeneities, we can calculate how the first vortex is pinned near the surface with the effect of increasing the nucleation barrier for the second one, as shown in Fig. 3. Hence, with this method, it is possible to accurately predict and design impurity density profiles able to protect the sample from vortex entry by pinning vortices in an area near the surface.
Real samples in general have rough surfaces. It is empirically known that in the presence of rough surfaces, the vortex entry barrier is altered [1]. A rough surface represents a challenging problem for the calculation of both superheating fields [14] and barriers [17] FIG. 2. Trapped vortex and corresponding energy profile for the minimum energy path in the case of a superconductor with inhomogeneities concentrated near the edge. The impurity layout σ is normally distributed, and it is decreasing to zero as we proceed towards the bulk. On top of it, we plot isolines of the order parameter modulus |ψ|. For the color notation of isolines, see the color bar in Fig For a quasi-two-dimensional superconductor, we can describe a rough edge as a sequence of geometrical defects, forming dents in the sample, as sketched in Fig. 4. The method is also straightforwardly applicable to the three-dimensional case.
The vortex entry occurs from a single dent, where the barrier is suppressed the most. Therefore, it is relevant to study how the geometrical properties of one dent affect the nucleation barrier ∆F n . As shown in Fig. 4 (b), we can characterize the dent by its depth l which can be comparable to or bigger than the characteristic length scales of the model (coherence length ξ and penetration In these three cases the profile is, going from left to right, comparable to the coherence length ξ, comparable to the penetration depth λ, and bigger than λ. (b) Model of a single dent: The depth is parametrized by l, while its sharpness by the angle φ. The black color indicates the superconductor (SC).
depth λ), and by an angle φ which determines its sharpness. Figure 5 compares the nucleation barriers ∆F n /F 0 Nucleation barrier computed for a dent with depth l as a function of the angle φ as shown in Fig. 4 (b). The system parameters are κ = 5 and H = 0.4Hc. The continuous lines are obtained by a quadratic regression (OLS) and agree with the input data. In particular, for the bigger dent, the nonlinear coefficient is negligible and the behavior is substantially linear.
for different dent depths l, as a function of the sharpness angle φ. For φ = π we have the limit of a perfectly straight edge. The green line shows the nucleation barrier for l = ξ. We notice that the energy barrier is unaltered by the presence of the dent with a comparable or smaller roughness. Hence, roughness profiles with depths of the order of a coherence length have no effect on the nucleation barrier, independently of the sharpness of the dents. As l increases, the barrier suppression becomes substantial, depending on the sharpness φ. For l = λ the dependence on φ is non-linear, whereas for l > λ (in our case l = 4λ), the barrier decreases linearly with φ.
In conclusion, we formulated a method to compute the minimum energy path in a gauge theory. We applied this method to solve the problem of vortex nucleation in the Ginzburg-Landau model of superconductivity. This is, a solution of the full nonlinear problem of the vortex entry barrier. The method allows calculating vortex entry barriers in the presence of surface roughness and impurities. We showed that previously developed models to estimate the vortex entry barrier, based on the London approximation, are, in general, inadequate to describe the nucleation process even for a perfectly flat surface. We find that the surface roughness at the scale of coherence length affects the barriers insignificantly, and thus for superconductors with short-scale surface roughness, the main mechanism for the surface barrier reduction is the modulation of the superfluid density at larger length scales near the surface. The method straightforwardly applies to three-dimensional configurations and to the geometries where demagnetization fields are important [27]. Such vortex entry barriers are of key importance in the design of quantum devices such as single-photon detectors and qubits, as well as superconducting rf cavities, transmission lines, and magnets. Finally, the method is straightforwardly applicable to other gauge theories, different boundary conditions [40] which includes microscopic models of superconductivity, and sphaleron identification in high-energy physics.
We thank Mats Barkman, Robert Vedin, and Oskar Palm for useful comments. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Quadro P6000 GPU used for this research. Published with the support from the Längman Culture Foundation. The effective free energy Ginzburg-Landau functional describing a superconductors reads: The complex field ψ is the superconducting order parameter, coupled to the electromagnetic field through the vector potential A. Here d is the effective thickness of the sample with a cross-section denoted by Ω, m is the mass of a superconductive carrier (i.e., the Cooper pair), q is the coupling constant with the gauge field, α and β are two temperature-dependent parameters, H is the uniform applied field while and µ 0 are, respectively, the Planck constant and the vacuum permeability. We consider the case H ê z .
To make the model dimensionless we scaled the variables as follows: and the spatial coordinates become x → λx and y → λy. ψ 0 = −α β is the uniform state order parameter, is the thermodynamic critical field, F 0 = µ 0 H 2 c λ 2 d is the characteristic energy value. We also introduce the Ginzburg-Landau parameter κ = λ ξ where ξ = 2 √ −αm is the coherence length. With this choice of definition for the coherence length [32], the critical coupling which separates type-I and type-II superconductors corresponds to κ = 1.
The resulting dimensionless Ginzburg-Landau model is expressed, up to constant terms, by the following free energy functional: where h = H/H c . Quantities in capital letters are measured in SI units, while variables in lower case are dimensionless. The Ginzburg-Landau model features local U (1) gauge symmetry. This means that a and the phase of ψ do not represent physically observable quantities. Instead, the observable fields of the system are the magnetic field b = ∇ × a and the superconductive current density j = − 2 κ Im ψ * −i 2 κ ∇ + a ψ , which are invariant to the local gauge transformation: In a type-II superconductor, when H < H c1 the equilibrium state is in the Meissner phase, no vortices are present in the bulk, and thus the magnetic field is completely screened by the superconducting current at surface. For H > H c1 , the system is in the Shubnikov phase in which vortices are present. However, due to the presence of the surface barrier the Meissner state can survive in a meta-stable way for higher fields than H c1 , i.e. up to the spontaneous nucleation field H n . In the absence of fluctuations, only when this field is exceeded the barrier is suppressed and vortices can nucleate freely from the boundaries. Moreover, if the applied external field is higher than the spontaneous escape field H e vortices cannot spontaneously escape from the sample without overcoming the surface barrier. This means that in the region H e < H < H n Shubnikov states are metastable as the presence of the surface energy barrier prevents changes in the number of vortices.

Minimum Energy Paths
Vortex nucleation in a magnetically superheated superconductor is an example of metastable state decay. Various approaches to the study of metastable decay have been proposed [41; 42], and they identify the transition rate as τ −1 = A exp − ∆E kBT , where A is a prefactor. ∆E is the energy barrier that has to be crossed for the transition to occur. In a system whose phase space is multidimensional, the barrier corresponds to a saddle point of the energy landscape, while in field theory, the infinite-dimensional saddle point is called sphaleron. Since the energy barrier appears as the argument of the exponential, the first step in the study of rare events is the sphaleron identification. A way to achieve this is the computation of the minimum energy path (MEP) between two states. The MEP is a path in the configuration space such that it crosses the minimum in the cotangent space of the path point by point. Considering a system described by a potential f (q) where q is the configuration of the system, the mathematical definition can be written like where the transition coordinate s has been introduced. The function f ⊥ is the potential function restricted to the cotangent space of the trajectory in point q(s).

Simplified string method
The most common family of methods to compute minimum energy paths are the so-called chain-of-states methods. These methods are based on the constrained optimization of a collection of systems {q n }, called frames, subjected to the same potential V (q) but a different fictitious force. In this way, it is possible to capture aspects of dynamics of a transition without involving real-time dynamics [33; 43]. These methods can be summarized as variants of the widely known Nudged Elastic Band method [44; 45]. However, nudged elastic band techniques require an additional undefined constant (string elasticity) for introducing fictitious forces. Besides, the chain can miss the saddle point and trying to prevent this requires additional considerations [46].
The approach we use is the simplified string method [31], generalized to a gauge field theory, in which we evolve the entire path q(s) in the configuration space towards the minimal energy one. The variable s continuously parametrizes the path which connects the initial state q(0) to the final state q(1). It is crucial to notice that, when talking of evolution, we refer to a pseudodynamics of the system and not real-time dynamics. To evolve the initial path to the minimal energy one, we apply a gradient descent algorithm according to: Here τ is the pseudo-time describing the evolution of the string towards the MEP. The quantity f [a, ψ] is the free energy defined in Eq. (5), q(s) = (a(r, s), ψ(r, s)) and ∇f = δf δa , δf δψ * . The term λ(q, s) ∂q/∂s has the function of maintaining a certain parametrization of the path q(s). Since ∂q/∂s is uniquely directed along the tangent direction, λ does not affect the trajectory of the string.

Algorithm
In the numerical implementation the path is discretized in a collection {q n } of N frames such that q n ≡ q n N −1 . At each iteration, the algorithm performs a minimization procedure and a reparametrization. In the minimization phase, a certain number of gradient descend steps are performed. In the reparametrization phase, the displaced frames are used to compute a new path by an interpolation method. Usually, the cubic spline interpolation is employed as it guarantees a smooth curve. The interpolation function is then used to generate a new collection of frames that are chosen to be equidistant with respect to the selected metric. There is a certain freedom in the choice of the metric, which we analyze shortly. To describe more in detail how the string method works, let us define q i,k n as the frame number n at iteration i, while k is the index of minimization steps. The simplified string method iteration works in the following way:

Minimization phase
Perform M minimization steps using a gradient descent method: where d i,k n is the descent direction and α i,k n is the step length. At the end of the minimization step the displaced framesq i n = q i,M n have been computed.

Reparametrization step
Let · be the norm induced by a chosen metric, the displacements between the frames are defined as ∆s n = q i n −q i n−1 . The parametrization is then computed as In this way the curve is parametrized as q n ≡ q(s n ) with s ∈ [0, 1]. Using the input set I = {(s n , q n )} the interpolating functionq i (s) is estimated. Finally, a new set of equidistant frames is generated by

Metric and Gauge Invariance
The application of the string method to sphalerons identification in gauge theories is associated with several technical challenges, mainly related to metric definition for a Hilbert space and gauge invariance preservation. For this reason, we have developed a variant that we called gauged string method, which we used to compute physically meaningful minimum energy paths while avoiding the use of global gauge fixing. If we consider a generic path in the phase space q(s), the system may evolve in s changing the gauge of the solution without affecting the observable quantities of the system. In other words, the string gets twisted in the gauge degree of freedom crossing zones of the configuration space where there is no movement in the physical fields.
This can be better understood thinking of the role of the metric, necessary to calculate the distances ∆s n = q i n −q i n−1 needed for the reparametrization step. In a straightforward extension of the prescription of the simplified string method to a Hilbert space, the metric for the Ginzburg-Landau free energy is the distance defined as: ∆s gd n = a n − a n−1 2 + ψ n − ψ n−1 where · is the standard L 2 norm. The problem arising with this choice is the use of the gauge dependent fields a and ψ. This in a way means that two configurations of the same physical system expressed in two different gauges are considered as two distinct systems. Hence, the string can move in the gauge degree of freedom. There can be a region of the string where the fields a and ψ are changing, but there is no variation in the observable quantities, such as b and j. The approach we use is to define a new gauge-invariant metric. With respect to the gauge transformation in Eq. (6) a possible choice, based on Eq.(12), reads: ∆s ph n = |ψ| n − |ψ| n−1 2 + 2 κ ∇θ n − a n − 2 κ ∇θ n−1 − a n−1 Eq. (13) entails some limitations. In fact, let us suppose to be interested in studying samples with particular geometry. Inside holes, or other domain exclusion, the order parameter ψ is not defined. Hence, Eq. (13) becomes ∆s ph n = a n − a n−1 , which is clearly not gauge invariant. Therefore, we can modify the definition of the metric in Eq. (13), to obtain a gauge invariant quantity which is related to the order parameter ψ and simultaneously defined both in the domain and geometrical exclusions. By using the observables b and j we define: In the development of our method, we compared the outcomes of choosing different metrics. The sphaleron does not depend on it, as in Fig. 6 displays. In fact, the height of the barrier is the same for s gd (continuous line) and s ph (dashed line), with the definition of Eq. (14), and it remains the same also for the metric in Eq. (13). If we focus on the line obtained using the metric s gd , we notice that the first part, highlighted in red, is an artifact that corresponds to pure gauge evolution.  (13) and gauge independent metric s ph defined in Eq. (14). The red highlighted zone at the beginning of the f (s gd ) path collapses to a single point in the f (s ph ). This means that the pseudo-motion in that zone is an artefact due to the gauge symmetry. Indeed, in the red zone the system moves along the nonphysical gauge degree of freedom.

Details of the numerical method developed
The energy landscape of the GL model is a system far more complex than the usual molecular dynamics problems addressed by the string method. The problem for 2D systems involves four different fields, it is strongly non-linear, and has a local gauge symmetry. For this reason, we need to introduce some expedients to apply the string method to this problem. One is the gaugeindependent metric described above, and the others concern the numerical implementation of the algorithm.
In our computation, the string q(s) is discretized as a vector of F frames q n . The number of frames we employ ranges between 80 and 120, depending on the transition we are considering. Each of these frames, q n = (a x , a y , Re{ψ}, Im{ψ}) represents a state of the system and is discretized on a mesh grid with a minimum of 400 × 400 lattice sites. We adopted a second-order finite difference scheme for the discretization of the functional. Then the gradient force is computed by simple analytical derivation of the function obtained by the discretization process.
For the minimization step, we applied the nonlinear conjugate gradient (NLCG) method with the Polak-Ribière-Polyak condition with automatic reset and exact line search. NLCG is necessary as the steepest gradient descent is too slow for this type of problem. Notice that, since the frames are displaced in the reparametrization step, the conjugacy is lost, and the conjugate direction needs to be reset at each iteration. In the final steps of the simulation, conjugacy can be switched off. By doing this, the algorithm applies the steepest gradient as descent direction. For what concerns the reparametrization step, we decided to use a linear interpolation step. Linear interpolation is computationally cheaper than the usual cubic interpolation, used within the string method, and more numerically stable. The only negative side of this choice is that the continuity up to the second derivative of q(s) is lost. However, the loss of analytical smoothness is counterbalanced by the use of an high number of frames, which results in a more dense discretization along s, generating a regular and well-behaved MEP.
The initial guess we used the ansatz for a single vortex ψ(r) = tanh −(r − r 0 ) 2 (x − x 0 + i(y − y 0 ))/|r − r 0 | where r 0 is the center of the vortex and it varies frame by frame. It starts from a point outside the superconductor at s = 0 and moves toward the center of the specimen for s = 1. Concerning the numerical extrema of the Hamiltonian (5), the current density perpendicular to the boundary may slightly deviate from zero due to the limited accuracy of the finite-difference scheme. This discrepancy tends to zero with increasing mesh density. Strictly speaking, for non-extremum cases, i.e., configurations along the MEP except initial, sphaleron (saddle point) and final states, this behavior is not guaranteed. Therefore, we monitored the perpendicular component of the current, and it turned out that this value is negligible. The computer implementation is developed using the General Purpose GPU (GPGPU) paradigm. In particular, the core of the tool is implemented in CUDA C++ language and run on a NVIDIA Quadro P6000 GPU.

Quantitative analysis of the planar edge problem
We simulated a 2D square system with side length L = 24λ. In such a vast domain, the interactions between the entering vortex and the other edges of the sample are negligible, allowing the study of the surface barrier of a flat edge, which is the standard case studied in the previous literature.
The minimum free energy path of the vortex nucleation process, in this case, presents only one sphaleron whose free energy is F sphaleron = max s F (s). Once we identify the sphaleron, we compute the nucleation and escape barriers according to: where Once these two quantities are computed for a wide region of the κ-H space, it is possible to study the stability limits of the Shubnikov and Meissner phases, defining the spontaneous nucleation and escape fields as ∆F e (H e ) = 0 .
Externally applying H n generates a complete barrier suppression, and allows vortex entery from the boundaries, filling the sample until inter-vortex repulsion force prevents further nucleations. When H e is applied to the superconductor, solutions with vortices are not stable anymore, and they spontaneously escape from the bulk. Notice that, in principle, H n and the superheating field H sh are different quantities. The former is computed studying the entry process of the first vortex, while the latter is the mathematical stability limit of the Meissner state. These two quantities are strictly the same only for type-II superconductors.
We have simulated this transition in a wide region of the κ-H region, and show the results in Fig. 8. From the computed ∆F n is possible to extrapolate the nucleation field H n . To do that, we interpolate the data points with the ansatz ∆F (fit) n (H)/F 0 = C 1 exp (−C 2 (H − H n )) where C 1 , C 2 and H n are the coefficients to be estimated for each value of κ. The calculated nucleation field is within a 2% error compared the superheating field computed in [10].
We defined the escape barrier ∆F e = F sphaleron − F Shubnikov as the barrier crossed in the event of the vortex exit from the bulk. The results for the barrier height in the κ − H diagram are shown in Fig. 9.
Studying the escape barrier as a function of the external field H, we can notice as the barrier is exactly zero only for null external fields. This is in agreement with experimental results showing that weak external magnetic fields are sufficient to keep a vortex inside a superconductor. This result fixes the escape field H e = 0 meaning that an applied field in the opposite direction is needed to expel a vortex from a sample.
The free energy F (s) is not the only interesting quantity which can be studied with this method. For example, it is well known that in the nucleation process the quantization of the magnetic flux does not hold near a surface [1], which can be understood via a mapping of the problem to a vortex-(image)antivortex pair. The variation of the total magnetic flux Φ(s) along the Minimum Energy Path can be calculated as displayed in Figure 10. Comparing the curves F (s), N (s) and Φ(s), we notice that  To assess the effects of nonlinearity, we consider the nucleation problem for the planar surface for which the linear model from Bean and Livingston was constructed [2]. The Bean-Livingston model adopts an image-charge method. In this approach the energy of the vortex can be seen as the sum of the interaction of the vortex at position x 0 with the decaying external field and the interaction with an image anti-vortex at position −x 0 as shown in 11.
x 0 FIG. 11. Sketch of the Bean-Livingston approach for the vortex nucleation process. In the London limit the free energy is obtained considering the interaction of the vortex with the combination of the decaying surface current and the imageantivortex. The image-antivortex is placed symmetrically with respect to the surface and its presence is necessary to impose the correct boundary conditions.
We can write the free energy as a function of the position in SI units as The first term is the interaction with the imageantivortex (with a 1 2 factor needed to avoid double counting), the second is the interaction of the vortex with the decaying magnetic field at the surface, the third term is the energy of the magnetic field of the vortex while the last term is the interaction with the external field.
To compare the Bean-Livingston picture with our results, we needed to extract a similar information from the computed minimum energy paths. This can be done by extrapolating, for each point of the minimum energy path, the position x 0 of the zero of the order parameter, which corresponds to the center of the vortex. In this way, we estimated the energy as a function of the position F (x 0 ). Even for this simplest geometry, the linear models is inaccurate when there is a non-negligible interaction between the vortex core and the edge. For low-κ the Bean-Livingston model fails to capture the mechanism of vortex nucleation as the vortex deformation near the boundaries plays a significant role in substantially decreasing the value of the nucleation barrier as can be observed in Fig. 12. Moreover, differently from the Bean-Livingston approach, The full Ginzburg-Landau model is able to describe also the transition states preceding the formation of a distinguishable vortex structure.
Vortex nucleation in a L-shape geometry A case of great interest for superconducting devices is a domain shaped as a L. In this case, as showed in Figure 14 which displays the sphaleron magnetic field configuration for vortex entry, the tip of the angle acts as nucleation center for the entering vortex. The L-shaped geometry is a well-studied benchmark for the numerical finding of critical fields [47]. The problem of sphaleron and energy barrier for this case was not solved. Differently from the flat boundary case, in this situation, the vortex enters from the concave corner. To gain further insight into the dependence of the barrier on the sample geometry, we continuously deform the L-shaped geometry by varying the curvature radius of the concave corner, as shown in Figure 14 b. For instance, curvature with a radius equal to zero corresponds to a sharp π/2 angle as in Figure 14 a, while r 1 eliminates concave corner and thus should yield results similar to the flat-boundary situation. The vortex entry barrier strongly depends on the geometry of the sample for high external magnetic fields Figure 14 c. In general, for H > 0, our numerical study shows that the vortex nucleation barrier increases with the curvature radius. By performing an exponential fit of the values of the nucleation barrier, it is possible Free energy of a vortex as a function of the distance x0 to the edge of the sample compared to the Meissner state free energy FM. The black line is calculated using the Bean-Livingston model, while the red one is obtained by extrapolating the position of the vortex core from the minimum energy path. An offset has been added to have the two curves coinciding for x0 deep in the bulk, since we are comparing a London model result with Ginzburg-Landau results. The value matches for long distances while near the boundary, nonlinear effects significantly diminish the value of the barrier. Moreover, the Bean-Livinston approximation fails to describe the steps of the process preceding the nucleation of a distinguishable vortex. This can be appreciated in the graph as the abrupt change from zero to about 2F0. The parameters of the simulation are κ = 1.3 and H = 0.80Hc. to derive the nucleation field H n as shown in Figure 13. The nucleation fields is significantly lower with respect to the flat edge.  Figure 14. The analytical estimate by [10] is plotted as reference. The gauged string method provide a reliable way to compute nucleation fields in arbitrary geometry. This quantity is accessible with simple experiments.