# Current-induced self-organisation of mixed superconducting states

Xaver S Brems, Sebastian Mühlbauer, Wilmer Y Córdoba-Camacho, Arkady

A Shanenko, Alexei Vagov, José Albino Aguiar, Robert Cubitt

► To cite this version:

Xaver S Brems, Sebastian Mühlbauer, Wilmer Y Córdoba-Camacho, Arkady A Shanenko, Alexei Vagov, et al.. Current-induced self-organisation of mixed superconducting states. Superconductor Science and Technology, 2022, 35 (3), pp.035003. 10.1088/1361-6668/ac455e. hal-05330203

HAL Id: hal-05330203

https://hal.science/hal-05330203v1

Submitted on 30 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Current-induced self-organisation of mixed superconducting states

To cite this article: Xaver S Brems et al 2022 *Supercond. Sci. Technol.* **35** 035003

View the [article online](#) for updates and enhancements.

## You may also like

- *Controlling soft vortex matter: edge effects on vortex configurations and partial vortices in a superconducting type-II/type-I bilayer*
W Y Córdoba-Camacho, R M da Silva, J Barba-Ortega et al.

- *Physical aging and rejuvenation in the vortex matter of superconducting Nb<sub>x</sub>Zr<sub>5-x</sub>*
Jagdish Chandra, Rashmi Singh and Meghmalhar Manekar

- *Finite-size effects in hyperuniform vortex matter*
Rocío Milagros Besana, Federico Elías, Joaquín Puig et al.
state (IS) of type I SCs, where Meissner regions coexist with normal state domains, the IMS, in contrast to the IS, can not be solely explained by the effect of a non-zero demagnetization factor. It can only exist in the presence of a partially attractive interaction between vortices, which results in an equilibrium vortex distance $a_0$ [2]. Bitter decoration (imaging of vortices via Fe particle decoration [11]) of pure Nb samples shows both laminar and tubular structure of the IMS [2] as also seen in the IS of type I SCs [12]. The behaviour of the IMS has been a subject of several theoretical and experimental studies in recent years [3–10, 13] where the temperature-field phase diagram in samples of different purity was mostly explored. Furthermore, a recent theoretical work [14] studied pinning effects in a two-dimensional cluster glass, as found in and motivated by systems with non-monotonic intervortex interactions such as, e.g. the IMS [2]. Despite these efforts, properties of the IMS in IT SCs are far from being fully understood.

This work reports results of combined small-angle neutron scattering (SANS) and transport measurements of vortex clusters in a single-crystal bulk sample of the archetypal IT SC Nb. We investigate the changes in the IMS vortex configurations induced by an applied current. Our results strongly indicate the emergence of a current-induced superstructure of parallel vortex stripes. It is well known that a current applied to a type II SC in the mixed state creates the Lorentz force acting on the vortices. It is balanced by the drag force resulting in a constant vortex velocity. In a typical type II SC these forces act similarly on all vortices, which then move with almost equal velocities and thus preserve their original Abrikosov lattice arrangement.

In contrast, our results reveal a totally different scenario for the vortex matter in the IMS. The initial configuration of isotropically distributed vortex clusters quickly rearranges itself by elongating in the direction perpendicular to the current flow. It eventually reaches a steady state of parallel vortex stripes. This state is robust—it is independent of the initial configuration and is conserved, when the current is rapidly switched off. The rearrangement dynamics suggests, that here—unlike in type II SCs—the Lorentz and drag forces act differently on different vortices leading to a considerable vortex velocity dispersion at the initial evolution stage.

We argue that the appearance of the stripe superstructure is a generic phenomenon closely related to the physics of the IMS and the infinite degeneracy of the Bogomolnyi point. To demonstrate this we perform numerical simulations using a time-dependent Ginzburg-Landau (TDGL) model with two components. It is the simplest approach to capture essential features of the long-range attractive and short range repulsive vortex-vortex interaction of a single band SC like Nb [2], and to qualitatively reproduce characteristics of its stationary IMS [15]. The numerical simulations reveal details of the cluster elongation and its relation with the asymmetric current distribution and the vortex velocity dispersion. The experimental results and theoretical analysis demonstrate, that we are dealing with a remarkable example of the dynamical pattern formation. This places IT SCs in line with a big group of systems where such self-organized phenomena take place, see, e.g. the

Rayleigh-Bénard convection or Turing reaction-diffusion patterns in chemical reactions and biological systems [16–19].

## 2. Experimental setup

For our study, we used a combined transport measurement and SANS setup, consisting of a low temperature cryostat mounted inside an electromagnet, installed on the SANS diffractometer D33 at the Institut Laue-Langevin [20]. A schematic drawing of the measurement setup can be seen in figure 1(a) with the orientation of the sample, current *I* and applied magnetic field $B_{app}$. The sample stick was equipped with normal conducting copper current leads in the upper part and NbSn superconducting current leads, spot welded to the sample in the lower part, enabling high currents with minimized ohmic heating effects in the vicinity of the sample. For experiments dealing with the effect of current on the IMS, Helium was allowed to condense in the sample space at $T = 4$ K, completely covering the sample in order to efficiently remove any heat created by vortex flow. Thin normal conducting copper voltage leads, required to record characteristic I-V curves, were attached to the sample using silver conducting paint. The distance between the voltage leads was $d \approx 8$ mm.

We used a thin strip of a single crystal Nb sample ($14 \times 1 \times 0.1 \text{ mm}^3$) in our study. The 100 crystal direction was measured to be perpendicular to the large sample face within a few degrees using a neutron Laue diffractometer. As indicated in figure 1(a), the large face of the sample was aligned perpendicular to the magnetic field direction, which results in a large demagnetizing factor. Small cadmium sheets were used to mask the current and voltage leads. The current *I* was applied perpendicular to the magnetic field *B* along the sample. The magnetic field was initially aligned parallel to the direction of the incident neutron beam but both cryostat and magnetic field could be rocked by the angles $\phi$ around a horizontal axis (x-axis) and $\omega$ around a vertical axis (y-axis) with both axes perpendicular to the neutron beam.

For the SANS measurements, the collimation was set to $12.8 \text{ m}$ with a sample aperture of $\approx 5 \times 2 \text{ mm}^2$. The scattered neutrons were detected using a position-sensitive 2D detector placed $12.8 \text{ m}$ behind the sample. A medium resolution setup using a neutron wavelength $\lambda = 10 \text{ Å}$ and square source aperture with $30 \times 30 \text{ mm}^2$ cross section was used for mapping the IMS phase diagram. For experiments dealing with the effect of current on the IMS we used a high resolution setup with neutron wavelength $\lambda = 14 \text{ Å}$ and a round source aperture with $d = 20 \text{ mm}$ diameter. The full width half maximum wavelength spread was $\Delta\lambda/\lambda = 10\%$ for both setups.

The niobium single crystal sample was prepared by spark erosion. It was cut from a Nb single crystal obtained from Heraeus previously used in other experiments on the IMS [6, 7]. The sample was left untreated after spark erosion cutting, since there are indications, that additional surface treatments (e.g. by means of electropolishing) increases the critical current [21]. The demagnetization factor is $D = 0.87 (B \perp$ to large sample face). From I-V measurements at room temperature and just above the transition temperature ($T_S = 9.5 \text{ K}$) a normal state
<figure><img src="image_2.png" /><figcaption>Figure 1. (a) Geometry of the SANS experiment (b) Average internal magnetic field $B_{int}$ and the colorplot of the IMS intensity. White crosses mark the applied magnetic field $B_{app}$ and temperature $T$. (c) Integrated intensities derived from the first order Bragg peaks. The dashed lines are guides to the eye.</figcaption></figure>

resistivity of $\rho_n(9.5) = (3.7 \pm 0.2) \times 10^{-10} \Omega m$ and residual resistivity ratio of $RRR \approx 390$ were deduced.

Numerical calculations were done using the computational facilities, equipped with NVIDIA GPUs, of the Laboratório de Supercondutividade e Materiais Avançados (SUPERLAB) do Departamento de Física da Universidade Federal de Pernambuco.

## 3. Results

### 3.1. Stationary IMS

We first focus on the stationary IMS without applied current. Figure 1(a) gives an example of a typical 2D SANS detector image of the scattering intensity profile, where the direct beam is excluded (central black circle). The image shows the first order Bragg peaks stemming from the ordered flux line lattice inside the mixed state domains and the IMS scattering around the blacked out direct beam resulting from the magnetic contrast between Meissner state domains and mixed state domains. We observe 12 Bragg spots resulting from two iso-scales domains of Nb ($\alpha \approx \beta \approx 65^{\circ}$, $\gamma \approx 50^{\circ}$) as reported previously [9]. Figure 1(b) shows the temperature dependence of the average internal magnetic field $B_{int}$ for a few selected values of the applied field $B_{app}$. The value of $B_{int}$ is obtained using the relation

$$B_{int} = \frac{\phi_0 \sqrt{3} q_{VL}^2}{8\pi^2}, \quad (1)$$

where $\phi_0$ is the elementary flux quantum and $q_{VL}$ is the average $|q|$-value of the first order Bragg peaks $^6$. The color density

map in figure 1(b) represents the IMS scattering obtained by summing the intensity around the direct beam (between the white circles in figure 2(a)). The data were collected during warming (W) after field cooling (FC) in an applied magnetic field $B_{app}$ from $T = 10$ K to $T = 4$ K and are corrected for a high temperature background ($T = 9.2$ K). In the data analysis we assume that the diamagnetic response after FC is vanishingly small using the fact that other samples obtained from the same single crystal demonstrate a negligible diamagnetic contribution in FC magnetization measurements [6].

At low temperatures, the internal field $B_{int}$ in figure 1(b) follows the same universal temperature dependence $B_{int}^*(T)$ for all values of the applied field $B_{app} < 800$ G. This behavior is known to be a hallmark of the IMS [6]. When $T$ increases further, $B_{int}(T)$ departs from $B_{int}^*(T)$ eventually approaching the corresponding value of $B_{app}$ in a way typical for the conventional mixed state. The crossover between the IMS and the conventional mixed state is also observed in the colour density plot in figure 1(b) representing the summed IMS scattering intensity around the direct beam, called IMS intensity. Starting at low temperatures, the IMS intensity peaks at an applied magnetic field of $B_{app} = 400$ G. It decreases for increasing and decreasing magnetic fields. For a given applied magnetic field $B_{app}$, the IMS intensity decreases with rising temperature and vanishes once the internal magnetic field $B_{int}$ leaves the common temperature dependence. Figure 1(c) shows the integrated intensities derived from the first order Bragg peaks. In agreement with [6], we see a linear increase with falling temperature (typical of the vortex lattice in the standard mixed state) and observe a downward bent curve (600 G) and an additional kink (400 G and 200 G) most pronounced at $B_{app} = 400$ G, indicating the transition to the IMS for applied magnetic fields $B_{app} < 800$ G. Inside the IMS the integrated intensity does not decrease with decreasing temperatures, as would be expected when approaching the Meissner state with vanishing internal field. We attribute this to considerable pinning effects, which lead to a macroscopic trapping of the magnetic flux. Note however, that these pinning effects do not hinder a microscopic rearrangement of the vortex lattice, when entering the IMS regime as seen from the change in internal magnetic field $B_{int}$ with $T$ shown in figure 1(b) and reported

<footnote>$^6$ In this equation we used the approximation of a perfect hexagonal lattice, since we were not able to resolve the exact position of the first order Bragg peaks for low applied magnetic fields and temperatures inside the IMS regime. Therefore it was not possible to extract the exact reciprocal unit cell, which has been shown to change as a function of applied field $B_{app}$ and temperature $T$ [9]. Equation (1) still gives a good approximation of the internal field $B_{int}$, since the observed angles at high applied magnetic fields and low temperatures are close to 60° and we only assume slight deviations for lower applied magnetic fields and high temperatures.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** SANS images from applied current measurements in a magnetic field of $B_{\text{app}} = 500$ G at $T = 4$ K: (a)–(c) FC approach with incremental current ramp from $I = 0$ A up to $I = 40$ A. (d) Subsequent current quench (QU) to $I = 0$ A. (e) ZFC approach without a current and (f) subsequent current ramp to $I = 40$ A. The measured flux flow voltage $V_{\text{ff}}$ is shown above each panel.</figcaption></figure>

elsewhere [6]. The temperature of the observed kink/the start of the downward bent part agrees sufficiently well with the temperature of deviation from $B_{\text{int}}^*(T)$ shown in figure 1(b). In summary we established the range of the IMS in good agreement with literature [6].

## 3.2. Current induced changes in the IMS

We will now investigate changes in the IMS induced by the applied current $I$. Previous to the neutron experiment we identified the regime of the vortex movement by measuring the $I-V$ curves of the sample for different values of $B_{\text{app}}$ at $T = 4$ K. The results are shown in figure 6 in appendix A.

Figure 2 shows examples of SANS images obtained with various applied currents in a magnetic field of $B_{\text{app}} = 500$ G at $T = 4$ K. Shown is the sum over rocking angles in the interval $\phi, \omega \in [-4^\circ, 5^\circ]$. The current is applied in x-direction perpendicular to the magnetic field, as indicated in panels (a) and (b). The flux flow voltage $V_{\text{ff}}$ is shown above the panels. All the data were corrected for a zero field cooled (ZFC) background at $T = 4$ K.

Figures 2(a)–(d) show measurements following a FC protocol: The sample was cooled in an applied magnetic field of $B_{\text{app}} = 500$ G with no applied current from $T = 10$ K to $T = 4$ K. We then performed rocking scans without applied current (figure 2(a)) and after applying an external current in incremental steps up to $I = 30$ A (figure 2(b)) and after further increasing to $I = 40$ A (figure 2(c)). While still being in an applied magnetic field $B_{\text{app}} = 500$ G and cold at $T = 4$ K we quenched the current to $I = 0$ A and subsequently performed another set of rocking scans shown in figure 2(d).

We also employ a different protocol where the sample is first cooled from $T = 10$ K to $T = 4$ K in the zero field (zero-field cooling—ZFC) and then the field is ramped up to $B_{\text{app}} = 500$ G. The SANS image obtained after this protocol for $I = 0$

is shown in figure 2(e) while figure 2(f) illustrates the changes after the current is further increased to $I = 40$ A.

We first consider the Bragg peaks stemming from the vortex lattice in the mixed state domains. For the FC protocol at $I = 0$ A (figure 2(a)) and up to $I = 20$ A (not shown) both the position and intensity of the Bragg peaks from two well ordered isosceles vortex lattice domains remain unchanged. When the current passes the critical value $I \gtrsim 30$ A and the vortices start to move resulting in $V_{\text{ff}} \neq 0$ the Bragg peaks get smeared out in the azimuthal direction and decrease in intensity due to the lattice disorder induced by the motion (figures 2(b) and (c)). The disorder is retained after the current is quenched to $I = 0$ A (figure 2(d)).

The ZFC—field ramp protocol does not create a well ordered vortex lattice as is indicated by the smeared Bragg peaks in figure 2(e). An external current of $I = 40$ A (figure 2(f)) restored the equivalent FC case (figure 2(c)) showing an ordered vortex lattice with azimuthally smeared out Bragg peaks. The Bragg peaks obtained after the FC and ZFC—field ramp protocol with external current $I = 40$ A are practically indistinguishable in the 2D SANS image. A further detailed analysis of the vortex lattice Bragg peaks can be found in appendices B and C (figures 7 and 8).

Details of the IMS structure are reflected by the scattering in the vicinity of the direct beam—between the white circles in figure 2(b). For the FC protocol and $I = 0$ A (figure 2(a)) up to $I = 20$ A (not shown) the IMS scattering is isotropic. However, when $I \gtrsim 30$ A (figures 2(b) and (c)) one observes a notable anisotropy in the scattering pattern which becomes elongated in the horizontal $q_x$-direction. The anisotropy is preserved after a rapid current quench to $I = 0$ A (figure 2(d)). After the ZFC—field ramp protocol the IMS scattering also demonstrates a slight anisotropy in the horizontal direction (figure 2(e)). A subsequent current ramp to $I = 40$ A results in the same IMS scattering pattern as obtained after the FC protocol (cf figures 2(c) and (f)) showing a clear anisotropic IMS scattering around the direct beam.

Further details of the IMS scattering are presented in figures 3(a)–(c) that plot the radially averaged IMS scattering intensity as a function of the azimuthal angle $\xi$ as defined in figure 1(a). The IMS intensity is evaluated between the white circles shown in the SANS image in figure 2(b)<fnref n="7" />. The detector pixel size and the small scattering angles limit the azimuthal bin size to $\Delta\xi = 20^\circ$. The corresponding values of $V_{\text{ff}}$ and $I$ are shown above each panel.

At $I = 0$ (figure 3(a)) up to $I = 20$ A (not shown) the IMS scattering is isotropic being independent of the azimuthal angle $\xi$. At $I \gtrsim 30$ A the scattering shown in figures 3(b) and (c) is angle-dependent having its maxima at $\xi = \pi/2$ and $3\pi/2$ (in the horizontal $q_x$-direction) and minima at $\xi = 0, \pi$ and $2\pi$ (in the vertical $q_y$-direction). When the current increases to $I = 40$ A the difference between the maximal and minimal values increases (see figure 3(c)).

<footnote n="7">⁷ To obtain the angle dependence shown in figures 3(a)–(c) we take the intensity in figures 2(a)–(c), respectively, represent it as function of the radial distance from the beam center and the azimuthal angle $\xi$ and then integrate the intensity over the radial distance inside the sector of interest.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Azimuthal averages ((a)–(c)) and rocking curves ((d)–(f)) of the IMS scattering in current experiments in a magnetic field of $B_{\text{app}} = 500$ G at $T = 4$ K. The solid line is a fit of two Gaussians. Shown are rocking scans around the x-axis in angle $\phi$ (curves (1) + (3)) and the y-axis in angle $\omega$ (curves (2) + (4)). The sectors for the azimuthal averaging and rocking scans are depicted in figures 2(b) and (a), respectively. The current $I$ and the voltage $V_{\text{ff}}$ are shown above each column.</figcaption></figure>

Finally, figures 3(d) and (e) show the rocking curves of the IMS scattering inside the sectors 1–4 shown in figure 2(a)<fnref n="8" />. The summed scattering intensities corresponding to sectors 1 and 3 are represented by green and red circles, respectively, and plotted as functions of the rocking angle $\phi$ around the x-axis. The summed scattering intensities inside sectors 2 and 4 are represented by orange and blue circles, respectively, and plotted as functions of the rocking angle $\omega$ around the y-axis. When the current is absent, $I = 0$ A, the IMS intensity dependence on the rocking angle is qualitatively similar for all four sectors (figure 3(d)). However, at $I = 30$ A figure 3(e) reveals a clear difference between the rocking scans of the vertical sectors (1 and 3) and rocking scans of the horizontal sectors (2 and 4). The intensity in the horizontal sectors increases and develops a sharper angle dependence, whereas the rocking curves of the vertical sectors flatten. This effect is getting more pronounced when the current increases (figure 3(f)).

## 4. Discussion

In the absence of transport current, the characteristics of the IMS such as depicted in figures 1(b) and (c) are in agreement with earlier experiments [6]. The hallmark of the IMS, a $B_{\text{app}}$ independent, but $T$ dependent internal magnetic field $B_{\text{int}}$ and the presence of very small angle scattering inside the IMS regime is clearly observed.

In our study we see, that the IMS still persists when vortices start moving, as evidenced by I–V-characteristics in different applied magnetic fields (see figure 6 in appendix A) and the neutron data. Figure 2 clearly demonstrates, that the hallmarks of the IMS—very small angle scattering and the constant internal field $B_{\text{int}}$—are preserved in the state of vortex movement.

Moreover, the applied current gives rise to an elongation of the IMS domains, that coincides with the onset of flux flow. The elongation is manifested by the transition from isotropic to anisotropic IMS scattering parallel to the applied current in x-direction at the onset of the vortex movement marked by a finite voltage at $I_c \approx 30$ A, as seen in figures 2(b) and 3(b). The anisotropy in the IMS scattering (figures 2(c) and 3(c)) is more pronounced with increasing current up to the maximal applied current $I = 40$ A. Furthermore, the rocking curve of the IMS scattering with respect to the y-axis in angle $\omega$ becomes more pronounced, while the rocking curve around the x-axis in angle $\phi$ flattens with increasing current (figures 3(d)–(f)). Both the anisotropy of the IMS scattering and the flattening of the rocking curve point to an emerging sheet-like superstructure that is orthogonal to the applied current and, thus, parallel to the flux lines motion.

At the same time, the current induces an overall increase of disorder of the vortex lattice. This fact is evident from the broadening of the vertical and horizontal Bragg spots of the vortex lattice with increasing current, seen from the increased azimuthal smearing (figure 2), the radial width (figure 8 in appendix C) and the rocking curve width (figure 7 in appendix B) of the Bragg peaks, as well as from the corresponding loss of their intensity. The latter cannot be attributed to an ohmic heating, since a temperature increase would also lead to a decrease in the total IMS intensity (figure 1(b)), not seen in figures 3(a)–(c).

A careful study of the rocking curve shape of the vortex lattice Bragg peaks as presented in appendix B allows for an estimation of the vortex bending due to the magnetic field induced by the applied current referred to as magnetic self field. The applied current splits into a bulk component $I_{\text{bulk}}$ and a non-dissipative surface component $I_{\text{surf}}$. Assuming the bulk current equals the applied current would result in a flat top shape of the rocking curves of vertical Bragg spots with a broadening of almost 30° due to the additional magnetic self field with a maximal value at the surface of $B_{\text{surf}} \approx 250$ G. In contrast, our measurements only show an additional broadening of $\Delta\sigma \approx 0.7 \pm 0.1^\circ$ at $I = 40$ A (rocking curves of the vertical Bragg spots with respect to the rocking curves of the horizontal Bragg spots). Our calculations therefore suggest a dissipationless critical current of $I_c = 36$ A flowing on the sample surface and the excess current of $I_{\text{bulk}} = 4$ A penetrating the bulk showing good agreement with the measured rocking curve. This splitting into bulk current and dissipationless surface current has also been shown previously [22–24] and is further strengthened by the shape of the I–V-curve, showing the absence of a discontinuity at the critical current $I_c$ (figure 6 in appendix A) in agreement with literature [22].

The current-driven IMS achieves a steady state, which represents a key result of this study. The steady state is robust,

<footnote n="8">⁸ To obtain the rocking curves we take the rocking scan intensity (not summed over angles $\phi$ and $\omega$ as is done in figures 2), and sum the intensity inside the sectors of interest for each rocking angle $\phi$ and $\omega$. The magnetic field was aligned via the standard method of using the Bragg spots of the vortex lattice. Rocking angles $\phi < -2^\circ$ were not achievable due to the size of the opening window of the electromagnet.</footnote>
remaining after the current is abruptly quenched to zero. This is illustrated in figure 2 that demonstrates no visible changes in the small-angle scattering pattern. Both, the elongated IMS superstructure and the disorder in the vortex lattice are retained, when the current is off, with the exception of the additional broadening associated with the current-induced field. Furthermore, the transport current suppresses the hysteretic behaviour of the IMS pattern with respect to FC-ZFC protocols. When the current is absent the ZFC protocol does not produce a well ordered lattice [6]. A slight anisotropy in the $q_x$-direction observed after applying the ZFC protocol is explained by flux entering the sample overcoming the Bean barrier [25] preferentially via the short sample dimension and therefore forming elongated IMS domains in the vertical y-direction in figure 1(a). However, a subsequent current ramp to $I = 40$ A drives the system to the state, which is achieved after applying the FC protocol for the same current. This fact is evident in both the Bragg peak and very small-angle IMS scattering shown in figure 2.

Our experimental results are compared with theoretical simulations performed using the model of two linearly coupled Ginzburg-Landau (GL) equations. It is the simplest model that allows one to study qualitative features of current-driven IMS vortex configurations (an interested reader can find justifications for the choice of the model and details of the calculations in appendix D). The model was originally derived for SCs with two separate pockets on the Fermi surface (bands) [26, 27] and applied to explain peculiar properties of the vortex matter in two-band SCs like $MgB_2$ [28]. However, vortex configurations given by this model [28, 29] are very similar to those calculated theoretically [3] and observed experimentally [2, 30] for low-$\kappa$ single-band BCS SCs. This is not surprising, given that both models reveal qualitatively similar inter-vortex interactions, in particular, attraction at large and repulsion at small distances [31, 32] and notable contribution of many-vortex interactions [4, 15, 33], which are both key for the IMS vortex pattern formation [3]. In this work we exploit this similarity to investigate the evolution of the IMS in the single band SC Nb.

Our calculations do not account for additional factors such as crystal anisotropy, pinning, sample geometry, and the stray field, which do affect details of the vortex arrangement inside the mixed state. For example, the crystal anisotropy gives rise to the isosceles vortex structure, evident in the Bragg peaks in SANS images. Also, the presence of pinning centers reveals itself in the stability of the IMS when the temperature is lowered. Nevertheless, these factors are not expected to play a major role in the current-induced formation of the vortex superstructure. Indeed, the final IMS state, obtained in the end of the time evolution, appears independent on the initial configuration, for example, on whether it is reached by following the ZFC or FC protocols (see figure 2). Our theoretical simulations also demonstrate that the appearance of such superstructures is described qualitatively well without considering these additional factors.

When the current in the sample exceeds the critical value $I > I_c$, vortices start moving in the perpendicular direction

to the current flow due to the Lorentz force, and the IMS configuration changes. Figure 4 shows snapshots of the time evolution of the initial IMS with random vortex clusters in figure 4(a). While moving, the clusters elongate in the direction parallel to their movement and merge. Finally, the time evolution creates a superstructure of vortex stripes elongated perpendicularly to the average current. Results of the simulations are fully consistent with the experimental observations. The type of the vortex structure inside the clusters depends on a point in the IT phase diagram that corresponds to the superconducting material [3]. When the material is close to type II ($\kappa \gtrsim \kappa_0$) vortices form a lattice, while in materials close to type I ($\kappa \lesssim \kappa_0$) the lattice melts becoming a liquid. However, for any initial IMS configuration, the current-driven evolution eventually yields qualitatively similar arrangements of vortex stripes.

The elongation of the vortex structures and subsequent formation of the stripe superstructure can also be traced in the Fourier transform (FT) of the spatial field profile, which corresponds directly to the measured intensity in a scattering experiment. The FT calculated from the initial field configuration in figure 4(a), obtained before the current is applied, is shown in figure 5(a). Correspondingly, figure 5(b) shows the results of the FT for the stripe configuration in figure 4(f). One sees that the initially almost isotropic FT profile changes to the visibly elongated structure. This fully agrees with the SANS images in figure 2. The absence of vortex lattice Bragg peaks in the FT is explained by the simplified numerical simulation neglecting the crystal anisotropy.

It is worth noting that several recent studies [34, 35] on type II superconducting films also demonstrated spontaneously emerging vortex configurations, elongated in the direction

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Numerical simulations of the IMS time evolution (the spatial distribution of the magnetic field is given). Panel (a) represents the initial configuration of the vortex clusters, panels (b)–(f) give snapshots at later times. The flowing current $I$ and the vortex motion with velocity $v_L$ are shown by arrows in panel (f). Details of the model and its microscopic parameters are described in appendix D.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The Fourier transform of the numerically obtained field distribution Panel (a) at the start of the current-induced evolution shown in figure 4(a), and Panel (b) when the stripes are formed as shown in figure 4(f). Appearance of stripes is reflected in the elongated central (small angle) part of the Fourier transform.</figcaption></figure>

perpendicular to the current flow, when the applied cur-
rent approached the pair-breaking limit. Those ‘vortex rivers’
appear similar to quasi-1D vortex chains and enable an
ultra-fast vortex motion inside them. In contrast, the stripes
observed in our study are much wider and the rearrangement
takes place at much lower currents only slightly exceeding the
critical value of vortex movement.

The rearrangement process into the elongated IMS super-
structure is a result of a combined action of three main factors:
the current-induced Lorentz force, the drag force and the inter-
vortex interaction potential. Their action on the IMS in IT SCs
demonstrates peculiar features, that distinguish these materi-
als from conventional type II SCs. In IT materials, the vortex
core size is comparable with the inter-vortex distances result-
ing in spatially non-monotonic vortex interactions with sizable
many-vortex contributions [4, 15]. As a result, vortices neces-
sarily clusterize and can no longer be viewed as isolated ‘ele-
mentary particles’ of the mixed state [3]. Consequently, the
IMS is characterized more by the collective vortex behavior
than by the properties of an individual vortex. The total action
of the applied current thus depends on a configuration of vor-
tex clusters. In principle, this can be seen already in type II
SCs where the current modifies the vortex-vortex interactions
such that it starts to depend on the vortex orientation relative to
the current flow [36]. The proximity to the degenerate Bogo-
molnyi point enhances this geometry dependence. The cal-
culations show that the applied current flows primarily along
boundaries between the mixed and Meissner state resulting in
the Lorentz force that pulls that boundary in the normal direc-
tion towards the Meissner phase (see figure 9 in appendix E).
The asymmetry of the current flow creates a non-zero cumu-
lative force which depends on the size and shape of a vortex
structure. This results in the dispersion of the vortex velocit-
ies, which in turn distorts the structure shape. Numerical cal-
culations demonstrate that vortex structures, that are large and
elongated in the motion direction, move faster (see appendix E
for details of the time evolution of vortex clusters of differ-
ent shape and size). It is intuitively explained by the larger
current density and thus larger Lorentz force at the bound-
ary of a larger cluster. Combined with the drag force and the
non-monotonic vortex interaction, that keeps the inter-vortex
distance, this configuration-dependent Lorentz force further

elongates vortex structures in the direction perpendicular to the
applied current eventually creating a superstructure of stripes
as shown in figure 4.

We finally discuss the limitations of our theoretical
approach before we highlight future perspectives of trans-
port phenomena in the IMS to model percolation and non-
trivial flow phenomena in two domain systems. Our model
of two coupled TDGL equations does not take into account
many factors such as pinning, crystal anisotropy, and phonon-
induced relaxations. It also uses a simplified model that
assumes that the average current is injected through the sys-
tem and then redistributed by the vortex configuration. In real
samples the current in the inner part can only flow in regions
with non-vanishing curl of B according to Ampere’s law.
Transport current is therefore constrained to the mixed state,
the surface and the interfaces between the Meissner and the
mixed state domains. In turn, current cannot exist in the Meiss-
ner regions where the magnetic field is zero. In order to pass
a current through the bulk of our sample, we therefore need
a connected mixed state bridging the contacts. For our experi-
ment, the volume fraction of the mixed state is estimated using
$f_{\text{MS}} = V_{\text{MS}}/V = B_{\text{app}}/B_{\text{int}}$ with the total sample volume $V$ and
the mixed state volume $V_{\text{MS}}$. This estimation holds since we
expect a negligible diamagnetic behaviour in our sample [6]
and therefore the applied magnetic field is equal to the internal
magnetic field averaged over the whole sample volume. In our
case this estimate yields $f_{\text{MS}} \approx 70\%$ at the field $B_{\text{app}} = 500$ G.
This is above the percolation threshold of $\rho_c = 44\%$ in 2D sys-
tems [37]. Therefore we can safely assume the connectivity of
the mixed state domains of our sample takes place at the onset
of the flux flow, where we made use of the quasi two dimen-
sional character of the vortex lattice<fnref n="9" />. This renders our sim-
plified theoretical model as a valid description of our experi-
mental findings.

We note, however, that at lower magnetic fields and the
associated lower volume filling of the IMS below the per-
colation threshold, one expects a much more complex situ-
ation. When the connectivity is broken, one needs to employ a
fully inhomogeneous theoretical description for a finite sample
that takes into account the influence of the stray fields out-
side the sample to capture the essential physics of the sys-
tem. In this case two orthogonal flows restricted to a single
phase of a two phase system might lead to interesting order-
ing phenomena. Further complications arise due to the influ-
ence of the surfaces of the sample, where (a) the vortices
are nucleated and destroyed, (b) pinning is significantly dif-
ferent as compared to the bulk and (c) the influence of the
dominant surface current is unclear. Besides the unknown
balance of surface versus bulk pinning, it is also unclear
whether the IMS structure is nucleated at the surface of the
sample or forms as a steady state deep inside the bulk of the
material.

<footnote n="9">⁹ At the same time this means that the Meissner state domains are not percolating as the Meissner state volume fraction is given as $1 - f_{\text{MS}} = 30\%$ (again assuming zero magnetization).</footnote>
## 5. Conclusion

This work studies the evolution of the IMS domains in the IT SC Nb under the influence of an external transport current using a combined SANS and transport measurement technique.

Our main result is an observation of a transition from isotropic to anisotropic IMS scattering, indicating, that the IMS rearranges itself into a stripe superstructure in the vortex movement regime. The stripe pattern is aligned perpendicular to the current direction along the motion of the vortices. A close examination of the rocking scans showed a splitting of the current into bulk and surface component with the latter being dominant. Most importantly, the absence of the hysteretic behaviour proves that the elongated superstructure is a steady state of the moving IMS. Numerical simulations of the time evolution of IMS vortex configurations using a model of two coupled Ginzburg-Landau (GL) equations qualitatively reproduced our experimental results and revealed details of the cluster elongation. Our findings highlight the importance of the IMS as a model system for universal domain physics and demonstrate, that we are dealing with a remarkable example of a self-organized pattern formation phenomenon. The current-induced movement of vortices in the IMS might act as a model system for the study of percolation and non-trivial (orthogonal) flow and self-ordering phenomena in two domain systems.

## Data availability statement

The data that support the findings of this study are available upon reasonable request from the authors.

## Acknowledgments

We express our gratitude to F Marchal and F Lapeyre for their support with sample preparation and M Bonnaud for support with the experiments. Further thanks are due to A Backs for fruitful discussions and his contributions to this project in the early stage. This work is based upon experiment EASY-568 [38] and experiment 5-31-2748 [39] performed at the instrument D33 at the Institut Laue Langevin (ILL), Grenoble, France. W Y C-C thanks the Brazilian agency Coordenação de Aperfeiçoamento de Pessoal de Nível Superior CAPES, Grant No. 88887.582633/2020-00. A V thanks the SUPER-LAB and his chair Prof. J Albino Aguiar for providing financial support for his visit to the SUPERLAB in March 2020. A A S thanks the Departamento de Física da Universidade Federal de Pernambuco for the hospitality during the years he worked there. J A A thanks the Brazilian agencies Conselho Nacional de Desenvolvimento Científico e Tecnológico CNPq, Grant No. 313810/2020-6 and Fundação de Amparo a Ciência e Tecnologia do Estado de Pernambuco FACEPE, Grant No. APQ-0936-1.05/15. X S B acknowledges the support by the Studienstiftung des deutschen Volkes. A A S and A V are grateful for support from the HSE University Basic Research Program. We also thank E Babaev for his interest in

our work, fruitful discussions, and drawing our attention to a number of results obtained using the two-component model.

## Appendix A. Transport measurement

We used a dedicated combined transport measurement and SANS setup presented in section 2 of the main text, which allowed us to study the I-V-characteristics of our sample prior to the experiment and also monitor the voltage response of our sample to an external current during the SANS experiment. During pre-characterization we measured I-V-curves at 20 different applied fields in the range of $B_{app} \in [500 \text{ G}, 1720 \text{ G}]$ at $T = 4 \text{ K}$. In contrast to the neutron experiment, we did not follow a FC measurement protocol for the pre-characterization measurements, since changing magnetic field $B_{app}$ implies heating up the sample above the transition temperature and therefore boiling away the condensed He in the cryostat. The data points of the I-V-curves of the pre-characterization were recorded on an average frequency of $\approx 0.1 \text{ Hz}$.

During the neutron experiment and the measurement of the I-V-curve the sample was FC. The I-V-data points during the neutron experiment represent an average of the voltage recorded over the whole duration of a rocking scan, which leads to an acquisition time of $\approx 2.5 \text{ h per point}$.

Figure 6 summarizes the I-V-characteristic of our sample at $T = 4 \text{ K}$ under liquid He at different applied magnetic fields $B_{app}$. Figure 6(a) shows examples of I-V-curves at $T = 4 \text{ K}$ for a few selected values of the applied field $B_{app}$ collected prior to the neutron experiment. As expected from literature [40], we see a decrease in the critical de-pinning current for vortex movement with increasing applied magnetic field $B_{app}$. Furthermore the slope of the I-V-curve in the linear regime related to the flux flow resistance $R_{ff}$ is increasing with increasing applied magnetic field$^{10}$.

For a normal type-II SC, exhibiting an Abrikosov lattice covering the whole sample for applied fields $B_{app} > B_{c1}$, the flux flow resistance has been shown to be proportional to the applied field $B_{app}$ according to (A1) derived in the Bardeen-Stephen model with the upper critical magnetic field $B_{c2}$ and the normal state resistance $R_n$ [42].

$$R_{ff} \propto R_n \frac{B_{app}}{B_{c2}}. \qquad (A1)$$

This should still hold for an IT SC in the pure mixed state. Inside the IMS regime, the sample splits into domains of mixed state and Meissner state. Assuming a homogeneous current distribution constraint to mixed state domains of the sample leads to a local decrease in the sample cross section and therefore an increased local resistance by a factor of $B_{int}/B_{app} = f_{MS}^{-1}$. The measured voltage is an average over

<footnote>$^{10}$ At this point we note that especially for low fields, the voltage response shows an extended regime in vicinity of the de-pinning current, where the voltage response is not ohmic, but rather shows a curved shape consistent with plastic de-pinning [41].</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) Transport measurements prior to the neutron experiment at T = 4 K as a function of different applied magnetic fields B<sub>app</sub>. (b) Comparing the I–V-curve measured prior to the experiment with the in-situ measurement in a magnetic field of B<sub>app</sub> = 500 G at T = 4 K. The insets show the corresponding 2D detector image. The dashed lines are guides to the eye.</figcaption></figure>

the whole sample and since only the mixed state contributes to the voltage build-up, we get an additional factor of B<sub>app</sub>/B<sub>int</sub> = f<sub>MS</sub>, which cancels the factor resulting from the local cross section decrease. This however only holds, if we have a sample with zero magnetization.

Figure 6(b) shows the I–V-curve measured prior to the neutron experiment (blue points) and the I–V-curve derived from the voltage measurement during the neutron experiment (orange points), both in an applied magnetic field of B<sub>app</sub> = 500 G. The 2D detector images corresponding to the applied current are shown in the insets. The dashed line is a guide to the eye. When comparing the two I–V-curves they agree well within errors. The critical de-pinning current for vortex movement is slightly smaller during the neutron experiment (I<sub>cexp</sub> ≈ 30 A vs I<sub>cpre</sub> ≈ 35 A). The data points match well for I > 35 A. The slight deviations can be explained by the different time scales on which the I–V-curves were recorded, as mentioned above.

## Appendix B. Flux line bending

We see from Ampere’s law (B1) that a transport current with a non-vanishing bulk current I<sub>bulk</sub> leads to curved flux lines due to its self field B<sub>self</sub>, whereas a pure surface current I<sub>surf</sub>, retains

the alignment parallel to the external magnetic field B<sub>app</sub> as schematically depicted in figures 7(a) and (b).

$$ \nabla \times \mathbf{B} = \mu_0 \mathbf{J}. \qquad (B1) $$

The curvature of the flux lines can be extracted from the width of rocking curves of the flux line lattice Bragg peaks as previously reported [23, 43]. The shape of the flux lines is affected by the self field resulting from the current distribution of a superconducting current-carrying wire. There are arguments for both, a constant r independent current distribution J(r) = const [23] and a square-root dependent current distribution J(r) ∝ √1/r [44], where r is the distance from the sample’s center. Here we limit the discussion to an assumed homogeneous current distribution over the cross section of the Nb strip resulting in a self field, that linearly decreases when approaching the center of the sample. The resulting maximum tilt angles of the flux lines are given by θ = ±B<sub>self</sub>/B<sub>app</sub> with the maximum value of the magnetic self-field B<sub>self</sub> on the surface of a strip with cross section w × t with w > t and current I, approximated by

$$ B_{\text{self}} = \frac{\mu_0 I}{2w}. \qquad (B2) $$

The resulting broadening of the rocking curve is given by

$$ \text{FWHM}(\phi) = \Delta\phi = 2 \cdot \frac{\mu_0 I}{2wB_{\text{app}}}. \qquad (B3) $$

Due to geometry, the flux lines are only bent in the y-direction, which means only the rocking curves in φ around the x-axis are affected by the magnetic self-field broadening.

Figures 7(c)–(f) give a detailed look on the rocking scans of the first order vortex lattice Bragg peaks in horizontal and vertical direction for different combinations of bulk and surface current. The sectors 1–4 used for the rocking scans are marked in the inset of figure 7(d).

Figures 7(c) and (d) compare the rocking curves of the vertical Bragg peaks (curves (1) and (3)) in rocking angle φ around the x-axis with the rocking curves of the horizontal Bragg peaks (curves (2) and (4)) in rocking angle ω around the y-axis for different applied currents I. Solid lines are fits to the data. For zero applied currents the rocking curves are fitted using a sum of two Gaussians, which captures the double peak-structure. In order to limit the number of free parameters, we made use of the symmetry of the experiment and constrained the relative centers, the intensity and the width of two corresponding rocking curves to be equal. To be more specific, the widths, centers and intensities of the peaks in curve (1) + (3) and (2) + (4) are constrained, respectively. The integrated intensity of an individual peak is given as the product of its weight and the total integrated intensity I<sub>0</sub> of the rocking curve. When the current increases, the individual peaks of the rocking curves are broadened becoming indistinguishable. Therefore, the rocking curves for I = 40 A are fitted using a single Gaussian with an individual set of parameters for each rocking curve. The essential fit parameters of the fits are summarized in table 1.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Scheme of a straight flux line due to the absence of a bulk current. (b) Scheme of a bent flux line caused by the self-field contribution of a non-vanishing bulk current. Panel (c) shows rocking curves and fits (solid lines) of the first order Bragg peaks around the x-axis in angle $\phi$ (curves (1) and (3)) and around the y-axis in angle $\omega$ (curves (2) and (4)) with no applied current while panel (d) depicts the corresponding rocking curves with an applied current of $I = 40$ A. The inset shows the sectors used for the rocking curves in panels (c) and (d). (e) Simulated rocking curves (dotted lines) for different combinations of bulk and surface current. The corresponding angular distribution is shown in the same color code at the bottom of the plot.</figcaption></figure>

<table><caption>**Table 1.** Total integrated intensity $I_0$, the weight of the principal peak, and Gaussian widths $\sigma_i$ of the rocking curve peaks. For zero applied current a sum of two Gaussians was fitted, while for $I = 40$ A a single Gaussian was fitted. $\sigma_1$ describes the width of the higher intensity peak in the double peak structure.</caption><thead><tr><th>Current<br>(A)</th><th>Curve</th><th><i>I</i><sub>0</sub><br>(arb. unit)</th><th>Weight<br>(%)</th><th>σ<sub>1</sub><br>(°)</th><th>σ<sub>2</sub><br>(°)</th></tr></thead><tbody><tr><td colspan="6">Rocking scans around the y-axis (rocking angle ω)</td></tr><tr><td>0</td><td>(1)+(3)</td><td>20.3 ± 0.4</td><td>52 ± 7</td><td>0.94 ± 0.09</td><td>0.64 ± 0.04</td></tr><tr><td rowspan="2">40</td><td>(1)</td><td>14.3 ± 1.3</td><td>100</td><td>2.61 ± 0.31</td><td>—</td></tr><tr><td>(3)</td><td>13.3 ± 6.0</td><td>100</td><td>2.50 ± 0.55</td><td>—</td></tr><tr><td colspan="6">Rocking scans around the x-axis (rocking angle φ)</td></tr><tr><td>0</td><td>(2)+(4)</td><td>9.8 ± 0.6</td><td>82 ± 14</td><td>1.14 ± 0.13</td><td>0.50 ± 0.12</td></tr><tr><td rowspan="2">40</td><td>(2)</td><td>12.9 ± 2.0</td><td>100</td><td>1.69 ± 0.15</td><td>—</td></tr><tr><td>(4)</td><td>12.8 ± 1.0</td><td>100</td><td>1.59 ± 0.09</td><td>—</td></tr></tbody></table>

We first describe changes in rocking curves of horizontal Bragg spots (curves (2) and (4)) after a current of $I = 40$ A is applied with respect to zero current. The rocking curves show an increase in the integrated intensity $I_0$ by a factor of $\approx 1.3 \pm 0.1$. We observe a slight increase in the rocking width ($\Delta\sigma \approx (0.5 \pm 0.1)°$). We see a shift of the respective peak centers $x_0$ in negative direction ($\Delta x_0 \approx -(0.6 \pm 0.1)°$).

The rocking curves of vertical Bragg peaks (curves (1) and (3)) show a decrease in integrated intensity $I_0$ by a factor of $\approx 0.7 \pm 0.1$ with respect to zero current. The rocking curve width $\sigma$ increases ($\Delta\sigma \approx (1.6 \pm 0.1)°$). In contrast to the rocking curves of the horizontal Bragg peaks, we observe a shift of $x_0$ in positive direction ($\Delta x_0 \approx (0.7 \pm 0.3)°$).

Figure 7(e) shows simulated rocking curves (dotted lines) for different combinations of bulk and surface current and the

rocking scan of Bragg spot (1) with applied current of $I = 40$ A. The solid lines are the fits. All curves are scaled such that their integral is equivalent to the integral of the 40A fit. The rocking curves were calculated by approximating the magnetic self field according to (B2) and convolving the resulting rectangular-shaped angular distribution with the fitted Gaussian of the 40 A rocking curve of horizontal Bragg spot (2)<fnref n="11" />. The corresponding angular distribution is shown in the same color code at the bottom of the plot. For low bulk currents, the shape of the calculated rocking curve is dominated by the

<footnote n="11">¹¹ Using the width of the horizontal Bragg peak, which we assume is independent of magnetic self field broadening, includes the disorder due to the vortex movement.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 8. Average radial width of all first order Bragg peaks as a function of current *I* in a magnetic field of *B* = 500 G at *T* = 4 K. The inset shows the radial average of an exemplary single Bragg peak, here in horizontal direction, for *I* = 0 A and *I* = 40 A. The resolution limit of the instrument is colored in grey.</figcaption></figure>

Gaussian fit. For high bulk currents the shape is dominated by the rectangular angular distribution and we see the flat top shape as expected from a rocking curve corresponding to a continuously bent flux line.

When comparing the calculated rocking curves with the measured rocking curve in $\phi$ of Bragg peak (1) at $I$ = 40 A, the experimental curve is best approximated by the simulated curve with a bulk current contribution of $I_{\text{bulk}}$ = 4 A, which results in a surface current of $I_{\text{surf}}$ = 36 A.

## Appendix C. Radial width of Bragg peaks

We can generate radial averages in $|q|$ to extract the radial width $\sigma_q$ of our Bragg peaks. The radial width $\sigma_q$ is inversely related to the radial correlation length and is therefore a measure of the size of the well ordered mixed state domains. Figure 8 shows the average radial width of all first order Bragg peaks from SANS measurements FC in a magnetic field of $B_{\text{app}}$ = 500 G at $T$ = 4 K under liquid He as a function of current $I$. The inset shows the radial average of an exemplary Bragg peak, here in horizontal direction, for $I$ = 0 A and $I$ = 40 A. The resolution limit of the instrument is colored in grey. The radial width of each pair of Bragg spots at a given current $I$ was fitted individually. Within errors, there was no difference between the radial width of horizontal and vertical Bragg spots. Therefore the average radial width of all first order Bragg peaks is shown. We see a slight increase of the average radial width with increasing current starting at $I$ = 30 A. This relates to a decrease in radial correlation length and therefore to a shrinking of the well ordered mixed state domains, for both horizontal and vertical Bragg spots. From the radial average of the horizontal Bragg peaks shown in the inset we can additionally clearly see the increasing intensity in horizontal direction at low $q$ values for increasing current as also seen in the azimuthal averaging of the IMS scattering and the corresponding rocking scans (see figure 3 in the main text).

## Appendix D. Theoretical model

The study of the IMS in a low-$\kappa$ IT SC requires an approach beyond the GL theory. At the same time, solving a full set of the microscopic BCS equations with irregular IMS vortex configurations is prohibitively expensive computationally. Stationary IMS configurations have been recently obtained using the extended GL (EGL) theory [3]. However, the corresponding time-dependent EGL formalism for the non-stationary IMS is not yet available.

This computational problem can be circumvented using a relatively simple model of two GL equations for the components $\Psi_{\nu}(\mathbf{r})$ ($\nu$ = 1, 2) coupled by the linear Josephson-like terms. The free energy functional for this model reads as

$$f = \sum_{\nu=1,2} \left( \frac{1}{2m_{\nu}} |\mathbf{D}\Psi_{\nu}|^2 + \alpha_{\nu} |\Psi_{\nu}|^2 + \frac{\beta_{\nu}}{2} |\Psi_{\nu}|^4 \right) - \Gamma \{ \Psi_1^* \Psi_2 + \Psi_1 \Psi_2^* \} + \frac{\mathbf{B}^2}{8\pi}, \quad (\text{D1})$$

where $\mathbf{D} = -i\hbar\nabla - 2e\mathbf{A}/c$, and the coefficients $\alpha_{\nu}$, $\beta_{\nu}$, $m_{\nu}$ and the coupling constant $\Gamma$ are to be chosen. This model was originally derived for a two-band SC with the pairing in two distinct pockets of the Fermi surface (bands) [26, 27] like MgB$_2$. It was employed to investigate peculiar properties of the vortex matter attributed to the existence of two superconducting gaps—the so-called ‘type 1.5’ superconductivity [28].

However, qualitative characteristics of the vortex matter obtained using this two-component model [28–30] close to the critical temperature turn out very similar to those calculated theoretically [3] and observed experimentally [2] for conventional single-band BCS SCs in the IT regime. In particular, the two-component model reveals the non-monotonic inter-vortex interaction, attractive at large and repulsive at small distances [15, 31], which is a defining feature of the IT (type II/1) superconductivity in single-band materials [32, 45–48]. Similarly, both models, with single [4, 49] and two [15, 33] components, demonstrate a notably increased role of the many-vortex interactions that are key to stabilize large vortex clusters [3]. Furthermore, both single- [32, 45, 47] and two-band [15] models have topologically equivalent phase diagrams of vortex configurations.

This similarity is natural in the light of a long-known fact [27] that solutions for both components $\Psi_{1,2}$ of the two-component model coincide with the standard single-component GL theory in the limit of $T \to T_c$ (see also discussions [50, 51]). The similarity between the models is closely related to the very origin of the IT superconductivity—the self-duality and infinite degeneracy of the BCS theory at $T \to T_c$ and $\kappa \to \kappa_0$ (the Bogomolnyi point). Lifting the degeneracy makes this point unfold into a finite IT domain in the phase diagram [32]. A qualitative picture of the IT superconductivity does not depend on details of the mechanism of the degeneracy removal, in particular, on whether this occurs due to the beyond-GL-theory contributions or because the GL theory acquires the ‘weak’ auxiliary component. In this work we
exploit this qualitative similarity of the single- and two-band models and investigate a non-stationary current-driven IMS in single-band superconducting Nb using the approach based on the model in (D1) augmented with the time dependence.

The equations of the two-component model are obtained from the extremum condition for the energy functional (D1) to which the time derivatives are added to account for the dynamics. It is convenient to write these equations using a system of units defined by the zero-temperature coherence length, uniform solution for the gap and the critical field, calculated for one of the equations with zero Josephson coupling [15]. This yields

$$ \eta D_t \psi_1 = \mathbf{D}^2 \psi_1 - (\chi_1 - |\psi_1|^2) \psi_1 - \gamma \psi_2, \quad (\text{D2a}) $$

$$ \eta D_t \psi_2 = \frac{1}{\alpha} \mathbf{D}^2 \psi_2 - (\chi_2 - |\psi_2|^2) \psi_2 - \frac{\beta_2}{\beta_1} \gamma \psi_1, \quad (\text{D2b}) $$

$$ \kappa_1^2 \nabla \times \nabla \times \mathbf{A} = \frac{1}{\kappa_1^2} \Re[\psi_1 \mathbf{D} \psi_1^*] + \frac{\alpha}{\kappa_2^2} \Re[\psi_2 \mathbf{D} \psi_2^*], \quad (\text{D2c}) $$

where $\mathbf{D} = -i\nabla - \mathbf{A}$ and $D_t = \partial_t - i\phi$ are the scaled gauge-invariant derivatives with $\mathbf{A}$ and $\phi$ being the vector and scalar potentials, respectively, $\gamma$ denotes the (scaled) interband coupling constant, $\kappa_{1,2}$ are the GL parameters for the components $\Psi_{1,2}$, calculated for each of the component (D2a) and (D2b) separately at $\Gamma = 0$, and $\alpha = \xi_2^2/\xi_1^2$ where $\xi_{1,2}$ are the GL coherence lengths calculated separately for both components. The GL parameters are related as $\kappa_2 = \alpha\kappa_1 \sqrt{\beta_1/\beta_2}$. Notice, that the GL parameter $\kappa$ for the entire system (at $T \to T_c$) differs from both $\kappa_1$ and $\kappa_2$, its expression can be found in [50]. Parameters $\chi_\nu = \tau - S_\nu$ with $\tau = 1 - T/T_c$ define the temperature dependence. Here constants $S_\nu$ are determined by the intraband coupling. They satisfy the condition $S_1 S_2 \alpha^2 = \gamma^2 \beta_1/\beta_2$ which ensures that the superconducting order parameter disappears at $\tau \to 0$, so that $T_c$ is the superconductivity transition temperature (one can find a detailed derivation of all parameters for a two-band model in [15]). Finally, $\eta$ is related to the system losses and determines the time scale. The influence of the normal current on the evolution of the IMS configurations is neglected.

In the calculations we take $S_1 = 0.043$, $S_2 = 0.188$, $\beta_1/\beta_2 = 0.92$, $\kappa_1 = 1.5$. We also assume $\alpha^2 = 0.65$ and $T = 0.74T_c$ which ensures that the system is in the IT domain, very close to the line of zero surface tension of the N-S domain wall, where the IT superconductivity is expected. These values are taken from the earlier work, where the equivalence of the two- and one-component models was established [15]. However, qualitative features of the IMS dynamics remain qualitatively similar in a wide range of parameters as long as the system is in the IT regime.

We consider a geometry where vortices are directed along the z axis and the corresponding magnetic field is then $\mathbf{B} = (0,0,B(x,y))$. We also assume that in the x-direction the sample has the finite width $a$ and that the current flows between the sample surfaces $x=0, a$. The length in the y-direction is taken much larger, $L \gg a$. A complete description of the current-induced vortex evolution requires solving (D2) for a

realistic finite-size sample together with the equation for the charge density distribution in the sample as well as in the contacts. In this work we follow a simplified approach where the current is 'injected' into the SC by applying a difference of the phase/field at the sample boundaries. Notice, that deep inside a sample the current is absent in the Meissner phase and therefore isolated vortex clusters do not move. However, a current flowing at the surface interacts with vortices creating the Lorentz force (figure 7(b)). We model this situation by injecting a small current by imposing the external linear potential $\phi(x,y) = gx$ which creates a phase difference between the boundaries $x=0$ and $x=a$, where $g$ controls the current value. This model neglects the current decay inside the sample and thus does not describe the vortex bending, discussed in the main text. It is nevertheless sufficient to capture general features of the evolution of the IMS superstructures.

Equation (D2) are solved with the SC-metal boundary conditions $\mathbf{D}_\perp \psi_j = -i\psi/b$ applied at $x=0$ and $x=a$. A particular value of the real parameter $b$ is not important, it is set to $b=20$ in the calculations$^{12}$. In the y-direction we assume the periodic boundary conditions. We note, that details of the model geometry and the boundary conditions have only little influence on the dynamics of the vortex matter far from the sample boundaries. (D2) are solved on a two-dimensional grid with the spacing $a_x = a_y = 0.25 \xi_1$, which is sufficient to describe vortices. In the calculations we first obtain a stationary vortex configuration by solving (D2) with $g=0$. Then, the potential is switched on ($g \neq 0$), and the evolution of the IMS begins. We allow the vortex matter to evolve until it achieves a quasi-stationary configuration (stripes).

## Appendix E. Details of the vortex cluster evolution

Results of the calculations are shown in figure 4 in the main text and in figure 9. When an IT SC is placed in the magnetic field and the current is absent, vortices inside the sample form the IMS of randomized vortex clusters (figure 4(a)). When the current flows the IMS changes, eventually forming stripes in the direction perpendicular to the current. We note that the evolution is qualitatively similar for all initial IMS configurations. Figures 9(a)–(l) illustrate main stages of the evolution by showing snapshots of the current density for a relatively large vortex cluster, where each row of the figure highlights a specific feature in the elongation process. To demonstrate the time evolution of the clusters shape, the observation frame in figures 4(a)–(h) moves downward and keep the cluster in focus so that one observes the relative motion of different clusters, all moving downwards (in the last row of figures 4(i)–(l) the frame does not move).

By looking at figure 9 we first note that the current distribution in voids of the Meissner phase inside the cluster forces them to shrink and disappear (cf figures 9(a) and (b)). The second important feature to note is that the total current-induced force increases with the boundary length such

<footnote n="12">12 An example of the effect of changing $b$ on the magnetization of a superconducting sample can be found for example in [52].</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The spatial dependence of the current density calculated for a vortex cluster at chosen time instants illustrating the time evolution (the time direction is shown by an arrow marked *t*), panel (a) is the initial state and panel (l) is the final. Current *I* flows horizontally rightwards, vortices move vertically downwards. In panels (a)–(h) the observation frame shifts downwards to keep the largest cluster in focus to demonstrate the time evolution of its shape and relative positions of neighboring smaller clusters. Areas of the larger current (red) correspond to a larger pulling force acting normal to the cluster boundary. Each panel row illustrates a different stage of the evolution, which produces vertical vortex stripes in the end.</figcaption></figure>

that larger clusters (or separate parts of a cluster) move faster. This is clearly seen in figures 9(a)–(d) where smaller semi-isolated structures at the rear move slower than the main cluster part and break off eventually. The evolution shown in figures 9(e)–(h) illustrates the same trend from the opposite perspective: here a large fast-moving cluster catches up with and then absorbs the smaller ones.

While moving, the cluster gradually elongates in the vertical direction which takes place due to a combined action of the Lorentz force, that pulls the cluster on its lower boundary, and the drag force, that slows down vortices at its rear. As a result, the cluster acquires the shape similar to that of a liquid droplet falling in the air. Vortices at the rear part of the cluster

are only weakly coupled to its main body and are eventually broken off and left behind. This process makes the neighboring clusters merge, forming a stripe in the direction of their motion (figures 9(i)–(l)).

The evolution of vortex clusters can be intuitively understood as a combination of three main factors: the Lorentz force, the drag force and the inter-vortex interaction potential. These are typical for the superconducting mixed state, however, their action on the IMS in IT SCs has a number of specific features. In a conventional type II SC the distance between vortex cores is much larger than their size so that each vortex can be considered separately with respect to the acting Lorentz and drag forces. Therefore the current-induced dynamics is practically
the same for all vortices (although even for type II SCs the vortex motion changes the inter-vortex interaction making it dependent on the vortex mutual position [36]).

In contrast, the size of the vortex core in IT SCs is comparable to the inter-vortex distance determined by the minimum of the non-monotonic vortex-vortex interaction potential [32]. In addition, many-vortex interactions play an increasingly important role for large vortex clusters [15]. As a consequence, vortices can no longer be regarded as separate ‘elementary particles’ of the mixed state. The properties of a vortex cluster therefore do not depend simply on its number of vortices but also on the cluster configuration—shape and size. The applied current is strongly distorted by the vortex configuration flowing mainly along cluster boundaries (see figure 9). Thus the Lorentz force acts mainly on the boundary pulling it in the normal direction. In addition the current profile is notably asymmetric, so that its density is larger at the down-facing boundaries and consequently leads to a larger net force in the downward direction. The drag force acting on a vortex inside a cluster is smaller than that for a separate vortex because the movement inside a cluster involves smaller changes in the field-condensate profile due to the comparable size of the vortex core to the intervortex distance. Then the drag force is largest for the boundary vortices (we note that our model does not take into account all mechanisms leading to the vortex drag [53, 54], however, this does not change qualitative conclusions). The dependence of the Lorentz and drag forces on the cluster configuration and size gives rise to the velocity dispersion for different clusters and cluster parts that is clearly visible in the numerical calculations in figure 9.

Finally, the unique vortex interactions in an IT material ensure the preferred mean intervortex distance, but not the vortex cluster shape. The Lorentz force, that pulls the cluster boundary downwards, the drag, that acts on vortices in the opposite direction and a ‘soft’ inter-vortex interaction result in elongated clusters in the direction of the movement, i.e. perpendicular to the current. The elongation increases the bypassing current flow and thus the Lorentz force at the boundary which, in turn, leads to a still faster elongation.

We conclude this discussion by noting that additional factors such as crystal anisotropy and pinning are not expected to change the superstructure formation mechanism qualitatively. For moving vortices, the static potential due to the pinning can be represented as a random force that does not change their averaged movement, which underlies the formation of superstructures. This argument applies on the large scales, for superstructure sizes much larger than the characteristic lengths associated with the pinning potential. It fails at short scales which define vortex arrangements inside those superstructures. A similar argument applies for the crystal anisotropy which affects only the short-ranged correlations between the vortex positions and not the shape and orientation of larger vortex clusters. The model must be amended only if the short range correlations of the vortex matter are of interest. Analysis of the role of those additional factors is therefore beyond the scope of this calculation, which aims to explain the mechanism for emerging superstructures in IT SCs, like Nb.

## ORCID iDs

Xaver S Brems https://orcid.org/0000-0002-0053-9116

Sebastian Mühlbauer https://orcid.org/0000-0002-3106-5785

Wilmer Y Córdoba-Camacho https://orcid.org/0000-0001-5813-1749

Arkady A Shanenko https://orcid.org/0000-0002-6031-5106

Alexei Vagov https://orcid.org/0000-0001-7446-7728

Robert Cubitt https://orcid.org/0000-0003-1778-3307

## References

[1] Brandt E H 1995 The flux-line lattice in superconductors *Rep. Prog. Phys.* **58** 1465

[2] Brandt E H and Das M P 2011 Attractive vortex interaction and the intermediate-mixed state of superconductors *J. Supercond. Nov. Magn.* **24** 57

[3] Vagov A, Wolf S, Croitoru M D and Shanenko A A 2020 Universal flux patterns and their interchange in superconductors between types I and II *Commun. Phys.* **3** 58

[4] Wolf S, Vagov A, Shanenko A A, Axt V M and Albino Aguiar J 2017 Vortex matter stabilized by many-body interactions *Phys. Rev. B* **96** 144515

[5] Ge J-Y, Gutierrez J, Lyashchenko A, Filipov V, Li J and Moshchalkov V V 2014 Direct visualization of vortex pattern transition in ZrB$_{12}$ with Ginzburg-Landau parameter close to the dual point *Phys. Rev. B* **90** 184511

[6] Backs A, Schulz M, Pipich V, Kleinhans M, Böni P and Mühlbauer S 2019 Universal behavior of the intermediate mixed state domain formation in superconducting niobium *Phys. Rev. B* **100** 064503

[7] Reimann T *et al* 2017 Domain formation in the type-II/1 superconductor niobium: interplay of pinning, geometry and attractive vortex-vortex interaction *Phys. Rev. B* **96** 144506

[8] Reimann T, Mühlbauer S, Schulz M, Betz B, Kaestner A, Pipich V, Böni P and Grünzweig C 2015 Visualizing the morphology of vortex lattice domains in a bulk type-II superconductor *Nat. Commun.* **6** 8813

[9] Mühlbauer S, Pfleiderer C, Böni P, Laver M, Forgan E M, Fort D, Keiderling U and Behr G 2009 Morphology of the superconducting vortex lattice in ultrapure niobium *Phys. Rev. Lett.* **102** 136408

[10] Laver M *et al* 2006 Spontaneous symmetry-breaking vortex lattice transitions in pure niobium *Phys. Rev. Lett.* **96** 167002

[11] Träuble H and Essmann U 1967 Die beobachtung magnetischer strukturen von supraleitern zweiter art *Physica Status Solidi* b **20** 95

[12] Hoberg J R and Prozorov R 2008 Current-driven transformations of the intermediate-state patterns in type-I superconductors *Phys. Rev. B* **78** 104511

[13] Ooi S, Tachiki M, Konomi T, Kubo T, Kikuchi A, Arisawa S, Ito H and Umemori K 2021 Observation of intermediate mixed state in high-purity cavity-grade Nb by magneto-optical imaging *Phys. Rev. B* **104** 064504

[14] Wang W, R Díaz-Méndez, Wallin M, Lidmar J and Babaev E 2021 Pinning effects in a two-dimensional cluster glass *Phys. Rev. B* **104** 144206

[15] da Silva R M, Milošević M V M V, Shanenko A A, Peeters F M and Albino Aguiar J 2015 Giant paramagnetic Meissner effect in multiband superconductors *Sci. Rep.* **5** 12695

[16] Cross M C and Hohenberg P C 1993 Pattern formation outside of equilibrium *Rev. Mod. Phys.* **65** 854
[17] Rabinovich M I, Ezersky A B and Weidman P D 2000 *The Dynamics of Patterns* (Singapore: World Scientific) (https://doi.org/10.1142/4207)

[18] Pismen L M 2006 *Patterns and Interfaces in Dissipative Dynamics* (Berlin: Springer) (https://doi.org/10.1007/3-540-30431-2)

[19] Hoyle R B 2006 *Pattern Formation* (New York: Cambridge University Press) (https://doi.org/10.1017/CBO 9780511616051)

[20] Dewhurst C D, Grillo I, Honecker D, Bonnaud M, Jacques M, Amrouni C, Perillo-Marcone A, Manzin G and Cubitt R 2016 The small-angle neutron scattering instrument D33 at the Institut Laue-Langevin *J. Appl. Crystallogr.* **49** 1

[21] Jones R and Rose-Innes A 1966 Distribution of current in a type-2 superconductor *Phys. Lett.* **22** 271

[22] Hocquet T, Mathieu P and Simon Y 1992 Continuum theory of the mixed-state and surface Joule effects in type-II superconductors *Phys. Rev.* B **46** 1061

[23] Pautrat A, Goupil C, Simon C, Charalambous D, Forgan E M, Lazard G, Mathieu P and Brûlet A 2003 Distribution of transport current in a type-II superconductor studied by small-angle neutron scattering *Phys. Rev. Lett.* **90** 087002

[24] Goupil C et al 2000 Small angle neutron scattering and vortex lattice dynamical phase diagram *Physica C* **341–348** 999

[25] Bean C P 1964 Magnetization of high-field superconductors *Rev. Mod. Phys.* **36** 31

[26] Tilley D R 1964 The Ginzburg-Landau equations for pure two band superconductors *Proc. Phys. Soc.* **84** 573

[27] Geilikman B T, Zaitsev R O and Kresin V Z 1967 Properties of superconductors having overlapping bands *Sov. Phys. Solid State* **9** 642

[28] Babaev E, Carlström J, Silaev M and Speight J 2017 Type-1.5 superconductivity *Superconductors at the Nanoscale* ed R Wördenweber *et al* (Berlin: De Gruyter) pp 133–64

[29] Babaev E and Silaev M 2012 Type-1.5 superconductivity in multiband and other multicomponent systems *J. Supercond.* *Nov. Magn.* **26** 2045–55

[30] Biswas P K et al 2020 Coexistence of type-I and type-II superconductivity signatures in ZrB$_{12}$ probed by muon spin rotation measurements *Phys. Rev.* B **102** 144523

[31] Babaev E, Carlström J and Speight M 2010 Type-1.5 superconducting state from an intrinsic proximity effect in two-band superconductors *Phys. Rev. Lett.* **105** 067003

[32] Vagov A, Shanenko A A, Milošević M V, Axt V M, Vinokur V M, Albino Aguiar J and Peeters F M 2016 Superconductivity between standard types: multiband versus single-band materials *Phys. Rev.* B **93** 174503

[33] Carlström J, Garaud J and Babaev E 2011 Length scales, collective modes and type-1.5 regimes in three-band superconductors *Phys. Rev.* B **84** 134518

[34] Dobrovolskiy O V, Vodolazov D Y, Porrati F, Sachser R, Bevz V M, Mikhailov M Y, Chumak A V and Huth M 2020 Ultra-fast vortex motion in a direct-write Nb-C superconductor *Nat. Commun.* **11** 3291

[35] Embon L et al 2017 Imaging of super-fast dynamics and flow instabilities of superconducting vortices *Nat. Commun.* **8** 85

[36] Kogan V G and Prozorov R 2020 Interaction between moving Abrikosov vortices in type-II superconductors *Phys. Rev.* B **102** 024506

[37] Scher H and Zallen R 1970 Critical density in percolation processes *J. Chem. Phys.* **53** 3759

[38] Cubitt R and Brems X 2020 The effects of current on the intermediate mixed state in niobium Institut Laue-Langevin (ILL) (https://doi.org/10.5291/ILL-DATA.EASY-568)

[39] Brems X, Cubitt R, Mühlbauer S and Oji U 2020 Self organisation of mixed superconducting states under the application of a current Institut Laue-Langevin (ILL) (https://doi.org/10.5291/ILL-DATA.5-31-2911)

[40] Huebener R P, Kampwirth R T and Seher A 1970 Electrical resistivity and peak effect in superconducting niobium *J. Low Temp. Phys.* **2** 113

[41] Okuma S, Motohashi A and Kawamura Y 2013 Scaling analysis on plastic depinning of vortex matter at low temperature in an amorphous $Mo_xGe_{1-x}$ film *Supercond. Sci. Technol.* **26** 025013

[42] Bardeen J and Stephen M J 1965 Theory of the motion of vortices in superconductors *Phys. Rev.* **140** A1197

[43] Kroeger D M and Schelten J 1976 Bending of flux lines by transport currents in type ii superconductors measured by neutron diffraction *J. Low Temp. Phys.* **25** 369

[44] Cubitt R, Campbell A S, Forgan E M, Dewhurst C D and Yang G 2009 Investigation of vortex structures in a current-carrying Nb wire *Supercond. Sci. Technol.* **22** 034014

[45] Jacobs A E 1971 Interaction of vortices in type-II superconductors near $T = T_c$ *Phys. Rev.* B **4** 3029

[46] Klein U 1987 Microscopic calculations on the vortex state of type II superconductors *J. Low Temp. Phys.* **69** 1

[47] Miranović P and Machida K 2003 Thermodynamics and magnetic field profiles in low-$\kappa$ type-II superconductors *Phys. Rev.* B **67** 092506

[48] Luk'yanchuk I 2001 Theory of superconductors with $\kappa$ close to $1/\sqrt{2}$ *Phys. Rev.* B **63** 174504

[49] Cordoba-Camacho W Y, Vagov A, Shanenko A A, Albino Aguiar J, Vasenko A S and Stolyarov V S 2021 Vortex interactions and clustering in thin super conductors *J. Phys. Chem. Lett.* **12** 4172

[50] Kogan V G and Schmalian J 2011 Ginzburg-Landau theory of two-band superconductors: absence of type-1.5 superconductivity *Phys. Rev.* B **83** 054515

[51] Shanenko A A, Milošević M V, Peeters F M and Vagov A V 2011 Extended Ginzburg-Landau formalism for two-band superconductors, *Phys. Rev. Lett.* **106** 047005

[52] Barba-Ortega J, Sardella E and Albino Aguiar J 2010 Superconducting boundary conditions for mesoscopic circular samples *Supercond. Sci. Technol.* **24** 015001

[53] Kopnin N B 1976 Forces acting on vortices moving in a pure type II superconductor *Sov. Phys. - JETP* **44** 7

[54] Blatter G, Feigel'man M V, Geshkenbein V B, Larkin A I and Vinokur V M 1994 Vortices in high-temperature superconductors *Rev. Mod. Phys.* **66** 1125