Ultra-fast vortex motion in a direct-write Nb-C superconductor

The ultra-fast dynamics of superconducting vortices harbors rich physics generic to nonequilibrium collective systems. The phenomenon of flux-flow instability (FFI), however, prevents its exploration and sets practical limits for the use of vortices in various applications. To suppress the FFI, a superconductor should exhibit a rarely achieved combination of properties: weak volume pinning, close-to-depairing critical current, and fast heat removal from heated electrons. Here, we demonstrate experimentally ultra-fast vortex motion at velocities of 10–15 km s−1 in a directly written Nb-C superconductor with a close-to-perfect edge barrier. The spatial evolution of the FFI is described using the edge-controlled FFI model, implying a chain of FFI nucleation points along the sample edge and their development into self-organized Josephson-like junctions (vortex rivers). In addition, our results offer insights into the applicability of widely used FFI models and suggest Nb-C to be a good candidate material for fast single-photon detectors.

T he dynamics of vortices at large transport currents is of major importance for the comprehension of vortex matter under far-from-equilibrium conditions and it sets practical limits for the use of superconductors in various applications [1][2][3][4][5][6][7][8][9] . The physics of current-driven vortex matter is getting especially interesting when the vortex velocity exceeds the velocity v ≈ 3 km s −1 of other possible excitations in the system, allowing for the Cherenkov-like generation of sound 10,11 and spin 12,13 waves by moving fluxons, which opens up novel routes to excite waves in magnon spintronics 14,15 . Furthermore, there is currently great interest in the interplay of Meissner currents and magnetic flux quanta with spin waves in the rapidly developing domain of magnon fluxonics 16,17 .
The maximal current a superconductor can carry without dissipation is determined by the pair-breaking (depairing) current I dep . However, a highly resistive state in real systems is usually attained at much smaller currents due to the presence of regions in which superconductivity breaks down long before I dep is reached. Namely, in a vortex-free state, the earlier breakdown of superconductivity is due to spatial variations of the order parameter caused by structural imperfectnesses and the sample geometry 18,19 . In the vortex state, fast-moving vortices are known to lead to a quench of the low-dissipative state at I * ≪ I dep as a consequence of the flux-flow instability (FFI) associated with the escape of quasiparticles (normal electrons) from the vortex cores 20,21 . Accordingly, to achieve I c ≲ I dep and high vortex velocities v ≳ 5 km s −1 , a high structural homogeneity and fast cooling of quasiparticles (governed by the quasiparticles' energy relaxation time τ ϵ and the escape time of nonequilibrium phonons to the substrate τ esc ) are both required. However, while short τ ϵ is inherent to disordered superconducting systems 22,23 , few of them have I c ≲ I dep in conjunction with weak volume pinning needed to maintain long-range order in the fast-moving vortex lattice. Variation in the local pinning forces induced by uncorrelated disorder (volume pinning) leads to a broader distribution of v and thereby prevents the exploration of vortex matter at high velocities [24][25][26][27] .
Recently, two approaches were used to demonstrate ultra-fast vortex motion at v ≳ 5 km s −1 . In the first case, a clean Pb bridge with both, an edge barrier for vortex entry and a high demagnetization factor (so-called geometrical barrier) was studied 6 . In the used geometry there was a strongly nonuniform current distribution both across and along the bridge due to a small Pearl length 2λ 2 /d ≪ w, where d and w are the film thickness and width, respectively. A weak pinning and a short electron-phonon relaxation time τ ep in Pb 28 allowed one to diminish nonequilibrium effects and achieve the regime with ultra-fast Abrikosov-Josephson vortices 6 . In the second case, an array of ferromagnetic Co nanostripes on top of a superconducting Nb film led to a dynamic ordering of flux quanta guided by the nanostripes and allowed to achieve a narrow distribution of their velocities 29 . In both of these approaches, specially designed, locally nonuniform structures were used. At the same time, a close-to-ideal uniform system where the fast heat removal from electrons rather than the finite width of the v distribution becomes the limiting factor for ultra-fast vortex dynamics was never investigated experimentally. Theoretically, however, it was recently predicted that dirty superconductors with weak volume pinning and strong edge barrier for vortex entry should also allow for ultra-fast vortex dynamics 30 . Extremely dirty superconductors are known to have a short electron-electron inelastic scattering time τ ee which leads to a decrease of τ ep 31 . This diminishes nonequilibrium effects and may lead to an increase of the critical velocity of vortices. One of the most important requirements for the observation of an edge-controlled FFI is a spatially homogenous edge in conjunction with a weak pinning in the superconductor's volume 30 . The presence of a strong edge barrier in such superconductors leads to a current gradient near the edge where vortices enter the superconductor and where FFI is actually nucleating.
Here, we demonstrate experimentally the phenomenon of edge-barrier-controlled FFI in direct-write superconductors with a close-to-perfect edge barrier and deduce vortex velocities up to 15 km s −1 from their current-voltage curves (I-V). The investigated system is the recently synthesized Nb-C superconductor fabricated by focused ion beam induced deposition (FIBID) 32 , with a very high resistivity ρ = 572 μΩcm. This implies a large effect of the inelastic electron-electron scattering with the characteristic times τ ee ≲ τ ep which speeds up the relaxation of disequilibrium. The Nb-C microstrips have a rather low depinning current and their critical current is controlled by the edge barrier for vortex entry. In contrast to ref. 6 , in our system λ 2 /d ≫ w, which means a negligible demagnetization factor (no geometrical barrier) and a uniform current distribution across the strip at zero magnetic field. The spatial evolution of the FFI is described in terms of the edge-barrier-controlled FFI model recently developed by one of the authors 30 , implying a chain of FFI nucleation points along the sample edge and their development into selforganized Josephson-like junctions (vortex rivers) evolving to normal domains which expand along the entire sample. In addition, our results offer insights into the applicability of widely used FFI models and render Nb-C to be a good candidate material for fast single-photon detectors.

Results
System under investigations. We study the vortex dynamics in a direct-write Nb-C superconducting microstrip fabricated by FIBID 32 . The microstrip is characterized by a transition temperature of T c = 5.6 K and close-to-depairing values of the zerofield critical current I c ≈ 0.7 − 0.74I dep above 0.5T c . The dimensions of the microstrip are: thickness d = 15 nm, width w = 1 μm, and length l = 6.6 μm, see Fig. 1 for the geometry. The microstrip is characterized by the coherence length at zero temperature ξ (0) ≈ 6.5 nm, the penetration depth λ(0) ≈ 1060 nm, and the Pearl length 2λ 2 (0)/d ≈ 150 μm, that is 2λ 2 (0)/d ≫ w ≫ ξ(0). The perpendicular-to-film-plane magnetic field with induction B = μ 0 H populates the microstrip with a lattice of Abrikosov vortices. The applied dc current exerts a Lorentz force on the vortices that causes their motion with velocity v across the microstrip. The associated voltage drop V along the microstrip is recorded as a function of the applied current I in the current-biased mode. The microstrip is capped with an insulating Nb-C layer fabricated by focused electron beam induced deposition (FEBID) 32,33 . Further details on the sample fabrication and its structural properties are given in the "Methods" section.
Current-voltage characteristics. The decreasing dependence of v * (B) below about 10 mT due to the decreasing vortex density (the so-called low-field crossover in the v * (B) dependence 34 ) is beyond our consideration, as we are especially interested in the regime of very high vortex velocities.
Influence of the edge barrier on the vortex dynamics. The magnetic field dependence of the critical current at 4.20 K is presented in Fig. 3c. At smaller fields, I c (B) decreases linearly with B, while at larger fields the decrease of I c becomes nonlinear and slower. This behavior can be explained by the presence of some threshold field B stop , which demarcates the Meissner (vortex free) and the mixed states of a superconducting stripe 35 . Namely, the dependence I c (B) in the Meissner state (B < B stop ) is linear and it is described by the expression I c (B) = I c (0 T)(1 − B/2B stop ), where B stop in the Ginzburg-Landau model 36 is given by πξðTÞwÞ. Here, B s is the field value at which the surface barrier for vortex entry is suppressed at I = 0, ξ is the superconducting coherence length, and w is the microstrip width. The definition of B stop following from I c (2B stop ) = 0 is illustrated in Fig. 3c. For 10 mT ≲ B ≲ 100 mT, the dependence of the critical current is described well by the dependence I c (B) = I c (0 T)B stop /2B, and I c (B) exhibits a linear decrease at low fields. At larger fields, B ≳ 100 mT, a further crossover at B * to a slower decrease of I c (B) as B −0.5 is observed. The totality of our    negative fields, the notch locally suppresses the edge barrier and thereby facilitates the entry of (anti)vortices. This leads to a small reduction of I c (B) up to larger field magnitudes at which the role of the volume pinning increases. At positive fields, when vortices enter the microstrip from the opposite side, the notch does not affect the vortex entry and this is why I c is not affected by the presence of the notch at B ≳ 15 mT. Remarkably, when vortices enter the microstrip via the edge with the notch, I * (B) at 20 mT ≲ B ≲ 100 mT decreases by up to about 10% in comparison with I * (B) when vortices enter from the opposite side, which is in line with the calculations 39 . Importantly, due to the nonlinear upturns of the I-V curves just before the instability jump, a decrease of I * by about 10% leads to a stronger decrease of the instability velocity v * . This provides a direct evidence of the decisive role of the edge barrier on the FFI, as will be detailed next.

Discussion
We first compare the experimental results with the widely used Larkin-Ovchinnikov (LO) FFI model 21,40 with the modifications introduced by Bezuglyj and Shklovskij (BS) 41 and Doettinger et al. 42 . Although edge-barrier effects are not considered in these models 21,[40][41][42] , it is still interesting to check what quasiparticle energy relaxation time τ ϵ values, related to the instability velocity, can be deduced from fitting of the experimental data to these models.
Within the framework of the LO theory 21,40 , the microscopic mechanism of FFI is the following. When the electric field induced by vortex motion raises the quasiparticle energy above the potential barrier associated with the order parameter around the vortex core, quasiparticles leave it and the core shrinks. The shrinkage of the vortex cores leads to a reduction of the viscous drag coefficient and a further avalanche-like acceleration of the vortex, eventually quenching the low-resistive state. The original LO theory was developed in the dirty limit near T c and in neglect of heating of the superconductor. To account for quasiparticle heating due to the finite heat-removal rate of the power dissipated in the sample, the LO theory was extended by BS 41 . In the BS generalization, the latter effect was considered in the framework of the kinetic equation LO approach, which assumes a nonthermal (non-Fermi-Dirac) electron distribution function, while Joule heating was taken into account using a thermal distribution function and the electron temperature T e was determined from the heat conductance equation. In contrast to the B-independent instability velocity v * in the LO model, a v * (B) variation is expected in the BS model 41 and takes the form: where h is the heat removal coefficient. While the magnetic field dependence v * (B) nicely fits Eq. (1) at B ≳ 50 mT, a notable deviation of v * (B) toward smaller values is observed in Fig. 3a at B ≲ 50 mT. This deviation will be commented in what follows. In all, the complete set of the instability parameters deduced from Fig. 2 nicely fits the BS scaling law, see Supplementary Fig. 1. However, if one associates τ ϵ with the electron-phonon scattering time τ ep in the LO model, the deduced τ ϵ is at least one order of magnitude smaller than one could expect from τ ϵ found in similar low-T c highly disordered superconductors [43][44][45] , see Supplementary Discussion. In the LO model modified by Doettinger et al. 42[,46 , the quasiparticle energy relaxation time can be found from the following equation: In Eq. (2), the term a= ffiffiffiffiffiffiffi ffi Dτ ϵ p , where a is the intervortex distance, has been added to incorporate the necessary condition of spatial homogeneity of the nonequilibrium quasiparticle distribution between vortices at relatively small magnetic fields. The calculation results by Eq. (2) are shown by solid lines in Fig. 3a where the energy relaxation time has been varied as the only fitting parameter. The best fits are achieved with τ ϵ = 16 ps which could be considered as a more accurate estimate for the energy relaxation time in the Nb-C-FIBID superconductor. We note that with this τ ϵ estimate, the quasiparticle diffusion length l ϵ ¼ ffiffiffiffiffiffiffi ffi Dτ ϵ p ≈ 28 nm is much smaller than the intervortex distance a at all used magnetic fields and, importantly, l ϵ ≲ 2ξ(T) with 2ξ (0.75T c ) ≈ 25 nm and 2ξ(0.9T c ) ≈ 38 nm.
The edge-barrier-controlled FFI scenario 30 is different from the FFI scenario of LO and BS. Indeed, LO and BS considered a moving periodic vortex lattice in an infinite superconductor in the Wigner-Seitz approximation and hence could not take into account the collective effects related to the transformation of the vortex lattice and edge-barrier effects. In contrast, in the edge- barrier-controlled FFI model 30 a nonuniform distribution of vortices is taken into account, as well as the local Joule heating and cooling (due to the time variation of the magnitude of the superconducting order parameter |Δ|) depending on the vortex position. The edge-barrier-controlled FFI model allows for studying a "local" instability and collective effects in the vortex dynamics relying upon the solution of a heat conductance equation for the electrons and a modified time-dependent Ginzburg-Landau equation for Δ(r, t). In this model, it was shown that, in the low-resistive state, there is a temperature gradient across the width of the microstrip with maximal local temperature near the edge where vortices enter the sample 30 . The higher temperature at the edge is caused by the larger current density in the near-edge area due to the presence of the edge barrier for vortex entry and, hence, the locally larger Joule dissipation. With increase of the current, there is a series of transformations of the moving vortex lattice. In Fig. 5, we show examples of the calculated I-V curves and snapshots of |Δ|(r) for the parameters of the superconductor as in ref. 30 . Similar transformations connected with reorientations of the moving vortex lattice in the insets 1-2 in Fig. 5b were experimentally observed 47 and theoretically analyzed 48 previously.
At currents just below I * , localized areas with strongly suppressed superconductivity and closely spaced vortices appear near the hottest edge (left edge in the insets in Fig. 5). Upon reaching I * , these areas begin to grow in the direction of the opposite edge and form a highly resistive Josephson SNS-like link (vortex river) along which vortices move 3,6,30,49 . These vortices are of the Abrikosov-Josephson type, as they are moving in areas with suppressed order parameter. Due to the increasing dissipation, vortex rivers evolve into normal domains which than expand along the microstrip. In consequence of this, a jump to the highly resistive state occurs at I * . In all, the simulation results demonstrate that transformation of the moving vortex array is a collective phenomenon, which involves correlated changes in the motion of many vortices with increase of the current and, at I * , results in the appearance of Josephson-like SNS links known as vortex rivers 3,6,49 .
In the edge-controlled FFI model 30 , the current I * increases linearly with the width of the strip, while V * does not depend on w as it does in the LO and BS models. This result holds at B ≫ B stop when a is much smaller than the microstrip width w and a becomes smaller than the width of the vortex-free region near the edge of the microstrip. This means that despite the nucleation of FFI points occurs near the edge where the local temperature and the current densities are maximal, far from the edge where the current density is uniform, the vortices should move at relatively high velocities. Otherwise the FFI will not develop across the whole microstrip and one has only origins of the vortex rivers, as it can be seen from Fig. 5 in 30 at I ≲ I * . The linear scaling of I * (w) with the microstrip width w is corroborated by the experimental observation in Fig. 6a, where the I-V curves for two microstrips with the widths w = 1 μm and 500 nm are shown at T = 4.2 K and B = 50 mT.
In the edge-barrier-controlled FFI model 30 , the energy relaxation time depends not only on the electron-phonon relaxation time τ ep (as in the LO model) but also on the escape time of nonequilibrium phonons to the substrate τ esc and the ratio of the electron and phonon heat capacities, C e and C p , respectively. At T ≃ T c and for a small deviation from equilibrium one has: where τ E ≃ τ ep /4.5 is the electron-phonon relaxation time renormalized due to fast electron-electron inelastic scattering. Here, τ ep is the electron-phonon relaxation time used in the LO model. Following the arguments of ref. 42 , one can claim that the instability occurs at the velocity v *~a /τ ϵ when the intervortex distance is a ≲ ffiffiffiffiffiffiffi ffi Dτ ϵ p . This condition leads to a dependence of v * (B), which was revealed in numerical calculations 30 . One important difference between the modified LO model 42 and the edgecontrolled FFI model is that in the latter 30 , a~B −1/2 only at relatively large magnetic fields, when the intervortex distance at Ĩ I c and I~I * is almost the same despite the change in the structure of the moving vortex lattice. At relatively small magnetic fields, a in the vortex rows is smaller than ð2Φ 0 =B ffiffi ffi 3 p Þ 1=2 at I~I * and, thus, the number of vortices is smaller than follows from the simple estimate nΦ 0 = BS, see Fig. 5a. Altogether, this leads to a weaker experimental dependence v * (B) than follows from the "global" instability model with v *~B−1/242 . Qualitatively, it is this behavior which is observed in the experiment, see Fig. 3a.
The large v * values observed in our system should be attributed not only to τ E < τ ep but, also, to a small τ esc in Eq. (3). Indeed, due to the insulating Nb-C-FEBID layer on top of the microstrip, there seems to be no phonon bottleneck which could exist due to an acoustic mismatch between a thin dirty superconductor and a substrate 44 . As an estimate, for our system we deduce τ esc4 d/u ≈ 24 ps, where u~2.5 km s −1 is the mean sound velocity. This value is larger than τ ε~1 6 ps deduced from the experimental data using the modified LO model. We have to stress that numerical coefficients in the LO model are strictly valid only rather close to T c (when Δ(T) ≪ k B T c , i.e., at T ≳ 0.9T c ) and in the case when τ ee ≫ τ ep and τ esc = 0. Therefore these coefficients may be different in our dirty system with τ ϵ~τesc and at temperatures further away from T c .
Finally, we would like to note that, unfortunately, there is no analytical relation between v * and τ ϵ in the edge-barriercontrolled FFI model 30 . Accordingly, a discussion of the relation between v * and τ ϵ has to remain on a qualitative level. From Eq. (3) it follows that a change of τ E , τ esc , and C e /C p leads to a change of the relaxation time τ ϵ . To illustrate this, in Fig. 6b we present a series of calculated I-V curves at different τ esc values, while the other parameters are kept fixed. Indeed, with increasing τ esc the critical velocity v *~E* decreases, but it decreases slower than τ À1 ϵ or τ À1=2 ϵ . Qualitatively, the same tendency is found if one increases the ratio C e /C p for a given τ esc value. Specifically, with an increase of τ esc /τ E by two orders of magnitude, E/E 0 decreases by only about a factor of three. In the inset of Fig. 6b, one can also see that with the increase of τ esc , the time-averaged temperature in the center of the superconducting microstrip increases, which indicates an increased contribution of Joule dissipation to the FFI. The increased temperature also affects v * because of the temperature dependence τ E~1 /T 3 and C e /C p~1 /T 2 in the used model 30 .
We would like to outline an applications-related aspect of the superconducting properties of the studied Nb-C-FIBID microstrip. Namely, the small diffusivity D ≈ 0.49 cm 2 s −1 and the low transition temperature T c = 5.6 K suggest that Nb-C-FIBID may be a candidate material for superconducting single-photon detectors (SSPDs). We refer to Table 1 for a comparison with parameters of some typical SSPDs and to ref. 31 for a further discussion. In this regard, it should be mentioned that for about a decade SSPDs were made of meandering nanostrips with widths in the range 50-150 nm as it was empirically found that the use of wider strips leads either to the loss of the single-photon nature of the response or to a decrease of the detection efficiency 50 . This observation was in line with a "geometric-hot-spot" detection model, in which the width of the supercurrent-carrying strip should be comparable with the diameter of the normal region where the superconducting state is suppressed due to the absorption of the photon.
Recently, a "photon-generated superconducting vortex model" was proposed 31,51 . It was revealed that the efficiency of the photon detection is not determined by the geometry, as long as the initial current density is uniform and close to the critical pair-breaking current I dep . It was emphasized that even several micron wide dirty superconducting stripes should be suitable to detect single near-infrared or optical photons if their critical current I c ≳ 0.7I dep 31 . The only requirement for the width of the strip is that it should be smaller than the Pearl length Λ = 2λ 2 /d that ensures the uniformity of the supercurrent across the superconductor width. Recently, this condition was satisfied in wide and short NbN 52 and MoSi 53 bridges, whose photon response was consistent with the vortex-assisted mechanism of initial dissipation 51 . In this way, given the superconducting properties of our samples, which drastically differ from much cleaner Nb-C films prepared by pulsed laser ablation in ref. 54 , Nb-C-FIBID appears to be a good candidate for fast single-photon detection. A further enhancement of the critical current in Nb-C-FIBID can be expected for tapered current leads 52,53 which should minimize the reduction of I c in consequence of undesired current-crowding effects 19 , and additional advantages of easy on-chip 55 or on-fiber 56 integration are provided by the direct-write nanofabrication technology. Furthermore, the ability to control the thickness of individual FIBID/ FEBID layers with an accuracy better than 1 nm 57,58 should allow for the fabrication of superconductor/insulator superlattices for studying quantum interference and commensurability effects 59 as well as photonic crystals with superconducting layers 60 .
To summarize, we have experimentally demonstrated ultra-fast vortex dynamics at velocities up to 15 km s −1 in a uniform superconducting microstrip fabricated by FIBID. A stable flux flow at such high velocities is a consequence of the combined effects of a strong edge barrier against a background of weak volume pinning, close-to-depairing critical currents, and fast quasiparticles relaxation in the investigated system. The distinctive feature of the direct-write Nb-C superconductor is a close-toperfect edge barrier which orders the vortex motion at large current values and allows for the description of the spatial evolution of the FFI relying upon the edge-barrier-controlled FFI model. The observed high vortex velocities in Nb-C-FIBID make accessible studies of far-from-equilibrium superconductivity 61 and vortex matter driven by large currents, opening prospects for Cherenkov-like generation of other excitations by the fast-moving vortex lattice in ferromagnet/superconductor hybrid structures. In addition, the small electron diffusion coefficient D ≈ 0.5 cm 2 s −1 , the low superconducting transition temperature T c = 5.6 K, and high I c values exceeding 70% of the depairing current render Nb-C-FIBID to be an interesting candidate material for fast singlephoton detectors.

Methods
Sample fabrication and its structural properties. Superconducting microstrips were fabricated by FIBID in a dual-beam scanning electron microscope (FEI Nova Nanolab 600). The substrates are Si (100, p-doped)/SiO 2 (200 nm) with lithographically defined Au/Cr contacts for electrical transport measurements 62 . FIBID was done at 30 kV/10 pA, 30 nm pitch and 200 ns dwell time employing Nb (NMe 2 ) 3 (N-t-Bu) as precursor gas. The as-deposited Nb-C-FIBID microstrips have well-defined smooth edges and an rms surface roughness of <0.3 nm, as deduced  NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-16987-y ARTICLE from atomic force microscopy scans in the range 1 × 1 μm. Right after the deposition, without breaking the vacuum, the microstrips were covered with a 10nm-thick insulating Nb-C layer prepared focused by FEBID 33,63 , see Fig. 1 for the geometry. While Nb-C-FEBID structures are amorphous, Nb-C-FIBID deposits have an fcc Nb-C polycrystalline structure, with grains about 15 nm in diameter 32 . The typical elemental composition in the Nb-C-FIBID microstrips is 43% at. C, 29% at. Nb, 15% at. Ga, and 13% at. N, as inferred from energy-dispersive X-ray spectroscopy on thicker replica of the fabricated structures. Experiments were done on a series of four samples. In the manuscript, we report typical data for one microstrip. An additional reference measurement has been made for a microstrip with an artificially fabricated edge defect. The defect (notch) was milled by focused Ga ion beam at a beam voltage of 30 kV and a beam current of 10 pA 64 .
Superconducting properties of the Nb-C-FIBID microstrip. The resistive properties of the microstrip are summarized in Fig. 7. The resistivity temperature dependence ρ(T) is shown in Fig. 7a, where the ρ(T) curve exhibits a transition from weak localization 65 to superconductivity at T c = 5.6 K. Here, the transition temperature T c is determined using the 50% resistance drop criterion, as illustrated in Fig. 7b. The resistivity at 7 K is ρ 7K = 572 μΩcm and the width of the superconducting transition, defined as the temperature difference between the 10 and 90% resistivity values at the transition, amounts to ΔT c ≈ 0.6 K. Application of a magnetic field B leads to a decrease of T c and a transition broadening, and we use the same 50% resistance drop criterion to deduce the temperature dependence of the upper critical field B c2 (T) shown in Fig. 7c. Near T c , the critical field slope dB c2 =dTj T c ¼ À2:24 T K −1 corresponds, in the dirty superconductor, to the electron diffusivity D ¼ À4k B =½πeðdB c2 =dTj T c Þ % 0:49 cm 2 s −1 . The coherence length and the penetration depth at zero temperature are estimated 52 as ξð0Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi _D=Δð0Þ p ¼ 6:5 nm and λð0Þ ¼ 1:05 Á 10 À3 ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ρ 7K =T c p % 1060 nm. By employing the 100 μV voltage drop criterion, from the I-V curves, we deduce the critical currents at zero field I c (0.75T c ) = 58 μA and I c (0.9T c ) = 16 μA. We assume that the temperature dependence of the depairing current can be described by the , which is justified for dirty superconductors 52,66,67 . Here, Δ(0) is the superconducting energy gap at zero temperature, e the electron charge, and R □ the sheet resistance. With the assumed BCS ratio Δ(0) ≈ 1.76k B T c , we obtain I dep (0) ≈ 268 μA. The calculated dependence I dep (T) is compared with the experimentally measured I c (T) in Fig. 7d. We note that I c varies between 0.7I dep ≲ I c ≲ 0.74I dep in the temperature range 0.5 < t < 1, where τ = T/T c is the reduced temperature.
Time-dependent Ginzburg-Landau simulations. To study the evolution of the superconducting order parameter, we numerically solve the modified TDGL equation 31 : p ÞÞ 2 = ð1 À T e =T c Þ, A is the vector potential, φ is the electrostatic potential, D is the diffusion coefficient, σ n = 2e 2 DN(0) is the normal-state conductivity with N(0) being the single-spin density of states at the Fermi level, and j Us s and j GL s are the superconducting current densities in the Usadel and Ginzburg-Landau models: where q s = ∇ϕ − 2eA/ℏc, ϕ is a phase of Δ = |Δ|e iϕ , and j GL s ¼ πσ n jΔj 2 4e_k B T c q s . It should be noted that at T e not very close to T c the Ginzburg-Landau expression for the superconducting current is not valid quantitatively and one needs to use the Usadel expression for j Us s . In this case, one should also modify the TDGL equation since the ordinary TDGL equation leads to divj GL s ¼ 0 in the stationary case, while one needs divj Us s ¼ 0. Accordingly, by adding the term divðj Us s À j GL s Þ in the TDGL equation we provide divj Us s ¼ 0. At T e → T c the modified TDGL equation reduces to the ordinary TDGL equation and divðj Us s À j GL s Þ goes to 0. The electron and phonon temperatures, T e and T p , respectively, are found from the solution of following equations: where E 0 ¼ 4Nð0Þðk B T c Þ 2 , E 0 E s ðT e ; jΔjÞ is the change in the energy of electrons due to the transition to the superconducting state, k s is the heat conductivity in the superconducting state: k s ¼ k n 1 À 6 π 2 ðk B T e Þ 3 Z jΔj 0 ϵ 2 e ϵ=k B T e dϵ ðe ϵ=k B T e þ 1Þ 2 ! ; k n ¼ 2Dπ 2 k 2 B Nð0ÞT e =3 is the heat conductivity in the normal state, the term jE describes Joule dissipation, and τ esc is the escape time of nonequilibrium phonons to the substrate. The parameter γ is defined as γ ¼ 8π 2 5 C e ðT c Þ C p ðT c Þ , where C e (T c ) and C p (T c ) are the heat capacities of electrons and phonons at T = T c , and the characteristic time τ 0 controls the strength of the electron-phonon and phonon-electron scattering 31 . It should be noted that the electron-photon scattering time enters the TDGL equation indirectly via the electron temperature T e whose dynamics is governed by τ e-ph~τ0 in the heat conductance equation. This is rather similar to the LO approach, where τ e-ph enters the kinetic equation for the electron distribution function f(E) (in our case this is the heat conductance equation for T e ) and f(E) enters the GL equation in the LO model 20,21 .
To find the electrostatic potential φ, we also solve the current continuity equation: divðj Us s þ j n Þ ¼ 0; where j n = −σ n ∇φ is the normal current density.
Values of the parameters γ = 9 and τ 0 = 925 ns used in the calculations are estimates for NbN. Their variation only leads to quantitative changes in the I-V curves.
At the edges where vortices enter and exit the microstrip, we use the boundary conditions j n | n = j s | n = 0 and ∂T e /∂n = 0, ∂|Δ|/∂n = 0 while at the edges along the current direction T e = T, |Δ| = 0, j s | n = 0, j n | n = I/wd. The latter boundary conditions model the contact of the superconducting strip with a normal reservoir being in equilibrium. This choice provides a way "to inject" the current into the superconducting microstrip in the modeling. The modeled length of the microstrip is L = 4w.
In the considered model, the penetration length of the electric field L E is about the coherence length ξ(T), which is a consequence of τ ee ≪ τ ep . If τ ee ≳ τ ep , then L E can be considerably larger than ξ(T). In general, L E stipulates the stability of the phase slip process in 1D superconducting wires at larger currents 68 . In the case of vortex rivers (phase slip lines with vortices) it should also lead to their stability at larger currents, providing a critical velocity of Abrikosov vortices close to the velocity of Josephson vortices, which could explain the experimentally observed v * ≳ 10 km s −1 . Within the framework of the considered model, a larger L E can be modeled by a smaller numerical coefficient at the time derivative ∂Δ/∂t. This simultaneously leads to a decrease of the relaxation time of |Δ|, which also leads to an increase of v * . For instance, a fivefold decrease of this coefficient (that corresponds to an increase of L E by a factor of % ffiffi ffi 5 p ) results in a twofold increase of V * and v * and a small decrease of I * at B = 0.1B 0 . One can also see that in this case vortex rivers are well formed at I = I * and Abrikosov vortices are closer to Abrikosov-Josephson vortices because of the stronger suppression of the order parameter along the vortex river, leading to higher instability velocities.