Instabilities of the normal state in current-biased narrow superconducting strips

TL;DR

This paper investigates the instabilities in the normal state of narrow superconducting strips near the critical temperature, focusing on high current biases and the resulting dynamical phase-slip states using the Ginzburg-Landau model.

Contribution

It provides a detailed analysis of the instability mechanisms and derives the reentrance current at the normal state transition in superconducting strips.

Findings

Identification of dynamical phase-slip states at high currents

Derivation of the reentrance current value at instability onset

Numerical calculation of the reentrance current based on the Ginzburg-Landau model

Abstract

We study the current-voltage characteristic of narrow superconducting strips in the gapless regime near the critical temperature in the framework of the Ginzburg-Landau model. Our focus is on its instabilities occurring at high current biases. The latter are consequences of dynamical states with periodic phase-slip events in space and time. We analyze their structure and derive the value of the reentrance current at the onset of the instability of the normal state. It is expressed in terms of the kinetic coefficient of the time-dependent Ginzburg-Landau equation and calculated numerically.

Figures

Tables

Table 1. Table 1: The 8 largest simultaneous eigenvalue pairs { ( 𝔲 ​ E ) c , Q c } subscript 𝔲 𝐸 𝑐 subscript 𝑄 𝑐 \{(\mathfrak{u}E)_{c},Q_{c}\} of L ^ ^ 𝐿 \hat{L} and L ^ † superscript ^ 𝐿 † \hat{L}^{\dagger} for N k = 24 subscript 𝑁 𝑘 24 N_{k}=24 Fourier components and corresponding solvability coefficients { λ 1 , λ 2 , λ 3 } subscript 𝜆 1 subscript 𝜆 2 subscript 𝜆 3 \{\lambda_{1},\lambda_{2},\lambda_{3}\} . Labels ν 𝜈 \nu correspond to the intersection points shown in Fig.  4 .

$\nu$	${\left(𝔲E\right)}_c$	$Q_c$	${\lambda }_1$	${\lambda }_2$	${\lambda }_3$
1	0.320	0.486	0.166	-0.503	-0.176
2	0.276	0.537	0.102	-0.230	-0.0787
3	0.21	0.335	-0.0478	0.128	0.030
4	0.189	0.359	0.0265	-0.0671	-0.0169
5	0.156	0.256	0.0133	-0.0345	-0.00570
6	0.144	0.27	0.0074	-0.0183	-0.00258
7	0.124	0.208	0.00348	-0.00922	-0.000833
8	0.116	0.217	0.00199	-0.00528	-0.000436

Table 2. Table 2: Similar to table I in the paper: The 8 largest simultaneous eigenvalue pairs { ( 𝔲 ​ E ) c , Q c } subscript 𝔲 𝐸 𝑐 subscript 𝑄 𝑐 \{(\mathfrak{u}E)_{c},Q_{c}\} of L ^ ^ 𝐿 \hat{L} and L ^ † superscript ^ 𝐿 † \hat{L}^{\dagger} for N k = 24 subscript 𝑁 𝑘 24 N_{k}=24 Fourier components and constants { λ 1 , λ 2 , λ 3 , λ 4 } subscript 𝜆 1 subscript 𝜆 2 subscript 𝜆 3 subscript 𝜆 4 \{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}\} . The corresponding coefficients Λ ​ ( 𝔲 ​ E ) Λ 𝔲 𝐸 \Lambda(\mathfrak{u}E) are given with the values of the eigenvectors (see text).

$\nu$	${\left(𝔲E\right)}_c$	$Q_c$	${\lambda }_1$	${\lambda }_2$	${\lambda }_3$	${\lambda }_4$
1	0.320	0.486	0.166	-0.503	-0.176	2.07
2	0.276	0.537	0.102	-0.230	-0.0787	1.97
3	0.21	0.335	-0.0478	0.128	0.030	2.93
4	0.189	0.359	0.0265	-0.0671	-0.0169	2.74
5	0.156	0.256	0.0133	-0.0345	-0.00570	3.83
6	0.144	0.27	0.0074	-0.0183	-0.00258	3.63
7	0.124	0.208	0.00348	-0.00922	-0.000833	4.74
8	0.116	0.217	0.00199	-0.00528	-0.000436	4.54

Equations

$$\mathfrak{u}\left({\partial }_t+i\mu \right)\psi ={\partial }_x^2\psi +\left(1-|\psi {|}^2\right)\psi ,$$

$$\mathfrak{u}={\mathfrak{u}}_0\frac{{\pi }^4}{14\zeta (3)}\frac{e^2\nu \mathfrak{D}}{\sigma },$$

$$j=j_{\mathrm{n}}+j_{\mathrm{s}}=-{\partial }_x\mu -ı\mathrm{Im}({\psi }^{\ast }{\partial }_x\psi ).$$

$$\mu =-Ex+\tilde{\mu },$$

$$\mathfrak{u}\left[{\partial }_t+ı(\tilde{\mu }+{\mu }_{\mathrm{ps}}-E\tilde{x})\right]\psi ={\partial }_x^2\psi +\left(1-|\psi {|}^2\right)\psi .$$

$${\mu }_{\mathrm{ps}}=E{x}_{\mathrm{ps}}{t}_{\mathrm{ps}}\sum _k\delta (t-t_{\mathrm{ps}}k)\times$$

$$\sum _{m}m\Theta \left[x-x_{\mathrm{ps}}\left(m-\frac{1}{2}\right)\right]\Theta \left[x_{\mathrm{ps}}\left(m+\frac{1}{2}\right)-x\right],$$

$$\psi =\sum _{k=-\infty }^{\infty }{A}_k{e}^{ıQkx},$$

$$j=E+Q\sum _{k}k|A_k{|}^2.$$

$$-ı\mathfrak{u}\left\{Ex+\frac{ı}{2}\sum _{k\ne 0}{k}^{-1}[A_1(2-k)A_{1-k}^{\ast }+A_1^{\ast }(2+k)A_{k+1}]e^{ıkQx}-\frac{ı}{2}\sum _{l\ne k\ne 1}\frac{k+l}{k-l}{A}_k^{\ast }{A}_l{e}^{-ı\left(k-l\right)Qx})\right\}\sum _{k=-\infty }^{\infty }{A}_k{e}^{ıkQx}$$

$$=(1-|A_1{|}^2-Q^2)A_1{e}^{ıQx}+\sum _{k\ne 1}\{(1-2|A_1{|}^2-k^2{Q}^2)A_k{e}^{ıkQx}-A_1^2{A}_k^{\ast }{e}^{-ı(k-2)Qx}$$

$$-\sum _{l\ne 1}\left[2A_1{A}_k{A}_l^{\ast }{e}^{ı(k-l+1)Qx}+A_1^{\ast }{A}_k{A}_l{e}^{ı(k+l-1)Qx}+\sum _{m\ne 1}{A}_k{A}_l{A}_m^{\ast }{e}^{ı(k+l-m)Qx}\right]\}.$$

$$\tilde{\mu }=-\frac{ı}{2}\left\{\sum _{k\ne 0}\frac{e^{ıkQx}}{k}\left[A_1(2-k)A_{1-k}^{\ast }+A_1^{\ast }{A}_{k+1}(2+k)\right]-\sum _{k\ne l\ne 1}\left(\frac{k+l}{k-l}\right)e^{-ı(k-l)Qx}{A}_k^{\ast }{A}_l\right\}.$$

$$j=j_{\mathrm{dp}}+E-E^2\gamma (\mathfrak{u}),$$

$$\gamma (\mathfrak{u})=-3\sqrt{3}{\mathfrak{u}}^2\frac{\partial }{\partial \mathfrak{u}}\left[\frac{\pi }{2\sqrt{2\mathfrak{u}}}\left\{\coth (\pi \sqrt{2\mathfrak{u}})-\frac{1}{\pi \sqrt{2\mathfrak{u}}}\right\}\left(1-\frac{6}{\mathfrak{u}}+\frac{12}{{\mathfrak{u}}^2}\right)+\frac{{\pi }^2}{\mathfrak{u}}\left(1-\frac{2}{\mathfrak{u}}\right)+\frac{4{\pi }^4}{15\mathfrak{u}}\right].$$

$$\gamma (\mathfrak{u})=\frac{18}{\sqrt{3}}{\mathfrak{u}}^2\left(\frac{61}{(1+2\mathfrak{u}{)}^2}+\frac{7}{(4+2\mathfrak{u}{)}^2}+0.0486531-0.0185438\mathfrak{u}\right).$$

$$\gamma (\mathfrak{u})=\sqrt{3}\left(3{\pi }^2+\frac{4{\pi }^4}{5}-\frac{3}{4}+\frac{3\pi }{8}\sqrt{2\mathfrak{u}}\right)-\frac{27\pi }{2}\sqrt{\frac{3}{2\mathfrak{u}}}-\frac{3\left(\sqrt{3}\left(4{\pi }^2-3\right)\right)}{\mathfrak{u}}+\sqrt{\frac{3}{2}}\frac{45\pi }{{\mathfrak{u}}^{3/2}}-\frac{27\sqrt{3}}{{\mathfrak{u}}^2}+\mathcal{O}\left(\mathfrak{u}{e}^{-2\pi \sqrt{2\mathfrak{u}}}\right).$$

$$\frac{\partial (\mathfrak{u}E)}{\partial Q}=0.$$

$$\hat{L}\mathbf{f}=0,$$

$${\hat{L}}_{k,l}=\{\begin{array}{ll}{Z}_{k-1}{\delta }_{l,1}+Z_{k-l}(1-{\delta }_{l,1})+(Q^2{k}^2-1){\delta }_{k,l}& ,k\ne 1\\ Z_{1-l}-(1-Q^2){\delta }_{l,1}& ,k=1\end{array}$$

$$Z_k=\{\begin{array}{ll}\frac{\mathfrak{u}E}{kQ}(-1{)}^k& ,k\ne 0\\ 0& ,k=0\end{array}.$$

$$\{f_k\}\leftrightarrow \{A_{k\ge 2},A_1,A_{k\le 0}\},$$

$$\delta {\hat{L}}_{k,l}=\delta Z_{k-1}{\delta }_{l,1}+\delta Z_{k-l}(1-{\delta }_{l,1})$$

$$\delta Z_k=\{\begin{array}{ll}\frac{\mathfrak{u}\left(E-E_c\right)}{k{Q}_c}{\left(-1\right)}^k& ,k\ne 0\\ 0& ,k=0\end{array}$$

$$\lambda (\mathfrak{u},E)f_k,$$

$$|\lambda (\mathfrak{u},E){|}^2\left(\mathfrak{u}{\lambda }_1+{\lambda }_2\right)+\sum _{k,k_1}{\tilde{f}}_k^{\ast }\left(\delta {\hat{L}}_{k,k_1}\right)f_{k_1}=0,$$

$${\lambda }_3\mathfrak{u}(E-E_c)=\sum _{k,k_1}\widetilde{f_k^{\ast }}(\delta {\hat{L}}_{k,k_1})f_{k_1},$$

$$j=E+|\lambda (\mathfrak{u},E){|}^2{\lambda }_{4}Q.$$

$$\sum _{k=2}^{\infty }k|A_k{|}^2+|A_1{|}^2-\sum _{k=1}^{\infty }k|A_{-k}{|}^2=0.$$

$$\left(\begin{array}[]{cc}\left[1-2|A_{1}|^{2}-\frac{\mathfrak{u}|A_{1}|^{2}\left(k+1\right)}{2\left(k-1\right)}-k^{2}Q^{2}\right]&-A_{1}^{2}\left(1+\frac{u\left(3-k\right)}{2\left(k-1\right)}\right)\\ -\left(A_{1}^{*}\right)^{2}\left(1-\frac{u\left(k+1\right)}{2\left(k-1\right)}\right)&\left[1-2|A_{1}|^{2}+\frac{u|A_{1}|^{2}\left(3-k\right)}{2\left(k-1\right)}-\left(k-2\right)^{2}Q^{2}\right]\end{array}\right)\left(\begin{array}[]{c}A_{k}\\ A^{*}_{-k+2}\end{array}\right)=\frac{\mathfrak{u}E\left(-1\right)^{k}}{\left(k-1\right)Q}\left(\begin{array}[]{c}-A_{1}\\ A_{1}^{*}\end{array}\right)$$

$$A_k=-\frac{1}{\mathcal{D}}\frac{\mathfrak{u}E{\left(-1\right)}^k{A}_1}{\left(k-1\right)Q}\left[1-3|A_1{|}^2-{\left(k-2\right)}^2{Q}^2\right],A_{-k+2}^{\ast }=\frac{1}{\mathcal{D}}\frac{\mathfrak{u}E{\left(-1\right)}^k{A}_1^{\ast }}{\left(k-1\right)Q}\left[1-3|A_1{|}^2-k^2{Q}^2\right],$$

$$\mathcal{D}={\left[1-2|A_1{|}^2-\left({\left(k-1\right)}^2+1\right)Q^2-\frac{\mathfrak{u}|A_1{|}^2}{2}\right]}^2-{\left[2\left(k-1\right)Q^2+\frac{\mathfrak{u}|A_1{|}^2}{k-1}\right]}^2-|A_1{|}^4\left[{\left(1-\frac{\mathfrak{u}}{2}\right)}^2-\frac{{\mathfrak{u}}^2}{{\left(k-1\right)}^2}\right].$$

$$\mathcal{D}=\frac{1}{9}(k-1{)}^2((k-1{)}^2+2\mathfrak{u})$$

Full text

reentrance current, phase-slips, critical current, time-dependent Ginzburg-Landau equation

I Introduction

Narrow superconducting strips are the subject of great interest for superconducting quantum electronic devices. Their dissipation-less state is very subtle and sensitive to thermal and quantum fluctuations, which can easily flip the superconducting strip into the resistive state, making them ideal candidates for very sensitive detectors. Various models have been proposed to explain the appearance of non-zero resistance in these strips and its temperature dependence in the region of low temperatures (for a review, see Refs. [Sahu et al., 2009; Bezryadin, 2012]).

The role of thermal fluctuations responsible for energy dissipation, when current flows through a one-dimensional superconductor, was considered for the first time in the seminal paper by Langer and AmbegaokarLanger and Ambegaokar (1967) over fifty years ago. Note, that a realistic “one-dimensional superconductor” is in fact a narrow strip with finite width $W$, much less than the Ginzburg-Landau coherence length $\xi(\tau)\propto(k_{\mathrm{B}}T_{\mathrm{c}}\tau)^{-1/2}$, where $\tau=1-T/T_{\mathrm{c}}$ is the reduced temperature and $T_{\mathrm{c}}$ the critical temperature. The energy dissipation in this system is related to phase-slip processes appearing in thin superconducting wires Tinkham (1996); Skocpol et al. (1974); Kramer and Baratoff (1977); Kramer and Rangel (1984); Rangel and Kramer (1989); Lau et al. (2001); Mooij and Nazarov (2006); McKay et al. (2008); Kimmel et al. (2017) or superfluids Glatz and Nattermann (2002); Scherpelz et al. (2014, 2015), i.e., the processes of vortices/flux quanta crossing the strip.

It is clear, that such events cannot be realized in the framework of a purely one-dimensional model. Indeed, as it was shown in Ref. [Langer and Ambegaokar, 1967], the minimal value of the order parameter magnitude remains finite and equal to $\left(2/3\right)^{1/2}\Delta_{\mathrm{BCS}}$ even when the density of current flowing through the one-dimensional superconductor reaches the “depairing” value, $j_{\mathrm{dp}}\propto e\nu\left(k_{\mathrm{B}}T\tau\right)^{2}\xi(\tau)$, and global superconductivity in the one-dimensional channel is only partially suppressed. Yet, the order parameter should become zero at least in one point of the strip in order to allow the system to perform a phase slip event. 111Note, that, in contrast to the one-dimensional superconducting strip, the order parameter in a higher-dimensional superconductor becomes zero when the current density approaches the “depairing” value.

The apparent paradox occurring in the one-dimensional case was resolved in Ref. [Ovchinnikov and Varlamov, 2015]. The authors demonstrated that the saddle point solutions of the static Ginzburg-Landau (GL) equations for the order parameter $\widetilde{\psi}$ in the presence of a fixed current density $j$, possessing at least one vortex, exist only for very weak current densities $j<j_{\mathrm{c}}=0.0312(W/\xi)j_{\mathrm{dp}}$. In the case $j>j_{\mathrm{c}}$, such phase-slip events are possible, which are random in space and time due to thermal fluctuations. When the current density exceeds the value $j_{\mathrm{dp}}$, the static scenario described above does not hold anymore. In this region ($j>j_{\mathrm{dp}}$) dynamical states are formed in the strip, where phase-slip events occur periodically in space and time 222Note that material defects and weak links also lead to phase-slip events Kramer and Watts-Tobin (1978); Watts-Tobin et al. (1981); Berdiyorov et al. (2014); Kimmel et al. (2017), but here we concentrate on the clean case..

In this paper we study the hysteretic structure of such dynamical states of a narrow superconducting strip and obtain the corresponding current-voltage ($I$–$V\,$) characteristics. In particular, we derive the value of the re-entrance current density $j_{\mathrm{r}}$, at the onset of the instability of the normal state, when the applied current decreases. Our consideration is valid in the gapless region, at temperatures slightly below the critical one, where the time-dependent Ginzburg-Landau equation holds. In this situation, we investigate the strong current regime $j>j_{\mathrm{dp}}$ of the superconducting strip being in its dynamical state up to second order perturbation theory in the electric field $E$.

It is important to note, that the considered situation is quite different from the problem of determining the critical current at which the superconductor becomes normal (see Refs. [Ivlev and Kopnin, 1984; Kopnin and Ivlev, 1984]). The authors of the cited papers, which are based on the works by Gorkov Gorkov (1970) and Kulik Kulik (1971), claim that the normal state remains stable for any finite value of $E$. However, they ignore the fact that the exponential growth of superconducting fluctuations in a time interval determined by $j\Delta t<\hbar\sigma/[e\xi(\tau)]$ leads to an instability of the normal state ($\sigma$ is the normal conductivity). As a result, the system enters a dynamical superconducting state at some finite electric field. We will show below that there are many values of the parameter $\mathfrak{u}E$ for which the normal state starts to be unstable even for infinitely small perturbations.

In the following we describe the model and show the analysis of the time-dependent Ginzburg-Landau equation near the critical point of the normal state. We derive the value of the reentrance current and order parameter values using first order perturbation theory by small deviation from the critical point. Details of the calculations can be found in the Appendices. We start with introducing the model in the following section and then analyze the $I$–$V\,$characteristics close to the deparing current and near the instability points of the normal state.

II Model

In this paper we approximate the narrow superconducting strip of width $W$ smaller than the superconducting coherence length by a one-dimensional (1D) system, described by the time-dependent Ginzburg-Landau equation (TDGLE). The TDGLE can be written in dimensionless variables, without accounting for thermal fluctuations and magnetic field (the latter does not appear in the 1D model) in the form

[TABLE]

where $\psi$ is the complex order parameter and $\mu$ is the scalar potential. The reduced relaxation rate $\mathfrak{u}$ controls the system’s evolution in time and is given by

[TABLE]

where $\nu$ is the density of states at the Fermi surface, $\mathfrak{D}$ is the effective diffusion constant, $\zeta$ is the Riemann-Zeta function, and $\mathfrak{u}_{0}$ is a numerical constant.Ovchinnikov (1997); Kopnin (2001)

The fixed total bias-current consists of the sum of normal ($j_{\mathrm{n}}$) and superconducting ($j_{\mathrm{s}}$) components and in its turn can be related to the space derivatives of the complex order parameter and the scalar potential:

[TABLE]

Time and distance in Eqs. (1)-(3) are measured in units of $t_{0}=8\pi\sigma\lambda^{2}/c^{2}$ and scaled superconducting coherence length, $\xi=\left(\frac{\pi\mathfrak{D}}{16k_{B}T_{\mathrm{c}}\tau}\right)^{1/2}$, respectively, with $\lambda$ being the London penetration depth. The electrical current density $j$ is measured in units of $j_{0}=c\Phi_{0}/(8\pi^{2}\lambda^{2}\xi\sqrt{2})$ ($\Phi_{0}$ is the flux quantum). In these units the depairing current density reads as $j_{\mathrm{dp}}=2/(3\sqrt{3})j_{0}\approx 0.385j_{0}$.

In what follows, it is useful to transform the scalar potential to the form

[TABLE]

where $E$ is the average electric field, which is equal to the average normal current in dimensionless units, and $\tilde{\mu}$ is the spatially and temporally fluctuating part of $\mu$.

Phase-slip events generate instabilities in the current-voltage characteristics, which will be the main subject of this work. In particular, these processes cause strong changes in the electric field and order parameter at their space-time coordinates. We take their effect into account explicitly by introducing a corresponding effective electric potential $\tilde{\mu}$ to the TDGLE.

A phase-slip process in a strip of finite width is related to the transfer of a magnetic vortex-antivortex pair across it. Each such event is accompanied by the suppression of the order parameter and, consequently, of the supercurrent. Due to the conservation of the total current, a sharp peak in the normal current appears at the time and space location of the phase-slip event.

In the large-current regime, $j>j_{\mathrm{dp}}$, these phase-slip events (ps) are periodic in time and space. The corresponding periods in time and space are denoted as $t_{\mathrm{ps}}$ and $x_{\mathrm{ps}}$, respectively. Outside a very narrow region in the time-space plane one can rewrite Eq. (1) in the form

[TABLE]

Here $\tilde{x}=x\mod(x_{\mathrm{ps}}/2)$ and the associated potential $\mu_{\mathrm{ps}}$

[TABLE]

where the quantization condition $Ex_{\mathrm{ps}}t_{\mathrm{ps}}=2\pi$ is implied. The shape of $\mu_{\mathrm{ps}}$ is illustrated in Fig. 1 (top) with smoothed step and delta functions. In the bottom panels $\mu_{\mathrm{ps}}$ and $j_{\mathrm{ps}}\propto\partial_{x}\mu_{\mathrm{ps}}$ are shown in the elementary space-time cell $\mathcal{A}_{\mathrm{ps}}\equiv]-x_{\mathrm{ps}}/2;x_{\mathrm{ps}}/2[\times]0;t_{\mathrm{ps}}[$. The latter illustrating the mentioned spikes in the normal current at the phase slip event location.

These additional terms do not contribute to the electric field inside $\mathcal{A}_{\mathrm{ps}}$ but only at its corners. The solution for the order parameter $\psi$ away from the corners of $\mathcal{A}_{\mathrm{ps}}$ follows from Eq. (5) with periodic boundary conditions. The size of the region with strong suppression of $\psi$ is of order $W\times Wt_{\mathrm{ps}}/x_{\mathrm{ps}}$, where $W$ is the width of the strip.

III Analysis of the $I$–$V\,$characteristics

Let us start with the analysis of Eqs. (1)-(4). Within the unit cell $\mathcal{A}_{\mathrm{ps}}$, one can make the Fourier-Ansatz for the complex order parameter

[TABLE]

where $Q$ should be found from local minimum conditions of $E$ for a given current density $j$.

Correspondingly, Eq. (3) acquires the form

[TABLE]

Plugging Ansatz (7) into Eq. (1) gives

[TABLE]

Finally, accounting for the charge conservation condition $\textrm{div}j=0$ and Eq. (3), one can find the expression for the fluctuating part of the scalar potential $\tilde{\mu}$ in terms of the introduced Fourier coefficients:

[TABLE]

In Eqs. (9)-(10) we explicitly separated the quantity $A_{1}$, since it is the dominant Fourier component in the vicinity of the critical point $j=j_{\mathrm{dp}}$ and plays an important role throughout the paper. Fourier coefficients with $|k|\gg 1$ quickly decay. All other coefficients in this region can be found in the framework of perturbation theory.

Eq. (9) enables us to obtain the $I$–$V\,$characteristics in the complete domain of dynamical resistive states. The above mentioned minimization with respect to $E$ allows us to finds the shape of the $I$–$V\,$characteristics, which turns out to be critically dependent on the value of the dynamic coefficient $\mathfrak{u}$ of the TDGLE.

One can expect to find three qualitatively different types of $I$–$V\,$characteristics, which are illustrated in Fig. 2. The first one (Fig. 2a) is reversible and could be realized when the dynamic coefficient $\mathfrak{u}$ is sufficiently small: $\mathfrak{u}<\mathfrak{u}_{c}^{(1)}$, where $\mathfrak{u}_{c}^{(1)}$ is the first critical value, which will be obtained below. In case $\mathfrak{u}>\mathfrak{u}_{c}^{(1)}$ the $I$–$V\,$characteristics becomes irreversible (Fig. 2b,c).

In the case of strong damping, when $\mathfrak{u}$ exceeds the second critical value $\mathfrak{u}_{c}^{(2)}$, the transition to a finite value of the order parameter happens at a smaller $E_{c}$ value, which is where we define the reentrance current $j_{r}$. Calculation of the order parameter value in the vicinity of the critical point $(\mathfrak{u}E)_{c}$ (see below) shows that in practice only the latter scenario of the $I$–$V\,$characteristics, shown in Fig. 2c), is realized.

III.1 Current density close to the depairing current

Next, we consider the case when the bias current density $j$ is close to $j_{\mathrm{dp}}$. Eq. (1) allows us to obtain the above mentioned value of $j_{\mathrm{dp}}$, which destroys superconductivity in the 1D channel, where the corresponding critical value of the order parameter is $\psi_{c}=\sqrt{2/3}$, and the value of wave-vector (for $E=0$) $Q_{c}=1/\sqrt{3}$. Close to this point, Eq. (9) is decomposed into two equations with $\left\{k,2-k\right\}$. A detailed analysis of the Fourier coefficients and the determination of the wave vector $Q$ is presented in Appendix A. As a result of these calculations we obtain the electric field dependence of the current density to second order, close to its depairing value as

[TABLE]

where the calculation of the coefficient $\gamma(\mathfrak{u})$ of the quadratic term as function of $\mathfrak{u}$ is a highly involved task, but can be performed exactly for any value of $\mathfrak{u}$ and be expressed as a full derivative:

[TABLE]

The explicit expression of the full derivative is given in Appendix A, Eq. (33), and $\mathcal{D}$ is defined there in Eqs. (26b) and (27). Note, that the difference in braces behaves as $\pi\sqrt{2\mathfrak{u}}/3$ for small $\mathfrak{u}$, such that all terms $\propto\mathfrak{u}^{-2}$ under the derivative cancel and the complete expression is non-singular at $\mathfrak{u}=0$. Therefore, in the limit of small $\mathfrak{u}$ we keep the first two terms of the sum in (12) and expand the remaining sum to first order. This gives

[TABLE]

For large values of $\mathfrak{u}$ ($\mathfrak{u}\gg 1$), we obtain from Eq. (12) (using the asymptotic expressions of the $\coth$ and $\sinh^{-2}$ terms)

[TABLE]

The $\mathfrak{u}$-dependence of $\gamma(\mathfrak{u})$ and the domains of validity of its approximations (13)-(14) are presented in Fig. 3. $\mathfrak{u}\sim 1$ separates the regions where small-$\mathfrak{u}$ and large-$\mathfrak{u}$ approximations work best, i.e., the relative deviations from the exact curve are both minimum at $\mathfrak{u}=1.061$, less than $10^{-3}$.

III.2 Vicinity of the critical points

We now consider the vicinity of critical points, $(\mathfrak{u}E)_{c}$, which are defined by the condition

[TABLE]

In this region we search for the solution to the non-linear problem (9) by its linear expansion over eigenfunctions. Therefore the linearized condition of Eq. (15) can be understood as the following eigenvalue problem

[TABLE]

where the form of the linear operator $\hat{L}$ follows from Eq. (9):

[TABLE]

with

[TABLE]

The eigenvector $\mathbf{f}$ of Eq. (16) is related to the Fourier coefficients at the critical points in the following way

[TABLE]

where $f_{1}\leftrightarrow A_{1}$.

The above eigenvalue problem has possibly an infinite number of solutions $\{(\mathfrak{u}E)_{c},Q_{c}\}$. Note, that the quantities $\mathfrak{u}$ and $E$ only appear in the form of a product at the critical point [in contrast to Eqs. (11) & (12)]. The linearized equation can be solved numerically and the largest eight critical points (defined by simultaneous eigenvalues of $\hat{L}$ and $\hat{L}^{\dagger}$) are listed in table 1 [the corresponding normalized eigenvectors of $\hat{L}$ and $\hat{L}^{\dagger}$ (here $\hat{L}^{\dagger}=\hat{L}^{\mathsf{T}}\neq\hat{L}$) are given in the supplementary information]. Below we discuss the numerical solution in more detail.

In order to get further insight into the behavior of the $I$–$V\,$characteristic near these critical points and to determine the type of instability point (first or second order transition), we introduce the operator $\delta\hat{L}$ as

[TABLE]

with

[TABLE]

[compare to Eq. (18)].

Near a critical point $\left(\mathfrak{u}E\right)_{c}$ the solution for the Fourier components in Eq. (9) in our linearized approximation can be written in the form

[TABLE]

with the proper permutation of $k$-indices as defined in Eq. (19) and where the coefficient $\lambda(\mathfrak{u},E)$ follows from the solvability condition

[TABLE]

see Appendix B for explicit expressions for all coefficients $\lambda_{i}$.

The expression for current density then takes the form

[TABLE]

(see Appendix B for definition of $\lambda_{4}$.)

We note, that in the critical region only the coefficient of the zero mode, $\lambda(\mathfrak{u},E)$, is dominant and all other coefficients are small, scaling with $|(\mathfrak{u}E)_{c}-\mathfrak{u}E|/(\mathfrak{u}E)_{c}$.

Altogether, we can now analyze the critical point in detail. Therefore we calculate the critical points $\{(\mathfrak{u}E)_{c},Q_{c}\}$ numerically for truncated Fourier series with $N_{k}$ components (indices $k\in\{-N_{k}/2+1,\ldots,N_{k}/2\}$). These are obtained as simultaneous solutions of the polynomial equations $Q^{N_{k}}\det(\hat{L})=0$ and $Q^{N_{k}+1}\partial_{Q}\det(\hat{L})=0$ of order $3N_{k}-2$ in $Q$. Figure 4 shows the solutions for $N_{k}=24$ (order $Q^{70}$), where the solid lines represent the solutions of the individual equations. Note, that solving the linearized equations for the truncated Fourier series leaves the largest critical values $\{(\mathfrak{u}E)_{c},Q_{c}\}$ invariant for sufficiently large $N_{k}$. For these solutions we can then obtain the eigenvectors of $\hat{L}$ and $\hat{L}^{\dagger}$, which allows us then to extract the behavior of the $I$–$V\,$characteristic near these critical points by evaluation of the parameters $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$. At those points new branches appear, which can bring the system out of the normal state. Using Eq. (23b) then defines the slope of the linearized $I$–$V\,$characteristic. The numerical calculation reveals that the critical value $\nu=3$, $(\mathfrak{u}E)_{c}^{(3)}\approx 0.21$, has locally a negative slope among the eight largest $(\mathfrak{u}E)_{c}$ values, indicating that a reentrance into the superconducting state can happen without a threshold (second order) at a current density $j_{\mathrm{r}}$. Fig. 5 shows the behavior of the order parameter $\Delta$, (a), and current $j$, (b), at the critical points as a function of $\mathfrak{u}E$. The connection of the critical points to the envelop (defined by the $I$–$V\,$curve for increasing current) is indicated by dotted lines (these cannot be realized physically). Practically, one can make a ‘hysteresis ‘loop” in the $j$-$(\mathfrak{u}E)$ diagram by starting in the superconducting state at zero current, following up to $j_{\mathrm{c}}$ upon increasing $j$, where the system becomes resistive and eventually jumps into to normal state (indicated by an arrow) following the normal $I$–$V\,$curve (blue, dashed). When decreasing $j$ from the normal state, one follows down the normal $I$–$V\,$-line till the critical point labeled $\ast$ ($\nu=3$), where the slope of $I$–$V\,$is negative, such that fluctuating superconducting regions can grow and one jumps back into the superconducting state at $j_{\mathrm{r}}$ (indicated by an arrow). Below this point the normal state is always unstable. At all other (larger) critical points we cannot follow the critical $I$–$V\,$(without threshold) as the slope is positive (first order). We note that the specific picture depends on the actual value of $\mathfrak{u}$ for the physical system under consideration, here we assume a value of $\mathfrak{u}$ of order one.

The numerical analysis demonstrates that probably an infinite set of solutions of the eigenvalue problem, Eq. (16) exists.

In analogy to finding the (global) extremum of a function on a finite support, where also the boundary values need to be checked, here we should also study the properties of the system close to the hypothetical “endpoints”, if those exist [besides the (local) critical points defined by (16)]. By “endpoints” we mean points of the surface $\left\{\hat{L}\psi=0\right\}$ in Hilbert space, where the value $\left(\mathfrak{u}E\right)$ reaches its maxima under the condition

[TABLE]

Our numerical evaluation of this condition reveal that such “endpoints” are irrelevant.

IV Conclusions

We have investigated the $I$–$V\,$characteristic of a superconducting strip in the region near the depairing current $j_{\mathrm{dp}}$ and studied the instability points of the normal state as function of the current. Interestingly, one finds a degeneracy in second order perturbation theory in the electric field $E$ by solving linearized equation at those critical points.

This degeneracy leads to the appearance of additional branches splitting off from the Ohmic behavior seen in Fig. 5. Numerically, we calculated the critical points and found that the largest electric field value, at which the normal state first becomes unstable upon decreasing the current, i.e., indicating the possibility of a transition into the superconducting state with finite value of the order parameter amplitude (see Fig. 5a), and a branch in the $I$–$V\,$appears is at $\mathfrak{u}E=0.3199$. However, the slope of the branch is positive, indicating a first order transition, which we can only follow when increasing the current (and eventually jumping back in the normal state). Therefore, when increasing the current for fixed $\mathfrak{u}E$ in the intervals $[0.21,0.275[$ or $]0.275,0.319]$, a transition into the superconducting state will happen at $\mathfrak{u}E=0.275$ or $\mathfrak{u}E=0.319$, respectively, before returning to the normal state at larger currents.

However, most importantly, we also found the smallest value for $\mathfrak{u}E$ when the normal state is always unstable to be equal to $0.2095$ (indicated by $\ast$ in Fig. 5), defining the reentrance current into the superconducting state. In contrast to the evaluation of the critical current, the evaluation of the re-entrance current is significantly more involved.

Acknowledgments

We are delighted to thank I.S. Aranson for interesting discussions. The research was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division.

Yu.N.O acknowledges Prof. Dr. Jeroen van den Brink for his hospitality in LIFW and the DFG for granting him a Mercator-Professoren-Stipendium. A. V. acknowledges financial support from the project CoExAn (HORIZON 2020, grant agreement 644076), from Italian MIUR through the PRIN 2015 program (Contract No. 2015C5SEJJ001).

Appendix A Solution of TDGLE near the depairing current

Close to the depairing current, Eq. (9) is decomposed into pairwise equations with $\left\{k,-k+2\right\}$.

For $k\neq$1

[TABLE]

Solving this system results in

[TABLE]

In our approximation, we obtain from Eq. (26b)

[TABLE]

when the electric field $E$ is much smaller than its critical value $E_{c}$.

For $k=1$, we obtain from Eq. (9), the following equation for the quantity $|A_{1}|^{2}$ in second order perturbation theory

[TABLE]

Using second order perturbation theory, we can set

[TABLE]

where $\{\alpha,\beta\}$ are constants. Inserting expression (26a) for the coefficients $A_{k}$ and expressions (29) into Eq. (28), we get the following relation for those constants:

[TABLE]

One important property of Eqs. (28) and (29) should be mentioned: in second order perturbation theory, defined by Eq. (30), $\alpha$ and $\beta$ appear only in combination $\alpha+\sqrt{2}\beta$. This implies that corrections to the quantities $|A_{1}|$ and $Q$ appear separately in perturbation theory only in order $\mathcal{O}\left(E^{4}\right)$

In the same approximation we obtain from Eq. (8)

[TABLE]

Inserting Eq. (30) into expression (31) yields an expression for the current $j$ in second order perturbation theory in $E$:

[TABLE]

The expression for $\gamma(\mathfrak{u})$ in (32b) can be evaluated explicitly as

[TABLE]

which can be written as the full derivative Eq. (12).

Appendix B Solvability and $I$–$V\,$at the critical points

The coefficients in the solvability equations, (23a) and (23b), are given by

[TABLE]

Here $\widetilde{f}_{k}$ are the components of the normalized eigenvector of the transposed operator $\hat{L}^{\dagger}=\hat{L}^{\mathsf{T}}$.

Next, we define the function

[TABLE]

where $f_{k}$ are the components of the eigenvector of $\hat{L}$ and $\widetilde{f}_{k}$ are the components of the normalized eigenvector of the operator $\hat{L}^{\dagger}$. Function $\delta\hat{L}_{k,k_{1}}$ is defined in Eq. (20). Eq. (9) in the vicinity of each critical point $\left\{\left(\mathfrak{u}E\right)_{c},Q_{c}\right\}$ can then be rewritten in the form

[TABLE]

where the coefficients are given in (34a)-(34c). Therefore:

[TABLE]

The functions $\Lambda(\mathfrak{u}E)$ corresponding to the eight largest critical value $(\mathfrak{u}E)_{c}$ are given by

[TABLE]

Note the signs of the coefficients in $\Lambda^{(3)}$. For completeness, we reproduce the corresponding eigenvectors $f_{k}$ and $\widetilde{f}_{k}$ in the supplementary information.

Appendix C Derivation of Fourier equations

Here we will obtain the equation system for the Fourier coefficients $A_{k}$. From Eq. (9) we obtain for $k_{0}\neq 1$

[TABLE]

From Eq. (9) we obtain a separate equation for $k_{0}=1$:

[TABLE]

Next, we introduce the operator $\hat{M}_{k_{0}}$ for $k_{0}\neq 1$ as

[TABLE]

With this, Eq. (1) can be written in the form

[TABLE]

where $\Phi_{k_{0}}$ is

[TABLE]

Note, that a free parameter $Q$ appears in Eqs. (37) and (39). The value of this parameter is found by the extremal condition for the electric field $E$ for a given current density $j$.

Appendix D Supplementary Information: Eigenvectors near critical points

Here we present detailed results of the numerical evaluation at critical points including the eigenvectors. All calculations are done for $N_{k}=24$ (Fourier) components of the eigenvectors of $\hat{L}$ and $\hat{L}^{\dagger}$, $f_{k}$ and $\widetilde{f}_{k}$, respectively. Here component indices range from $k=-N_{k}/2+1\ldots N_{k}/2$.

The results for the eight largest critical $(\mathfrak{u}E)_{c}$ are reproduced in table 2, including the critical structural constant $Q_{c}$, and the four parameters $\lambda_{i}$ ($i=1,2,3,4$); see manuscript. The bold printed critical value 3 corresponds to the reentrance point.

The normalize eigenvectors and functions $\Lambda(\mathfrak{u}E)$ (see manuscript) are listed below, where bracketed superscripts correspond to the index $\nu$ in table 2. The $k=1$component is printed in bold.

[TABLE]

[TABLE]

Note, the eigenvectors are only defined up to an overall $(-1)$ factor.