Metastability and bifurcation in superconducting nanorings

TL;DR

This paper introduces a numerical approach to analyze metastability and bifurcation phenomena in superconducting nanorings, providing insights into critical currents, activation energies, and their dependence on temperature and current.

Contribution

It presents a novel numerical method based on the Usadel equation to study stationary points, metastability, and bifurcation in superconducting nanorings, including a correction to Bardeen's formula.

Findings

Computed critical current as a function of temperature.

Mapped activation energy landscape for nanorings.

Identified discrepancies with experimental data at high currents.

Abstract

We describe an approach, based on direct numerical solution of the Usadel equation, to finding stationary points of the free energy of superconducting nanorings. We consider both uniform (equilibrium) solutions and the critical droplets that mediate activated transitions between them. For the uniform solutions, we compute the critical current as a function of the temperature, thus obtaining a correction factor to Bardeen's 1962 interpolation formula. For the droplets, we present a metastability chart that shows the activation energy as a function of the temperature and current. A comparison of the activation energy for a ring to experimental results for a wire connected to superconducting leads reveals a discrepancy at large currents. We discuss possible reasons for it. We also discuss the nature of the bifurcation point at which the droplet merges with the uniform solution.

Tables

Table 1. Table 1: The coupling strengths and the numerically obtained values of the critical current at T = 0 𝑇 0 T=0 for several values of K D subscript 𝐾 𝐷 K_{D} . The current is in units of 𝕀 0 subscript 𝕀 0 \mathbb{I}_{0} , eq. ( 30 ).

$K_D$	5	10	20	40
$\lambda N\left(0\right)$	0.280	0.234	0.202	0.177
$I_c\left(0\right)/{𝕀}_0$	0.07752	0.07796	0.07819	0.07831

Equations

$$I_c(T)\approx I_c(0)[1-(T/T_c{)}^2{]}^{3/2}$$

$$\Delta F=F(\mathrm{droplet})-F(\mathrm{uniform})$$

$$\mathcal{F}=\frac{|\Delta {|}^2}{\lambda }+2\pi TN(0)\sum _{\omega >0}\left\{-{\Delta }^{\ast }F-\Delta F^{\ast }+\frac{1}{2}D\left[\nabla F^{\ast }\nabla F+(\nabla G{)}^2\right]+2\omega (1-G)\right\}.$$

$$\Delta (x)=2\pi T\lambda N(0)\sum _{\omega >0}F(x,\omega ).$$

$$\frac{1}{\lambda N(0)}=\ln \frac{T}{T_c}+\sum _{k=0}^{K-1}\frac{1}{k+\frac{1}{2}}$$

$$-\Delta -\frac{D}{2}\left({\nabla }^{2}F-\frac{F}{G}{\nabla }^{2}G\right)+\frac{\omega F}{G}=0,$$

$$J=-i\pi TN(0)D\sum _{\omega >0}(F^{\ast }\nabla F-F\nabla F^{\ast }).$$

$$F(x,\omega )=e^{iqx}\tilde{F}(x,\omega )=e^{iqx}[R(x,\omega )+iI(x,\omega )],$$

$$\nabla F^{\ast }\nabla F=(\nabla -iq){\tilde{F}}^{\ast }(\nabla +iq)\tilde{F},$$

$$|\tilde{F}{|}^2+G^2=1.$$

$$P:F(x,\omega )\to F^{\ast }(-x,\omega ).$$

$$F(x,\omega )=-F^{\ast }(-x,\omega ),$$

$$F(x,\omega )=f(\omega )e^{iQx},$$

$$R(0,\omega )=I^{\prime }(0,\omega )=0,$$

$$F(x,\omega )=-F^{\ast }(L-x,\omega ).$$

$$\tilde{F}(x,\omega )={\tilde{F}}^{\ast }(L-x,\omega ).$$

$$R^{\prime }(L/2,\omega )=I(L/2,\omega )=0.$$

$$\tilde{F}(\omega ,x)=f(\omega )e^{i(Q-q)x},$$

$$Q=Q_m=q+\frac{\pi (2m+1)}{L},$$

$$W(\omega )=\frac{1}{2\pi }{\int }_{-L/2}^{L/2}{\partial }_x\arg F(x,\omega )dx,$$

$$Q=q+\frac{\pi }{L}.$$

$$R(x,\omega )$$

$$I(x,\omega )$$

$$2\pi W=qL+2{\int }_0^{L/2}{\partial }_x\arctan \frac{I(x,\omega )}{R(x,\omega )}dx$$

$$=qL-\pi \text{mbox}sgn\frac{I(0,\omega )}{R(0^{+},\omega )}.$$

$$R(x,\omega )>0,x>0.$$

$$\frac{1}{\lambda N(0)}=\sum _{k=0}^{K_D-1}\frac{1}{k+\frac{1}{2}}.$$

$$\ln \frac{T}{T_c}=-\sum _{k=K_D}^{K-1}\frac{1}{k+\frac{1}{2}},$$

$${\bar{\xi }}^2=\frac{\hslash D}{4\pi T_c}.$$

$${\mathbb{F}}_0=(2\pi T_c{)}^{2}N(0)A\bar{\xi },$$

$${\mathbb{I}}_0=2e\times 2\pi T_{c}N(0)AD{\bar{\xi }}^{-1}=\frac{4e{\mathbb{F}}_0}{\hslash },$$

Full text

I Introduction and summary

Destruction of superconductivity in thin wires at high currents is a classic topic. Bardeen’s 1962 review Bardeen summarizes the state of the art at the time. In particular, it presents an interpolation formula

[TABLE]

for the critical (depairing) current as a function of the temperature. An important development subsequent to Bardeen’s article has been Little’s work Little , which emphasized the role of large thermal fluctuations (phase slips) as a cause for transition to the normal state at currents below $I_{c}$. Experimental studies of this switching transition have developed fast in recent years Sahu&al ; Li&al ; Aref&al .

On the theoretical side, a study of thermal phase slips begins with identifying the saddle-point of the free energy (the critical droplet) that determines the activation barrier. A number of models have been used for this purpose. The original computation LA of Langer and Ambegaokar (LA) was in the context of the Ginzburg-Landau (GL) theory. The GL theory is well motivated microscopically but is limited to a vicinity of the critical temperature. One alternative may be to use a discrete model—essentially, a chain of Josephson junctions. Relevance of such models to nanowires has been discussed in Matveev&al ; Goswami&Chakravarty ; Kh:2016 .

Here, we proceed with a continuum description and present results for activation energies obtained by a direct numerical solution of the one-dimensional Usadel equation Usadel . The latter arises as the dirty-limit reduction of the Eilenberger-Larkin-Ovchinnikov theory Eilenberger ; LO . It is well motivated microscopically, applies at any temperature, and contains the GL theory as a limit. For an infinite (very long) wire, the critical droplets have been considered on the basis of the Usadel equation in Semenov&al . Here, we will be interested in solutions for a wire of a finite length.

Usadel’s equation is second-order in spatial derivatives and requires boundary conditions. Let us say a few words about those. Bardeen’s formula (1) pertains to a uniform superconducting state: the current density and the gap are the same everywhere along the wire (which is assumed here to be in the $x$ direction). If the wire is connected by leads into an external circuit, such a uniform state can only be an idealization, the more so the shorter the wire is. The true equilibrium in this case is necessarily $x$-dependent. Such $x$-dependent solutions have been found for very short wires (bridges of length $L\ll\xi$, where $\xi$ is the coherence length), both in the GL theory AL and on the basis of the Usadel equation KO . In principle, the numerical approach we describe here can be used to find such solutions also for longer wires, provided one is willing to do some modeling of physics in the leads.†††In the context of the GL theory, transition to longer wires has been recently considered in Marychev&Vodolazov .

Another aspect of the boundary conditions, which is particularly relevant to a study of phase slips, is whether the leads prevent rapid changes in the boundary values of the phase of the order parameter (as bulk superconductors do) or allow such changes to occur easily (the case, for instance, for normal contacts). In the first case, the critical droplet does not have to have the same value of the current as the equilibrium state from which it originates, while in the second case it typically does.‡‡‡An intuitive picture of how bulk superconducting leads allow for fluctuations of the current in the wire can be obtained by viewing the leads as small impedances connecting the ends of the wire to the ground. An impedance will be of order $Z=({\cal L}/C)^{1/2}$, where ${\cal L}$ is the kinetic inductance (which is small for a large superconductor), and $C$ is the capacitance of (a portion of) the lead relative to the ground. The two impedances shunt the ends of the wire, allowing the current in it to fluctuate. For long wires, the two scenarios are not that different (as shown by McCumber McCumber in the case of the LA saddle point), but for shorter wires there is a genuine difference, reflecting the difference in the experimental setup. We refer to the first scenario as nucleation at a fixed winding number, and to the second as nucleation at a fixed current.

The simplest system in which nucleation of the droplet occurs at a fixed winding number, rather than a fixed current, is a superconducting ring. One facet of that simplicity is that, for a ring of a uniform cross section, there is a uniform equilibrium state regardless of the length. This fortuitously circumvents the problem of finding the $x$-dependent equilibrium characteristic of short wires in the presence of leads. In particular, results for a ring of a finite length can be directly compared to the predictions of Bardeen’s formula (1), and a correction factor to (1) can be obtained. With this in mind, we have chosen the ring geometry for the present study. To treat the winding number as a continuous variable, we consider twisted boundary conditions, such as would arise in the presence of a fractional magnetic flux through the ring.

After some preliminaries in Sec. II, we discuss details of the boundary conditions in Sec. III and then describe the numerical results in Sec. IV. Here, we present in advance the main result: the metastability chart (Fig. 1), which shows level contours of the activation free energy

[TABLE]

for the critical droplet mediating transition from the uniform equilibrium with given values of the temperature and current to the state with one fewer unit of winding. When $\Delta F$ is expressed as a multiple of the free energy unit (29), the metastability chart depends on only two parameters: the length $L$ of the ring and the effective electron-electron coupling $\lambda N(0)$. For the typical weak-coupling case $\lambda N(0)<0.3$, the dependence of $\Delta F$ on the coupling is rather weak.

For a uniform ring, the properties of the equilibrium state do not depend on $L$, but those of the critical droplet do. One might expect that, as $L$ becomes large compared to a suitably defined coherence length $\xi(T)$, $\Delta F$ approaches the constant value corresponding to an infinite wire, with finite-size corrections of order $\xi(T)/L$. That is true over most of the chart but not in a narrow band of currents near the critical. The reason is that the size of the droplet grows for currents near $I_{c}(T)$, so even for a small $\xi(T)/L$ there is a range of currents for which finite-size effects are important.§§§For the LA droplet in the GL theory, this effect can be deduced already from the expressions presented in LA . In what follows, we will often quote the length in units of the diffusion length $\bar{\xi}$, defined by eq. (28) below. To estimate the significance of finite-size effects on $\Delta F$, we can compare the results for $L=82$, which is the value used for the chart of Fig. 1, to those for about twice the length, $L=162$. We concentrate on temperatures $0.25<T/T_{c}<0.5$ and currents $I/I_{c}(T)\sim 0.9$, where thermal switching transitions are typically observed. More specifically, we consider the “observability line” $\Delta F=0.03$ (which lies in the middle of the second metastability band in Fig. 1). We have found that, for these ranges of the parameters, increasing the length to $L=162$ leads to a decrease in $\Delta F$ by about 20%.

We should remark that, at a finite $L$, the critical (maximum) current $I_{c}(T)$, represented by the thick upper curve in Fig. 1, does not coincide with the bifurcation current $I_{\rm bif}(T)$, at which the droplet merges with the uniform solution and the activation barrier disappears. The two would be strictly the same only in the $L\to\infty$ limit. The separation between $I_{\rm bif}(T)$ and $I_{c}(T)$ is possible because the value of the current does not define the uniform solution uniquely: the winding number does. As we increase the winding from zero, the current first increases until it reaches the maximum, $I_{c}(T)$, where $\Delta F$ is still nonzero, and then decreases down to $I_{\rm bif}(T)$, where $\Delta F$ finally vanishes.¶¶¶Very recently, this effect has been observed experimentally Petkovic&al . The value of $L$ chosen for Fig. 1, however, is already large enough for the difference between $I_{c}$ and $I_{\rm bif}$ to be unnoticeable on the scale of the plot. We return to discussion of this point in Sec. IV, where results for a smaller length are presented.

In Sec. IV.4, we carry out a comparison of results for $\Delta F$ obtained on the basis of our numerical solutions to the experimental results for one of the samples of Ref. Aref&al . The comparison uses only experimentally measured quantities and has no free parameters. We find a very good agreement at low currents, but a significant discrepancy at large ones, where switching transitions are observed. The experimentally determined $\Delta F$ scales as $(1-I/I_{c})^{3/2}$ near the critical current, while the one obtained numerically is almost linear. As a consequence, the numerical result there is significantly larger than $\Delta F$ deduced from experiment. We mention a couple of possible reasons for this discrepancy in Sec. IV.4. Here, we remark only that the resolution of the problem would be much aided by experiments on rings thin enough for the “premature” switching ($I<I_{c}$, $\Delta F\neq 0$) to have a chance to be observed. The experiments Petkovic&al on thicker rings have already accessed the “deterministic” switching regime $\Delta F=0$.

For a uniform ring, the equilibrium solution continues to exist even after the winding number $W$ is increased beyond the bifurcation point $W_{\rm bif}$, although it becomes absolutely unstable. One may ask if this property is generic and will hold also for a wire connected to leads. In general, we consider that unlikely: as we discuss in Sec. V, the persistence of the equilibrium solution for a ring can be seen as a consequence of the translational invariance. (More generally, is can be seen as a consequence of there being different “versions” of the critical droplet that are degenerate in free energy.) It would nevertheless be interesting to see if there are any distinct physical consequences of such absolutely unstable solutions: on the one hand, fabricating a suitable ring sample may not be out of question; on the other, for a wire connected to leads, there is an approximate translation symmetry in the middle, so the same consequences may show up in that case as well.

II Preliminaries

Usadel’s equation Usadel can be obtained by variation of a certain free energy functional, which is essentially that of an $O(3)$ nonlinear sigma model. The corresponding free energy density is

[TABLE]

The sigma-model variables are the complex $F(x,\omega)\equiv F$ and the real $G(x,\omega)\equiv G$, related by the nonlinear constraint $|F|^{2}+G^{2}=1$. Variation with respect to $\Delta$ produces the self-consistency condition

[TABLE]

In a thin superconductor of uniform cross section, $F$ and $G$ depend on only one spatial coordinate, $x$ (so in particular $\nabla=\partial_{x}$). In addition, they depend on the Matsubara frequency $\omega$ that runs over the values $\omega_{k}=2\pi T(k+\frac{1}{2})$, where $k$ is an integer. Note that the summations over $\omega$ in (3) and (4) are restricted to $k\geq 0$. If one wishes to pursue the sigma-model context, one can think of $F$ and $G$ as the in-plane and out-of-plane components of a (fictitious) magnetization, and of the superconducting gap $\Delta$ as the in-plane component of the self-consistent fictitious magnetic field.

The parameters appearing in (3) are as follows: $\lambda$ is the electron-electron coupling, $T$ is the temperature, $N(0)$ is the density of states at the Fermi level (for a single spin projection), and $D$ is the electron diffusion constant. Following Eilenberger , we will use the substitution

[TABLE]

to relate $\lambda N(0)$ to the critical temperature. The sum in (5) is cut off just below some integer $K$, which corresponds to the maximum (Debye) frequency in units of $2\pi T$.

Variation of the free energy with respect to $F$ (with $G$ expressed as $G=(1-|F|^{2})^{1/2}$) reproduces the Usadel equation

[TABLE]

which should be solved together with (4).

The current density, in units of twice the electron charge, is Usadel

[TABLE]

Although this nominally depends on $x$, for static solutions (the only case considered here) in samples of uniform cross section, charge conservation ensures that $J$ is $x$-independent.

It is convenient to define a new field, $\tilde{F}(x,\omega)$, by removing “most of the winding” from $F$ (we will say precisely how much below), as follows:

[TABLE]

where, for now, $q$ remains unspecified. For future use, we have separated the real and imaginary parts of $\tilde{F}$. In terms of $\tilde{F}$, the gradient term in (3) becomes

[TABLE]

and the nonlinear constraint becomes

[TABLE]

In what follows, we consider $\tilde{F}$, rather than $F$, as the main dynamical variable and $q$ as a parameter, on which the free energy depends explicitly via (9).

III Boundary conditions

We begin by considering the theory on a ring of length $L$, with periodic conditions on $F(x,\omega)$ for all $\omega$. Later, we will allow for twisted boundary conditions, such as those associated with a magnetic flux through the ring.

For our purposes, it is not sufficient to simply state that the solution is periodic. For one thing, because of the perfect translational invariance, any $x$-dependent solution will have a translational zero mode (a zero eigenvalue of the Hessian matrix), which will render the numerical method we use inapplicable. The same applies to the zero mode generated by the symmetry with respect to global phase rotations. We wish to eliminate these zero modes by placing the “core” of the solution at a particular point (say, $x=0$) and also fixing its overall phase.

We first note that the free energy is invariant under the discrete transformation

[TABLE]

To eliminate the zero modes, we concentrate on solutions that have definite parity under $P$. It is then sufficient to consider only $P$-odd solutions, namely those that satisfy

[TABLE]

because any $P$-even solution can be turned into a $P$-odd one by multiplication by $i$. The condition (12) breaks the translational symmetry to translation by $L/2$, and the symmetry with respect to global phase rotations to multiplication by $-1$. These are discrete symmetries, which do not produce zero modes.

As an example of how (12) fixes the overall phase, consider one of the uniform solutions, which are of the form

[TABLE]

where $Q$ is a parameter. For this to satisfy (12), $f(\omega)$ must be purely imaginary.

In Sec. II, we have defined a new field $\tilde{F}$, obtained from $F$ by taking out some of the winding, cf. eq. (8). The amount of winding taken out is represented by $e^{iqx}$, where $q$ has so far been left undetermined. Regardless of the value of $q$, the condition (12) translates into the identical condition for $\tilde{F}$ or, in terms of the real and imaginary parts of $\tilde{F}$, into the condition that $R$ is odd, and $I$ is even about $x=0$. In other words,

[TABLE]

where the prime denotes derivative with respect to $x$.

As long as $F(x,\omega)$ is periodic, the condition (12) implies that it is odd also under $P$ combined with a translation by the total length $L$, that is

[TABLE]

We now choose $q$ so that $\tilde{F}$, on the contrary, is even under this combination:

[TABLE]

In terms of the real and imaginary parts of $\tilde{F}$, this becomes

[TABLE]

For a periodic $F$, the restriction this imposes on $q$ is $e^{iqL}=-1$. We observe, however, that we now have enough boundary conditions to formulate a boundary problem on $[0,L/2]$ directly for $\tilde{F}$ (or, equivalently, for its real and imaginary parts), so we can switch to more general twisted $F$ simply by abandoning the restriction, i.e., by allowing for arbitrary values of $q$. In this way, $q$ becomes a continuous parameter, which enters the calculation only through the expression (9) for the gradient term.

Next, we note that, for a given $q$, the boundary conditions (14) and (17) still do not determine the solution uniquely. The reason is that the configuration space is comprised of sectors corresponding to different winding numbers, and one can in principle look for a solution to the boundary problem in each of these sectors individually. For example, for the uniform solutions (13),

[TABLE]

and the boundary conditions require only that $Q$ belongs to the discrete set

[TABLE]

where $m$ is any integer. The uniform states corresponding to different $Q_{m}$ in general have different energies, and their activated decays are described by different critical droplets. We want to make sure that we are solving for the right droplet, the one that describes the decay of a particular uniform state. In other words, we want to be able to restrict to one particular sector out of many.

Loosely speaking, the different sectors can be viewed as neighborhoods of different uniform solutions, those with different $Q$ from the sequence (19). More precisely, separation into sectors can be done with the help of the winding number

[TABLE]

which is well defined for any configuration for which the absolute value $|F(x,\omega)|$ does not vanish anywhere. As indicated in (20), a priori one can imagine the winding number to be a function of $\omega$, but here we restrict attention to solutions for which it is the same for all $\omega$, $W(\omega)=W$.

For a uniform solution, $W=QL/2\pi$. We can scan over all uniform solutions by fixing $m$ in (19) and varying $q$, and we now wish to concentrate on the sector with $m=0$, for which

[TABLE]

Recalling that by virtue of (14) $f$ is purely imaginary, we find the real and imaginary parts of (18) to be

[TABLE]

These expressions reflect the characteristic property of the $m=0$ sector: the real part of $\tilde{F}$ has no zeroes except the one at the origin, prescribed by (14), and in addition, at $x$ just above zero, the signs of the real and imaginary parts are opposite. Indeed, for a given $q$, any field that has these properties and satisfies the boundary conditions (14) and (17) will have the same winding number as the uniform solution (22)–(23). This can be seen from the following relation, which applies whenever $R$ has no zeroes at $x>0$:

[TABLE]

[TABLE]

Unlike the relative sign of $R$ and $I$, their overall sign is unimportant: there is a companion set of solutions, with the same energy and current, that differ from the ones we find by an overall sign. As already noted after (12), this is a remnant of the original invariance with respect to global phase rotations. So, in what follows we identify the $m=0$ sector by the condition

[TABLE]

Eq. (24) then shows that a solution (uniform or not) will be in the $m=0$ sector if $I(0,\omega)<0$ and in the $m=-1$ sector if $I(0,\omega)>0$.

Let us note that the critical droplet responsible for fluctuations out of the $m=0$ sector can itself be in either $m=0$ or $m=-1$ sector, depending on the sign of $q$. To see that, let us first consider the special case $q=0$. For the uniform solution, eq. (21) now gives $Q=\pi/L$, so the winding number is $W=1/2$. The uniform solution with one fewer unit of winding (corresponding to $m=-1$) has $Q=-\pi/L$ and $W=-1/2$. These solutions have exactly opposite currents and equal energies. If we think of the twisted boundary conditions as being a result of a magnetic flux through the ring, this special situation corresponds to the ring biased by half the flux quantum. The critical droplet in this case has $I(x,\omega)=0$ identically, and the current on it is zero; because of the zero of $F$ at $x=0$, the winding number is not well defined. One can say that for $q=0$ the droplet lies at the boundary between the $m=0$ and $m=-1$ sectors.

Recall that, in our computation, $q$ is a parameter that can be chosen arbitrarily. As we change it from positive to negative, $I(0,\omega)$ changes from negative to positive, and the droplet moves from the $m=0$ sector to $m=-1$. The underlying physics is that, at $q>0$, the $m=0$ uniform solution is metastable, and the droplet describes its activated decay to the $m=-1$ solution, which has a lower free energy. At $q=0$, the $m=0$ and $m=-1$ uniform solutions become degenerate in energy, and at $q<0$ the situation is reversed: for sufficiently small $|q|$, the $m=0$ solution is absolutely stable, and the droplet describes the decay of the $m=-1$ solution to it. Alternatively, even for $q<0$, the droplet can be interpreted as a fluctuation out of the $m=0$ sector. In that interpretation, its activation energy $\Delta F$ will always be counted from the free energy of the $m=0$ uniform solution. This is how it has been presented on the chart of Fig. 1, where the states corresponding to $q=0$ are shown by the thick dashed line.

The sign condition (25), together with the boundary conditions (14) and (17), completes our formulation of the boundary problem for $\tilde{F}$. Note that this boundary problem produces only half of the configuration, namely, the part corresponding to $0\leq x\leq L/2$. The other half can be obtained by reflecting $R(x,\omega)$ as odd and $I(x,\omega)$ as even about $x=0$.

IV Numerical results

Before we proceed to the results, let us discuss some aspects of the substitution formula (5) that are important for numerical work. Because $\lambda N(0)$ is constant, the cutoff integer $K$ in (5) must increase with decreasing temperature. On can express this by saying that, as $T$ is lowered, the system builds up a “synthetic” frequency dimension, represented by a chain of values of $k$. Physically, this dimension is associated with motion of electrons within a pair.

Suppose we take the smallest value of $K$, which we will call $K_{D}$, to correspond to $T=T_{c}$. Then, from (5),

[TABLE]

Thus, we can use $K_{D}$ instead of $\lambda N(0)$ as a parameter determining the strength of the interaction. A typical value is $K_{D}=10$, which corresponds to $\lambda N(0)=0.234$. As long as $K_{D}$ is not much smaller than this, the results depend comparatively weakly on it. As a consequence of (5), once $K_{D}$ is specified, the size $K$ of the frequency dimension is related to $T/T_{c}$ as follows:

[TABLE]

One limitation of this procedure is that, since $K$ can only increase in discrete (integer) steps, $T/T_{c}$ is likewise limited to a discrete set of values. The effect is significant mostly for $T$ near $T_{c}$. Indeed, the smallest $K>K_{D}$ one can possibly choose is $K=K_{D}+1$. Then, from (27), the largest $T/T_{c}$ not equal to unity is $T/T_{c}=\exp[-1/(K_{D}+\frac{1}{2})]$. For $K_{D}=10$, this translates into the cutoff at $T/T_{c}=0.91$, seen in the chart of Fig. 1.

Of course, one can choose $K=K_{D}$ to correspond not to $T=T_{c}$ but to a slightly lower temperature. Since (26) will no longer apply, this will result in a slightly different coupling, but will allow one to obtain a data point in the missing range of temperatures near $T_{c}$. Given that the region near $T_{c}$ is not our main interest here, we will not pursue that.

IV.1 Units of length, energy, and current

A natural unit of length in the present problem is the characteristic diffusion length $\bar{\xi}$ defined by

[TABLE]

In this subsection only, we have restored $\hbar$. For orientation, for $T_{c}=k_{B}\times(5~{}\mbox{K})$ and $D=1.2\times 10^{-4}$ m${}^{2}/$s (values appropriate for amorphous MoGe wires Bezryadin:book ), $\bar{\xi}=3.8$ nm. In what follows, we will often quote the length $L$ in units of $\bar{\xi}$, and the winding number parameter $q$ in units of $\bar{\xi}^{-1}$.

A natural unit of the free energy is

[TABLE]

where $A$ is the cross-sectional area of the sample. This amount represents a condensation energy of order $2\pi T_{c}$ for each electron in a length $\bar{\xi}$ of the wire in an energy layer about $2\pi T_{c}$ thick near the Fermi surface. Similarly, a natural unit of the electric current is

[TABLE]

where $e$ is the electron charge. In what follows, we present results for the free energy and the current in units of $\mathbb{F}_{0}$ and $\mathbb{I}_{0}$, respectively.

IV.2 Uniform solutions

Uniform solutions are those of the form (13). For these, our main interest is in quantifying the accuracy of Bardeen’s formula (1). When we do not wish to refer to a specific length of the wire, $L$, we will use $Q$ rather than $q$ as a parameter. If a value of $L$ is available, the two can be related by (21). As per discussion at the end of Sec. III, either is now considered a continuous parameter.

For a uniform solution, we do not need to solve the full boundary problem: once $Q$ is chosen, the $x$-dependence of the solution is known, and the Usadel equation becomes an equation for the amplitude $f(\omega)$. It reads

[TABLE]

where $\omega$ takes the values

[TABLE]

Dividing (31) by $2\pi T_{c}$, we see that the full set of parameters on which the solution depends can be chosen as follows: $Q\bar{\xi}$, where $\bar{\xi}$ is the diffusion length (28); $K_{D}$, which determines $\lambda N(0)$; and $K$, the upper cutoff in the sum. From the latter two, the ratio $T/T_{c}$ can be calculated via (27). In what follows, we present the results as functions of that ratio, rather than $K$ itself.

We solve (31) numerically by the multidimensional Newton-Raphson method. For given $K_{D}$ and $T/T_{c}$, we compute the current as a function of the one remaining parameter, $Q\bar{\xi}$. The current reaches a maximum at some ($T$-dependent) critical value $Q_{c}(T)$; the value at the maximum is designated as the critical current, $I_{c}(T)$. The zero-temperature limit of $I_{c}(T)$ may be of independent interest, and we present results for it, for a few values of $K_{D}$ in Table 1. The table also lists the values of the coupling strength $\lambda N(0)$, computed from $K_{D}$ via (26). One may note the weakness of the dependence of $I_{c}$ on the coupling, especially for smaller couplings.

The ratio of the numerically computed $I_{c}(T)/I_{c}(0)$ to the interpolating function $[1-(T/T_{c})^{2}]^{3/2}$ proposed by Bardeen Bardeen is shown in Fig. 2. This ratio constitutes a correction factor to Bardeen’s formula (1). The weak dependence of the result on $K_{D}$ means that, for weak coupling, the correction factor is essentially universal.

One may observe that the correction factor retains a significant dependence on $T$ even at the lowest temperatures: it grows by almost 2% between $T=0$ and $T/T_{c}=0.1$. This dependence, however, is almost entirely due to that in Bardeen’s expression $[1-(T/T_{c})^{2}]^{3/2}$. The numerically obtained $I_{c}(T)$ is essentially flat for $T/T_{c}<0.1$ and can be well approximated there by the corresponding value of $I_{c}(0)$ (as found in Table 1). This saturation of $I_{c}(T)$ at $T/T_{c}\approx 0.1$ is visible already on the chart of Fig. 1 and is probably the main qualitative feature distinguishing the numerical result from Bardeen’s formula.

Next, we consider what happens to the uniform solution when $q$ is increased past $q_{c}$, i.e., past the value at which the current reaches the maximum. As a point of comparison, we recall results for the theory of a single degree of freedom (DOF), that is, a “particle” with a potential in the form of a cubic parabola, as used for instance in studies of Josephson junctions. The minimum of the potential corresponds to the ground state, and the maximum to the critical droplet. In that case, the maximum current is, at the same time, a bifurcation point at which the droplet merges with the ground state.

Turning to the present case, note that, for a uniform solution, the current is proportional to the derivative of free energy density with respect to $q$, $d{\cal F}/dq$. Thus, $q=q_{c}$, at which the current is maximal, is an inflection point of ${\cal F}$. For an infinite wire, the standard convexity argument then guarantees that at $q=q_{c}$ the uniform solution becomes unstable to spinodal decomposition. Finite-size effects can delay onset of the instability, so for a wire of a finite length we expect that it occurs, if at all, at some $q=q_{\rm bif}>q_{c}$. (Very recently, this effect has been observed experimentally Petkovic&al .)

As for a single DOF, the change in the stability properties (now at $q=q_{\rm bif}$) is accompanied by a bifurcation—a merging or splitting up of two or more solutions. In our case, the solution merging with the uniform state is the critical droplet; the process is illustrated in Fig. 3. Note that, unlike for a singe DOF, the uniform solution continues to exist at $q>q_{\rm bif}$, even though it becomes absolutely unstable. As we discuss in Sec. V, this can be attributed to the fact that there are two “variants” of the critical droplet, both merging with the uniform solution at $q=q_{\rm bif}$. They are related by a discrete transformation which is a symmetry of the theory and so have the same free energy and current.

IV.3 Critical droplet

The critical droplet is the saddle point of the free energy that sits at the top of the potential barrier separating the uniform solution (13) with $Q=q+\pi/L$, as given by (21), from the uniform solution with one fewer unit of winding, i.e., $Q=q-\pi/L$.

The droplet is $x$-dependent, so finding it requires solving the full eq. (6). Numerically, we proceed by discretizing (6) on a one-dimensional spatial grid of $N$ equally spaced sites (typically, $N=41$ or $N=81$) and applying the multidimensional Newton-Raphson (NR) method to the resulting nonlinear equations for $N\times K$ complex unknowns $F(x_{j},\omega_{k})$. As usual, the success of the NR method depends on having a good initial approximation. We adopt the following protocol. We start at a value of $K$ close to $K_{D}$, which corresponds to $T$ close to $T_{c}$. There, we use the LA solution LA of the GL theory as the initial approximation for the gap $\Delta(x)$, from which we reconstruct an approximation for $F(x,\omega)$ via

[TABLE]

(the usual approximation for transitioning from the Usadel equation to the GL theory). After the NR method finds the exact solution at this value of $K$, we compute the exact $\Delta(x)$ and increase $K$ to a larger value, using (33) to populate the missing modes of $F$ (those with $k$ between the old value of $K$ and the new one). That forms the initial approximation at the new $K$. In this way, increasing $K$ in steps, we move to progressively lower values of the temperature.

Results of the computations have already been presented in Figs. 1 and 3.

IV.4 Comparison to experiment

Comparing our results to experiment requires, as the input data, the $T=0$ value of the critical current $I_{c}(0)$, the critical temperature $T_{c}$, and the length $L$ of the sample. The weak dependence of $I_{c}(0)$ on the coupling strength (see Table 1) allows us to determine the unit $\mathbb{I}_{0}$ of current by using the approximate relation $I_{c}(0)=0.078\mathbb{I}_{0}$. The unit $\mathbb{F}_{0}$ of free energy is then computed from (30). Note that this does not require separate knowledge of the density of states and the cross-sectional area. With the value of $\mathbb{F}_{0}$ in hand, we can refer to a chart such as that of Fig. 1 to find the activation barrier in physical units.

As an illustration, let us carry out this program for sample B of ref. Aref&al (the longest wire described there). The parameters as determined in Aref&al are $I_{c}(0)=12.11$ $\mu$A, $T_{c}=5.48$ K, and $L=221$ nm (for $T_{c}$, we use the value referred to in Aref&al as $T_{c}^{\prime}$). Following the steps outlined above, we obtain $\mathbb{F}_{0}=159$ meV. Next, we use numerical results for $L=62$ (in units of $\bar{\xi}$), a somewhat smaller value than that used for Fig. 1 but matching more closely the physical length. For illustration, we take $T/T_{c}=0.2$. The resulting activation barrier is plotted as a function of the equilibrium current in Fig. 4. The dashed line is the activation barrier obtained in Aref&al by fitting the experimental data. At small currents, the agreement is very good. Of main interest, however, is the region of large currents (say, those within 10% of the critical), where the barrier becomes small enough for the switching transition to be observable. There, the numerical result is almost a straight line and significantly overestimates the observed value. In fact, for $\Delta F=6$ meV (the value obtained numerically for $I/I_{c}=0.9$), the Boltzmann exponent at the temperature in question is over 60; in all likelihood, such an exponent would render the transition unobservable.

We can offer two possible explanations for the discrepancy at large currents. One is that, the reasoning given in the Introduction notwithstanding, nucleation of a critical droplet in a wire connected to superconducting leads does not occur, even for large leads, exactly at fixed winding, so the results obtained for a ring are not immediately applicable. The second possibility is that, at large currents, in addition to the uniform equilibrium states considered here, there are other, possibly nonuniform, states available, activation from which proceeds more easily. For instance, for a periodic structure equivalent to a chain of Josephson junctions, the $b=3/2$ scaling of $\Delta F$ at $I\to I_{c}$ applies even for rings of relatively short lengths Kh:2016 . As the example of an array of phase-slip centers Skocpol&al indicates, such a nonuniform state may exist even in a morphologically uniform wire.

V Consequences of the translational invariance

In addition to the symmetry under the discrete transformation (11), the free energy is invariant, at an arbitrary fixed $q$, under the transformation

[TABLE]

This can be seen as a composition of (11) with a translation by $L/2$ and a global phase rotation. In terms of the real and imaginary parts of $\tilde{F}$, it corresponds to

[TABLE]

The boundary conditions (14) and (17) are also invariant under $P^{\prime}$. Indeed, this symmetry can be viewed as a discrete remnant of the full transitional invariance that the theory had before we imposed (14) and (17). As such, it is characteristic of the rather special situation presented by a uniform ring and is not available generically.

The uniform solution (22)–(23) is invariant under $P^{\prime}$, but the critical droplet is not. Thus, there must be a twin droplet, with the exact same values of the free energy and current as the original, obtainable from it by an application of $P^{\prime}$. The only physical difference between the two is that the “core” of the original droplet is at $x=0$, while the core of the twin has been shifted to $x=L/2$.

One consequence of the existence of the twin is that, at the bifurcation point, there are in fact not two but three solutions merging. Thus, for instance, in Fig. 3, the dashed line should now be taken to mean two solutions: the original droplet and the twin. An index theorem can then be used to relate the properties of the solutions before and after the bifurcation. The argument is parallel to that used in Kuznetsov&Tinyakov to study bifurcation of a different type of solution, the periodic instanton periodic .

The requisite index theorem is the Morse equality Milnor , which states that, for a smooth function $f$ with non-degenerate critical points on a compact smooth manifold $M$,

[TABLE]

where the sum is over all the values of the index (the number of negative modes of $f$ at a critical point), $N_{s}$ is the number of critical points of index $s$, and $\chi(M)$ is the Euler characteristic of $M$. In our case, $M$ is the product of $N\times K$ spheres of the form (10), and $f$ is the discretized version of the free energy with the self-consistency condition (4) and the expression (9) substituted in. As a result, $f$ depends on $q$ as a parameter.

The critical points are, in our case, the various static solutions described earlier. Assuming that no bifurcations except the one in question occur at $q=q_{\rm bif}$, we can compute the change in the right-hand side of (37) between $q<q_{\rm bif}$ and $q>q_{\rm bif}$. That change must be zero, as $\chi(M)$ is a topological invariant. At $q<q_{\rm bif}$, there are three solutions: the uniform solution of index 0, and the droplet with its twin, each of index 1. Their total contribution to (37) is $-1$. This is different from the case of a single DOF (the cubic parabola), where only two solutions are merging, and the total index is zero. In that case, the two solutions can simply disappear at $q=q_{\rm bif}$. In our case, the minimal structure needed at $q>q_{\rm bif}$ to preserve the total index is a single solution with one negative mode. That is indeed what we have seen numerically: the uniform solution remains but becomes absolutely unstable.

It is natural to ask where the instability that the uniform solution acquires at $q>q_{\rm bif}$ leads. Numerically, this can be answered by displacing a little along the negative mode and following a relaxation algorithm. To better describe the results, let us first recall that, in addition to the now unstable uniform solution, which we will call $F_{0}$, there are, for the same $q$, uniform solutions with other winding numbers. For the solution with one unit of winding fewer than $F_{0}$, the real and imaginary parts of $\tilde{F}$ are

[TABLE]

The difference with (22)–(23) is that $R$ and $I$ now have the same sign for $x>0$. There are in fact two solutions of this form: one with $g>0$, and the other with $g<0$; let us call them $F_{+}$ and $F_{-}$. The transformation (35)–(36) maps one into the other. Numerically, we have found that the instability of $F_{0}$ at $q>q_{\rm bif}$ develops into $F_{+}$ or $F_{-}$, depending on the direction of the initial displacement.∥∥∥ It makes sense to consider two fields that differ only by an overall sign as physically equivalent. Then, we would be talking about a loop, which starts at $F_{+}\sim F_{-}$, goes up to $F_{0}$, and then down back to $F_{-}\sim F_{+}$.

Finally, we remark that, as $q$ is increased further, beyond $q_{\rm bif}$, the uniform solution acquires additional negative modes and so is expected to go through additional bifurcation points. We have not studied those in any detail.