Nonlinear dynamics of Josephson vortices in merging superfluid rings

TL;DR

This paper investigates the complex nonlinear dynamics of Josephson vortices during the merging of two superfluid rings with different vorticities, revealing 3D vortex evolution, instabilities, and the influence of initial conditions.

Contribution

It provides new insights into 3D vortex behavior and the effects of initial population imbalance and trap shape on superfluid ring mergers.

Findings

Josephson vortices emerge in the barrier during ring merging.

Kelvin-Helmholtz instability can develop at the interface.

Final vortex states depend on initial conditions and trap geometry.

Abstract

We consider merger of two parallel toroidal atomic Bose-Einstein condensates with different vorticities in a three-dimensional (3D) trap. In the tunnel-coupling regime, Josephson vortices (rotational fluxons) emerge in the barrier between the superflows. When the barrier is gradually eliminated, we observe essentially three-dimensional evolution of quantum vortices, which may include the development of the Kelvin-Helmholtz instability at the interface between the rings, in the framework of a weakly dissipative Gross-Pitaevskii equation. An initially more populated ring, carrying a persistent current, can drag an initially non-rotating less populated one into the same vortex state. The final state of the condensate crucially depends on an initial population imbalance in the double-ring set, as well as on the shape of the 3D trapping potential, oblate or prolate. In the prolate (axially elongated) configuration, robust 3D hybrid structures may appear as a result of the merger of persistent currents corresponding to different vorticities.

Equations

$$(i-\gamma )\hslash \frac{\partial \psi }{\partial t}=-\frac{{\hslash }^2}{2M}{\nabla }^2\psi +V_{\text{ext}}(\mathbf{r},t)\psi +g|\Psi {|}^2\psi -\mu \psi ,$$

$$V_{\mathrm{ext}}(\rho ,z,t)=\frac{1}{2}M{\omega }_r^2(\rho -{\rho }_0{)}^2+\frac{1}{2}M{\omega }_z^2{z}^2+V_{\text{b}}(z,t),$$

$$V_{\text{b}}(z,t)=U_{\text{b}}(t)\exp \left(-\frac{1}{2}\frac{(z-z_0{)}^2}{a^2}\right),$$

$$U_{\text{b}}(t)=\left\{\begin{array}[]{c}(1-t/t_{d})u_{b},~{}\text{at~{}}~{}t<t_{d},\\ 0,~{}~{}\text{at~{}}~{}t>t_{d}.\end{array}\right.$$

$$(i-\gamma )\frac{\partial \psi }{\partial t}=-\frac{1}{2}{\nabla }^2\psi +V_{\text{ext}}\psi -\mu \psi +g|\psi {|}^2\psi ,$$

$$V_{\text{ext}}=\frac{1}{2}(\rho -{\rho }_0{)}^2+\frac{1}{2}{A}^2{z}^2+V_{\text{b}},$$

$$A={\omega }_z/{\omega }_r.$$

$$\Psi (\mathbf{r})=|{\Psi }_0(x,y,z)|e^{iS(z)\varphi },$$

$$⟨{\Psi }_0|{\Psi }_0⟩=N\equiv N_1+N_2,$$

$$N_1$$

$$N_2$$

$$P=(N_1-N_2)/(N_1+N_2).$$

$$E=\int \left[\frac{1}{2}|\nabla \psi {|}^2+V_{\text{ext}}(\mathbf{r})|\psi {|}^2+\frac{g}{2}|\psi {|}^4\right]d\mathbf{r},$$

$$\mathbf{L}=-\frac{i}{2}\int \left\{{\psi }^{\ast }\left[\mathbf{r}\times \nabla \psi \right]-\psi \left[\mathbf{r}\times \nabla {\psi }^{\ast }\right]\right\}d\mathbf{r},$$

Full text

I Introduction

One of the most remarkable manifestations of quantum properties of superconductors A. V. Ustinov (1998) and superfluids Kaurov and Kuklov (2005, 2006) is formation of Josephson vortices (JVs), alias fluxons, in long Josephson junctions. In particular, much interest to JVs in atomic Bose-Einstein condensates (BECs) has been drawn since circulating atomic supercurrents, counterpropagating in two parallel 1D BEC, were theoretically investigated in Refs. Kaurov and Kuklov (2005, 2006). Multidimensional JVs in spin-orbit-coupled BECs were considered in Ref. Gallemí, A. and Guilleumas, M. and Mayol, R. and Mateo, A. Muñoz (2016). Rotational fluxons in BEC loaded into coplanar double-ring traps were studied in Ref. Brand et al. (2009). Further, it was demonstrated in Ref. Su et al. (2013) that atomic BECs, confined by a dual-ring trap, support JVs as topologically stable defects in the filed of the relative phase. Stationary JVs were predicted in an array of linearly-coupled one-dimensional Bose-Einstein condensate in Ref. Christian Baals, Herwig Ott, Joachim Brand, Antonio Muñoz Mateo (2018). Very recently, JVs were investigated in a closed-loop stack of 1D and 2D atomic BECs Gil Granados et al. (2019).

In our recent work Oliinyk et al. (2019) it was demonstrated that the azimuthal structure of the tunneling flow implies formation of Josephson vortices (rotational fluxons) with zero net (integral) current through the junction for hybrid states, built of persistent currents with different topological charges in coupled rings (in particular, these include the case of opposite topological charges – this is the configuration which was defined as hybrid states in Ref. R. Driben, Y. Kartashov, B. A. Malomed, T. Meier, and L. Torner (2014), where they were supported by an effective nonlinear potential). It turns out that the fluxons’ cores rotate and bend, following the action of the quench, i.e. formation of tunnel junction with chemical potential difference. Further, it was found that, as the barrier decreases, and the effective coupling between the rings respectively increases, the JVs accumulate more and more energy. This scenario suggests the question if all quantum vortices, which are present in the settings, viz., the fluxons and circular currents, are always unstable, and what the final state of the merging counter-propagating superflows is in the double-ring system.

In this connection, it is relevant to mention that the Kelvin-Helmholtz instability (KHI), which is well known in classical hydrodynamics, leads to growth of roll-up perturbations at the interface of two streams with different velocities W. R. Peltier and C. P. Caulfield (2003). KHI plays a crucial role in destroying laminar flows, driving them into a turbulent state. In classical viscous liquids and gases, adjacent streams with initially different velocities evolve into a state with the averaged velocity. For superfluids the dynamics of KHI is fundamentally different M. Tsubota, M. Kobayashi, H. Takeuchi (2013). In particular, a pair of weakly interacting persistent currents in toroidal Bose-Einstein condensates (BECs), carrying quantum vortices, may undergo transitions between different angular-momentum states via phase slips, at which the angular momentum changes by unitary quanta. In sharp contrast to classical streams, there is no state with an “intermediate velocity” for two interacting persistent currents initially created with different vorticities. The quantum nature of the transition between different angular-momentum states leads to far-reaching consequences for the development or suppression of the KHI in ultracold atomic gases.

Dynamical KHI at the interface in phase-separated two-component Bose-Einstein condensates was investigated in Refs. Suzuki et al. (2010); Takeuchi et al. (2010). In recent theoretical work A.W. Baggaley, N.G. Parker (2018) two merging counter-propagating streams have been investigated in the framework of conservative 2D (two-dimensional) Gross-Pitaevskii equation (GPE). Planar counterflow channels with periodic boundary conditions mimic two narrow (in the radial direction) ring-shaped condensates which merge when the separating barrier is switched off. However, the description of KHI in realistic atomic BECs calls for a more detailed analysis. On the one hand, vortex lines, which drive KHI in the double-ring system, are oriented along the orthogonal axis, which requires to introduce a full 3D model for the consideration of the vortex dynamics. Further, weak dissipative effects must be taken into account to analyze the relaxation to the stationary state after the transition induced by KHI.

In the present work we address the following issues: (i) How nonlinear dynamics of quantum vortices drives the relaxation of two merging persistent currents, and (ii) identification of a final state of two merging ring-shaped condensates with different initial vorticities. The analysis reveals two noteworthy effects: (i) In pancake-shaped configurations (those compressed in the axial direction), a more populated ring can impose its quantum state onto the originally less populated one. In particular, if the ring carries topological charge $1$, the persistent current merges with a non-rotating ring. The final state of the resulting toroidal condensate is either one with zero angular momentum, or a single-charged quantized flow, depending on the initial difference in populations of the two rings. (ii) Instead of the development of KHI at the interface of the merging persistent currents, we observe formation of robust hybrid vortex structures in the potential trap sufficiently elongated in the axial direction.

The rest of the paper is organized as follows. The model is formulated in Section II. Results of the systematic numerical analysis are summarized in Section III, separately for the pancake-shaped and elongated shapes of the trapping potential. The paper is concluded by Section IV.

II The model

In modeling nonequilibrium dynamics, such as quantum turbulence M. Tsubota, M. Kobayashi, H. Takeuchi (2013) or nucleation of vortices A. I. Yakimenko, Yu. M. Bidasyuk, O.O. Prikhodko, S.I. Vilchinskii, E.A. Ostrovskaya, and Yu. S. Kivshar (2013), dissipative effects, even if they are weak, are of crucial importance for relaxation to equilibrium states. In particular, the dissipation drives the drift of the vortex core to the edge of the BEC cloud Rooney et al. (2010); Moulder et al. (2012); Yakimenko et al. (2013); Allen et al. (2013); Bidasyuk et al. (2015). Such effects naturally arise in a trapped condensate due to its interaction with a thermal component, and can be captured phenomenologically by the dissipative GPE derived by Choi et al. S. Choi, S. A. Morgan, and K. Burnett (1998); N. P. Proukakis and B. Jackson (2008). Close to the thermodynamic equilibrium, the weakly dissipative GPE is written as

[TABLE]

where $g=4\pi a_{s}\hbar^{2}/M$ is the coupling strength, $M$ is the atomic mass ($M=3.819\times 10^{-26}$ kg for 23Na atoms), $a_{s}$ is the $s$-wave scattering length (positive $a_{s}=2.75$ nm, corresponding to the self-repulsion in the same atomic species, is used below), $\mu$ is the chemical potential of the equilibrium state, and $\gamma\ll 1$ is a phenomenological dissipative parameter. This form of the dissipative GPE has been used extensively in previous studies of vortex dynamics (see, e.g., M. Tsubota, M. Kobayashi, H. Takeuchi (2013); R. Carretero-Gonz´alez, N. Whitaker, P. G. Kevrekidis, and D. J. Frantzeskakis (2008); K. Kasamatsu, M. Tsubota, and M. Ueda (2003); A. I. Yakimenko, Yu. M. Bidasyuk, O.O. Prikhodko, S.I. Vilchinskii, E.A. Ostrovskaya, and Yu. S. Kivshar (2013)). In what follows below, we assume $\gamma$ to be spatially uniform, and set $\gamma=0.03$ as in Refs. M. Tsubota, M. Kobayashi, H. Takeuchi (2013); Yakimenko et al. (2013). Actually, we have verified that results reported below do not essentially depend on a specific value of $\gamma\ll 1$.

We consider a toroidal condensate, split by a blue-detuned sheet beam in upper and lower weakly coupled rings-shaped components. The respective total trapping potential is

[TABLE]

where $\rho\equiv\sqrt{x^{2}+y^{2}}$ and the sheet potential is

[TABLE]

with the time-dependent strength,

[TABLE]

Here the switching time is $t_{d}=0.015$ s, and $z_{0}$ is a possible shift of the barrier shift along the $z$-axis.

For numerical simulations of the 3D GPE we use scaled time, $t\rightarrow t\omega_{r}$, length $\mathbf{r}\rightarrow\mathbf{r}/l_{r}$, chemical potential $\mu\rightarrow\mu/(\hbar\omega_{r})$, external potential $V_{\text{ext}}\rightarrow V_{\text{ext}}/(\hbar\omega_{r})$, and wave function $\psi\rightarrow\psi\cdot l_{r}^{3/2}$, which casts GPE (1) in the following form:

[TABLE]

where the scaled positive nonlinearity strength is $g=4\pi a_{s}/l_{r}$, and the scaled trapping potential is

[TABLE]

where the aspect ratio of the toroidal trap is

[TABLE]

It turns out that dynamics of the quantum vortices, observed after the merger of the rings, crucially depends on $A$. In this work we first concentrate on the pancake-shaped trapping potential with typical values of the trapping frequencies K. C. Wright, R. B. Blakestad, C. J. Lobb, W. D. Phillips, G. K. Campbell (2013); F. Jendrzejewski, S. Eckel, N. Murray, C. Lanier, M. Edwards, C. J. Lobb, and G. K. Campbell (2014): $\omega_{r}=2\pi\times 123$ Hz and $\omega_{z}=2\pi\times 600$ Hz, hence $A=4.88$, the oscillator length of the longitudinal trapping potential is $l_{r}=\sqrt{\hbar/\left(M\omega_{r}\right)}=1.84$ $\mathrm{\mu}$m, $\rho_{0}=19.23$ $\mathrm{\mu}$m, and $g=1.88\cdot 10^{-2}$. Subsequently, we address the evolution of the merging rings for an elongated trap with $A<1$. Scaled parameters of the potentials in Eqs. (3) and (4) are fixed to be $a=0.3,u_{b}=80$. In the text following below, we use the same notation for the scaled wave function $\psi$, spatial coordinates $\left(x,y,z\right)$, and time $t$ as above, as it will produce no confusion.

A steady-state solution of Eq. (5) was found numerically by dint of the imaginary-time-propagation method, applied with $\gamma=0$. To obtain stationary rings, with different vortex phase profiles in the upper and lower rings, the evolution in imaginary time was initiated by the following ansatz:

[TABLE]

where $\phi$ is the polar angle in cylindrical coordinates, and different integer topological charges $m_{1}$ and $m_{2}$ are imprinted in the top and bottom rings: $S(z)=m_{1}$ for $z<z_{0}$ and $S(z)=m_{2}$ for $z\geq z_{0}$, cf. a similar procedure used in the model with a “peanut-shaped” nonlinear trapping potential in Ref. R. Driben, Y. Kartashov, B. A. Malomed, T. Meier, and L. Torner (2014). The imaginary-time propagation converges to steady states with required accuracy for an arbitrary input amplitude $|\Psi_{0}(x,y,z)|$ in Eq. (8) with a fixed norm:

[TABLE]

where $N_{1}$ and $N_{2}$ are scaled numbers of atoms in the bottom and top rings, respectively:

[TABLE]

.

Note that, by shifting center $z_{0}$ of the splitting barrier (3), it is easy to prepare an initial state with dominant population in the ring with topological charge $m_{1}$ ($N_{1}>N_{2}$ for $z_{0}>0$) or $m_{2}$ ($N_{2}>N_{1}$ for $z_{0}<0$), the respective asymmetry parameter being

[TABLE]

It is relevant to mention two-component models with incoherent nonlinear interaction between the components, which conserve the norms in each component separately. Unlike the present setting, such systems readily admit stationary states with different vorticities and different chemical potentials in the components. In particular, systems of this type give rise to stable states with “hidden vorticity”, i.e., ones with opposite vorticities and equal norms in the two components, the total angular moment of the states being zero, as predicted in BEC M. Brtka, A. Gammal, and B. A. Malomed (2010); L. Wen, Y. Qiao, Y. Xu, and L. Mao (2013); A. I. Yakimenko, K. O. Isaieva, S. I. Vilchinskii, and M. Weyrauch (2013); S. Ishino, M. Tsubota, and H. Takeuchi (2013); S. Ishino, M. Tsuboto, and H. Takeuchi (2015); M. Hoashi, Y. Nakamura, and Y. Yamanaka (2016); S.-L. Xu, Z.-Q. Wang, J.-R. He, L. Xue, and M. R. Belic (2018); Li et al. (2018); Chen et al. (2019) and optics D. Leykam, B. Malomed, and A. S. Desyatnikov (2013); J. R. Salgueiro (2016); T. Mayteevarunyoo, B. Malomed, and D. Skryabin (2018).

III The evolution of merging persistent currents

The dynamics of merging toroidal condensates in real time was simulated by means of the split-step fast-Fourier-transform method. The barrier strength $U_{\text{b}}(t)$ in potential (2) decreasing in time until extinction at $t=t_{d}$, two rings tend to merge into a single toroidal condensate at $t>t_{d}$.

We considered different regimes of the dissipative evolution, controlled by parameter $\gamma$ in Eq. (1). In the conservative limit, with $\gamma=0$ and a constant trapping potential, the total norm, which is defined by Eqs. (9) and (10), energy

[TABLE]

and the $z$-component of the total angular momentum,

[TABLE]

are conserved in our simulations with high accuracy. As seen in Fig. 1 in non-dissipative case ($\gamma=0$), the total energy decays at $0<t<t_{d}$, whilst the potential barrier between the two rings gradually vanishes, and remains constant at $t>t_{d}$.

The system described by Eq. (5) with $\gamma>0$ conserves neither the energy nor the number of particles. The chemical potential $\mu(t)$ of the equilibrium state in our dynamical simulations was adjusted at each time step so that the number of condensed particles slowly decays with time: $N(t)=N(0)$exp$(-t/t_{0})$ , where $t_{0}=10$ s corresponds to the $1/e$ lifetime of the BEC reported in the experiment S. Beattie, S. Moulder, R. J. Fletcher, and Z. Hadzibabic (2013). In our simulations we assume that total number of atoms in the initial state is $N(0)=6\cdot 10^{5}$.

The goal of the numerical analysis is two-fold. First, we aim to trace the 3D dynamics of quantum vortices, in the course of the relaxation following the vanishing of the separating barrier. The second objective is to identify final states of the dissipative evolution for different values of the initial asymmetry parameter (11).

III.1 Pancake-shaped toroidal condensates

Figures 2 and 3 provide typical examples of merging pancake-shaped toroidal condensates [with aspect ratio $A=4.88$, see Eq. (7)], in the cases when the initial atomic population in the bottom ring, with vorticity $m_{1}=1$, is, respectively, weakly or strongly dominant [$P=0.06$ and $P=0.49$, see Eq. (11)] over the top one, with $m=0$. The lower ring, carrying the persistent current in the initial state, contains a single-charged vertical vortex line threading the internal toroidal hole [the red dashed lines in Fig. 2 (a) and 3] pha .

The fact that the initial states in the bottom and top rings have different vorticities, $m_{1}\neq m_{2}$, implies that the superflow tunneling across the separating potential creates $|m_{1}-m_{2}|$ JV lines (rotational fluxons) located in the low-density region between the weakly interacting rings Oliinyk et al. (2019). In the course of the evolution, the JV, which is initially oriented horizontally between the rings, splits in vortices oriented vertically along the $z$ axis, see Fig. 2(b,c) and 3(b,c)]. Thus, Figs. 2 and 3 demonstrate that the initial states, with the horizontally-oriented JV and vertically-oriented vortex line, evolve into series of vortices and antivortices, with axes oriented parallel to the $z$-axis. Note that a similar snake instability and decay of the JV into vortex dipoles was also observed in spin-orbit-coupled BECs Gallemí, A. and Guilleumas, M. and Mayol, R. and Mateo, A. Muñoz (2016). Subsequent evolution of the vortices depends on two factors: (i) interaction between adjacent vortices (attraction for a vortex-vortex pair and repulsion for an vortex-antivortex one), and (ii) drift of the vortex line caused by the dissipation. Under the action of the weak dissipation, the vortex lines spiral out from the bulk of the atomic cloud towards low-density regions at the surface of the condensate. This may imply slow drift of the vortex line either to the external surface of the toroidal condensate, where the vortex decays into elementary excitations, or trapping of the vortex line in the central toroidal hole, which adds a positive or negative unit to the topological charge of the persistent current.

One can see from the Fig. 4, that for the initial state with weak asymmetry ($P=0.06$, $z_{0}=0.02$), the evolution of which is displayed in Fig. 2, the angular momentum per particle, $L_{p}=L_{z}/N$, vanishes as a result of the relaxation [the total angular moment is defined in Eq. (13)]. This observation agrees with that outlined in the caption to Fig. 2, viz., that, under the action of the dissipation, the vortex line drifts in Figs. 2(c,d) to the peripheral (external) low-density area and decays, while the antivortex drifts towards the central $z$ axis and eventually annihilates with the central vortex (shown by dashed red line). In this case the system relaxes to the irrotational ground state, as seen in Fig. 2(d).

In the case of strong asymmetry ($P=0.49$, $z_{0}=0.175$) the value of $L_{p}$ evolves towards $L_{p}=1$, see Fig. 4 (b). As seen in Fig. 3(c,d) and is outlined in the caption to the figure, the fluxon bends and splits into a vortex-antivortex dipole, which escape from the condensate. The final state contains the central vortex (designated by the red dashed line) corresponding to the single-charged ($m=1$) persistent current, as seen in Fig. 3(d).

Systematic simulations demonstrate that there is a threshold value of the initial asymmetry parameter, $P_{\text{cr}}$ [defined as per Eq. (11)], in the pancake-shaped double-ring system, which determines the final state of the toroidal condensate after the relaxation (see the inset in Fig. 4). After long-time dissipative evolution, the toroidal condensate relaxes to the state with zero angular momentum at $P<P_{\text{cr}}$, while at $P>P_{\text{cr}}$ the dominant component (with $m_{1}=1$) imposes its quantum state onto the less populated ring (with $m_{2}=0$), so that the entire condensate relaxes to the state with $m=1$. Summarizing results of the simulations, we conclude that the critical initial asymmetry takes values in interval $0.20<P_{\text{cr}}<0.22$.

The conclusion that the more populated ring, carrying the persistent current with $m=1$, can drag an initially zero-vorticity one into the same vortex state with $m=1$ does not contradict the energy balance (in the dissipationless system, with $\gamma=0$). The potential energy of the condensate, trapped by oscillatory potential in $z$-direction, decreases when the horizontal repulsive separating barrier is eliminated. Thus the total energy of the initial double-ring system is greater than the energy of the final single-ring condensate state, both for $L_{p}=0$ and $L_{p}=1$. As is known, in toroidal condensates a persistent current with topological charge $0\leq m\leq m_{\text{max}}$ corresponds to the stable state with a local energy minimum, where the maximum topological charge $m_{\text{max}}$ admitting the stability depends on the condensate parameters (see, e.g., Moulder et al. (2012); Yakimenko et al. (2015)). Thus, by reducing the potential barrier between the rings with charges $m_{1}=1$ and $m_{2}=0$, one triggers the relaxation of the appearing single toroidal condensate either to the zero-vorticity ($L_{p}=0$) ground state, or to the stable state with $L_{p}=1$, depending on the initial value of imbalance $P$, defined as per Eq. (11).

We have performed a similar analysis for the double-ring system with initial topological charges $(m_{1}=2,m_{2}=0)$. Figure 5 displays the evolution of the angular momentum per particle for different initial values of the imbalance, $P$, in this case. Note that, at $P>P_{\text{cr}}\approx 0.29$, the system relaxes to the state with double-charged ($m=2$) persistent current. Note that the final state with $m=1$ is not observed in the present case, even for the symmetric input with $P=0$ and $L_{p}=1$, which is closest to a possible outcome with $m=1$. The point is that, due to the symmetry of the system, two fluxons featured by the initial configuration split into vertical vortices which either escape towards the external periphery (which leads to the final state with $m=0$), or both vortex cores are trapped together in the central toroidal hole, which corresponds to the final state with $m=2$. Detailed consideration of spontaneous symmetry breaking in the course of the relaxation process after the merger of superflows with different vorticities may be a relevant extension of the present work.

It is interesting to compare our results for the 3D merging rings (where the quantized vortices drive the relaxation dynamics) and the merger of two concentric 2D BECs, studied in Refs. Kanai et al. (2018, 2019), which reported the emergence of a spiral dark soliton when one condensate carried a nonzero initial angular momentum. The spiral dark soliton enables the transfer of the angular momentum between the condensates, allowing the merged condensate to rotate even in the absence of the overall vorticity.

III.2 Elongated toroidal condensates

Our next objective is to compare the evolution of the rings elongated in the vertical ($z$) direction, with small values of aspect ratio (7), and those supported by the pancake-shaped toroidal traps with large values of $A,$ which were considered above. Thus, we aim to figure out the role of the trapping-potential geometry in the JV dynamics.

Note that a strongly elongated condensate is related to the 2D setup considered in Ref. A.W. Baggaley, N.G. Parker (2018), thus one may expect formation of a vortex lattice similar to one observed in the 2D single-component setting Ohta et al. (2010); A.W. Baggaley, N.G. Parker (2018), or a vortex sheet observed in two-component K. Kasamatsu, M. Tsubota (2009) condensates. The lattice of 2D vortices in the 3D setup corresponds to fluxons at an interface of superflows with different vorticities.

The simulations clearly demonstrate that 3D hybrid states, built of two axially separated parts with different vorticities, emerge as a result of the evolution of the merging elongated rings with different vorticities. A typical example is displayed in Fig. 6, for unequal vorticities $m_{1,2}=\pm 1$, aspect ration $A=0.16$ [see Eq. (7)], and zero initial imbalance between the two parts, $P=0$ [see Eq. (11)]. Strictly speaking, the emerging hybrid states are transient ones, as they will eventually transform into usual ones. However, they remain robust in the course of long evolution, and do not decay even being strongly perturbed; instead, perturbed hybrids relax back to the unperturbed shape. velocity of the vertical drift of the fluxons is determined by the dissipation rate parameter $\gamma$ and gradient of the condensate in the $z$-direction: the greater the dissipation, the faster fluxons drift to the region with lower condensate density. Since the condensate’s density distribution is nearly uniform halfway between two elongated condensates, which inhibits degradation of the hybrid. A typical example of hybrid decay is shown in in Supplemental video sup , where the lifetime of hybrid with particles imbalance $P=0.1$ is about $0.7s$. The lifetime of the hybrids in our simulations depends on the imbalance $P$ and for the symmetric setup $P\approx 0$ a perturbed hybrid can survive during several seconds, which is comparable with duration of typical experiment.

Families of the long-lived hybrids may be well characterized by the respective dependencies of the total number of particles on the chemical potential, and of the total energy on the number of particles, as shown in Fig. 7.

These findings are in sharp contrast with the results of 2D simulations reported in Ref. A.W. Baggaley, N.G. Parker (2018), where unstable evolution of the vortex street led to the onset of the KHI-driven dynamics. Note that the decay of the vortex lattice and development of KHI was demonstrated in Ref. A.W. Baggaley, N.G. Parker (2018) for $|m_{1}-m_{2}|=20$, hence the respective vortex street contained $40$ vortex cores. As pointed out in A.W. Baggaley, N.G. Parker (2018), in the absence of perturbations the vortex streets survived for a long time; however, under the action of noisy perturbations (whose magnitude was less than one percent of the background wave function) the 2D vortex lattice rolled-up, which was interpreted in as the development of KHI A.W. Baggaley, N.G. Parker (2018).

Thus a natural question arises if KHI actually develops in full 3D system for $|m_{1}-m_{2}|~{}\sim 1$. In our simulations, the lattice of Josephson vortices is not affected by KHI even in the case of relatively strong noisy perturbations, with amplitudes up to $\sim 5\%$ of the unperturbed state.

As shown in Fig. 6, following the merger of the rings with topological charges $m_{1}=1$ and $m_{2}=-1$, vortex lines are horizontally oriented, in sharp contrast to the vertically oriented vortices observed in Figs. 2 and 3. Results of our 3D simulation are in qualitative agreement with findings reported in Ref. A.W. Baggaley, N.G. Parker (2018), which were obtained in the framework of the simplified 2D model for merging flows with vorticities $m_{1}=+1$ and $m_{2}=-1$.

In our simulations the development of the snake-type instability is strongly suppressed for elongated condensates with aspect ratio $A<1$, see Eq. (7). Further, it was found that for the strongly elongated toroidal traps with $A<0.6$ JVs keep the radial orientation in the course of the relaxation. To gain deeper insight into the physical nature of the suppression of the snake instability, it is instructive to perform a simple energetic analysis. As is known, energy of a straight vortex line in BEC is proportional to its length, and grows with the condensate density. Thus, the vortex line oriented in the radial direction has the maximum energy in the $z=0$ plane, while the vertically oriented line attains the maximum energy if it is located in the peak-density region, at distance $\rho_{0}$ from the axis $z$ [see Eq. (6)]. Therefore, the maximum energies of the vertical and horizontal vortex line relate to each other as the radial, $R_{\text{TF}}=\sqrt{8\mu/(M\omega_{r}^{2})}$, and vertical, $Z_{\text{TF}}=\sqrt{8\mu/(M\omega_{z}^{2})}$, Thomas-Fermi widths of the condensate, i.e., as aspect ratio (7) of the toroidal trap: $A=R_{\text{TF}}/Z_{\text{TF}}=\omega_{z}/\omega_{r}$. From these straightforward estimates, it follows that the fluxons are most likely to keep the horizontal orientation when elongated rings merge. In our numerical simulations we observed that this property of elongated condensates does not depend on the number of atoms, but is determined by the aspect ratio of the trap. The robust dynamics of JVs in the elongated trap is in sharp contrast with the evolution of fluxons in the pancake-shaped setting, which, as demonstrated above, bend and eventually split into vertically oriented vortices and antivortices.

The formation of $|m_{1}-m_{2}|$ fluxons at the interface between flows with topological charges $m_{1}$ and $m_{2}$ necessarily follows from the azimuthal periodicity of the condensate wave function. For elongated condensates the horizontal orientation of the JVs is energetically preferable even when the barrier between the rings vanishes, and they are located in the peak-density regions. In the same time, vertically oriented vortex lines remain trapped in the potential well formed by the central toroidal hole. Eventually, this scenario gives rise to formation of long-lived transient states in the form of the above-mentioned hybrids, which are qualitatively similar to those supported by means of a completely different mechanism in Ref. R. Driben, Y. Kartashov, B. A. Malomed, T. Meier, and L. Torner (2014). The hybrids are topological structures with $(m_{1},m_{2})$ vertically oriented vortex lines and $|m_{1}-m_{2}|$ horizontally oriented JVs.

Furthermore, as illustrated in Fig. 8 (see also Supplemental movie sup ), the rotational fluxons can be pinned and completely stabilized if the amplitude of the barrier is reduced to some nonzero value, maintaining a residual potential well with the bottom point at $z=z_{0}$. As the result, one obtains completely stable stationary (rather than transient) 3D hybrid vortex complexes with any value of the angular momentum per particle in the range of $m_{1}<L_{p}<m_{2}$. This value can be readily tuned by choosing position of the barrier, $z_{0}$, in the initial state, and even dynamically tuned by slow motion of the sheet beam after formation of the stable hybrid as illustrated in Fig. 8 for $t>0.3$ s. The stabilized JVs are not only topological modes of fundamental interest in physics of superfluidity and superconductivity, but also may be used for design of topologically protected qubits, cf. Ref. S. Das Sarma, M. Freedman, and C. Nayak (2005).

IV Conclusion

We have demonstrated that the evolution of merging persistent currents in parallel BEC rings is driven by essentially three-dimensional nonlinear dynamics of JVs (Josephson vortices, alias rotational fluxons). Quite different dynamical scenarios are reported for the fluxons, which give rise to different final states of the condensate, depending on the imbalance in initial populations of the two rings, and on the aspect ratio of the toroidal trap. In sharp contrast to classical flows with viscosity, in the considered weakly dissipative quantum system one can observe transitions into states with higher-order values of the angular momentum. It turns out that for elongated merging condensates the vorticity of the final state is always determined by the more populated ring component, while for pancake-shaped rings there is a threshold value of imbalance, such that below it the final vorticity is zero, while above the threshold the dominant component imparts its vorticity to the final state.

In the course of the merger of the rings carrying different initial vorticities, which are strongly elongated in the axial direction, we observed formation of dynamical hybrids, which are long-lived transient states in 3D, with different vorticities, $m_{1}$ and $m_{2}$, in their two axially separated parts, which share a common vertically oriented axis threading the merged double-ring system and $|m_{1}-m_{2}|$ radially oriented JVs (rotational fluxons). Under the action of weak dissipation, which drives slow vertical drift of the fluxons, the dynamical hybrid suffers slow degradation into the state with the vorticity imposed by the initially dominant component. The JVs can be pinned and stabilized by a residual nonvanishing potential barrier separating the parallel rings, which thus leads to the formation of completely stable 3D hybrid vortex complexes.

In addition to revealing the dynamics of the merger of superfluid rings with different vorticities, and possibilities of the formation of transient and permanent hybrids incorporating different vorticities in their axially separated parts, these findings may also help to understand the correspondence between rotational fluxons in long Bose-Josephson junctions and the formation of vortex lattices at the interface of counter-propagating superflows.

Acknowledgment

The work of B.A.M. was partly supported by the Israel Science Foundation through grant No. 1287/17.