Our experimental findings raise many fundamental questions: What are the mechanisms of vortex confinement in the channels? Why does the branching instability occur and what controls the number of vortex stems? What are the mechanisms which determine the terminal velocities of vortices? How does the structure of a vortex evolve at the superfast velocities observed in our experiments? Given the lack of theory to describe vortices under such extreme conditions, we only limit ourselves to a qualitative discussion of essential effects which may help understand our SOT observations.

The stationary pattern of vortex channels shown in Fig. 2 seems counterintuitive, since vortices repel each other and should therefore disperse over the film. Moreover, each stem grows into a tree through a series of subsequent bifurcations but the branches of different trees do not merge. In order to keep vortices within each channel a mechanism for dynamic alignment of fast moving vortices must be present. One such mechanism is that a rapidly moving vortex leaves behind a wake of reduced order parameter which attracts the following vortex. As a result, a confined chain of vortices in a self-induced channel of reduced superfluid density can be formed (Supplementary Note 1), as it was observed previously in numerical TDGL simulations15, 16 at high currents, J ~ J d , and T ≈ T c . As discussed below, this mechanism apparently becomes dominant at velocities substantially higher than those accessible in our experiment.

Vortex alignment may also result from a weak quasiparticle overheating: the power P = ηv 2 generated by each vortex produces a channel of enhanced temperature along the moving vortex chain. The resulting bell-shaped temperature distribution T(y) across the channel causes a restoring force f r = −s *dT/dy which stabilizes the vortex chain against buckling distortions (Supplementary Note 4), where s *(T) is the transport entropy that defines thermoelectric effects in superconductors7. The thermal confinement resulting in long-range alignment of vortices is particularly effective at large vortex spacing a ≫ ξ relevant to a ≃ 1–2 μm. The observed spacing along the stem varies from 0.6 to 2 μm (Fig. 4a) much larger than ξ = 46 nm, which suggests that the thermal confinement is dominant. As shown in Supplementary Note 4, this mechanism can be effective, even at weak overheating, without causing any thermo-magnetic instabilities.

Once a confined channel is formed, why would it bifurcate? A chain of aligned vortices can become unstable and buckle as repelling vortices get closer to each other. In a thin film, the long-range surface currents produced by vortices result in the repulsion force \({f_{\rm{m}}} = \phi _0^2{\rm{/}}{\mu _0}\pi {a^2}\) between vortices spaced by \(a \gg 2{\lambda ^2}{\rm{/}}d \simeq 250\) nm8. As the vortex spacing drops below a critical value a c , the net repulsion force pushing a vortex displaced by u perpendicular to the channel, \({f_ \bot } \sim u\phi _0^2{\rm{/}}{\mu _0}{a^3}\), exceeds the restoring force f r = −ku, leading to chain bifurcation, where the spring constant k(a) is determined by the confinement mechanism. The thermal confinement leads to the critical spacing that decreases with the voltage V across the channel as a c ∝ V −1/2 (Supplementary Note 4), in qualitative agreement with the experimental data shown in Fig. 4d. As described in Supplementary Note 5, transverse displacements of fast vortices at J ≫ J c can be enhanced by the weak effect of disorder which may cause premature bifurcation of vortex channels.

To gain an insight into the structure and the viscous dynamics and channeling of fast vortices, we performed numerical simulations of TDGL equations for the bridge geometry of Fig. 1 and the material parameters of our Pb films (see Methods section). We limit ourselves to a minimal model of superfast vortices driven by strong current densities J ≫ J c for which disorder and heating was neglected. The simulated Cooper pair density Δ 2(x, y) shown in Fig. 5a reproduces the main features of the SOT images at I ≲ I c , namely vortices displaced to the right edge and a pronounced vortex-free region along the left edge. Notice that in the absence of disorder, the stationary vortices in Fig. 5a form an ordered structure within a smooth confining potential of the geometrical barrier, in contrast to Fig. 2f which shows a disordered vortex configuration determined by pinning in the right-hand side of the sample where J < J c . At I = I c the calculated current density J(x) shown in Fig. 5b reaches the depairing limit J d at the left edge of the constriction and vanishes at the opposite edge.

Fig. 5 TDGL simulations of stationary and fast moving vortices at the experimentally accessible velocities. a Calculated Cooper-pair density Δ 2(x, y) of a stationary vortex configuration at applied current density and magnetic field corresponding to the experimental conditions in Fig. 2f. b Corresponding distribution of the supercurrent density |J(x, y)/J d | in the sample showing edge currents in the constriction reaching J d at the verge of vortex penetration. The black arrows point to the local direction of the current. c Time-average of the Cooper-pair density over 5 × 104τ GL at I = 1.05I c , revealing branching vortex trajectories coexisting with adjacent stationary vortices. d Snapshot of moving vortices in c with arrows denoting the relative displacement of each vortex following an entry of a new vortex into the sample. e,f Same as c,d but at highest applied current before an additional stem is formed. g Experimental vortex velocity along the stem for B a = 2.7 mT and indicated applied currents with the TDGL data from d and f in normalized units (scaled to v GL = ξ/τ GL ). The animation of the vortex flow dynamics corresponding to e,f is presented in Supplementary Videos 5 and 6. The scale bar in a is 20 ξ Full size image

At I > I c vortices start penetrating through the left edge and move along a network of preferable paths forming a branching tree with an overall shape determined by the bridge geometry. The vortex chains are curved on larger scales (Fig. 5e) due to the lensing effect of the current distribution in the constriction, which tends to orient the vortex chains perpendicular to the local current J(x, y). The calculated vortex flow pattern is similar to the SOT image in Fig. 2j and also exhibits the coexistence of moving and stationary vortices as observed in Fig. 2 where bulk pinning further hampers the motion of remote vortices. The Copper pair density averaged over the simulated period of time shows a non-uniform distribution along the vortex channels (Fig. 5c, e) with distinct bright spots of reduced Cooper pair density indicating that the vortex velocity varies non-monotonically along the channels. The bright spots describe the regions where the vortices slow down or even stop momentarily, giving rise to vortex crowding (Supplementary Movies 5–8). Similar features of B z (x, y) along the channels are observed in Fig. 2i–p. Such vortex “traffic jams” can be understood as follows. A vortex penetrating from the left edge slows down as it moves along the channel since the driving current J(x, y) decreases across the strip. The subsequent penetrating vortices move along the same trajectory, causing jamming in the regions where vortices slow down. The resulting mutual repulsion of vortices either pushes them further along the channel, where vortices speed up due to attraction to the right edge of the strip, or causes bifurcation of the channel into branches.

The fact that there is a single vortex entry point and several exit points necessitates that the penetration frequency per stem is higher than the exit frequency per channel. Figure 5d, f shows snapshots of vortex motion at different applied currents, with arrows proportional to the instantaneous velocities of vortices right after penetration of a new vortex. We find that the periodically-entering vortices take alternating routes at the bifurcation points due to interactions with other vortices, which slow down further after the bifurcation (Supplementary Movies 5–8).

Figure 5g presents the calculated vortex velocity v(x) along the stems in Fig. 5d, f juxtaposed with the experimental data. The decreasing vortex velocity follows the drop in J(x) from the left edge of the constriction, v(x) = ϕ 0 J(x)/η, showing no significant velocity dependence of η(v) at these values of v. This allows us to extrapolate the experimental v(x) to x ≃ ξ at the entrance point, where the velocity is maximal, \(v\left( \xi \right) \simeq {J_{\rm{d}}}{{\rm{\phi }}_0}{\rm{/}}\eta = 24\) km/s, assuming a constant η(v).

Using the current density J(x) obtained from the simulations and v(x) extracted from our experimental data, we obtain the vortex viscous drag coefficient η = 2.6 × 10−8 kgm−1 s−1. This value of η is of the order of \({\eta _0} \simeq \phi _0^2{\rm{/}}2\pi {\xi ^2}{\rho _{\rm{n}}} = {10^{ - 8}}\) kgm−1 s−1 of the Bardeen-Stephen model, which indicates no excessive changes in the structure of the Abrikosov vortex core even at velocities of the order of 10 km/s. This conclusion is corroborated by our TDGL simulations in Fig. 5 which reproduce the channel bifurcations due to vortex repulsion and the variation of B z (x, y) along the channels due to disorder and interactions induced variations in vortex velocity. The totality of our SOT and TDGL results indicate that vortices maintain their integrity as stable topological defects even at the observed extreme velocities for which the magnetic field of a moving vortex does not deviate substantially from that of a stationary Abrikosov vortex. In particular, we have observed no evidence of the transition of Abrikosov vortices into Josephson-like phase slip lines48, 49 extending across the bridge.

Now we turn to the numerical study of even faster vortices, beyond our experimentally accessible range of parameters, for which a significant change in the internal vortex structure is expected. For instance, nonequilibrium effects can give rise to a velocity dependence of η(v) and to the Larkin–Ovchinnikov (LO) instability caused by diffusion of quasiparticles from the vortex core50. The LO instability results in jumps in the vortex velocity above \(J >{J_{{\rm{LO}}}} \simeq {\eta _0}{v_0}{\rm{/}}2{\phi _0}\) for which the force balance ϕ 0 J = η(v)v at v > v 0 is not satisfied because \(\eta \left( v \right) = {\eta _0}{\rm{/}}\left( {1 + {v^2}{\rm{/}}v_0^2} \right)\) decreases with v 50. The LO or overheating instabilities51, 52, have been observed on various superconductors with v 0 ranging from 1 to 10 km/s17,18,19.

Our TDGL calculations at twice higher current and field as compared to those shown in Fig. 5, reveal three different types of vortices described in Fig. 6. Far from the constriction region, J is lower and the moving vortices (red dot in Fig. 6a) have a regular, nearly isotropic shape with no wake of reduced order parameter. Closer to the constriction, a chain of vortices (marked by a black dot in Fig. 6a) is confined in a channel of reduced order parameter. These faster-moving vortices are slipstreaming one another because their velocity v exceeds a/τ Δ , where \({\tau _{\rm{\Delta }}} = \pi \hbar \sqrt {1 + 4{\tau ^2}_{\rm in}{{\rm{\Delta }}^2}} {\rm{/}}\hbar^{2}{/}8{k_{\rm{B}}}\left( {{T_{\rm{c}}} - T} \right)\) is a recovery time of the superconducting order parameter in the wake of the moving vortex and τ in is an electron-phonon inelastic scattering time. Our TDGL simulations show that these vortices, moving in channels, have elongated cores along the direction of motion, and their drag coefficient can be approximated by the LO dependence \({\eta _{{\rm{LO}}}}\left( v \right) = {\eta _0}{\rm{/}}\left( {1 + {v^2}{\rm{/}}v_0^2} \right) + {\eta _{\rm{i}}}\) with η i ≈ 0.25η 0 and \({v_0} \approx \xi {\rm{/}}{\tau _{\rm{\Delta }}} \approx 20\) km/s for our sample parameters (see Methods section). These anisotropic slipstreamed vortices can undergo a kinematic transition to conventional vortices upon stem bifurcation which leads to additional vortex slowdown, as marked by a blue dot in Fig. 6a.

Fig. 6 Different morphologies of ultra-fast vortices at velocities significantly higher than in our experiment. a,b A snapshot a and time-averaged Cooper-pair density Δ 2(x, y) b as in Fig. 5, but for twice higher applied field and twice the current. Three vortex phases are found with distinctly different core structure, level of quasiparticle tailgating, velocities and resulting kinematic trajectories (see text and Supplementary Movies 7 and 8), namely the extremely fast Abrikosov-Josephson vortices (marked by green dot), the ultrafast slipstreamed vortices (black dot), and conventional Abrikosov moving vortices (red dot). c Spatial profiles of vortex velocities v(x) (scaled to v GL = ξ/τ GL , see Methods) for the three main vortex phases, and for one detected branch of vortices going through an in-motion transition (dynamic transition from slipstreamed vortices to conventional Abrikosov vortices, identified by a blue dot in a). The scale bar in a is 20 ξ Full size image

The most radical change in the structure of moving vortices occurs in the narrowest part of the constriction, where J is maximal. Here a channel with a significant reduction of the mean superfluid density appears, in which ultrafast vortices (green dot in Fig. 6a) are moving with velocities that are 3–5 times higher than the speed of slipstreamed vortices, as shown in Fig. 6c. The ultrafast vortices in the central channel can be regarded as Josephson or mixed Abrikosov-Josephson vortices similar to vortices at grain boundaries53, 54, in high-critical-current planar junctions55, or S/S′/S weak links56. The TDGL results shown in Fig. 6c suggest that Josephson-like vortices in these channels can move with velocities as high as ~100 km/s, because the viscous drag coefficient η(v) in the channel is reduced to just a few percent of η 0 due to strongly elongated and overlapping vortex cores. Spatial modulation of the order parameter between these vortices is rather weak and effectively the channel behaves as a self-induced Josephson junction, which appears without materials weak links. Similar flux channels in thin films were previously interpreted in terms of phase slip lines48, 49. In the case of strong suppression of the order parameter and weak repulsion of Josephson vortices which extend over lengths exceeding Λ =2λ2/d, the channel does not bifurcate as shown in Fig. 6a and the magnetic field along the channel is nearly uniform. This feature of Josephson-like vortices appears inconsistent with our SOT observations of vortex channels which always bifurcate and show noticeable variations of B z (x, y) along the channels. The SOT results thus indicate an essential effect of intervortex repulsions and weak suppression of the order parameter along the channel, consistent with the dynamics of Abrikosov vortices shown in Fig. 5.

Another interesting SOT observation shown in Fig. 2n–p is the nucleation of additional stems of vortices as current increases. This effect can be understood as follows. The first stem appears at I = I c as the local current density J s at the edge of the constriction reaches J d . As I increases above I c , vortices start penetrating at the narrowest part of the constriction in such a way that a counterflow of circulating currents produced by a chain of vortices moving in the central channel maintains the current density J s (y) < J d everywhere along the curved edge of the film except for the vortex entry point. This condition defines the spacing a(I) between the vortices in the chain. However, above a certain current I > I 1 , a single chain of vortices can no longer maintain J s (y) < J d along the rest of the constriction edge, leading to nucleation of an additional stem as seen in Fig. 2n. Our calculation presented in Supplementary Note 6 and Supplementary Figs 5 and 6 show that the second stem appears at the current \({I_1} = {I_{\rm{c}}}\left[ {1 + {{\left( {5\sqrt 3 \xi {\rm{/}}d} \right)}^{1/2}}\lambda {\rm{/}}R} \right]\) which depends on the radius of curvature R of the constriction. For ξ = 46 nm, d = 75 nm, λ=96 nm, and R = 2 μm, we obtain I 1 = 1.11I c in good agreement with the observed I 1 = 1.09I c at B a = 2.7 mT in Fig. 3b. In addition, the edge roughness can affect the location and the dynamics of stem evolution, favoring stem nucleation at points of local edge protuberances (Supplementary Note 6). We have incorporated the actual details of the edge shape of our sample derived from the SEM image into our TDGL simulations resulting in the observed asymmetry between the vortex channels in the upper and lower parts of Fig. 6a, b.

As the magnetic field increases, the width of the vortex-free region near the edges and vortex velocities decrease. Figure 2 shows how dissipative vortex structures evolve from a few mesoscopic chains and branches sustaining extremely high vortex velocities at low field (Fig. 2j, n) to a multi-chain structure with much lower vortex velocities at higher fields (Fig. 2l, p). Remarkably, the vortex channeling is preserved even at high fields that would usually be associated with the conventional flux flow of an Abrikosov lattice. The dynamic structure revealed in Fig. 2p, in which vortices move in parallel channels, appears consistent with the predictions of the moving Bragg glass theory57, thus providing microscopic evidence with a single vortex resolution for the existence of this dynamic phase.

In conclusion, this work uncovers the rich physics of ultrafast vortices in superconducting films and offers a broad outlook for further experimental and theoretical investigations. By proper sample design and improved heat removal it should be possible to reach even higher velocities for investigation of non-equilibrium instabilities17, 18, 50,51,52, 58. Our detection of vortices moving at velocities of up to 20 km/s, significantly faster than previously reported, strengthens the recently renewed appeal of vortex-based cryogenic electronics59. The observed frequencies of penetration of vortices in excess of 10 GHz may be pushed to the much technologically desired THz gap in the case of dynamic Abrikosov-Josephson vortex phases. This work shows that the SOT technique can address some outstanding problems of nonequilibrium superconductivity and ultrafast vortices in type II superconductors as well as dynamics of the intermediate state in type I superconductors on the nanoscale. These issues can also be essential for further development of superconducting electronics, opening new challenges for theories and experiments in the yet unexplored range of very high electromagnetic fields and currents.