Current-driven production of vortex-antivortex pairs in planar Josephson junction arrays and phase cracks in long-range order

Proliferation of topological defects like vortices and dislocations plays a key role in the physics of systems with long-range order, particularly, superconductivity and superfluidity in thin films, plasticity of solids, and melting of atomic monolayers. Topological defects are characterized by their topological charge reflecting fundamental symmetries and conservation laws of the system. Conservation of topological charge manifests itself in extreme stability of static topological defects because destruction of a single defect requires overcoming a huge energy barrier proportional to the system size. However, the stability of driven topological defects remains largely unexplored. Here we address this issue and investigate numerically a dynamic instability of moving vortices in planar arrays of Josephson junctions. We show that a single vortex driven by sufficiently strong current becomes unstable and destroys superconductivity by triggering a chain reaction of self-replicating vortex-antivortex pairs forming linear of branching expanding patterns. This process can be described in terms of propagating phase cracks in long-range order with far-reaching implications for dynamic systems of interacting spins and atoms hosting magnetic vortices and dislocations.

Topological defects such as vortices or dislocations determine many key properties of systems with long-range order, including superconductivity, superfluidity, magnetism, liquid crystals, and plasticity of solids 1,2 . These defects are characterized by a topological charge, for example, by an integer winding number n in the phase θ of the complex order parameter in exp( ) θ Ψ = ∆ for quantized vortices in superconductors. Destruction or creation of a single topological defect in the bulk requires overcoming a macroscopic energy barrier, resulting in the conservation of topological charge. In this work we address the question whether such topologically-protected stability of a vortex can be broken if a rapidly moving, two-dimensional vortex (2D) is driven by a strong force.
The physics of superfast vortices has attracted much attention in light of the development of superconducting qubits and digital memory [3][4][5] , sources of THz radiation 6 , or radio-frequency superconducting cavities for particle accelerators 7 . Theoretical and experimental investigations have uncovered a wealth of dynamic vortex phases [8][9][10][11][12][13][14][15][16][17][18][19] , while new imaging tools have made it possible to probe vortices at nanometer scales 20 and reveal hypersonic vortices moving much faster than the velocity of superfluid condensate 21 . The physics of ultrafast dislocations and their instabilities have also attracted much attention [22][23][24] . It has been usually assumed that a driven vortex preserves its identity as a topological defect, no matter how fast it moves, because instability of a vortex would violate the conservation of its topological charge. Yet it has been shown recently that a vortex driven by strong current in a long Josephson junction 25,26 or 1D Josephson junction array (JJA) 27 can become unstable due to Cherenkov radiation producing a cascade of expanding vortex-antivortex (V-AV) pairs which destroy the global superconducting phase coherence.
The issue of stability of driven topological defects becomes particularly intriguing in 2D systems in which long-range order can be destroyed by the Berezinskii-Kosterlitz-Thouless (BKT) transition resulting from thermally-activated unbinding of V-AV pairs 28,29 or a superconductor-insulator transition in JJAs 30,31 which model ultrathin granular films and monolayers [32][33][34] . So far theories of the BKT transition under current drive have assumed that current only affects thermally-activated or quantum V-AV pair production 31,35 , but the topological defects remain intact. Here we show that long-range order in driven 2D systems can be destroyed by a rapidly moving topological defect which causes a cascade of self-replicating pairs of defects of opposite polarity. For instance, a driven vortex in a planar JJA triggers a cascade of V-AV pairs which destroy the global phase coherence in the way similar to the propagation of cracks or elastic twins is solids 36 .
We performed numerical simulations of coupled sine-Gordon equations describing a particular case of planar and t ( ) ij β are phase differences between i-th and j-th junctions along the x and y axis, respectively, I I / c γ = , and I is a transport current applied uniformly along the y-axes, as shown in Fig. 1 quantifies inductive coupling of neighboring junctions, where I d and ξ are the depairing current and the coherence length in superconducting links of length s between the junctions. Here we take into account only the kinetic inductance of the superconducting condensate, disregarding long-range charge and inductive coupling between the grains through the electromagnetic stray fields, and focus on classical JJAs with negligible quantum charge effects 30,31 . Generic equations (1-3) can model dynamics of driven topological defects in many systems with long-range order, including artificial JJAs and granular superconducting films 30,31 , magnetic vortices and domain walls described by a classical XY model 37 , commensurate-incommensurate transitions and domain walls in charge density waves [38][39][40][41] , or dislocations in crystals 36 . The 2D systems are special as they exhibit a continuous BKT transition due to unbinding of topological defects in equilibrium 28,29 .

Results
We addressed the behavior of driven topological defects by simulating Eqs (1-3) for a 50 × 100 array with periodic boundary conditions along y and open boundary conditions along x at the edges of the array. Vortices were revealed by the circulation matrix ij ξ which identifies the vortex core in a cell where the sum of circulations along the plaquette is close to 2π. Because ξ ij only contains the circulation of the phase differences on the Josephson junctions but does not include phase gradients in the superconducting links due to circulating current in the vortex, ξ ij does not vanish outside the plaquette where the core is centered but decreases with the distance from that plaquette. In the subsequent figures the core boundary is outlined by the contour along which the sum of ij ξ equals π. To visualize the field of a moving vortex, we also calculated the z-component of the magnetic field H x y z ( , , ) z at the hight z above the array using the discrete Biot-Savart law: ij are the components of currents at position x′, y′ at the midpoint of each junction. We used a moving frame that keeps the vortex in the middle of the array after the vortex was initially introduced as explained in the Methods. The simulations were continued until a stationary state was reached, and the vortex velocity was obtained from the updating rate of frames. Steady state behavior. Our calculations show that at small currents γ γ < p the vortex is immobilized by intrinsic pinning in a JJA, in agreement with previous works 30,31 . If γ exceeds the depinning current γ p , the vortex moves and emits Cherenkov radiation along with bremsstrahlung caused by hopping of a vortex in the discrete JJA. Cherenkov radiation with the wave vectors k appears if the velocity of the vortex is smaller than the Swihart velocity at k 0 = but exceeds the phase velocity ϕ v k k ( , ) x y of small-amplitude waves in the JJA at finite k x and k y , where v k k ( , ) x y ϕ decreases as k x and k y increase 25,30,31 . For instance, Fig. 1

shows a snapshot of the field map
The radiated wake in the Cherenkov cone behind the vortex and the interference of waves reflected from the edges are apparent, the radiation appears at currents well below I c . The amplitude of Cherenkov wake diminishes as the damping parameter η increases, and radiation behind moving vortices was detected at η up to 0.5. As the inductance parameter ε increases radiation becomes more intense.
The stationary velocity of the vortex v( , , ) η ε γ was calculated as a function of current γ for underdamped arrays with η < 1 for which the effect of radiation is most pronounced. In our simulations the driving current γ = kt was gradually increased with a small ramp rate k 5 10 6 = ⋅ − to avoid artifacts due to transient instabilities. Shown in Fig. 2a is a mean vortex velocity v(γ) as a function of current calculated for 0 3 η = . and different values of ε. Likewise, Fig. 2b shows v(γ) calculated for 0 3 ε = . and different values of η. In both cases v 0 = below the depinning current γ ε ( ) p . The velocity has a jump at γ γ = p and then increases smoothly up to a second critical current γ s at which the second jump in v(γ) occurs.
In the region of currents p s γ γ γ < < the vortex moves with a steady-state velocity controlled by the balance between the driving Lorentz force and the drag forces caused by ohmic currents in the junctions and radiation losses, as illustrated by Fig. 1c. Here the threshold current γ s calculated for different values of η is a nonmonotonic function of the inductance parameter ε, as shown in the inset in Fig. 2. The velocity of the vortex v(γ) defines the dc electric field-current characteristics of the array in the flux flow state, where n n is the areal density of vortices and φ 0 is the magnetic flux quantum.

Types of instabilities.
As γ exceeds γ s the steady-state vortex becomes unstable, causing striking dynamic patterns shown in Fig. 3. The instability proliferates from the Cherenkov wake behind a moving vortex where a critical nucleus is formed by the junctions being in the unstable state with the phase differences confined between π/2 and 3π/2, similar to that of a driven vortex in a long JJ 25,26 . Manifestations of this instability are different for weak ( 1) η  and moderate η  ( 1) damping, yet in any case the dynamic structures represented in Fig. 3 result from continuous production of V-AV pairs and the subsequent separation of vortices and antivortices which accumulate at the opposite sides of an expanding domain. The pattern shown in Fig. 3 can be regarded as an expanding dipole of growing positive and negative topological charges with a fixed net topological charge of the initial vortex. This process bears a remarkable similarity with the crack propagation in crystalline solids where the relevant topological defects are edge dislocations 36 . Indeed, a chain of vortices and antivortices produces a staircase phase difference x ( ) β ∆ on the Josephson junctions along the chain as shown in Fig. 4. For the spatially separated V and AV cores marked by the yellow and magenta contours in Fig. 3a, the phase jumps on the junctions connecting the aligned V-AV pairs first increase with the distance from the last vortex, reaches the maximum in the middle of the chain and then decreases to zero at the other end of the chain 25 . As the vortex pattern expands, the staircase distribution of phase differences along the V-AV chain evolves into a growing dome-like x t ( , ) β ∆ which can be regarded as a "phase crack" in the long-range superconducting order.
In underdamped arrays represented by Fig. 3a,b the expanding vortex dipole is mostly caused by Cherenkov and bremsstrahlung radiation. At γ γ > s the amplitude of the radiation wake behind the moving vortex exceeds a threshold of generation a V-AV pair. Here the critical wake first produces an expanding V-AV pair which in turn generates enough radiation to produce two more V-AV pairs, triggering a subsequent V-AV chain reaction and continuous pair production. Our simulations performed in a wide range of ε and η have shown that in a moder- ( 1)  this radiation mechanism produces V-AV pairs aligned along the direction of motion of the initial vortex. Despite a rather complex 2D magnetic field pattern shown in Fig. 3a, the multiquantum vortex dipole essentially results from a 1D Cherenkov instability similar to that was observed in previous simulations of the Cherenkov vortex instability of vortices in nonlocal Josephson junctions 25,26 or 1D junction arrays and dislocations in the Frenkel-Kontorova model 27 . As η increases, the amplitude of radiation wake decreases, so the instability occurs at higher currents γ s which are still smaller than the critical current of the Josephson junctions γ = ( 1). Shown in Fig. 3a are the color map of the magnetic field and the circulations ij ξ at the initial stage of the V-AV pair production, where the red and blue regions correspond to vortices and antivortices, respectively, and the yellow and purple contours outline the boundaries of vortex cores. Figure 3a is characteristic of V-AV pairs production in a moderately underdamped array with η = . 0 5 where the V-AV pairs propagate along the trajectory of the initial vortex. One can see waves radiated from the central zone where the V-AV pairs are generated. As the ). In this case the vortices structure grows in the vertical direction and the radiation wake is not visible. on the row of junctions along the initial V-AV pair propagation. Each ascending step in β x corresponds to an antivortex and the decending step to a vortex, so that the net phase difference 2π between the left and right ends remains fixed. Calculation was done for η = . Here the red region γ γ ε < ( ) p corresponds to the pinned vortex, and the yellow region corresponds to a stable vortex moving with a constant mean velocity (see Fig. 1). The uniform motion of the vortex is unstable below the green line ( , ) s η η γ ε = . In turn, the unstable region is divided into two regions corresponding to two different mechanisms of generation of V-AV pairs represented in Fig. 3: the magenta region corresponds to the radiation-assisted generation of V-AV pairs, and the cyan region corresponds to the doublet V-AV pair emission. Here the color maps shown in Fig. 3a,c were calculated for the parameters corresponding to the two stars in Fig. 5. The transition between the line and doublet V-AV pair production occurs in a crossover region where the slope of the green boundary line η γ ( ) s increases sharply. As the parameter ε decreases, the depinning and threshold lines shift to the right.

Discussion
Our simulations of JJAs revealed rich dynamics of driven vortices controlled by interplay of ohmic losses, Cherenkov radiation and bremsstrahlung. Emission of V-AV pairs by a moving vortex at s γ γ > can be described in terms of propagating phase cracks in the order parameter, which can be applicable to driven topological defects in a wide class of systems, including vortices in JJAs, granular superconducting films or superfluid films, magnetic vortices and domain walls or dislocations in crystals. Dynamic states of propagating phase cracks can be controlled by current and the damping parameter, giving rise to either 1D self-replication of V-AV pairs and branching phase cracks at η 1  or generation of doublets of V-AV pairs in the transverse direction producing beads of macrovortices at  1 η , as illustrated in Fig. 3. The above radiation mechanism of vortex splitting resulting in the subsequent V-AV pair production in the 2D JJAs is markedly different from the well-known proliferation of pairs of Abrikosov vortices and antivortices from the edges of a wide, current-carrying superconducting film where a transition of a chain of expanding Abrikosov V-AV pairs into a phase slip can occur 21,42,43 . The Cherenkov V-AV pair production triggered by a moving vortex in single long JJ in wide thin film or monolayers 25,26 can be strongly affected by pinned Abrikosov vortices in the electrodes. Interaction of Josephson and Abrikosov vortices 44,45 could in principle be used to tune the V-AV pair production in long Josephson junctions. In large planar junctions, the V-AV pair production results in expanding Josephson vortex loops 46 .
In conclusion, this work addresses a fundamental question of what happens to long-range order if a driven topological defect becomes unstable. We show that the long-range order in driven 2D systems can be destroyed Figure 5. Phase diagram. The η γ − phase diagram was calculated from the data shown in Fig. 2b at ε = 1. Here the read region corresponds to a pinned vortex and the yellow region corresponds to a stable moving vortex. The stars show the values of the parameters for which the color maps in Fig. 3 were calculated.
ScieNtific REPORTS | (2018) 8:15460 | DOI:10.1038/s41598-018-33467-y by a single topological defect triggering a cascade of self-replicating pairs of defects of opposite polarity. This process can be initiated be either a quenched defect which is already present in the system (for example, a trapped vortex in a superconductor or a residual dislocation in a crystal) or appears due to thermal fluctuations. The resulting V-AV pair production can manifest itself in jumps on the V-I characteristics in underdamped granular superconductors at I I c  , the jumps being unrelated to Joule heating 47 . Moreover, one can expect new effects of current on the vortex-charge duality and superconducting-insulator transition in Josephson arrays and granular films 30,31 . Propagation of phase cracks could be one more mechanism of jumps on V-I characteristics observed in granular thin films at the onset of insulator transition 48,49 which so far have been attributed to electron overheating 50,51 . The V-AV pair production triggered by a moving vortex can greatly enhance the output signal in single photon detectors based on arrays of Josephson junctions 52,53 . Our results also suggest a mechanism of initial stage of origination and propagation of mechanical cracks or twins initiated by dislocations in graphene and other atomic monolayers under strong tensile stress.

Methods
The 4th order Runge-Kutta method has been used to solve the set of equations (1-3) numerically. The implementation of the boundary conditions has been imposed through the circulations ( ij ξ ). For periodic boundary conditions, the circulation of each plaquette at the top boundary of the array has been forced to be the same as the circulation at the corresponding plaquette at the bottom boundary, and vice versa. For open boundary conditions, all the circulations at the left and right boundaries have been set equal to zero. Unlike the XY model, where there is a freedom of choosing the phase of an element, in our JJA model we deal with phase differences between the junctions and it is not trivial to establish a vortex configuration. In order to do so, we construct an initial environment that includes a weak link in the direction of the current which triggers the appearance of vortices. At position = i 25 and j 0 = , we include a factor χ in the sine term: The value of χ = .
0 01 was chosen in our simulations. If the current γ is high enough then a cascade of equispaced vortices starts to emerge from the weak link. The distance between vortices (or the generation rate) depends on the current, and we choose a value low enough to produce a clear isolated vortex. Once the vortex is far from the weak link, the configuration of αs and βs is stored in order to be used as the initial configuration of successive runs.
To avoid the vortex reaching the right end of our system, we used a moving frame. The full configuration of the array is shifted by a single plaquette when the center of the vortex reaches a certain position (equal to the lattice spacing). The updating rate gives us the vortex velocity.
The way the velocity at Fig. 2 is calculated has some limitations if we have more than one vortex. Frames are updated according to the position of the maximum circulation plaquette of the whole array, which may not corresponds to the original vortex after the first V-AV emission, i.e., for γ γ > s . Furthermore, when new vortices appear they reach the left boundary in a short time with the corresponding leak of topological charge.