Resistance of 2D superconducting films

We consider the problem of finite resistance $R$ in superconducting films with geometry of a strip of width $W$ near zero temperature. The resistance is generated by vortex configurations of the phase field. In the first type of process, quantum phase slip, the vortex worldline in 2+1 dimensional space-time is space-like (i.e., the superconducting phase winds in time and space). In the second type, vortex tunneling, the worldline is time-like (i.e., the phase winds in the two spatial directions) and connects opposite edges of the film. For moderately disordered samples, processes of second type favor a train of vortices, each of which tunnels only across a fraction of the sample. Optimization with respect to the number of vortices yields a tunneling distance of the order of the coherence length $\xi$, and the train of vortices becomes equivalent to a quantum phase slip. Based on this theory, we find the resistance $\ln R \sim -g W/\xi$, where $g$ is the dimensionless normal-state conductance. Incorporation of quantum fluctuations indicates a quantum phase transition to an insulating state for $g \lesssim 1$.

We consider the problem of finite resistance R in superconducting films with geometry of a strip of width W near zero temperature. The resistance is generated by vortex configurations of the phase field. In the first type of process, quantum phase slip, the vortex worldline in 2+1 dimensional spacetime is space-like (i.e., the superconducting phase winds in time and space). In the second type, vortex tunneling, the worldline is time-like (i.e., the phase winds in the two spatial directions) and connects opposite edges of the film. For moderately disordered samples, processes of second type favor a train of vortices, each of which tunnels only across a fraction of the sample. Optimization with respect to the number of vortices yields a tunneling distance of the order of the coherence length ξ, and the train of vortices becomes equivalent to a quantum phase slip. Based on this theory, we find the resistance ln R ∼ −gW/ξ, where g is the dimensionless normal-state conductance. Incorporation of quantum fluctuations indicates a quantum phase transition to an insulating state for g 1.
Introduction. Since its original discovery by H. Kamerlingh Onnes a defining feature of superconductivity is its vanishing electrical resistivity in the thermodynamic limit. Yet, experimental samples are finite and can therefore be expected to display a non-zero, albeit small, resistance R. This naturally provokes questions about the resistance of superconductors 1,2 -particularly in twodimensional (2D) films, where resistivity and resistance have the same physical dimension: What is the dependence of R on the system size and disorder strength? Is its scaling with increasing system size "dual" to the exponentially vanishing zero-temperature conductance of an insulator?
An approximate duality of this kind can be formalized 4 using the particle-vortex duality within the bosonic description of the superconductor-insulator transition (SIT), [5][6][7] when the normal state resistance is tuned to The linearly decaying asymptotes of ln R are consistent with experiment. 3 the quantum of resistance R Q ≡ h/4e 2 . Such theories are optimized for samples, where Cooper pairs on superconducting granules Anderson-delocalize at the SIT. In this paper, instead, we concentrate on the different experimental situation of homogeneous films for which the impurity induced reduction and ultimate annulment of T c (as defined by the onset of a spectral gap) is well described by fermionic theories 8,9 and the study of resistivity below T c constitutes a separate, subsequent question: At finite temperature, but infinite system size, resistance is established by vortex proliferation above a renormalized Berezinskii-Kosterlitz-Thouless temperature. 10,11 On the other hand, near zero temperature, but at finite system size, the resistance as a function of size and disorder strength is unknown and is the subject of this paper. We focus on a 2D superconducting strip of width W and length L > W , Fig. 1a.
We begin by briefly reviewing the literature on the resistance of superconductors. Josephson's equation relates the voltage V = (d∆φ/dt) /2e to the time derivative of the phase difference ∆φ = φ(x = L/2) − φ(x = −L/2). The presence of a bias current I lifts the degeneracy of quantum states with integer difference in ∆φ/2π, see Fig. 1b. Near zero temperature, quantum tunneling between adjacent minima dominates the decay of the phase difference. The voltage is then given by the effective tunneling rates 1/τ ±2π for the change of the phase by ±2π, Theoretically, the decay of metastable vacua is described by instanton-like field configurations [12][13][14][15] i.e., imaginarytime solutions of the semiclassical equations of motion which connect the adjacent minima. For the decay of the supercurrent, the general form of such solutions is not known. However, certain trial solutions given by vortex configurations are believed to be good approximations. For the present strip geometry, those are the following: (i) vortex configurations which are y-independent and swirl in the x − τ plane (τ is the imaginary time), i.e., the worldline of the core is space-like; (ii) vortices which wind in the x−y plane and tunnel across the system perpendicularly to the current (a time-like worldline). We will refer to case (i) as quantum phase slips and (ii) as vortex tunneling.
Most of the literature on the resistance of superconductors is devoted to thermally activated decay of the supercurrent [16][17][18] and to vortex physics in the presence of a magnetic field [19][20][21] or for current bias close to the critical current. We restrict ourselves to summarizing the work in the quantum regime for quasi-1D and 2D systems at small current and without external field. Supercurrent decay in homogeneous 1D wires was considered theoretically in Refs. [22][23][24][25][26][27][28]. There is a consensus that the semiclassically dominant field configuration is a dipole of phase slips at distance δτ ∝ 1/I in time direction, while the role of electromagnetic fields led to a debate 29,30 between the authors of Refs. [23] and [25]. In addition, it was shown that inhomogeneities in the wire can be crucially important. 31,32 For 2D samples, studies of voltage generation in ultraclean samples focus on the theory of vortex tunneling, [33][34][35][36] partly in terms of a very elegant mapping 37,38 to 2+1 dimensional electrodynamics. At the same time, the role of (i) finite width (ii) electromagnetic gauge potentials were disregarded in these works. Further, the influence of the Magnus force and the size of the vortex mass remained controversial (for review, see Ref. [39]). Quantum tunneling in the presence of inhomogeneities (pinning centers) was addressed in Refs. [33,[40][41][42].
Experimental evidence 43 for quantum tunneling of vortices in 2D superconductors at finite current bias, in particular in the context of dark photon counts, 44,45 additionally augmented the interest in this research field. New experimental tools, such as SQUID-on-tip microscopy, 46 allow accessing both vortex motion and energy dissipation (local heating). At the same time, the direct measurement of the extremely small resistance of large samples of 2D superconductors is experimentally challenging; most recent studies 47 demonstrate the technical subtleties and advances in the careful filtering of external radiation. A remedy is to study higher-resistance samples, e.g., not too far from the SIT or with a smaller width W . In the regime where a finite zero-temperature saturation value of R is measurable, experimental data 3 is consistent with − ln R ∝ W .
In this work, we consider moderately disordered homogeneous films with 1/τ el ∆ 2 /E F under infinitesimal current bias. We will refer to the case 1/τ el ∆ (1/τ el ∆) as clean (dirty) limit. Here, E F is the Fermi energy, τ el the elastic scattering time of electrons, and ∆ the spectral gap of the superconductor. As discussed below, for systems considered here the Magnus force is unimportant and the vortex mass is finite. We find that the combined tunneling of several vortices dominates over single-vortex events and demonstrate that, for the optimal number of involved vortices, the tunneling process has the same contribution as quantum phase slips. This allows us to determine the linear-response resistance of superconducting strips, which is exponentially small but finite. Estimating the preexponential factor and pushing the theory to the border of its applicability range, we also show that the SIT due to the antagonistic interplay of energy versus entropy (here due to quantum fluctuations) is captured at sufficiently strong disorder 1/τ el ∼ E F .
Quantum phase slips. We first consider the voltage generation by quantum phase slips. In perfectly coherent systems, the decay out of a "false" vacuum is technically encoded in "bounce solutions": 14,25 a phase slip takes the system to a lower vacuum where it dwells for a time δτ and subsequently an anti-phase-slip takes it back to the point of departure. The action of a bounce is (Φ 0 is the flux quantum) Each phase slip described by Eq. (2) is weighted by its core action S core ∼ K 1D ∼ (W/ξ)E F min(τ el , 1/∆), where the factor (W/ξ) represents the length of the space-like worldline of the core. Here, we introduced the dimensionless stiffness of 1+1 dimensional (i.e., yindependent) phase fluctuations K 1D and the coherence The action is maximal at the typical dwell time δτ typ = K 1D /(Φ 0 I), and a steepest-descent evaluation 48 of the instanton contribution to the partition sum 49 Z instanton ∼ (L/ξ) dδτ exp[−S I (δτ )] yields the decay rate 25 This non-linear current-voltage characteristics is a manifestation of a perfectly quantum-coherent bounce. However, at infinitesimal bias current, δτ typ exceeds the relaxation time τ coh which is always finite in a finite system connected to the metallic leads or in the presence of an external bath. The broadening of levels inside the core described by τ coh leads to a finite rate for quantum phase slips and anti-phase-slips 1/τ ±2π ∝ e −S ±I (τ coh ) L/(τ coh ξ). Then, Eq. (1) leads to a linear current-voltage relation, where the prefactor is determined by matching the nonlinear resistance at I K 1D /(Φ 0 τ coh ) mentioned above. We conclude this consideration with three remarks. First, we comment on τ coh which is finite in any realistic situation and may result from a variety of systemdependent origins, such as phonons, external radiation, and noise in the leads. Different sources generally imply different temperature and system-size dependence of the relaxation time which we do not study here, since the main factor in − ln R is the tunneling action, while τ coh enters in the form of a logarithmic prefactor. Second, spatial fluctuation of the worldline of the phase-slip core around the y−independent line have been disregarded up to now. These fluctuations lead to an additional preexponential factor, which we estimate at the end of the paper. Third, magnetic screening effects (see Ref. [29] and Supplement 50 ) are unimportant in the above consideration 30 so long as τ coh ∆ exp(λ M /W ), where λ M is the Pearl length (2D analog of the London penetration length), which is usually macroscopic. For longer relaxation times, the result for − ln R gets modified according Vortex dynamics. When the width of the sample exceeds the superconducting coherence length, vortices in the x − y plane of the superconducting film become welldefined topological excitations, and the supercurrent acquires an additional decay channel related to vortex tunneling across the sample. This effect therefore relies on the dynamics of vortices, described by the action: where x (τ ) = (x(τ ), y(τ )) is the time-dependent position of a given vortex. The first term represents the kinetic energy for a point particle of mass m, the term proportional to α describes the Magnus force, which is somewhat similar to the Lorentz force in a magnetic field, and the last term, proportional to η, describes dissipation. The values of these parameters drastically depend on whether the vortex core is featureless, or whether it contains a quasi-continuum of bound states. In the first case, dissipation is absent (η = 0) and the parameter α is topologically quantized (it is given by the superfluid density). 51 The mass diverges logarithmically if the superfluid is neutral, but the logarithm is cut inside the core for charged superconductors, where the mass is therefore minute, m ∼ m electron λ TF /ξ (here, λ TF is the Thomas-Fermi screening length). 52- 56 Here, we concentrate on the more realistic second case, in which the Caroli-deGennes-Matricon-type 57 bound states inside the vortex core cannot be neglected. In s-wave superconductors, the spacing ω 0 of these subgap states 58 can be estimated to be ω 0 ∼ ν −1 ξ −2 , where ν is the metallic density of states. Typically, in the clean (dirty) limit, For a moving vortex that explores different microscopic realizations of disorder, this leads to a quasi-continuum of states, which allows for energy dissipation, yielding η ∼ nω 0 τ el , where n is the density of electrons in the normal state. 40,59 At the same time, the Magnus force acquires a second topological contribution resulting from the spectral flow of bound states, which is equal in magnitude but opposite to the hydrodynamic contribution discussed above, so that α effectively vanishes. 60,61 Finally, the vortex mass is given by the total mass of particles trapped inside the vortex, m ∼ m electron ξ 2 /λ 2 F . 61 Single vortex tunneling. Let us now inspect tunneling events involving a single vortex in the system. We remind the reader of the logarithmic attraction of a vortex-antivortex pair, which physically reflects the principle of minimizing the region of costly superflow between defects of opposite winding. It is, perhaps, less known that a single vortex in a finite superconducting strip is attracted to the boundaries. Technically, this results from the implementation of no-current-outflux boundary conditions by means of mirror charges (see Fig. 2a,c,e and Supplement 50 ) and leads to a potential V (y) = 2J ln[W cos(πy/W )/ξ], that should be included in the total action of the vortex in the strip geometry. Here, J is the 2D superconducting stiffness: J ∼ E F (J ∼ E F ∆τ el ) in the clean (dirty) limit at zero temperature.
We begin by considering a 'kink solution', y kink (τ ), of the equation of motion, i.e., the unbinding of a single vortex from lower boundary and subsequent tunneling across the system, see Fig. 2 for illustration. The tunnel-ing action is S kink (W ) = S pot (W ) + S cond (W ), where (5) is the contribution from the potential barrier. Here, we have introduced the dimensionless function f (x) that satisfies f (x) ln(1/x) in the limit x 1 and f (π/2) = 0. The second contribution, S cond = E cond T , is the condensation energy E cond ∼ J times the total length of the worldline T = T (W ). Here, is the time for a vortex to nucleate. Thus, for narrow strips, S kink (W → 2ξ) J/∆ ∼ K 1D | W ∼ξ approaches the same value as the core action of a quantum phase slip, while for wide strips, both S pot (W ) and S cond (W ) are of the order of JW/v ∼ (W/λ F )E F min(τ el , 1/∆), which is (ξ/λ F ) times larger than the core action of a phase slip in Eq. (2).
A finite current bias lifts the degeneracy of the potential minima at y = ±W/2 by an additional term V I = −Φ 0 Iy/W . As mentioned, the supercurrent decay is governed by a bounce solution, 14,49 which we approximate by a kink and an antikink at temporal distance δτ : Using this Ansatz, we obtain the following δτ -dependent action, see Supplement 50 , At variance with the 1D bounce solution, Eq. (2), the attraction between kink and antikink, encoded in S int (δτ ), can not be determined exactly. It asymptotically vanishes if the vortex disappears at the opposite boundary, 62 S int (δτ ) → 0 for δτ > T and monotonically increases in the preceding regime 0 < T −δτ T , where S int (δτ ) −E cond (T − δτ ), which shows that the dominant reason for attraction is the string tension of the vortex worldline.
Formally, the instanton contribution to the partition sum can be evaluated using the steepest-descent method analogously to the bounce action (2). However, the properties of S int (δτ ) imply that the typical dwell time δτ diverges as a function of I → 0, exceeding the coherence time τ coh T . Consequently, evaluation of Eq. (1) similarly to (3) but using Eq. (6) instead of Eq. (2) yields The pre-exponential factor A vortex depends on the nature of relaxation and should be evaluated for a given microscopic model of relaxation. We conclude the discussion of a single-vortex tunneling with three comments. First, the logarithmic vortex interaction is cut beyond the Pearl length 50 λ M so that Eq. (5) is valid only for samples W < λ M . Second, the effect of dissipation on the tunneling has been disregarded here, which is valid if the tunneling time T is small compared to the dissipation time m/η ∼ 1/ω 2 0 τ el . This is equivalent to the condition W min(v F τ el , ξ)/(ω 0 τ el ). We will see in the next section that the resistance is dominated by vortices tunneling across an effective distance d opt W , which satisfies both the bounds imposed by screening and by dissipation, even for wide junctions. Third, in deriving the effective tunneling action, we disregarded the mesoscopic fluctuations of the superconducting gap (and, hence, stiffness J). 11,63 As we demonstrate in Ref.
[50] for E F τ el 1, these fluctuations do not affect the exponential factor in Eq. (7).
Multi-vortex tunneling. A simple tunneling event that involves more than one vortex is a dipole dissociation, 34 in which a quantum fluctuation creates a dipole of vortices inside the strip and subsequently the dipole constituents tunnel towards opposite edges, leading to an overall phase slip of 2π. Here, we consider a generalization of this event where n dipoles are nucleated in the strip and subsequently tunnel to either side of the strip, see Fig. 3b. We also consider a generalization of the edge unbinding event discussed above, Fig. 2, where the tunneling of a vortex from one edge to the other is assisted by n dipoles in the bulk, see Fig. 3a.
As a natural Ansatz, we consider (see Supplement 50 ) the total action S kink, is the kink action for a single vortex tunneling in a strip of size d introduced before Eq. (5), d i is the tunneling distance of the ith vortex, and N = 2n + 1 in the case of Fig. 3a and N = 2n in case of Fig. 3b. The minimization of the action prescribes that all dipoles nucleate at the same moment in time and same x-position, and that the nucleation points are furthermore equally spaced in y-direction, i.e., each vortex travels the same distance d N = W/N , so that the kink action is S kink,N = N S kink (d N ). We remark that this is an accurate tunneling action for the optimum event; yet, the Ansatz that we used is heuristic and employing it for the calculation of fluctuation determinant is not exact, although para-metrically correct.
Increasing N in S kink,N has two effects. On one hand, N S pot (d N ) decreases with the number of dipoles. On the other hand, the contribution of the string tension N E cond T (d N ) increases with N (even though T (d N ), of course, decreases). The first effect prevails, see Fig. 3c, in which the solid curve is obtained on the basis of the vortex motion at length scales large compared to ξ, and the dashed curve is the extrapolation of T (d N ) 1/∆ as d N approaches ξ (shaded region). The action is monotonically decreasing with increasing number of vortices, leading to an optimum number N opt ∼ W/ξ and d N ∼ 2ξ. The resistance produced by multi-vortex configurations, follows along the lines of the calculation for the singlevortex tunneling event (the preexponential factor A, which may depend on relaxation details, is considered separately). Clearly, as d N ∼ ξ, vortices cease being welldefined excitations, and the multi-vortex tunneling event cannot be distinguished from a quantum phase slip. It is, therefore, a reassuring consistency check that the exponent in the resistance due to N co-tunneling vortices, , is in accordance with the exponent of the resistance due to quantum phase slips, Eq. (3).
Conclusion. In summary, we have investigated the resistance of superconducting strips with moderate disorder 1/τ el ∆ 2 /E F near zero temperature. The obtained resistance R is dominated by cotunneling of N ∼ W/ξ > 1 vortices, a process that is equivalent to a quasi-1D quantum phase slip. This yields ln R ∼ −(W/ξ)E F min(1/∆, τ el ), see Eqs. (3) and (8), which qualitatively agrees with the low-T saturation values of resistance presented in Fig. 1 of Ref. [3]. Our result for exponential dependence of resistance on W is reminiscent of exponential suppression of conductance with the system size on the insulating side of SIT.
While we concentrated on the exponent of the superconducting decay rate, we conclude with a discussion of the preexponential factor. Its W -dependence stems from quantum fluctuations of the tunneling trajectory, only, if the mechanism of dissipation is local. Using the above-mentioned Ansatz for multi-vortex events, we estimate that each of the W/ξ vortices contributes to A in Eq. (8) a fluctuation determinant of the order of the inverse kink action, 50 leading to a preexponential factor [S kink (2ξ)] −W/ξ . (A similar factor is expected to arise from spatial fluctuations about the straight y-independent phase-slip world-line.) By reexponentiation, this quantum-fluctuation-induced prefactor can be viewed as an entropic contribution to the tunneling action, leading to where, in the dirty limit, S kink (2ξ) ∼ E F τ el ≡ g is given by the normal-state conductance, see Fig. 1c. Pushing our theory to the limit of its applicability, we observe a sign change of the action at g ∼ 1, reminiscent of SIT. We

S-II. VORTEX INTERACTION PARAMETERS AND SCREENING
For the convenience of the reader, this section summarizes basic properties of vortex interactions in the finite strip.

A. Long-wavelength action
We consider the following effective long-wavelength action of the superconducting plate, see Fig. 1

S =
supplemented with the Maxwell action of electromagnetic fields (A 0 , A) living in (3+1) dimensional space-time. As the phase field φ may contain singularities it is important to keep in mind that ∇φ is in general not a gradient field. The notation used in Eq. (S1) is suggestive and mathematically correct only away from the singularities. In our calculations, we use v 0 c and additionally assume the following realistic hierarchy of length scales L {W, λ M } ξ {λ TF , λ 0 }. Note that, in typical experimental samples, the Pearl-London length, which is about four to five orders of magnitude larger than λ 0 1nm, can be smaller, larger or comparable to the system's width.
Since we are considering a charged superconductor, the smooth part of the Goldstone mode φ can be reabsorbed in a redefinition of the electromagnetic fields (A 0 , A) ("Anderson-Higgs-mechanism"). However, singular configurations of φ, such as vortices, and an externally applied current bias can not be "eaten" by a gauge transformation of the vector potential. Differentiation of the action (S1) with respect to the gauge potentials (A 0 , A) implies the 2D charge and current densities, ρ and j tot = j + j 0êx , respectively, where and is the externally injected, fixed background current density. Here, we have separated the phase into two contributions where the time independent field φ I enters j 0 and ψ represents the fluctuating part of the phase field. The finite geometry implies the boundary conditions We will be interested in the limit of linear response in I.

B. Vortex configurations.
As mentioned in the introduction, the finite voltage in superconductors relies on quantum phase slips and tunneling of vortex excitations. Here, we provide further details leading to the technical definition of the problem. Generally, in a 2D system characterized by coordinates (x 1 , x 2 ), a vortex of winding n w at (x 1,0 , x 2,0 ) is defined by the singular field configuration ψ nw,x1,0,x2,0 (x 1 , x 2 ) = n w arctan This mapping represents the principal branch of the phase field. In our convention, the branch cut is chosen on the line x 1 = x 1,0 , x 2 < x 2,0 . In the limit λ M W vortex field lines in a charged superconductor are the same as in a neutral superfluid. A vortex field complying with the no out flux boundary conditions (S4), can be constructed by an infinite series of mirror vortices and has the form, 18 (see Sec. S-III B) This configuration is shown in Fig. 2.

C. Vortex interactions
The definining property readily leads to the single vortex potential presented in the main text, i.e. a logarithmic attraction for a vortex -antivortex pair which is at distance x . The trigonometric function in the potential of a single vortex in a finite strip can be derived from Eq. (S6) or as follows: Consider a sequence of vortices and anti-vortices in an infinite plane at positions y (i) (S9) For d = W/2 and y = y 0 + d/4 this sequence corresponds to the sequence of mirror charges of a single vortex in a strip of width W and at position y 0 . The total potential energy is (S10) This obviously diverges linearly with the number of dipoles. The energy per dipole is The constant term is y and d independent and thus physically irrelevant, we therefore do not discuss possible regularization schemes of this logarithmically divergent term. Clearly, this reproduces the potential for a single vortex in a finite strip as given in the main text, We also use Eq. (S12) for the calculation of multivortex tunneling actions.

S-III. SCREENING OF VORTICES VIA LONDON ELECTROMAGNETISM
We here briefly summarize the effect of screening of vortices and phase slips and determine that none of these effects is important in the parameter regime of interest. We largely follow 54,55 and omit retardation effects in view of the small parameter v 0 /c.
Since the action of electromagnetic fields is Gaussian, integration of electromagnetic fields is equivalent to the saddle point treatment dictated by the Maxwell equations (we use Lorenz gauge), We introduced the notationδ(y) = θ(y + W/2)θ(W/2 − y)/W , where θ(x) is the Heaviside function. Since the background current j 0 is kept fixed, it should not be further screened and does not enter these Maxwell's equations. Integration of electromagnetic fields leads to the effective action S = S 0 + S I of the matter field 54,55 The electromagnetic fields entering Eq. (S14a) are solutions to Eq. (S13) for a given field configuration ψ. The regular part of ψ can be removed by partial integration combined with the continuity equation (gauge invariance).

A. Screening of the phase-slip
We consider a phase slip in the phase field ψ at spacetime (x 0 , τ 0 ). The gradient field is (τ = v 0 τ ) with x T = (x,τ ). In Fourier space we find where q = (q x , qτ ). We 4D Fourier transform Eq. (S13) usingδ(y) = δ(y) and introduceÃ 1D (q ) = A(q x , y = 0, z = 0, qτ ), use the integral with r x ∼ ξ the extension of the asymmetric vortex core in x−direction to obtain Here, γ M,E = W/λ M,E . This leads to We observe that screening is important for − ln |q x r x | > γ −1 M,E and thus only for large systems L > ξe γ −1 M . Anticipating, that the saddle point configuration of the dipole occurs at x + = x − we obtain The symbol means equality in the limit |∆τ |/r x → ∞. We repeat that magnetic screening should only be kept provided γ −1 M = λ M /W < ln[L/ξ]. So long as γ M ln(δτ v 0 /ξ) 1, i.e. δτ ∆ e λ M /W can be omitted, while electric screening effects can always be omitted. This concludes the derivation underlying our statement of the main text: Magnetic screening effects of quantum phase slips are negligible in the parameter regime of W/ξ, and lead to the double logarithm quoted in the main text.

B. Screening of a vortex
We first discuss the screening of a static vortex configuration in the bulk of the system Eq. (S13) become Fourier transformation and the definition of A 2D (q ) = A(q , 0) leads to and thus to In the case of a vortex, ∇ µ ψ(q ) = µν 2πiq ν /q 2 , we can write after inverse Fourier transformation Here, γ EM is the Euler-Mascheroni constant, H 0 the Struve function of order zero, Y 0 the zeroth Bessel function of the second kind and the symbol means asymptotic equality for both a 1 and a 1. Since only the derivative of F enters into the action, we omit the second term, which has vanishing derivative in both limits.
It is also possible to calculate the resummation of image charges of a vortex at (0, y 0 ) such that j satisfies the boundary conditions. We define F (M ) (x, y) via Here, the summation has been evaluated by a contour integral. We introduced the notationx = x/πW,ȳ = y/πW,λ = λ M /πW and sinh(x)π/2W cosh(x) + cos(ȳ +ȳ 0 ) .
Note that g(x)| y0=±W/2 = 0, and thus j y | y0=±W/2 = 0, in accordance with the boundary conditions. We define to obtain In the screeningless limit λ M /W → ∞ this can be regarded as the derivation of Eq. (S6).

S-IV. VORTEX TUNNELING
In this section we present details on vortex tunneling.
The tunneling time in the inner part of a wide strip (W ξ) can be calculated as Note that there is an additional contribution of order 1/∆ which corresponds to the time to nucleate a vortex and thus is the lower bound for T entering S cond = JT int he main text. This concludes the derivation of the results presented around Eq. (6) of the main text.

Asymptotic evaluation of auxiliary integral
We here present the asymptotic evaluation of the auxiliary integral I α . Note that the sum I 0 (ξ) + I 0 (ξ) is plotted in Fig. 3 c) of the main text. where Here, we twice integrated by parts in order to reduce I α (x) to a simpler integral, I 1,α (x), and an integral which is small and determined by u ∼ 1, I 2,α (x), which we evaluated approximately in the limit x 1. Here erfi is the imaginary error function. Expanding the result up to next to leading order in 1/ α 2 − ln(x), we obtain . (S36)

B. Bounce solution and variational action
We now consider edge unbinding for the potential, Eq. (S12) supplemented with a tilt, V I = −βy, where β = Φ 0 I/W at small β.
As an Ansatz, we consider two kinks located at positions ±δτ /2 (S37) For δτ > T , the two kink solutions do not overlap and the action is However, if 0 < (T −δτ )/2 ≡ τ f T the overlap leads to non-linear interaction between the kinks. Specifically, the β = 0 part of the potential contribution to the action is Here, we used the τ → −τ symmetry of the Ansatz and we have expanded the action in small δy kink (τ ) = y kink (τ − δτ /2) + (W/2 − ξ). The second line vanishes, because y kink solves the equation of motions. This leads to Here, we have used the approximate solution Eq. (S34) near the turning point. In total this leads to a δτ dependent bounce action . S1: Graphical representation of the variational Ansatz (S42) for the action in case a) and n = 3. Each arrow contributes s(d), Eq. (S42c), to the total action, where d is the distance in real space between the points of nucleation and annihilation of the vortices. The major approximation is that the motion along the arrows is linear with constant speed v(d).
Here, the second last line is the correction to the condensation energy paid at during the time of the bounce T + δτ , which is shorter than the time of two kinks. This concludes the derivation of Eq. (6) of the main text.

S-V. MULTIVORTEX CONFIGURATIONS
In this section we present details on the co-tunneling events of multiple vortices.

A. Variational Ansatz and solution
We consider processes which involve n dipoles in the bulk of the system. We estimate the contribution of these processes by means of the following variational Ansatz for the action, see also Fig. S1. For the edge unbinding, i.e. processes of type Fig. 3 a) we write and for the case of Fig. 3 (S42b) and we remind the reader that for d ξ, (the general expression is presented in Fig. 3  c). Following the graphical representation of Fig 3 and Fig. S1, we label initial positions X i = (X i , Y i ) and initial times τ i with indices i ordered from top to bottom, i.e. Y i > Y i+1 . The positions X i,i+1 = (X i,i+1 , Y i,i+1 ) and times τ i,i+1 represent the points in space time, at which recombined vortex pairs annihilate. The tunneling event is dominated by contributions which fulfill Y i,i+1 > Y i+1,i+2 . Furthermore, for case a) we define X n+1 = (X n+1 , −W/2) and for case b) X n,n+1 = (X n,n+1 , −W/2). In both cases a) and b) X 0,1 = (X 0,1 , W/2).
Of course, processes with vortices and antivortices interchanged or processes of the type a) with a net vortex transfer from top to bottom are also present. Their contribution equals the contribution of one of the processes discussed here and are thus not discussed separately.
The calculation involves two assumptions. The first is the approximation of linear motion with constant velocity v(d) in an interval of distance d (e.g. d = |X 1 − X 1,2 |). In fact, writing T = d/v(d) leads to v(d) which is nearly consistent with the maximal vortex speed as imposed by energy conservation (up to a factor of 2). The second approximation is that each space time position (X i,i+1 , τ i,i+1 ) depends only on the two starting points of adjacent vortices, i.e. on (X i , τ i ) and (X i+1 , τ i+1 ). This second approximation has no influence on the saddle point action, but it does affect the fluctuation determinant. The parameter regime of validity of these approximations shall be discussed below.

Saddle-point configurations
We vary the Ansatz (S42) with respect to the starting times and positions of the vortices. It is convenient to use a lighter notation S kink (d) = s(d). This yields (here Recall that in case a) [case b)] Y 01 = W/2 and Y n+1 = −W/2 [Y 01 = W/2 and Y n,n+1 = −W/2] are fixed and thus not variational parameters. By consequence, ∂ X 1 Y 01 ≡ 0 in both cases, and additionally in case b) ∂ X n Y n,n+1 ≡ 0 while in case a) the terms proportional to ∂ Yn+1 X n,n+1 should be omitted.
Clearly, a solution of these equations is given by (X k−1,k − X k ) = (X k−1 − X k−1,k ) ∀k, X i = X j ∀i, j, and equal distances in y-direction In case a) d N = W/(2n + 1), and in case b) d N = W/(2n). By consequence of the linear motion all starting times are the same, τ i = τ j ∀i, j. This is the physical solution we perturb about.

Determination of space time positions of vortex annihilation
For the derivation of the fluctuation determinant, it will be necessary to determine the dependence X i,i+1 (X i , τ i ; X i+1 , τ i+1 ). Here we present this dependence perturbing weakly about the saddle-point solution.
The Ansatz of linear motion implies for the motion of vortices on downward oriented arrows of Fig. S1 x (S45a) while the motion on upwards oriented arrows is We determine the moment of annihilation τ k,k+1 by equating the two expressions. This leads to τ k,k+1 = ∆X k + v(∆X k )τ k + ∆X k+1 v(∆X k+1 )τ k+1 v(∆X k ) + v(∆X k+1 ) , (S46) where ∆X k = |X k − X k,k+1 | and ∆X k+1 = |X k,k+1 − X k+1 |. Inserting this solution into the Ansatz for x (τ ) in the quasi classical equations of motion yields

Fluctuation determinant (Gaussian approximation)
We here present the calculation of the fluctuation determinant. Thus, space time fluctuations of initial positions are treated at Gaussian level.
(S48d) an analogous result for the ∂ X n,ν (X n,n+1 ) µ | SP , ∂ τn (X n,n+1 ) µ | SP in the case b). For case a), we remind the reader that Y n+1 is not a variational parameter.
By means of the previous expressions one can determine the correction to the action determining Gaussian fluctuations around the saddle point solution. Direct calculation of second derivatives of the action Eq. (S42) show that fluctuations in x, y and τ directions do not and are a consequence of translation invariance in x-and τ -direction, respectively. Following the standard procedure, zero mode integrals have to be performed exactly and to be kept out of the fluctuation determinant.
The fluctuation determinant is incorporated (up to a constant) into the entropic contribution to the action δS H which in case a) is In the concluding section of the main text, we present a simplified estimate, which matches the large n limit of the calculations above, i.e.
δS H ∼ n ln[s (d N )(2ξ) 2 ] = n ln ∂ 2 ∂(d N /2ξ) 2 s(d N ) (S51) We can use that s(d N ) is actually a function of (d N /2ξ), so that modulo unimportant factors of order unity where d N ∼ 2ξ at the optimum.

S-VI. EFFECT OF MESOSCOPIC FLUCTUATIONS ON VORTEX TUNNELING
In this section, we estimate the effect of randomness in the superconducting stiffness J (corresponding to the disorder-induced mesoscopic fluctuations of the gap) on the resistance generated vortex tunneling events. The effect of gap fluctuations can be modeled by a random stiffness where ρ ∼ 1/[E F τ el (E F τ el − C)] and C ∼ 1 is small for moderate disorder. 63 The tunneling action for a single-vortex tunneling in the presence of such randomness is approximately given by Taylor expansion S = S kink + W/2−ξ −W/2+ξ dy 2mV (y)ρ(x )/2, and thus the "local resistance" is given by where R 0 is the result without fluctuations. Then, it readily follows that R(x) is white-noise, log-normal distributed, dy V (y).

(S56)
The total resistance is given by the integral of R(x) over x, i.e., by the average local resistance. Both typical and average resistance behave as ln R ∝ −W .
In multi-vortex tunneling events, the tunneling distances d i depend on the disorder configuration and are