Dynamic metastable vortex states in interacting vortex lines

The electron transport in current-biased superconducting nano-bridges is determined by the motion of the quantum vortex confined in the internal disorder landscape. Here we consider a simple case of a single or two neighbouring linear defects crossing a nano-bridge. The strong anharmonicity of the vortex motion along the defect leads, upon RF-excitation, to fractional Shapiro steps. In the case of two defects, the vortex motion becomes correlated, characterized by metastable states that can be locked to a resonant RF-drive. The lock-unlock process causes sudden voltage jumps and drops in the voltage-current characteristics observed in experiments. We analyze the parameters promoting these metastable dynamic states and discuss their potential applications in quantum devices.


I. INTRODUCTION
Quantum vortices are famous topological objectslines of 2π-phase singularities in the many-body wave function of coherent quantum condensates.In superconductors, where condensed particles are electrically charged Cooper pairs, the phase gradients generate vortex currents circulating around singularities, and producing a magnetic flux.The vortices strongly influence the characteristics of superconductors, limiting their critical currents and fields.In fact, externally applied currents and fields interact with the vortex, forcing it to move.In the vortex centers -cores -the superconductivity is suppressed, and the normal state is recovered.The vortex motion is therefore dissipative, often triggering the transition to the normal state of the entire system.
In their motion inside superconductors, vortices interact with various local and extended defects, as well as with other vortices and obstacles.The collection of defects along with other moving and pinned vortices form a potential landscape in which a given vortex evolves.This landscape is generally dynamic and intricate, comprising local minima and saddle points.Consequently, the formation of various metastable states can occur, with their characteristic energies and stability subject to perturbation by external magnetic fields, DC or AC currents.
In this work, we focus on the over-critical behaviour of current-biased superconducting nano-bridges (see Fig. 1a as an example) aiming to explain the spectacular results of recent experimental findings.In most of the cases, a nano-bridge behaves like a Josephson junction: when the critical current I c is reached, the bridge transits to the normal state.However, this transition is not always abrupt: the voltage V (I > I c ) across the bridge increases progressively, often over several orders of magnitude, before the device reaches a fully normal state.In this progressive transition, the differential resistance dV /dI(I) The model system is a rectangular superconducting sample representing the central part of a typical nanobridge, as shown in Fig. 1.The model bridge has a length of L = 60ξ and a width of W = 40ξ, where ξ represents the Ginzburg-Landau (GL) coherence length.Within the TDGL framework, the temporal and spatial evolution of the complex superconducting order parameter ψ(t, r) can be expressed as [11] (see Methods): where A is the vector potential due to magnetic field, J S and J N are superconducting and normal components of the electric current.The GL parameter κ = λ/ξ (λ is the field penetration depth) is taken equal to κ = 4, meaning that the bridge is in the type-II regime (see Methods).
The defects, as depicted in Fig. 1b,c, are introduced by spatially varying the parameter ǫ(r), which is associated with the local critical temperature T c (r) and the global sample temperature T : The superconducting part of the bridge is described by ǫ=1, while the defects are characterized by a locally reduced critical temperature T c (r), and are described by a lower ǫ(r).For instance, the linear (grey) defect of width 1×ξ crossing the bridge is characterized by ǫ = 0.5.
It represents an extended structural defect -a grain boundary crossing the real sample or an artificial weaklink (possible experimental realizations are discussed in Sec.III).At a given temperature T , this defect is superconducting, but its local critical temperature is 3/4 of the critical temperature T c in the rest of the sample.The two point defects at the edges, presented by black rectangles 2ξ × 5ξ, are characterized by ǫ = 0, that corresponds to a fully suppressed superconductivity.These edge defects appear at the ends of the grain boundaries as a result of a damage caused during the nano-bridge fabrication processes.
When simulating the Shapiro step experiments, the microwave illumination is added as an AC-current of amplitude I AC and frequency f AC .The total transport current through the bridge is therefore: The state of the bridge is determined by calculating the voltage V between y = 0 and y = L boundaries for each value of transport current (see Methods).By averaging this voltage over the sample width and time one gets the DC-voltage V measured in experiments.

B. Single linear defect
As a starting point, we consider a single linear defect, as shown in Fig. 1b, that simulates a grain boundary crossing the bridge.Additionally, two point defects are introduced at the ends of the linear defect, representing suppressed superconductivity in the locations where the grain boundary reaches the sample edges.
The results of calculations are presented in Fig. 2. When the transport current I tr is well below a critical value I c , the order parameter in the bridge is steady.It is depleted at the two local edge defects and at the linear defect, Fig. 2a, following the imposed ǫ(r).At I tr I c , there is already one vortex and one anti-vortex pinned at the two edge defects.The entire bridge remains in the superconducting state, with V =0, as expected.Figs.2c-e are snapshots of the temporal evolution of the order parameter amplitude when a constant I tr = 0.10 > I c is applied.Under this condition, one vortex and one antivortex simultaneously enter the bridge, Fig. 2c.They accelerate towards each other under the action of the Lorenz force and experience mutual attraction, Fig. 2d, and annihilate, Fig. 2e.The process is periodic, with the period and details of the vortex-antivortex dynamics depending on the TDGL parameters.Moving vortices dissipate energy and generate an instantaneous voltage V (t) proportional to the relative vortex velocity.Fig. 2f illustrates the evolution of V (t), with points c-e correspond to snapshots in Figs.2c-e.At the moment (c), the voltage FIG. 3. Transport properties of the nano-bridge with one linear defect.Solid lines -normalized V (I DC ) characteristics calculated for different values of the AC component I AC of the total transport current.The right vertical axis displays the numbers of Shapiro plateaus.Dashed line -V (I) characteristic of a SNS Josephson junction calculated within the RSJ model [13,14].To fit the curve into the plot window, its vertical scale was divided by a factor of 10.
rapidly rises as vortices accelerate due to their interaction with the edges and the transport current.In (d), V crosses a local minimum as the vortex velocity drops in a region where interaction with the edge is already sufficiently small, and the transport current is reduced on the scale of ∼ λ.A sharp increase in the vortex velocity due to the vortex-antivortex attraction just before annihilation produces a peak in V (t) at moment (e).The Fourier spectrum of V (t) is presented in Fig. 2g.It contains the fundamental frequency of an amplitude V 1 and several harmonics with comparable amplitudes.Both V (t) and spectrum indicate the strong anharmonicity of the vortex motion.It is important to note that the fundamental frequency is not fixed but grows with I [12] as the increasing Lorentz force pushes vortices to move and annihilate faster.
By repeating the calculations for different I DC and time-averaging V (t) , we retrieve V (I DC ) dependencies measured in experiments.The dark green curve in Fig. 3 is the result of these calculations for the range of I DC around the transition from the non-dissipative to the dissipative state.The plots are presented in reduced coordinates V /µ 0 vs I DC /(J 0 W ) (see Methods).The shape of the curve resembles the V (I) characteristics of an ordinary Superconductor -Normal metal -Superconductor (SNS) Josephson junction.The latter, represented by the dashed line in Fig. 3, was calculated using the Restively Shunted Junction (RSJ) model.Both curves exhibit a non-dissipative branch at low currents, a rise at some critical current, and a smooth increase at higher currents.However, the resemblance is lim-ited.First, in the present case, the SNS junction does not exist, and the intrinsic critical (depairing) current of the bridge ∼ J 0 W is much higher than the calculated value I c ≃ 0.072J 0 W . Second, in SNS junctions, the V (I DC ) curve asymptotically approaches the normal branch V = R N I DC as I DC increases, while the bridge "resistance" V /I DC × (J 0 W/µ 0 ) remains much lower than its normal state resistance R N (this is why the vertical scale of the RSJ curve was reduced to fit it in the plot window).Third, the RSJ model fails in reproducing an almost linear rise of V (I DC ) for I > I c .The three deviations originate from the fact that, contrary to SNS junctions in which the voltage appears as a result of the suppression of the proximity-induced superconducting correlations in the N-part, in the nano-bridges it is due to individual vortex motion inside a still superconducting device.This difference is essential, leading to unique transport properties that we focus on in this work.
The SNS-like behaviour of the bridge is further evidenced by simulating its response to microwave illumination.When an AC-current is added, the oscillating voltage V (t) can be locked to the frequency f AC of this external drive, resulting in plateaus of constant voltage on V (I DC ) curve, as shown in Fig. 3.This effect resembles the well-known Shapiro steps observed in ordinary Josephson junctions under microwave illumination, where the n th voltage plateau is defined by the locking condition f 1 = n • f AC (where n is an integer) and the second Josephson relation V (t) = hf 1 /2e, where f 1 is the fundamental (Josephson) frequency.
In addition to the integer Shapiro plateaus, the fractional ones are also revealed.These plateaus appear in Fig. 3 at voltages satisfying the condition n k • hf AC = 2e V (t) , where n and k are integers.The presence of fractional plateaus is directly linked to a high anharmonicity of V (t) oscillations in Fig. 2f.The Fourier spectrum of V (t), presented in Fig. 2g, is indeed characterized by high amplitudes V k of k th voltage harmonics (at a frequency f k ) which are comparable to the amplitude V 1 at its fundamental frequency f 1 .This enables an efficient locking of these harmonics to the AC-drive when In Fig. 4, the instantaneous voltage V (t) and its spectrum for the Shapiro plateau 1/3 are shown.When the frequency f AC is slightly detuned from f k /n, a lowfrequency envelope of a beat frequency |f AC − f k /n| is observed.This effect allows for the detection of higher harmonics experimentally, even when their magnitude is small and imperceptible in V (I DC ) curves in the Shapiro step experiment.Consequently, by tunning the amplitude and frequency of the AC excitation, it becomes possible to induce fractional Shapiro step when f AC is not only a multiple of f 1 but a multiple of higher f k harmonics.All these features reveal a rich spectral characteristic of the considered system.It should be mentioned that when the amplitude I AC becomes comparable to I DC , the AC-excitation cannot be considered as a perturbation anymore.Instead, one should think of a complex dynamical system whose spectrum (amplitudes V k and frequencies f k ) depends on both components of I tr .
Concluding this section, it is important to recall pioneering works [15,16] that predicted similarities between the transport properties of Josephson junctions and those of nano-bridges crossed by vortices (or phaseslips).These predictions were later confirmed in several experiments, where both integer [5,6,8,17] and fractional [4] Shapiro steps were observed.This analogy was also explored, both experimentally and theoretically, in the case of vortices jumping between pinning sites [12,[18][19][20][21].

C. Two neighbouring linear defects
In experimentally studied nano-bridges, the disorder is rarely represented by only one grain boundary.Most non-epitaxial superconducting films exhibit granularity on a scale of 20-200 nm, which can be significantly shorter than the nano-bridge width W .For instance, nanomeanders studied in [22] were elaborated out of thin YBa 2 Cu 3 O 7−δ films.They possess a specific morphology [23] and form a network of grain boundaries.Statistically, several such boundaries can cross the bridge.The vortex motion in these networks is much more complex than in a single grain boundary studied above.While the vortex cores are confined within the grain boundaries [9], the vortex currents extend far beyond; they circulate on a scale of λ or, in ultrathin films, on an even larger scale of the Pearl penetration depth [24].This results in a mutual interaction between vortices present in different grain boundaries, affecting their collective motion.As a step towards accounting for this complexity, we consider now two linear defects (grain boundaries) characterized by ǫ=0.5.The defects are separated by a distance l = 5ξ ∼ λ, Fig. 1c, thus inducing an interline vortex-vortex interaction.
The calculated V (I DC ) characteristics in the case of two identical linear defects is presented as a dark green solid line in Fig. 5.The shape of this curve is almost identical to that obtained in the case of a single linear defect.As in the previous case, at DC-currents just above the critical one I DC I c , one vortex-antivortex pair enters the nano-bridge and moves along one of the two linear defects, thus generating a non-zero voltage.The only difference with the single defect case is that after vortex-antivortex annihilation in one line, a new vortexantivortex pair enters the other line, and the process repeats.Further increase of the DC-current leads to an acceleration of vortices and, consequently, to an increase of voltage V .At a high enough DC-current, the system enters a new state in which the second vortex-antivortex pair enters into the second line before the first pair annihilates in the first one.This moment is witnessed by a slight inflexion of V (I DC ) curve at I DC ≃0.084.In this state, there are two vortex-antivortex pairs in the nanobridge at the same time.Due to the mutual repulsion of vortices of the same sign, they try to position themselves as far from each other as possible, while remaining inside linear defects.This leads to a lateral x-shift of the vortex positions in neighbouring lines, as shown in Fig. 6a.This dynamic vortex pattern is reminiscent of the static Abrikosov vortex lattice.As time advances, the vortexantivortex pair in the bottom line annihilates, the one in the top line advances towards the center, and a new one enters the bottom line.
By adding a low AC-current, one gets the Shapiro plateaus that also look very similar to the single defect case (brown line in Fig. 5).Though, at higher AC-currents new features appear.These are large 2/3 Shapiro plateaus on V (I DC ) curves with a rapid voltage raise on their left side and a voltage drop on their right side (hand-added smooth dashed lines help to appreciate the amplitude of the effect).Unlike other plateaus, the width of the 2/3 plateau rapidly grows with the ACcurrent (compare the curves at I AC =0.02, 0.03, and 0.05).
To understand the origin of this phenomenon, let us consider the dynamics of the system close to the voltage drops.In the specific case of I AC = 0.03, this occurs at I drop DC ≃ 0.076, as indicated by the arrow on the V (I DC ) curve in Fig. 5.The calculations show that just above I drop DC , the vortex-antivortex motion in the two defects is sequential, as presented in Fig. 6a, while just below I drop DC (that is on the plateau) it is synchronous: Vortexantivortex pairs enter the defects simultaneously, move in parallel to each other (see the snapshot Fig. 6b), and annihilate at the same time.This leads to high peak-topeak voltage spikes in V (t), as those visible on the left side of Fig. 6c.
The synchronous configuration is not stable itself.Indeed, when a vortex in one line is located under a vortex in the other, the projection of a vortex-vortex repulsion force on a x-axis is zero, and any x-shift of their position gives rise to the x-axis component of vortex-vortex repulsion which drives the system out of this unstable balance towards a more stable checkerboard configuration, Fig. 6a.Thus, the metastable configuration of Fig. 6b is stabilized by the external AC-drive that works as a periodic force; if the amplitude of this force (proportional to I AC ) is sufficient, the configuration is stabilized, in some range of external parameters, giving rise to a plateau on V (I DC ) curve.When the DC-current is slightly increased above I drop DC , the Lorenz force increases, and the system jumps down to the stable configuration of Fig. 6a.The corresponding evolution of V (t) is presented in Fig. 6c.One can observe that after a few periods of high peak-to-peak voltage oscillations, the system transits to oscillations with a nearly twice lower peak-to-peak voltage (compare left and right parts of Fig. 6c).This change is due to the fact that in the metastable state, the vortex-antivortex annihilation takes place simultaneously in the two lines, while in the stable configuration the process is sequential.The system is no more locked to the 2/3 Shapiro step in the stable configuration.A movie illustrating the oscillatory dynamics of this transition is provided in the Supplementary Material.The motion of vortices in the two close linear defects can be seen as a system of two coupled identical anharmonic oscillators.In this representation, the two oscillation patterns of Fig. 6 can be seen as two modes, one of which is low in energy (E 0 ) and therefore stable, while the other, at higher energy E 1 , is metastable.Each of these modes depends on DC current I DC , and the evolution of the lowest mode corresponds to V (I DC ) curve at I AC = 0. Another mode can only be achieved with an external excitation, in a certain range of pumping powers and frequencies.
The calculated Fourier spectra of V (t) in the states E 0 and E 1 are presented in Fig. 7.In the metastable configuration E 1 , the AC-drive locks to the third harmonic of the system as 2f AC = 3f 1 .The Josephson frequency is f 1 = (2/3)f AC and, consequently, the DC-voltage measured in the experiment is V (t) = (2/3)hf AC /2e.This voltage remains constant as long as the system is locked to the drive, resulting in the unusual 2/3 Shapiro plateau in Fig. 5. Immediately after the drop, the drive locks to the second harmonic as In principle, it could be the usual 1/2 Shapiro plateau, due to anharmonicity.Though, when the AC-current increases and the width of the unusual 2/3 plateau rapidly grows, the plateau 1/2 shrinks and disappears (compare the curves at I AC = 0.02, 0.03 and 0.05 in Fig. 5).Note that as I DC is further increased above I drop DC , the Lorentz force rises pushing vortices to move faster, the corresponding frequencies grow, and the lock to the fixed frequency of the AC-drive is lost.This roller-coaster ride between different metastable, stable locked and unlocked states is reflected in voltage spectra and as a consequence in a V (I DC ) curve.
Till now we have considered a very idealistic case where the two coupled linear defects were identical.This situation could be realized in artificial stacks of SNS junctions [25], periodic pinning arrays [20,26] but not in nano-bridges made of films in which the intrinsic pining landscape is aperiodic and the inter-grain coupling varies from one grain boundary to the other.To account for this diversity, we also studied asymmetric linear defects.In Fig. 8, we show the V (I DC ) characteristics for the case of two linear defects located as in Fig. 1c, but characterized by different ǫ parameters: ǫ = 0.5 and ǫ = 0.42.The curves differ significantly from the previous case, even without AC-excitation (green curve).The critical current is lower, and for 0.074 < I DC < 0.0805, the voltage appears exclusively due to the vortex motion in the ǫ = 0.42 line; the ǫ = 0.5 line contains no vortices.Above I DC ≃ 0.0805, the vortices start to penetrate the second line as well, and at high enough currents, I DC 0.085, their motion becomes mutually synchronized, similarly to the previous case displayed in Fig. 6a.Remarkably, in the intermediate current region, 0.0805 < I DC < 0.085, the two anharmonic oscillators have very different spectral fingerprints and, as a result, there is no clear synchronization of the vortex motion in the two lines; in this region, the V (I DC ) characteristics demonstrates a bump with several local maxima and minima.When AC-excitation is added, integer and fractional Shapiro steps are observed, the latter stemming from the anharmonic nature of vortex motion.The transitions to/from metastable modes are also observed, although their number is larger, their shape more complex and intricate than in the case of identical defects.Clearly, the vortex dynamics in the presence of asymmetric defects leads to a greater variety of collective motion modes.

III. DISCUSSION
The evolution of Shapiro features in Figs.5,8 with increasing I AC is not trivial.At low AC-excitation, I AC ≪ I DC , conventional Shapiro plateaus are narrow, and no signatures of metastable states are seen.In this regime, the AC-component acts as a probe that locks, at a fixed f AC , onto the spectrum of the vortex motion, solely determined by the main driving (Lorentz) force ∼ I DC .As I DC increases, the vortices move faster, f 1 and f k increase.At some I DC , a given f k gets close enough to nf AC , and the motion locks to f AC ; f 1 remains fixed in some range of I DC .As I DC further increases, the locking effect is lost.This results in a series of integer and fractional Shapiro plateaus visible on V (I DC ) curve at I AC =0.02.
When I AC is increased and becomes comparable with I DC , two phenomena appear.The first one is the wellknown enlargement of Shapiro plateaus, due to a stronger locking effect at higher AC-currents.The second one is related to the perturbation of the vortex motion spectrum by the oscillatory force ∼ I AC , whose amplitude becomes comparable to the Lorentz force due to DCcurrent.The combined action of I DC and I AC enables the existence of metastable states.They can be locked to f AC , resulting in jump-plateau-drop features as observed in Fig. 5.The same phenomenon takes place in Fig. 8, where many more voltage bumps and drops are observed (some jumps to metastable states are indicated by black arrows) as compared to Fig. 5.The lift of degeneracy, resulting in a more rich and complex metastable state spectrum, is certainly behind these differences.Finally, the DC-current range where the feature appears rapidly extends with increasing I AC .
In the limit of a dense, on the scale of W , network of defects, one would expect a huge number of apparently chaotically arranged voltage jump-bump-drops to appear on V (I DC ) curves, reflecting a vast number of accessed vortex motion modes and the complexity of the related spectra.The term "chaotic" is justified here due to a high sensitivity of the accessed metastables configurations to external parameters such as I DC , I AC , f AC , the disorder landscape, etc.Indeed, after unlocking from one metastable state, the system can jump down to a more stable configuration or lock up to another metastable state, from the available set.As a result, the position and shape of bumps-drops on V (I DC ) curves would appear arbitrary (see Fig. 8), while they are deterministic.
The revealed voltage drops correspond to a negative dynamic resistance dV /dI(I DC ).The latter has been experimentally observed in periodic pinning arrays subject to a specific external magnetic field [27][28][29], where a complex collective dynamics of vortices led to multiple phase transitions in their collective motion, with no need for additional AC-drive, resulting in various features in the V (I) characteristics [27][28][29].Another system is a perforated Nb film put in an external magnetic field, where the negative dynamic resistance can appear due to the Ratchet effect under AC-drive [3].More recently, both Shapiro steps and negative dynamic resistance were observed in a MoN strips with an artificial cut [30].The authors attributed the negative dynamic resistance to the chaotic aperiodic vortex motion at high AC-excitation amplitude.
The ability to use AC-excitation both as a pump and as a probe opens up interesting possibilities for realization, spectroscopy and control of metastable states in superconducting weak-links.The obtained results demonstrate the potential for designing artificial disorder landscapes to achieve desired responses to AC-amplitude and/or frequency.The general nature of weak-links suggests that there would be multiple ways of experimental realization of these functionalities.One of straightforward routes is to engineer superconducting films with a controlled disorder by using Focused Ion Beam approaches [3] or to deposit superconducting materials onto faceted structures [31].By carefully designing the spatial distribution of defects or grain boundaries, one could tailor the response of the weak-links to both DC-and AC-excitation.Another avenue is to overlap superconducting weak-links by ferromagnetic strips that can locally suppress the superconducting order parameter due to the inverse proximity effect [32][33][34][35][36].This can introduce additional complexity in the vortex dynamics and lead to novel effects under microwave excitation.

IV. CONCLUSION
In this work, we numerically studied transport properties of current-carrying superconducting nano-bridges subject to microwave illumination.The granularity of experimentally measured devices was accounted for by introducing one or two linear defects (simulating grain boundaries) which were directed perpendicularly to the applied current.We revealed a rich and complex dynamics of the vortex motion along these defects.Its strong anharmonicity enabled us to lock the spectrum of the system to an external periodic drive, and to obtain both integer and fractional Shapiro plateaus in DC voltagecurrent characteristics.In the case of two close linear defects, the inter-vortex coupling leads to the appearance of collective modes of correlated motion, with multiple stable and metastable states.These transitions are revealed in the current-voltage characteristic as regions of negative differential resistance dV /dI(I DC ).By playing with the external drive amplitude and frequency it becomes possible to pump the system to higher-resistance metastable modes and stabilize it there, in a finite range of DC transport currents.A step out of this range leads to a relaxation to a lower-resistance modes.The ability to control and stabilize different modes of the vortex motion opens up new possibilities for designing superconducting devices with tunable transport properties and novel functionalities.

V. METHODS
Within the TDGL framework, the temporal and spatial evolution of the complex superconducting order parameter ψ(t, r) can be expressed as [11]: where ψ is in units of ψ 0 = |a| b , with a and b being phenomenological parameters of the GL theory.The parameter u=1 is taken since we focus only on vortex motion but not on its nucleation dynamics.The coordinates r = (x, y) are in units of ξ.The scalar potential µ is measured in units of µ 0 = 2eτGL , where τ GL = 4πσλ 2 c 2 denotes the GL relaxation time, and λ is the London penetration depth.The parameter σ corresponds to the normal state conductivity of the material.The variable t is measured in units of τ GL , while the vector potential A is in units of H c2 ξ, with H c2 = c 2eξ 2 representing the upper critical field.The parameter ǫ(r) is associated with the local critical temperature T c (r) through Eq.(2); it enables spatially modulating the strength of the order parameter.
In the second GL equation, the total current J has superconducting (J S ) and a normal (J N ) components; it can be expressed in units of J 0 = cΦ0 8π 2 λ 2 ξ as: As TDGL equations are invariant under a gauge transformation, we use the zero scalar potential µ = 0 gauge to eliminate the scalar potential from both equations.For simplicity, we set the order parameter equal to zero ψ = 0 on boundaries y = 0, L. To apply external transport current I tr , we use boundary conditions for the vector potential on boundaries x = 0, W as ∇ × A = H I , where H I = 2πI tr /c represents the magnetic field induced by the transport current.On the other boundaries, we set ∇× A = 0. Additionally, we impose the superconductorvacuum boundary condition n • (∇ − ıA)ψ = 0 on boundaries x = 0, W , where n is the normal vector to the boundaries.
The state of the bridge is determined by calculating the voltage V between y = 0 and y = L boundaries for each value of transport current.In the chosen gauge, the electric field is written as E = −∂ t A. The corresponding instantaneous voltage drop V y1,y2 between two arbitrary points y 1 , y 2 in y-direction can be calculated as V y1,y2 (x, t) = − y2 y1 E y (x, y, t) dy = y2 y1 ∂ t A y (x, y, t) dy.(6) By averaging this voltage over the sample width and time we get the DC-voltage V measured in experiments.To avoid the voltage drops at y = 0, L boundaries, we calculate the voltage inside the bridge where the order parameter is fully restored (ψ = 1), as indicated by the blue dashed lines in Fig. 1.When simulating the Shapiro step experiments, we consider the microwave illumination as an additional time-dependent transport current of amplitude I AC and frequency f AC .The total transport current through the bridge is given by Eq.( 3).
In the model, all non-equilibrium quasiparticle processes are omitted, and for all considered frequencies, the microwave illumination acts on vortices only as an additional periodic Lorenz force.
The system of Eqs.(1), with the above-described boundary conditions, was solved using the commonly used link-variable method [11,[37][38][39] on the finitedifference grid.Spatial derivatives were approximated using the central difference method, and for time integration, the forward Euler method was employed [40].In all calculations of Shapiro steps, we set f AC τ GL = 0.03.

FIG. 1 .
FIG. 1. Defects in superconducting nano-bridges.a Scanning Electron Microscope image of a typical nano-bridge studied in experiments [10].b and c Two sample geometries studied theoretically, representing the central (narrowest) part of the real device.They contain one (b) or two (c) linear defects (grey regions).Edge defects situated at the ends of linear defects are indicated by black rectangles.The direction of the transport current is showed by arrows.The voltage is calculated between the two blue dashed lines.Further details are provided in the text.

FIG. 2 .
FIG. 2. Vortex dynamics in single linear defect with no AC current applied (I AC = 0).a Static map of the order parameter amplitude |ψ(r)| at low transport currents I DC ≪ Ic. b Static |ψ(r)| map at I DC IC indicates the presence of one vortex and one anti-vortex at the edge defects, ready to enter.d-e Snapshots of |ψ(r)| at different moments of vortex propagation for a fixed I DC = 0.10 > Ic. f Periodic temporal evolution of the instantaneous voltage V (t) at the same conditions.The dots c, d and e on the graph correspond to snapshots c, d and e.The period of V (t) oscillations provides the fundamental frequency f1 of the process.g Fourier spectrum of V (t).

cFIG. 4 .
FIG. 4. Evolution of the instantaneous voltage V (t)for the Shapiro plateau 1/3 of Fig.3.a V (t) at constant I DC and I AC for detuned frequency f AC = 0.625 f3.The lowfrequency envelope due to the beat effect is visible.b V (t) at the resonance f AC = f3.c Frequency spectrum of V (t) in the case (b).

FIG. 5 .
FIG. 5. V (I DC ) characteristics for different values of I AC in the case of two identical linear defects (displayed in fig.1c).The right vertical axis displays the numbers of Shapiro steps.Dashed lines are used as eye-guides (see in the text).

FIG. 6 .
FIG. 6. Vortex dynamics in the case of two identical linear defects of fig.1c. a and b Snapshots of the order parameter amplitude in the stable (a) and metastable (b) states near the transition (see in the text).c Evolution of V (t) at the transition from the metastable to the stable state at I AC = 0.03.The initial DC-current I DC = 0.076 switches to I DC = 0.077 at the moment t=500.

FIG. 7 .
FIG.7.Schematic energy diagram of the states considered in Fig.6.Stable (E0) and metastable (E1) states are presented along with their frequency spectra.

FIG. 8 .
FIG. 8. V (I DC ) characteristics for different I AC in the case of two different linear defects, ǫ = 0.5 and ǫ = 0.42.The right vertical axis displays the numbers of corresponding Shapiro steps.