Dynamical instability of the electric transport in superconductors

We develop a nonlinear theory of the electronic transport in superconductors in the framework of the time-dependent Ginzburg-Landau (TDGL) equation. We utilize self-consistent Gaussian approximation and reveal the conditions under which the current-voltage V(I) dependence (I–V characteristics) acquires an S-shape form leading to switching instabilities. We demonstrate that in two-dimensions the emergence of such an instability is a hallmark of the Berezinskii-Kosterlitz-Thouless (BKT) transition that we have detected by transport measurements of titanium nitride (TiN) films. Our theoretical findings compare favorably with our experimental results.

A panoply of nonlinear and instability effects in electronic transport in superconductors includes instabilities stemming from overheating and vortex contraction effects and switching instabilities triggered by vortex depinning-like instabilities 1 . The signatures of these instabilities are the related S-shaped I(V) dependences. The common view assign the emergence of S-like I(V) curves to the classical overheating of an electron gas by the external field [1][2][3][4] . If a bottleneck of the relaxation process is the energy transfer from electronic system to the thermal bath, the effective temperature of electrons T e becomes higher than the temperature of phonons T ph . A steep temperature dependence of the resistance R(T) is viewed as a signal that the system in question has a tendency towards an instability 2 . In particular, thermal instabilities are likely to develop near the superconducting transition temperature T c . The reason is the sharp R(T) dependence in the critical vicinity of T c primarily due to divergent Aslamazov-Larkin contribution ∼ T T 1/ ln( / ) c . The applicability of the above model to superconducting films was examined in work 5 and it was revealed, in particular, that the power-balance equation 3,4 correctly describes the behavior of the I-V curves of thin superconducting films at T > T c . At these temperatures, however, bistability is absent in films and develops only at T < T c .
At the same time, the early classical works by Gorkov 6 and Masker, Marcelja and Parks 7 suggested the existence of the intrinsic I-V, purely superconducting instabilities caused by the superconducting (SC) fluctuations rather than the overheating effects. In one dimension the analysis presented by Tucker and Halperin 8 utilized the dynamical Hartree -Fock approximation and indicated the possibility for instability. Similar instability was found in 9 for 1D TDGL. Experiments on quasi-1D and 2D superconducting systems revealed instability effects [10][11][12][13] , where voltage jumps were found in the current driven settings [14][15][16][17] . Moreover, experiments on systems with Cooper pair chanel revealed that while the high-resistive branches of the I(V) curves may fairly well fit the overheating behavior, the latter fails to describe the low-resistance branches 18 . This calls for a thorough understanding of the role of the intrinsic SC fluctuations-induced instabilities not related to overheating effects and the associated with them emergence of hot spots and/or moving heating domains. Our paper takes up on the task and expands the self-consistent approach following refs 19-23 that provides a quantitative description of transport instabilities in superconductors.

Thermal Fluctuations and Electric Field in TDGL Equation
The model. The GL free energy has a standard form 23 with Ψ being the order parameter, D being the dimension of the system and A being either the cross section for the wire (for the 1D case) or the thickness of the film. In the GL potential term, T c0 is the mean-field critical temperature. The GL coefficients define as usual the superconducting coherence length ξ 2 = ℏ 2 /(2m * αT c0 ) and the London penetration depth λ 2 = bc 2 m * /(16πe 2 αT c0 ). The relaxation dynamics of the superconducting order parameter in the presence of electric field E is described by the gauge-invariant TDGL equation 22 with the Langevin white noise: Here the order parameter relaxation time is given by , where the inverse diffusion constant γ/2, controls the time scale of the dissipation dynamical processes. The scalar potential for constant homogeneous electric field applied along the x axis is ϕ = −Ex. The thermal forces, which induce the thermodynamical fluctuations, satisfy the fluctuation-dissipation theorem The electric current density includes two components, J = J n + J s , where J n = σ n E is the current density contributed by the Ohmic normal part and J s is fluctuation supercurrent density given by s Characteristic scales and dimensionless variables. We will be measuring the distances in the units of the coherence length ξ, the time in the units of the GL " relaxation" time 22,23 τ GL = γξ 2 /2. The fluctuation strength is characterized by the parameter where D is dimensionality of the system; A is the thickness for D = 2 or the cross-section area for D = 1. The order parameter is normalized by Ψ 2 = (2αT c0 /b)ψ 2 and electric field by The TDGL Eq. (2), written in dimensionless units reads, where t ≡ T/T c0 , = + τ τ ∂ ∂ D i y  and ζ = ζ (2αT c0 ) 3/2 /b 1/2 , the white noise correlation takes a dimensionless form: Finally, the dimensionless current density j s = J s /J GL , with J GL = cH c2 ξ/2πλ 2 as the unit of the current density, is

The Self -Consistent Approximation Calculation of the I-V Curve
In what follows we will use the Hartree -Fock type self-consistent Gaussian approximation (SCGA) 25-28 used in the past to calculate fluctuation contribution to magnetization 29 , Nernst effect 25 and conductivity above T c 30 .
Dynamical Gaussian approximation. The TDGL in the presence of the Langevin white noise, Eq. (7), is nonlinear, so cannot generally be solved. Since we will need only the thermal averages of quadratic in ψ quantities, like the superfluid density and the electric current, a sufficiently simple and accurate approximation (similar in nature to the Hartree-Fock approximation in the fermionic models) is the gaussian approximation 25,28,29 . The nonlinear |ψ| 2 ψ term in the TDGL Eq. (7) is approximated by a linear one 2〈|ψ| 2 〉ψ (there are two possible contractions between ψ * , ψ in |ψ| 2 ψ, see discussion of this point in [19][20][21]: For stationary homogeneous processes considered here, the superfluid density 〈|ψ| 2 〉 is just a constant. Now it takes a form, 2 where the excitations energy gap 23 is, The solution therefore can be written via the Green's function, Then the superfluid density, using the noise correlator, Eq. (8), can be expressed via the Green's function as, and is a function of the parameter ε which is determined self consistently by Eq. (13).
Green's function for a homogeneous constant electric field. To calculate the response of the system, one needs the Green's function in the presence of electric field: The invariance with respect to the time translations is already taken into account by setting τ = τ 1 − τ 2 . Using these expressions, the superfluid density of Eq. (15) takes a form, The integrand in Eq. (17) is divergent as 1/τ when τ → 0 when D > 1. The cutoff τ cut is thus required to account for the inherent UV divergence of the Ginzburg-Landau theory and it will be addressed below.
Finally the gap equation assumes the form After (numerical) solution for the energy gap ε, we turn to calculation of the supercurrent. While the upper limit of the integration in Eq. (18) is safe (both terms in exponent are positive), the lower limit (UV) depends on dimensionality.
In the paper [ 25 ], it was shown that τ cut in time dependent Ginzburg Landau and the energy cutoff Λ in static Ginzburg Landau theory are related by where γ E is Euler constant and Λ is the energy cutoff 25,29 . Our calculation show that taking value τ cut from 0.1 to 10, the physical quantities is essentially unchanged and is taken as τ cut = 1 in what follows.
The electric current density. The dimensionless supercurrent density along the electric field direction x, defined by Eq. (9), expressed via the Green's functions is x s 2 1 2 12 Performing the integrals, one obtains, Returning to the physical units, the total electric current density reads where E GL was defined in Eq. (6) and ω is the dimensionless fluctuation stress parameter. The gap equation determining the dimensionless energy gap ε in this units is The dynamical instability point. The dynamical instability transition temperature on the phase diagram, T * , see Fig. 1, defined as a maximal temperature at which the instability appears. Mathematically is determined by vanishing of the first two derivatives,   Now the dynamical instability transition temperature T * is determined numerically taking into account the gap Eq. (18).

Comparison with Experiments and Discussion
The results are first applied to a one dimensional superconductors -metallic wires and then for several qualitatively different types of 2D superconductors (as explained above, it is very difficult to observe the instability phenomenon in purely 3D materials, although in layered high T c cuprates close to T c the fluctuations become nearly 3D and the phenomenon was observed in magnetic field 17 .

I-V curves of 1D Sn nanowires.
We start with 1D nanowires. Granular superconducting Pb and Sn nanowires of quite regular cross -section and length have been produced by electro -deposition in nanoporous membranes 10 . It is important to note that the series of experiments of ref. 10 on Pb and Sn nanowires is the only one (known to us) in which both the current and the voltage drives were employed. This allows a qualitative understanding of the important difference between the dynamical behaviour two. We focus on the voltage drive I-V curves of Sn.
The I-V curves, measured using the voltage drive at three temperatures, are shown by dotted lines in Fig. 1. The voltage drive employed clearly demonstrates the non -monotonic character below the onset T c ≈ T c0 = 3.8K slightly above the bulk temperature of Sn (T c = 3.72 K). An experiment on the same sample (see Fig. 3b in ref. 10 ) demonstrates the voltage jumps over unstable domains of the dynamical phase diagrams. The jumps are more pronounced in Pb, see Fig. 3a of ref. 10 . This is consistent with the existence of the dynamical instability and was observed in numerous experiments (see 2D examples below).
The experimental data are fitted by Eqs (22,23) for D = 1, see solid curves. The normal-state conductivity is given, σ n = 3.6⋅10 4 (Ω*m) −1 , nanowires are 50 μm long with 55 nm in diameter. Measured material parameters are: coherence length 31 ξ = 210 nm, penetration depth λ = 420 nm and the normal conductivity was obtained from the red doted line in Fig. 1. The values of fitting parameters are: the conductivity ratio k = σ n /σ GL = 0.08 and the fluctuation strength parameter ω = 0.004, the value ω estimated from Eq. (5) is ω −  10 2 . This experiment has been already attempted to discuss in the framework of TDGL equations neglecting thermal fluctuations in ref. 9 , but only a qualitative explanation of S-shaped V(I) curves was achieved. Our approach yields the theoretical curves providing a fair agreement with the experiment and demonstrate that thermal superconducting fluctuations play the prime role in forming the instability.

I-V curves of thin quasi-2D TiN films below
Tc. The overheating model [1][2][3][4] predicts that bistability occurs in V(I) at  T T c . In experiment on thin films where the superconducting transition occurs in two stages 32 , the bistability develops at  T T BKT 5 . Thus the naive thermal balance theory does not describe the switching behavior of a 2D system near T BKT and in what follows we check if the experimental data can be reasonably fit by formulas of the inherent superconducting instability derived above. Before doing that, however, we present qualitative consideration indicating that in 2D the instability emerges at ≤ T T BKT . Now we turn to the analysis of experimental data. The measurements were taken on the titanium nitride (TiN) film having the thickness d = 7 nm, the sheet resistance in the normal state = R 320  Ω (σ n = 4.3⋅10 4 (Ω cm) −1 ) and the superconducting critical temperature T c = 3.06 K. The film was formed on the Si/SiO 2 substrate by the atomic layer deposition. Temperature dependence of the resistance of this film has been examined in 33 where it was shown that the film is quasi-2D in the considered temperature range. The sample was patterned into bridges 50 μm wide and 250 μm long. Transport measurements are carried out using low-frequency ac and dc techniques in a four-probe configuration. The current drive was applied in a wide temperature range.
We identify the BKT transition to superconducting state from the analysis of V(I) curves (see Fig. 2a) by the jump in the exponent α in the ∼ α V I -dependence, α = 1 ↔ α = 3. Upon further cooling down, the bystabilily of V(I) curves with respect to applied current (see arrows if Fig. 2b) develops and we observe voltage jumps between the low-resistive and high-resistive states. The I-V dependencies calculated using Eqs (22,23) for D = 2, for different temperature are shown in Fig. 2b as solid curves. The experimental data are fitted best with T c0 = 3.065 K and k = 0.08 which is close to k = σ n /σ GL = 0.05, where we take σ GL = 2πc 2 ξ 2 /Dλ 2 (see Eq. 10) and D is the diffusion constant D = 0.7cm 2 /sec 5,33 ; ξ = 8.2 nm; λ d = 0.6 mkm. The fitting parameter is ω = 4⋅10 −4 ÷ 8⋅10 −4 (see inset in Fig. 2b) which is in agreement with ω ≈ 5⋅10 −4 estimated from Eq. (5). One further sees that in the low voltage region fitting is expectedly worse since in the nonlinear response regime the superconducting fluctuations are practically absent. The values of the parameters ω and k were determined from fitting of experimental V(I) at currents I > I c , i.e. in the regions of instability and linear response. The linear resistance at currents I > I c is practically independent on temperature. Since this linear response is determined by σ n , we assume that k does not depend on temperature as well. The parameter ω that characterizes the fluctuation strength (5) demonstrates the maximum around  T T BKT , which is natural since fluctuations are expected to enhance around the transition temperature 24 . The observed instability gets very quickly suppressed by the magnetic field (see Supplementary). A simple qualitative picture one can think of is as follows. If the magnetic field suppresses superconductivity, it should suppress the superconducting fluctuations as well.
Layered materials. Instability in the form of the voltage jumps was observed recently in FeSeTe thin film on Pb(MgNb)TiO substrate 11 . Only the current drive was used, so that the S-shaped I-V curved cannot be determined. Only the voltage jumps were observed close to T c . The thickness of FeSeTe thin films is 200 nm. The layer distance L z = 0.55 nm 34 . Here, the normal-state conductivity is taken to be σ n = 1.3⋅10 4 (Ω*cm) −1 .
The calculated I-V curves of the 2D FeSeTe thin film with different temperature are shown in Fig. 3 as solid curves. The experimental data of FeSeTe in a current driving setup from ref. 11 are fitted best by the following values of parameters: T c0 = 7K, k = σ n /σ GL = 0.07 and the fluctuation strength parameter ω = 0.01.
As well as in TiN film (Fig. 2b), the theoretical I-V curves for FeSeTe predict a re-entrant behavior for  ≈ T T 6 BKT K (at  T 6 K the low-voltage part of the V(I) curve is V ∝ I α , where α ≈ 3 at 6 K and increases with cooling down 11 ). This re-enrance of V(I) is hard to observe directly in the current driving experimental setup. Experiments show that the current driven measurements lead to the "jump" in the I-V curves and the voltage driven measurements lead to the re-entrant S-Shaped I-V curve in the superconducting nanowires at low temperature 9 .
Other layered materials. The instability in the ultra-thin granular YBa 2 Cu 3 O 7−δ nanobridges was clearly observed in a series of works in ref. 35 . Unfortunately a 2D or a 3D model cannot quantitatively describe these I-V curves since the fluctuations in this layered material and the temperature range can be described by a more complicated Lawrence -Doniach model. The generalization is possible but was not attempted in the present work. Also, the "jump" I-V curves in a current driven setup was also reported in BSCCO 14 that is clearly 2D. Unfortunately, I-V curve at the zero magnetic field was measured only at one temperature (76 K for T c = 85.2 K).

Conclusion
In this paper, I-V curve of a D-dimensional superconductor including the thermal fluctuations effects is calculated in arbitrary dimension using the dynamical self consistent gaussian approximation method. It is shown how the thermal fluctuations generate the S-shaped instability of I-V curves. The results are applied to analyse the transport data on various materials that possess sufficiently strong fluctuations in 1D or 2D. While it is found that the unstable region can exist also in 3D, the S-Shaped I-V curve in realistic materials show only in 1D superconductors.