The critical current of disordered superconductors near 0 K

An increasing current through a superconductor can result in a discontinuous increase in the differential resistance at the critical current. This critical current is typically associated either with breaking of Cooper-pairs or with the onset of collective motion of vortices. Here we measure the current–voltage characteristics of superconducting films at low temperatures and high magnetic fields. Using heat-balance considerations we demonstrate that the current–voltage characteristics are well explained by electron overheating enhanced by the thermal decoupling of the electrons from the host phonons. By solving the heat-balance equation we are able to accurately predict the critical currents in a variety of experimental conditions. The heat-balance approach is universal and applies to diverse situations from critical currents to climate change. One disadvantage of the universality of this approach is its insensitivity to the details of the system, which limits our ability to draw conclusions regarding the initial departure from equilibrium. Increasing the critical current of superconductors has been a central scientific effort, but the fundamental understanding of critical currents near 0 K is lacking. Here, Doron et al. report that in disordered superconductors the critical current near 0 K is well explained by a thermal bi-stability where electrons thermally decouple from phonons in a discontinuous manner.

O ne of the central properties superconductors is their critical current (I c ) 1-3 the maximal current (I) they are able to maintain. The sudden onset of resistance (R) at I c is usually associated with one of two mechanisms: de-pairing, which occurs when the kinetic energy of a Cooper-pair exceeds its binding energy (the superconducting gap) 4,5 , or de-pinning of vortices, when the I-induced Lorentz force acting on vortices exceeds their de-pinning force, setting them in motion 4,6,7 . Typically, in type-II superconductors under the application of magnetic field (B), de-pinning occurs at lower I rendering the depairing current a theoretical upper-bound 4,8 .
Due to the practical significance of I c the bulk of the scientific effort was centered around increasing its value at finite temperatures (T's) rather than on its fundamental, T = 0, value. In a recent publication, I c 's of superconducting amorphous indium oxide films (a:InO) have been studied at low T's and high B's near the high critical field of superconductivity, B c2 (~13 T) 9 . The authors found that I c~| B − B c2 | α , with α ≈ 1.6 that is close to the mean-field value of 3/2 indicating, as they pointed out, that I c is a result of the combined action of de-pairing and de-pinning where the increasing I initially suppresses the order parameter (by pairbreaking), helping the Lorentz force to overcome the pinning. While by using this approach they were able to suggest a resolution to the ubiquitous linear B c2 (T) as T → 0 10,11 , their theory is not yet refined enough to offer a quantitative prediction to the value of I c itself.
Our purpose in this article is to suggest that a different mechanism, which is Joule self-heating, is behind I c , this mechanism inevitably becomes more dominant as T → 0. Selfheating occurs when the power dissipated by the measurement I exceeds the rate of heat removal from the electrons. To analyze this process we model our experiment as being comprised of four independent subsystems (Fig. 1a) that are thermally coupled via lumped thermal resistors (R's): The electrons, the host a:InO phonons, the substrate phonons and the liquid helium mixture (in which our sample is immersed in our dilution refrigerator). While our system as a whole is driven out of thermal equilibrium by our measurement I, it maintains a steady-state where we assume that we can treat each subsystem as being at local equilibrium albeit at different T's represented by T el , T ph , T sub and T 0 as indicated in Fig. 1a.
Our electronic subsystem is thermally linked to its phonons viã R elÀph , mediated by electron-phonon coupling [12][13][14][15] . The a:InO  phonons are, in turn, linked to the substrate's phonons via  acoustic transfer at the interface between the different solids 16 ,  which transfer their heat to the helium mixture through a  thermal-boundary resistance at the interface known as Kapitza  resistance,R K  17-19 . Under steady state conditions the power (P) flowing across each boundaryR is equal to the Joule heating P ≡ I ⋅ V delivered to the electronic subsystem. A finite P flowing through theR's results in a T-difference between each pair of subsystems. If one of theR's is significantly larger than the others it will constitute a thermal bottleneck, impeding the cooling process, and the largest T difference will develop across it. A straightforward analysis, given in Supplementary Note 6, reveals that the thermal bottleneck is between the electrons and the phonons (R elÀph ) and henceforth we assume that all other subsystems are in equilibrium with each other (the functional form of Equation (1) is general and also applies to the otherR's illustrated in Fig. 1a, therefore we can write the terms of Equation (1) as ifR elÀph is the thermal bottleneck and not lose generality). The T-differences across R elÀph is determined by a heat-balance equation: where Ω is the sample's volume and Γ and β are parameters characterizingR elÀph 12,14,15 . It turns out that Equation (1) can lead to dramatic behavior. If a rise in T el that results from an increase in I causes a sufficiently steep increase in R(T el ), which is certainly the case in our type-II superconductor (in our experiment RðT el Þ % R 0 e ÀT 0 =T el ), then below a critical T ph value the heat-balance equation acquires two stable solutions for T el (I): a low T el solution where T el ≳ T ph and a high T el solution where T el ≫ T ph 12,20 . The jump at I c is simply a manifestation of the system switching discontinuously between the two stable T el solutions, and the sudden increase in V at I c results from V = I c R(low T el → high T el ). We present the results of a systematic study of I c in superconducting a:InO at 0.5 > T > 0.01 K and B's 12 ≥ B ≥ 9 T (where B c2 ≃ 13T), for samples of various thicknesses in both perpendicular B (B ⊥ ) and in-plane B (B || ). We demonstrate that the V-jump at I c results from the behavior expected from the heat-balance Equation (1). Furthermore, the value of I c can be accurately determined using only measurements done at I → 0 and at I ≫ I c . We also show that I c is not consistent with either de-pairing or de-pinning mechanisms, nor with their combined action.
Results I c at low T and high B. In Fig. 1b, c we depict several low-T current-voltage characteristics (I-V's) typical of our study obtained from the 280 nm sample. In Fig. 1b we plot dV/dI vs. I measured at B ⊥ = 12 T and at T = 60mK. The measured dV/dI clearly separates to two regimes, a low-resistive (LR) state at low |I|'s and a high resistive (HR) state at high |I|'s. The transition between these states occurs discontinuously at two different I c value; we define the measured I where the HR → LR transition occurs as I H!L c and the measured I where the LR → HR transition occurs as I L!H c (as marked in the Fig. 1b). Typically our measured I-V's did not show a pronounced hysteresis as I H!L c % I L!H c (see "Discussion" section) therefore throughout the manuscript we refer to both as I c . In the inset of Fig. 1b we re-plot the data of the main graph on a semi-logarithmic scale. Plotted this way it can be seen that in the LR state, even at low I's, there is a small but finite R. This is because a thin film type-II superconductor in the presence of high disorder and B reaches R = 0 only at T = 0 6,21 . In Fig. 1c we plot dV/dI vs. I on a semi-logarithmic scale measured at 20mK and at 6 B-values between B = 9.5-12 T (B c2 for this sample, plotted in Supplementary Fig. 2, was~13 T). All curves exhibit a jump of several orders of magnitudes in dV/dI at a well-defined, and B-dependent I c : Increasing B towards B c2 results in a decrease in I c . Similar B dependence of I c is observed in all of our samples.
Heat-balance analysis. We next demonstrate, using the heatbalance approach (Equation (1)), that the I-V's are well described in terms of electron self-heating [12][13][14][15][22][23][24] . Inspecting Equation (1) we see that the only unknown variable is T el , which we need to obtain independently. For that, we assume that all deviations from Ohm's law are due to an increase in T el and not from other non-linear effects (we shall review the flaws of this assumption in the discussion). Under this assumption, we convert the raw I -V's obtained from our 280 nm film at B ⊥ = 12 T at several T's plotted in Fig. 2a to effective T el and plot the results in Fig. 2b 20,23,25 (see Supplementary Note 3 for a detailed description of the heatbalance analysis).
Finally we plot, in Fig. 2c, P þ ΓΩT β ph vs. ΓΩT β el alongside the fit to Equation (1) (dashed black line), which our data follow for more than 4 decades, and from which we extract the parameters β = 5.1 and ΓΩ = 1.48 × 10 −5 W ⋅ K −β (the values of β for our samples at various B's are given in Supplementary Table 1). The systematic deviations at low P's are addressed in Supplementary Note 7.
Calculation of I c from the heat-balance analysis. Encouraged by the excellent fit of our data to the heat-balance theory, we further provide a quantitative test for its validity by using it to predict the values of I c for our superconductors. This is achieved by numerically solving Equation (1) to obtain the theoretical lower and upper bounds of the bi-stability I-interval, I min c and I max c , for our B and T range as demonstrated in Fig. 3.
In Fig. 3a we plot the measured dV/dI vs. I for the 280 nm thick sample at T ph = 60 mK and B ⊥ = 12 T. In Fig. 3b we plot both sides of Equation (1), the magenta curve is the right-hand-side (ΓΩðT I extracted from the data of (a) using the zero-bias R(T) (inset of (c)) as an electron thermometer. c Fitting the data to Equation (1). By collapsing the different isotherms of (a) such that P þ ΓΩT where the black and magenta curves intersect twice as they become tangent at T el~7 0 mK. Above I max c the low-T el stable solution and the unstable high T el solution coincide and vanish leaving the system only with the high-T el stable solution.
In Fig. 3c we plot the theoretical J min c which is the critical current density corresponding to I min c (squares) for samples with various thicknesses together with our measured J c (crosses and ||'s representing B ⊥ and B || respectively where we have used the measured I H!L c ). For all samples and B values there is a remarkable quantitative agreement between theory and experiment. We emphasize that the measured value of J c was not used in the heat-balance analysis and so our accurate prediction of J c is a good test for the validity of this theoretical framework. Note that this result shows that although the transition can theoretically occur anywhere between I min c to I max c , in practice it occurs at I min c . In the discussion section and in Supplementary Note 8 we show that this is consistent with a switching wave that propagates through our superconductors 22 .

Discussion
Similarities with insulators -The heat-balance approach we used throughout this article is a general concept that can account for thermal bi-stabilities in various systems such as superconductors 12,22,[26][27][28][29] , insulators 20,25 , and even in earth's T 30,31 . Here its use was inspired by earlier studies of the B-driven insulating phase of a:InO 32 . There, the discontinuities in the I-V's were attributed to bi-stable T el assuming thatR elÀph dominates the electrons cooling rate at low T's 20,25,33 . In Fig. 4a and b we plot V vs. I of one of our superconducting samples alongside I vs. V obtained from the B-driven insulating phase of a more disordered a:InO sample. The color-coding indicates the measurement T. We draw attention to the qualitative similarity between both measurements, and to the fact that the values of the parameters β and Γ do not vary significantly between superconducting and insulating samples (see Supplementary Table 1).
The Ohmic assumption: In our analysis, we assumed that all deviations from Ohmic transport are due to heating and other mechanisms leading to non-linearity, while present, are less effective and do not influence our main results. For example, we do not take into account intrinsic non-linearities that are known to exist in type-II superconductors at finite T and B 5,34-36 . Our analysis, therefore, fails to quantitatively account for the onset of nonlinearity at I < I c limiting its range of applicability to I ≥ I c and I ≪ I c . In Supplementary Fig. 5, we demonstrate these deviations. A similar discrepancy was also reported in the heat-balance description of the insulating phase in a:InO 20,25 . We note that self-heating also applies in the presence of intrinsically non-linear effects and a complete description of our I-V's awaits a theory that integrates both self-heating and intrinsic non-linearities.
Other mechanisms for I c : The main result of this work is that, at low T's, I c is a result of thermal bi-stability. While we clearly demonstrated this by accurately predicting I c under various conditions, it is also important to show that the other mechanisms for I c are not relevant in our experiments. We focus here on the role played by vortices at a finite B 4,6,7 and refer the reader to Supplementary Note 2 where we show that the mechanism of depairing of Cooper-pairs is unlikely.
To examine whether vortex de-pinning can be the mechanism causing our I c 's we oriented B in the sample's plane (B || ) and conducted two measurements of I c : one where B || is aligned parallel to the source-drain I (I SD ) and one where B || was at an angle of φ ≈ 45 ∘ from I SD (φ is defined in the inset of Fig. 4c). Because the coherence length of our films ξ~5 nm 37 is smaller than the film thickness vortices penetrating the plane of the sample experience a Lorentz force / I SD sinðφÞ. In Fig. 4c we plot dV/dI vs. I of the 26 nm thick sample at T = 13 mK and at B || = 11 T where the dashed black line and the continuous red line correspond to φ ≈ 45 ∘ and φ ≈ 0 ∘ , respectively. It is apparent that the entire dV/dI curves, and in particular I c , are completely independent of φ demonstrating that I c is not due to collective depinning of vortices (similar insensitivity of transport properties to φ was reported in high T c superconductors [38][39][40][41] . While these results were still interpreted in terms of vortex motion, the different theoretical models rely heavily on the large anisotropy in high T c 's. The contrast between the φ dependence of high T c 's and an amorphous MoGe alloy, which is a conventional type II superconductor, is demonstrated in ref. 41 ).
Our I c results are not different from those recently presented in 9 . These authors offered an interpretation very different from ours. They claim that I c is a result of a combination of de-pairing and de-pinning. Their main experimental evidence are that J c ∝ |B − B c2 | α with α~1.6 which is similar to the mean-field depairing value of 3/2 and that J c is comparable to the de-pairing J c (smaller by a factor of 4 according to their calculation). We do not intend to counter their claims. We do think, on the other hand, that our analysis better describes the data for three reasons: (1.) contrary to their results, the value of α is actually nonuniversal (see Supplementary Fig. 3). (2.) the de-pairing J c is actually 10-15 times larger than their measured J c and 10-400 times greater than in our measurements (see Supplementary Note 2). (3.) unlike their model the heat-balance analysis provides a good quantitative prediction to J c . In the supplementary material of ref. 9 Sacépé et al. discuss the possibility of a thermal bi-stability and provide several arguments against this interpretation. In Supplementary Note 7, we respond to these arguments.
Lack of hysteresis due to a propagating switching current:  Table 3). However, the limited hysteresis and the results displayed in Fig. 3c indicate that both I H!L c and I L!H c occur near I min c . A common mechanism driving the transition between LR and HR states in bi-stable conductors (such as our samples at I min c <I max c ) is a propagating switching wave 22 . It turns out that such a switching wave is consistent with our observation that In this scenario, we consider a bi-stable material in the LR state (similar arguments apply to the HR state) where, due to disorder, there is local nucleation of high T el domains embedded in the low T el medium. In ref. 22 it is shown that above some minimum propagating current (I p ) these hot T el domains expand and that for I < I p these domains shrink. I p is extracted using the equal area condition SðT 3 ; IÞ Effects of the contacts: In the heat-balance analysis presented above we assumed that heating is a result of the non-vanishing R of our type-II superconductor at finite B and T. A different possibles source of heating that can potentially lead to the destruction of superconductivity at high I is dissipation that originates at the contacts due to their finite R. To reduce contact R we prepared our samples with a large overlap area of 333μm on 50 μm between the source and drain Ti/au contacts and the a:InO (see "Methods"). We find this possibility unlikely for two reasons: First, in our heat-balance analysis we obtain I c relying strictly on very low bias 4-terminal resistance measurements that are independent of contact resistances and were measured at |I| ≪ I c where heating is irrelevant. The accuracy of our analysis, as displayed in Fig. 3c, and its correspondence with the 4-terminal R makes it unlikely that the effect we present is caused by heating at the contacts. We emphasize that because the electron thermal length in our samples at 10mK is L T ≈ 0.2 μm (approximated using a free electron gas approximation and a typical charge density of n = 5 × 10 20 cm −3 ), while our samples are 1mm long, the onset of the finite 4-terminal R cannot be due to electron heating from the contacts. Second, because the contact R is not typically strongly B dependent, one would equally not expect I c to be strongly B dependent, which it clearly is, see Fig. 1c. Due to these arguments we find it more likely that the heating we report above is in the bulk of our finite R type-II superconductor.
In summary, we have showed that the I c 's of superconducting a:InO films measured at low T's and high B's are well described by thermal bi-stabilities originating from a model of heat-balance between electrons and phonons (Equation (1)). Using this model we predicted quantitatively I c for samples of different thicknesses for both B ⊥ and B || .

Methods
Sample fabrication. a:InO was deposited in an Oxygen rich environment of 3 × 10 −5 Torr by e-gun evaporation of high purity In 2 O 3 pellets onto a Si/SiO 2 substrate (a boron doped silicon wafer with ρ < 5 mΩ ⋅ cm with a 580 nm thick oxide layer). The sample thickness was measured in situ during evaporation using a crystal monitor and verified later by atomic force microscopy. The contacts of the samples are Ti/Au, prepared via optical lithography prior to the In 2 O 3 evaporation.
In Fig. 5a we present an illustration of the measured samples. The samples are Hall-bar shaped where the distance between source and drain contacts is 1 mm and the width of each sample is 1/3 mm. Adjacent V contacts are located 0.8 squares apart (267 μm). Our study was performed using four such a:InO films of thicknesses 26, 57, 100, and 280 nm. Each sample was thermally annealed post deposition to a room-T resistivity (ρ) of 4 ± 0.2 mΩ ⋅ cm, which places them in the relatively low-disorder range of a:InO.
Measurement setup. The samples were measured in an Oxford instruments Kelvinox dilution refrigerator with a base T of 10 mK, equipped with a z-axis magnet. In order to apply B's in both perpendicular and in-plane orientations we mounted our samples on a probe with a rotating head. While measuring, all lines were filtered using room-T RC filters with a cutoff frequency of 200 KHz.
The transport method we use to measure the zero bias R (defined as lim I!0 V I ) is a 4-terminal method. In this measurement, the sample is probed using an input ac I (low frequency of~10Hz) and we measure the resulting ac V drop between a pair of contacts along the I path. The 4-terminal configuration is described in Fig. 5b where during a zero bias R measurement we fix I in dc ¼ 0 and set a different I ac for samples of different thicknesses while maintaining a current density of J ≈ 0.1 A ⋅ cm −2 . To measure a non-Ohmic response we use the same configuration as in Fig. 5b and measure how the differential resistance dV dI dV ac dI ac varies as a function of I I in dc .

Data availability
The data that support the findings of this study are available in Mendeley with the identifier https://doi.org/10.17632/24kvcdtkjn.  Fig. 5 Sample geometry and schematics of transport measurements. a Schematics of the sample geometry. b Schematics of the 4-terminal measurement scheme. A probing I is driven between the sample's source and drain and the resulting V drop is measured. An ac I is produced by dropping an oscillating V (from an EG&G 7265 lock-in amplifier) on a large resistor R ac . To produce a dc I we output a dc V (from a Yokogawa 7651 programmable DC source) and drop it over a large resistor R dc . I from the drain is amplified using a variable gain Femto DDPCA 300 inverting I-amplifier. The measured V drop is amplified using a home made instrumentation amplifier. The outputs of both I and V-amplifiers are measured using an EG&G 7265 lock-in amplifier (for ac) and an HP 34401A multimeter (for dc). All lines are filtered using external low pass filters (LPF).