Fragmentation of Fast Josephson Vortices and Breakdown of Ordered States by Moving Topological Defects

Topological defects such as vortices, dislocations or domain walls define many important effects in superconductivity, superfluidity, magnetism, liquid crystals, and plasticity of solids. Here we address the breakdown of the topologically-protected stability of such defects driven by strong external forces. We focus on Josephson vortices that appear at planar weak links of suppressed superconductivity which have attracted much attention for electronic applications, new sources of THz radiation, and low-dissipative computing. Our numerical simulations show that a rapidly moving vortex driven by a constant current becomes unstable with respect to generation of vortex-antivortex pairs caused by Cherenkov radiation. As a result, vortices and antivortices become spatially separated and accumulate continuously on the opposite sides of an expanding dissipative domain. This effect is most pronounced in thin film edge Josephson junctions at low temperatures where a single vortex can switch the whole junction into a resistive state at currents well below the Josephson critical current. Our work gives a new insight into instability of a moving topological defect which destroys global long-range order in a way that is remarkably similar to the crack propagation in solids.

Scientific RepoRts | 5:17821 | DOI: 10.1038/srep17821 THz radiation 9 , or polycrystalline superconducting resonator cavities for particle accelerators 10 , and have broader implications for other systems with long-range order.
We start with a standard theory of a Josephson vortex in a long junction described by the sine-Gordon equation for the phase difference of the order parameter θ(x, t) = ϕ 1 − ϕ 2 between two bulk electrodes 3,4 : Here the prime and the overdot denote partial derivatives with respect to the dimensionless coordinate x/λ J and time ω J t, ω J = (2πJ c /φ 0 C) 1/2 is the Josephson plasma frequency, J c is the tunneling critical current density, C is the specific capacitance of the junction, λ J = (φ 0 /4πμ 0 λJ c ) 1/2 is the Josephson penetration depth, λ is the London penetration depth, η = 1/ω J RC is the damping constant due to the ohmic quasiparticle resistance R, and β = J/J c is the driving parameter controlled by a uniform transport current density J.
The sine-Gordon equation has been one of the most widely used equations to describe topological defects in charge and spin density waves 11 , commensurate-incommensurate transitions [12][13][14] , magnetic domain walls 15 , dislocations in crystals 16,17 , kinks on DNA molecules 18 moving with a constant velocity v, where c s = ω J λ J is the Swihart velocity of propagation of electromagnetic waves along the junction 3 . As v increases, the vortex shrinks at η ≪ 1 and expands at η > 1 4 .
Unlike the sine-Gordon equation, the nonlocal equation (2) at η = 0 is not Lorentz-invariant, so a uniformly moving vortex can radiate Cherenkov waves δθ(x, t) ∝ exp(ikx − iω k t) with the phase velocities ω k /k smaller than v 23,29 . The condition of Cherenkov radiation at η = 0 is given by: 0 0 2 , and G(k) is the Fourier image of G(x). Here G(k) decreases as 1/k at k > Λ −1 so equation (4) is satisfied for k > k c where the maximum wavelength 2π/k c increases with v 30 . To address the effect of Cherenkov radiation on the moving vortex, we performed numerical simulations of equation (2) for SIS junctions of different geometries.
Shown in Fig. 1 are the numerical results for a planar bulk junction at η = 0.05 and the large ratio λ J /λ = 10 usually described by the sine-Gordon equation (1). Yet the more general integral equation (2) reveals the effects which are not captured by equation (1), particularly a trailing tail of Cherenkov Scientific RepoRts | 5:17821 | DOI: 10.1038/srep17821 radiation behind a vortex moving with a constant velocity 29 . Moreover, as the amplitude and the wavelength of radiation increase with v, the vortex becomes unstable at β > β s , the instability is triggered at the highest maximum of Cherenkov wave where θ m reaches a critical value θ c ≈ 8.65-8.84, depending on η, λ/Λ , and the junction geometry 30 . Here θ c is confined within the interval 5π/2 < θ c < 3π in which a uniform state of a Josephson junction is unstable 3,4 . As the velocity increases, the domain where 5π/2 < θ(x − vt) < 3π behind the moving vortex widens and eventually becomes unstable as its length exceeds a critical value. This suggests a qualitative picture of the vortex instability caused by the appearance of a trailing critical nucleus being in the unstable π-junction state 3,4 caused by strong Cherenkov radiation. The latter appears entirely due to the Josephson nonlocality described by equation (2), which has no steady-state vortex solutions at J > J s where J s can be well below J c at which the whole junction switches into a resistive state.
The dynamic solutions of equation (2) at β > β s change strikingly. Our simulations have shown that the instability originates at the highest maximum θ = θ m of the trailing Cherenkov wave which starts growing and eventually turning into an expanding vortex-antivortex pair 30 , as shown in Fig. 1. As the size of this pair grows, it generates enough Cherenkov radiation to produce two more vortex-antivortex pairs which in turn produce new pairs. Continuous generation of vortex-antivortex pairs results in an expanding dissipative domain in which vortices accumulate at the left side, antivortices accumulate at the right side, while dissociated vortices and antivortices pass through each other in the middle 30 . As a result, θ(x, t) evolves into a growing "phase pile" with the maximum θ m (t) increasing approximately linear with time and the edges propagating with a speed which can be both smaller and larger than c s , the phase difference θ(∞) − θ(− ∞) = 2π between the edges remains fixed. We observed the phase pile dynamic state for different junction geometries and η ranging from 10 −3 to 0.5 30 . For instance, Figs 2 and 3 show the 3D images of the initial stage of dynamic separation of vortices and antivortices calculated for a bulk junction and a thin-film edge junction. Here the local magnetic field B(x, t) oscillates strongly at the moving domain edges but becomes rather smooth away from them, as shown in Fig. 4. In the most part of the phase pile overlapping vortices are indistinguishable, yet the net flux φ = φ 0 of this evolving multiquanta magnetic dipole remains quantized. Shown in Fig. 5 are the steady-state vortex velocities v(β) calculated for different junction geometries. The instability corresponds to the endpoints of the v(β) curves which have two distinct parts. At small  β η the velocity v(β) increases sharply with a slope limited by a weak quasiparticle viscous drag. At larger  β η the increase of v(β) with β slows down, as the vortex velocities are mostly limited by radiation friction 29 and depend weakly on the form of dissipative terms in equation (2). For a low-J c junction with λ J /λ = 10, the effect of Cherenkov radiation on v(β) is weak, but for a high-J c bulk junction with λ λ / = 10 J and η ≪ 1, radiation friction dominates at practically all β, significantly reducing both v(β) and β s . For thin film edge junctions, the critical splitting current density J s gets reduced down to J s ≈ 0.4J c at η = 10 −3 , as shown in Fig. 5. In the extreme nonlocal limit described by equations (2) . Dynamics of θ(x, t) in the nonlocal limit at J > J s is similar to that is shown in Figs 1-3, except that the edges of phase pile can propagate with "superlu- 30 . Once vortex-antivortex pairs start replicating, the speed of leading vortices at the edges gradually increases from v s to a limiting value v ∞ , for instance, from v s ≈ 0.72lω J to v ∞ ≈ 1.12lω J for an edge junction with l = Λ /2 and η = 0.1 30 .
The effects reported here are most pronounced in underdamped SIS junctions between s-wave superconductors at low temperatures for which the viscous drag coefficient η ∝ exp(− Δ /T) due to thermally-activated quasiparticles 3 is small. Here η ≪ 1 also implies that a moving vortex does not generate additional quasiparticles because the induced Josephson voltage θ = ′ / V v eL 2 m is smaller than Δ /e, where θ ′ m is the maximum phase gradient. These conditions are satisfied for the parameters used in    This equation shows that the power P is independent of J c and is greatly reduced in the underdamped limit at low temperatures as the quasiparticle resistance R of SIS junctions becomes exponentially large at T ≪ T c . To estimate P, it is convenient to write equation (5) in the form ηε ω 2 is a characteristic line energy of Abrikosov vortex 31 . For an edge junction in a Nb film with t = 1 nm, λ = 40 nm, ε 0 ~ 10 4 kelvin/nm, and ω J = 100 GHz much smaller than  ∆/ .  2 4 THz 10 , equation (5) yields P ~ 0.16 nW at η = 10 −2 . Local overheating δT = PY K caused by vortex dissipation is further reduced in thin film junctions for which the energy transfer to the substrate due to ballistic phonons is much more effective than diffusive phonon heat transport in thick samples, where Y K is the Kapitza interface thermal resistance 32 . Such weak overheating caused by a moving vortex cannot result in thermal bistability and hysteric switching due to hotspot formation 32 .
Proliferation of vortex-antivortex pairs triggered by a moving Josephson vortex can be essential for the physics and applications of weak link superconducting structures where the formation of expanding phase pile patterns can switch the entire junction into a normal state at currents well below the Josephson critical current, > ( . − . )  J J J 0 4 0 7 s c . Such dynamic vortex instability can result in hysteretic jumps on the V-I curves which appear similar to those produced by heating effects 4,9 , yet this instability is affected by neither cooling conditions nor the nonequilibrium kinetics of quasiparticles. Indeed, heating is most pronounced in overdamped junctions with η > 1 in which Cherenkov radiation is suppressed. By contrast, the Cherenkov instability is characteristic of the weakly-dissipative underdamped limit η ≪ 1, although Fig. 5 shows that this instability in thin film edge junctions can persist up to η = 0.5. Therefore, the crucial initial stage of the phase pile formation at η ≪ 1 is unaffected by heating which may become more essential at the final stages of the transition of the entire junction into the normal state. At η ~ 1 the Cherenkov instability may be masked by heating effects, particularly in bulk junctions for which heat transfer to the coolant is less efficient than in thin films.
It should be emphasized that the instability reported here does not require special junctions with J c ~ J d . In fact, even for the seemingly conventional bulk junction with λ J = 10λ shown on the top panel of Our results can be essential for other topological defects such as crystal dislocations or magnetic domain walls described by the generic nonlocal equation (2) in which the integral term results from a common procedure of reduction of coupled evolution equations for several relevant fields to a single equation. For Josephson junctions, such coupled fields are θ and B, but for domain walls in ferromagnets, the nonlocality can result from long-range magnetic dipolar interactions 35 . For dislocations, the nonlocality and Cherenkov radiation of sound waves in equation (2) come from the discreteness of the crystal lattice 17 and long-range strain fields 16 , although the dynamic terms in the Peierls equation 36,37 are more complex than those in equation (2). Dynamic instabilities of dislocations have been observed in the lattice Frenkel-Kontorova models 17 in which sonic radiation can also result from periodic acceleration and deceleration of a dislocation moving in a crystal Peierls-Nabarro potential 16 . The latter effect becomes more pronounced as the dislocation core shrinks at higher velocities and becomes pinned more effectively by the lattice. By contrast, the instability reported here results entirely from Cherenkov radiation, the condition (4) can be satisfied for any system in which G(k) in equation (4) decreases with k. This instability can thus have broader implications: for instance, the phase pile dynamics of Josephson vortices appears similar to a microcrack propagation caused by a continuous pileup of subsonic dislocations with antiparallel Burgers vectors at the opposite tips of a growing crack described by equations (2) and (3) 16 .
Our results give a new insight into breakdown of a global long-range order which has been usually associated with either thermally-activated proliferation of topological defects (like in the Berezinskii-Kosterletz-Thouless transition) or static arrays of quenched topological defects pinned by the materials disorder 2 . Here we point out a different mechanism in which a long-range order is destroyed as a single topological defect driven by a strong external force becomes unstable and triggers a cascade of expanding pairs of topological defects of opposite polarity.

Methods
We have developed an efficient MATLAB numerical code to solve the main integro-differential equation (2) using the method of lines 38 . By discretizing the integral term in equation (2) it was reduced to a set of coupled nonlinear ordinary differential equations in time which were solved by the multistep, variable order Adams-Bashforth-Moulton method 39 . We have checked our numerical results using a slower iterative method to make sure that the logarithmic singularity of G(x − u) is handled properly, the absolute and relative error tolerances were kept below 10 −6 . The length L b of computational box x 1 < x < x 1 + L b along the x-axis (either co-moving with the vortex or expanding with the phase pile) was taken large enough to assure no artifacts coming from possible reflected waves at x = x 1 and x = x 1 + L b . We set 6 and made sure that changing L b does not affect the results, where L b was typically taken at least three times larger than the spatial extent of θ(x, t), be it a single vortex or expanding phase pile. The steady state phase distribution θ(x − vt) in a uniformly moving vortex at a given β was computed by solving the full dynamic equation (2) using the single-vortex solution calculated at a smaller preceding value of β as an initial condition. The code then run until the velocity of the vortex stabilizes to the accuracy better than 0.1%.