Critical transport and vortex dynamics in a thin atomic Josephson junction

TL;DR

This paper investigates how vortex ring nucleation causes dissipation in a thin atomic Josephson junction between Fermi superfluids, linking microscopic vortex dynamics to macroscopic transport properties.

Contribution

It provides a detailed simulation-based analysis of vortex ring nucleation and its role in dissipation, aligning with experimental observations and elucidating the microscopic mechanisms.

Findings

Identification of critical population imbalance and maximum Josephson current for dissipation onset.

Demonstration that quantum phase slips via vortex rings cause resistive current.

Quantitative agreement with experimental data on vortex dynamics and transport.

Abstract

We study the onset of dissipation in an atomic Josephson junction between Fermi superfluids in the molecular Bose-Einstein condensation limit of strong attraction. Our simulations identify the critical population imbalance and the maximum Josephson current delimiting dissipationless and dissipative transport, in quantitative agreement with recent experiments. We unambiguously link dissipation to vortex ring nucleation and dynamics, demonstrating that quantum phase slips are responsible for the observed resistive current. Our work directly connects microscopic features with macroscopic dissipative transport, providing a comprehensive description of vortex ring dynamics in three-dimensional inhomogeneous constricted superfluids at zero and finite temperatures.

Equations

$$V_{\mathrm{trap}}(x,y,z)=\frac{1}{2}M({{\omega }_x}^2{x}^2+{{\omega }_y}^2{y}^2+{{\omega }_z}^2{z}^2)+V_0\cdot e^{\frac{-2x^2}{w^2}}$$

$$i\hslash \frac{\partial \psi (\vec{r},t)}{\partial t}=\left(-\frac{{\hslash }^2}{2M}{\nabla }^2+V_{\mathrm{trap}}+g|\psi (\vec{r},t){|}^2\right)\psi (\vec{r},t)$$

$$\mu \psi (\vec{r})=\left(-\frac{{\hslash }^2}{2M}{\nabla }^2+V_{\mathrm{trap}}+g|\psi (\vec{r}){|}^2\right){\psi }_0(\vec{r})$$

$$\mu =\frac{1}{2}\hslash \overline{\omega }{\left(\frac{15N_{\mathrm{BEC}}{a}_M}{\overline{l}}\right)}^{2/5}$$

$$\psi (\vec{r},t)={\psi }_L(t)\cdot {\eta }_L(\vec{r})+{\psi }_R(t)\cdot {\eta }_R(\vec{r})$$

$$\int {\eta }_i\cdot {\eta }_{j}d\vec{r}={\delta }_{i,j}$$

$$z_{BEC}=\frac{N_L-N_R}{N}$$

$$E=\frac{UN}{2}{z}^2-2K\sqrt{1-z^2}\cos \left(\Delta \varphi \right).$$

$$E=\frac{UN}{2}{z}^2-2K\cos \left(\Delta \varphi \right).$$

$$z_{cr}=\sqrt{\frac{8K}{UN}}.$$

$$2K=\frac{\left(E_a-E_s\right)}{N}.$$

$$U=g\int {\eta }_L^{4}d\vec{r}$$

$$U_{NL}=2\frac{\partial \mu }{\partial N}$$

$${\omega }_J=\frac{\sqrt{E_J{E}_c}}{\hslash }=\frac{\sqrt{E_{J}UN}}{\hslash }.$$

$$i\hslash \frac{\partial \psi }{\partial t}=\left[-\frac{{\hslash }^2{\nabla }^2}{2M}+V_{\mathrm{trap}}+g(|\psi {|}^2+2n_{\mathrm{th}})\right]\psi .$$

$$\mu (T){\psi }_0=\left(-\frac{{\hslash }^2}{2M}{\nabla }^2+V_{\mathrm{trap}}+g(|{\psi }_0{|}^2+2n_{\mathrm{th}}^0)\right){\psi }_0$$

$$\frac{\partial f}{\partial t}+\frac{\vec{p}}{M}\cdot {\nabla }_{\vec{r}}f-{\nabla }_{\vec{r}}{V}_{eff}\cdot {\nabla }_{\vec{p}}f=0$$

$$n_{th}=\frac{1}{(2\pi \hslash {)}^3}\int d\vec{p}f(\vec{p},\vec{r},t).$$

$$I={\int }_{R_{TF}^z}{\int }_{R_{TF}^y}{j}_x(x=0,y,z)dydz,$$

$$j=\frac{\hslash }{2iM}({\psi }^{\ast }\nabla \psi -\psi \nabla {\psi }^{\ast })$$

$$I={\int }_{R_{\mathrm{TF}}^z}{\int }_{R_{\mathrm{TF}}^y}|\psi (0,y,z){|}^2\cdot v_x(x=0,y,z)dydz$$

$$⟨v_x⟩=\frac{\int v_x\cdot |\psi (x=0,y,z){|}^{2}dydz}{\int |\psi (0,y,z){|}^{2}dydz},$$

$$I={\rho }_x⟨v_x⟩,$$

$$I(\phi )=I_c\sin (\phi )+J_1\sin (2\phi )\equiv I_c\cdot (\sin (\phi )+\bar{g}\cdot \sin (2\phi ))$$

$$\frac{I_{\max }}{I_c}=\frac{(\sqrt{1+32{\bar{g}}^2}+3{)}^{3/2}(\sqrt{1+32{\bar{g}}^2}-1{)}^{1/2}}{32|\bar{g}|}$$

$$⟨{\mu }_L-{\mu }_R⟩=h⟨\frac{dn}{dt}⟩$$

$$\mathbf{v}={\mathbf{v}}_f+{\mathbf{v}}_{\omega }.$$

$${\mathbf{v}}_f(O,t)=\mathbf{v}(O,t)-{\mathbf{v}}_{\omega }(O,t),$$

$${\mathbf{v}}_{\omega }(O,t)=-\frac{\kappa }{4\pi }\sum _{i=1}^N{\oint }_{{}_{{\mathcal{C}}_i(t)}}\frac{{\mathbf{s}}^{\prime }(\zeta ,t)\times \mathbf{s}(\zeta ,t)}{|\mathbf{s}(\zeta ,t){|}^3}d\zeta ,$$

$$E=\int \left[e_{int}(\vec{r},t)+e_k(\vec{r},t)+e_V(\vec{r},r)\right]d^3\vec{r}$$

Full text

Interest is growing in model systems that allow for investigating the interplay between resistive and dissipationless quantum transport phenomena. In this context, ultracold gases in tailored optical potentials represent an ideal framework, owing to the real-time control over the relevant parameters in experiments [1, 2], combined with the ability for ab initio modelling [3, 4]. A paradigmatic example is the study of the dynamics between two atomic superfluids weakly coupled through a thin tunnelling barrier. This realizes a Josephson junction [5, 6], which represents a minimal platform to observe both coherent quantum transport [7, 6], and its breakdown driven by dissipative microscopic mechanisms [8, 9].

The coherent dynamics of atomic Josephson junctions (JJs) [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20] is governed by the competition between the charging energy $E_{C}$ and the Josephson tunneling energy $E_{J}$ [10, 11]. $E_{C}$ relates the chemical potential difference across the tunnelling barrier to the relative population imbalance between the reservoirs, and depends on interparticle interactions. $E_{J}$ promotes the delocalization of the superfluid across the two reservoirs and sets the maximum coherent flow through the weak link. When $E_{J}$ dominates, superfluid current and relative phase oscillate in quadrature at the Josephson plasma frequency. In the opposite regime, and in the absence of dissipation [11, 12], the system may enter the Macroscopic Quantum Self-Trapping (MQST) regime. This is characterized by high-frequency coherent oscillations of the population imbalance around a non-zero value, driven by a monotonically increasing relative phase [10, 13, 15, 16, 17, 18]. Even without thermally induced decay of the population imbalance [17, 12, 21], the stability of MQST depends on whether vortices nucleated inside the barrier annihilate therein [22, 23], or penetrate into the superfluid reservoirs. Recent experiments with inhomogeneous three-dimensional Fermi superfluids [24, 25] revealed the intimate connection between phase slippage and dissipation arising from vortices created within the barrier and shed into the superfluid. Similar effects have been studied in ring-shaped bosonic condensates [26, 27, 28, 29], mesoscopic structures [30, 31] and lower-dimensional geometries [32, 33, 34]. While vortices crossing the weak link are known to yield a finite resistance [27, 30, 25], the relation between microscopic vortex nucleation, dynamics and macroscopic dissipative flow is still poorly understood.

In this work we demonstrate the connection between resistive superfluid currents and vortex ring (VR) dynamics in an atomic JJ of fermionic superfluids. We obtain the critical population imbalance and the maximum coherent current delimiting the boundary between dissipationless and dissipative transport even at finite temperatures. We find excellent agreement with recent measurements [24, 25], thus clarifying their interpretation. Trap asymmetry is shown to foster the emergence of elliptical VRs exhibiting Kelvin wave excitations, while thermal fluctuations reduce the VR lifetime.

Methodology. Our numerical simulations are based on the experimental parameters of Ref. [25]. We consider two molecular Bose-Einstein condensates (BEC) of about $10^{5}$ atom pairs of 6Li, weakly coupled through a thin optical barrier, at $1/k_{F}a\sim 4.6$ (where $k_{F}$ is the Fermi wavevector and $a$ the interatomic scattering length). The harmonic trapping potential is asymmetric with approximately (1:12:10) ratio along the $x$, $y$ and $z$-axis, respectively. The Gaussian barrier bisects the gas along the weakest ($x$) direction (with $\nu_{x}=15$ Hz), with a $1/e^{2}$ waist, $w\sim 4\xi$, where $\xi$ is the superfluid coherence length [25]. The superfluid transport through the barrier is triggered by an initial non-zero population imbalance $z_{0}=z_{\mathrm{BEC}}(0)$ between the two reservoirs. Here, $z_{\mathrm{BEC}}(t)=(N_{R}(t)-N_{L}(t))/N_{\mathrm{BEC}}$, with $N_{L}\,(N_{R})$ the BEC number in the left (right) reservoir, and $N_{\mathrm{BEC}}=N_{L}+N_{R}$ the total condensate number. The imbalance corresponds to a chemical potential difference $\Delta\mu=\mu_{L}-\mu_{R}=E_{C}z_{0}N_{\rm BEC}/2$. Dynamics in the $T=0$ limit are simulated via the time-dependent Gross-Pitaevskii equation (GPE), extended to non-zero temperatures $T\lesssim 0.4\,T_{c}$ (where $T_{c}$ is the BEC critical temperature), via its coupling to a collisionless Boltzmann equation [3, 37, 38]. We stress that the standard two-mode model [10, 13] that captures both Josephson and MQST dynamics of previous experiments [18, 16] is out of its validity range due to the considered values of the ratio $V_{0}/\mu$ and to the thinness of the junctions [38]. Although dissipative effects can be phenomenologically modeled by damped two-mode [35, 21, 36] and RSJ-circuital models [7], such approaches provide limited insight into the microscopic dissipative processes.

Dynamical Regimes and Phase Diagram. We study $z_{\mathrm{BEC}}(t)$, varying both the initial population imbalance $z_{0}$ and the barrier height $V_{0}$. At each value of $V_{0}$ we observe two distinct dynamical regimes. For $z_{0}$ smaller than a critical value $z_{\mathrm{cr}}$, $z_{\mathrm{BEC}}(t)$ exhibits sinusoidal plasma oscillations (Josephson regime). For $z_{0}\geqslant z_{\mathrm{cr}}$, we instead observe an initial rapid decay of $z_{\mathrm{BEC}}(t)$ (dissipative regime), followed by plasma oscillations with amplitude smaller than $z_{\mathrm{cr}}$. We validate our numerics by comparing $z_{\mathrm{BEC}}(t)$ with experiments under the same conditions, finding excellent agreement [Fig. 1(a), insets].

Combining calculated and newly extracted experimental $z_{\mathrm{cr}}$ values, we map out the phase diagram delimiting Josephson and dissipative regimes as a function of the normalized barrier height $V_{0}/\mu(T)$ [Fig. 1(a)], where $\mu(T)$ is the chemical potential including the thermal mean-field contribution [37, 38]. Increasing $V_{0}/\mu(T)$, the onset of dissipation appears at smaller $z_{\mathrm{cr}}$. This reproduces the observed boundary within experimental uncertainty up to $T\approx 0.3\,T_{c}$ upon keeping the condensate number equal to the $T=0$ case. Our findings can also be interpreted in terms of the critical current $I_{\mathrm{max}}$ across the junction, defined as the maximum value of $I=\dot{z}_{\mathrm{BEC}}\,N_{\mathrm{BEC}}/2$ at $z_{0}=z_{\mathrm{cr}}$ [Fig. 1(b)]. Numerically, $\mathinner{\!\left\lvert I\right\rvert}_{\mathrm{max}}$ is well approximated by $z_{\mathrm{cr}}\,\omega_{J}\,N_{\mathrm{BEC}}/2$, where $\omega_{J}$ is the Josephson plasma frequency. The corresponding $\mathinner{\!\left\lvert I\right\rvert}_{\mathrm{max}}$ from the experimentally determined $z_{\mathrm{cr}}$ and $\omega_{J}$ is in excellent agreement with the theoretical prediction. The overall trend of $\mathinner{\!\left\lvert I\right\rvert}_{\mathrm{max}}$ against $V_{0}/\mu(T)$ is also quantitatively captured by extending to inhomogeneous systems an analytical model, originally developed for two homogeneous BECs weakly coupled through a rectangular barrier [20, 38, 59].

Vortex Ring Nucleation and Evolution. The onset of the dissipative regime for $z_{0}\geq z_{\mathrm{cr}}$ in [24, 25] has been linked to the appearance of topological defects in the superfluid. Here, we fully characterize such dynamical features at $T=0$, by computing the superfluid velocity $\mathbf{v}=(\hbar/M)\nabla\phi$, where $M$ is the atom pair mass and $\phi$ the condensate phase (see later for thermal effects). Given the symmetry of our junction, we consider the $x$-component of the superfluid velocity, weighted over the transverse density in the $x=0$ plane, $\langle v_{x}\rangle$ [38]; we identify three distinct dynamical stages (I, II and III, see Fig. 2) in the nucleation process of the first VR (the emerging pattern applies to subsequent VRs). In stage I, following the Josephson-Anderson relation $M\dot{\bf{v}}=-\nabla\mu$ [60, 61, 62], the chemical potential gradient $\nabla\mu$ drives an unidirectional, accelerated superfluid flow across the junction. When $-\langle v_{x}\rangle$ reaches a critical value exceeding the mean sound speed $\langle c\rangle$ (Fig. 2(a)), a VR is nucleated, associated with a relative-phase jump of $2\pi$ [38]. It originates outside the Thomas-Fermi (TF) surface on the central radial plane, $x=0$, where, the superfluid velocity is maximum, due to the flow constriction, and the local speed of sound is minimum, since density vanishes. After its nucleation [stage II], the VR moves axially very slowly away from its nucleation region $x_{\mathrm{V\!R}}=0$ (Fig. 2(b)) with its mean radius $R_{\mathrm{V\!R}}$ rapidly shrinking (Fig. 2(c)) due to the strong radial density inhomogeneity in the barrier region, until it is comparable to the transversal TF radius of the BEC, and enters the bulk superfluid [22]. During such evolution, $-\langle v_{x}\rangle$ exhibits a rapid decrease, possibly even changing sign. Then, in stage III, the VR gradually leaves the barrier region with the axial velocity $-\langle v_{x}\rangle$ re-accelerated by $\nabla\mu$ (Fig. 2(a)), until some time later ($\Delta t_{\mathrm{slip}}\simeq h/\Delta\mu$), when it has already travelled a considerable distance from the barrier edge, another VR is nucleated at the trap centre (see, e.g. Fig. 2(e)). Note that early on in stage III, before the VR exits the barrier region (i.e. before reaching the point of maximum transversal TF radius), $R_{\mathrm{V\!R}}$ continues decreasing due to the strong background density gradient.

For a deeper insight into the underlying superfluid dynamics, we decompose at $\bf{x=0}$ the total axial superfluid velocity $v_{x}=v_{f}+v_{\omega}$, where $v_{f}$ is the main flow velocity (which is slowly varying compared to the timescale of the early vortex dynamics) and $v_{\omega}$ is the VR-generated swirling velocity [38]; we also initially neglect compressibility effects (addressed in the next paragraph). By the end of stage II, the shrinking VR has just left the trap centre, and so the vortex contribution $v_{\omega}$ evaluated at $\bf{x=0}$ (where the local superfluid velocity $\langle v_{x}\rangle$ shown in Fig. 2(a) is calculated) tends to 0. This leads to a drop of $\langle v_{x}\rangle$ with amplitude $\Delta\langle v_{x}\rangle\sim\kappa/R_{\mathrm{V\!R}}$, corresponding to the change in the axial superfluid velocity at the trap centre due to the lost vortex contribution, where $\kappa$ is the quantum of circulation [38]. This sawtooth-like profile of $\langle v_{x}\rangle$ (Fig. 2(a)) is typical of phase slippage phenomena seen in superfluid helium [60, 61, 64, 62], with the less abrupt drop found here stemming from the initial persistence of the VR within the barrier region.

The drop $\Delta\langle v_{x}\rangle$ can even overcome the generating flow velocity, leading to flow reversal (i.e. backflow) in the post-nucleation dynamics, in agreement with Biot-Savart calculations [38]. The amplitude of each subsequent velocity drop is reduced due to the overall decay of $z_{\mathrm{BEC}}(t)$.

To connect the dissipation with the microscopic VR nucleation and dynamics, and phonon emission, we decompose the temporal evolution of the BEC total kinetic energy in its incompressible $E_{k}^{i}$, compressible $E_{k}^{c}$, and quantum pressure $E_{q}$ contributions [63, 38]. $E_{k}^{i}$ and $E_{k}^{c}$ correspond respectively to the kinetic energy of the flow (both potential flow driven by $\nabla\mu$ and vortex generated swirls) and to the sound wave energy in the superflow. $E_{q}$ accounts for the energy arising from density inhomogeneities due to the trapping potential [38]. When each VR enters the TF surface (end of stage II), and while still propagating within the barrier’s region of increasing density, sound waves are emitted and $E_{k}^{c}$ increases at the expenses of $E_{k}^{i}$ (Fig. 2(d)) [38]. The dissipation of Josephson oscillations [60, 61, 62, 25] thus stems from two effects: the incompressible kinetic energy transferred from the axial flow to the vortex swirling flow and the phonon-emission occurring during vortex nucleation and propagation within the barrier region.

We further quantify both those effects by considering the effect of $z_{0}$ on the velocity $v_{\mathrm{V\!R}}$ and the incompressible kinetic energy $E_{k,V\!R}^{i}$ of the first VR nucleated; this is shown in Fig. 3 for $z_{0}\in[0.13,\,0.37]$ and $V_{0}/\mu=0.8$. We find that increasing $z_{0}$ leads to a decreasing $v_{\mathrm{V\!R}}$ and to a monotonic increase of $E_{k,V\!R}^{i}$ [Fig. 3(a)] [38]. Calculating the fraction of the total kinetic energy flowing through the junction until the nucleation of the first VR which is dissipated in sound ($\epsilon_{c}$) or transferred to the first VR ($\epsilon_{i}$) [38], we find that both sources of dissipation increase as $z_{0}$ gets larger [see Fig. 3(b)], and can cumulatively account for a significant fraction of the total energy. Surprisingly, the acoustic dissipation $\epsilon_{c}$ is always larger than the incompressible contribution $\epsilon_{i}$.

Increasing $z_{0}$ leads to more nucleated vortices $N_{\mathrm{V\!R}}$ (Fig. 3(b)), due to the larger time-averaged chemical potential difference [60], consistently with [25] and with earlier studies of controlled vortex generation [68, 69]. Similarly, the VR lifetime increases by increasing $z_{0}$ [Fig. 3(b) (inset)]. The VR survival during its propagation in the superfluid bulk is thus determined by two competing effects: On the one hand, the VR tends to expand [70] to conserve its incompressible kinetic energy as it moves towards lower-density regions with decreasing transverse size. On the other hand, the radial trapping asymmetry ($\omega_{y}\neq\omega_{z}$) leads to elliptical VR profiles with oscillating aspect ratio, corresponding to a $m=2$ Kelvin wave excitation on a circular VR [71]. This wobbling motion induces dissipation of the VR incompressible kinetic energy via the emission of phonon-like excitations, reducing the VR radius [72, 73]: When $R_{\mathrm{V\!R}}\sim\xi$, the VR loses its circulation and annihilates in a rarefaction pulse [22].

This picture remains qualitatively correct for the probed $T\lesssim 0.4\,T_{c}$ (red symbols in Figs. 2-3), for which no thermal VR activation occurs ($k_{B}T<0.8V_{0}$). For a fixed BEC number, superfluid flow, VR generation and early dynamics are not noticeably affected by the thermal cloud, whose main effect is to add an extra potential to the BEC [65, 66, 67, 38, 37]. Over longer timescales, dissipation due to relative BEC-thermal motion becomes relevant, decreasing the VR lifetime (Fig. 3(b)).

To connect with Ref. [24, 25], we implement in our simulations the same protocol by which vortices were observed in time of flight after gradually removing the barrier with a 40 ms linear ramp. The dynamics of the $4^{\rm th}$ VR generated in the same conditions as in Fig. 2 [right VR in Fig. 2(e)] is shown in Fig. 4, including or excluding the barrier removal procedure. Upon removing the barrier (orange curve), the VR propagates for longer time and longer distance [Fig. 4(a)]. This facilitates the direct observation of Kelvin-wave oscillations [visible in Fig. 4(b)], whose period is consistent with the dispersion relation $\omega(k)\sim\kappa k^{2}/(4\pi)[\ln(2/(k\xi)-0.5772]$ [71, 38]. The longer lifetime can be attributed to the larger kinetic energy of VRs nucleated during the gradual barrier removal process. As the VR approaches the edge of the condensate, it breaks up into two anti-parallel vortex lines [Fig. 4(b), final snapshot] [22, 74, 70]. Critically, thermal fluctuations destabilize the VR, causing it to drift off-axis, and reach the transversal boundary asymmetrically [Fig. 4(c)]. There, it reconnects with its image and forms a ‘vortex handle’ [78, 75, 76]. This could explain why a single vortex line is typically detected in each experimental run after barrier removal [24, 25].

Conclusions. We have studied the complex interplay between coherent and dissipative dynamics in a thin atomic Josephson junction. We have shown that resistive currents are directly connected with nucleations of vortex rings and their propagation into the superfluid bulk. In particular, dissipation originates from two irreversible effects: phonon emission when vortex rings are nucleated, and incompressible kinetic energy transfer from the superfluid flow to the swirling one of the nucleated vortex rings. The detailed understanding of the connection between vortex-ring dynamics and dissipation is valuable for advancing our comprehension of the complex superfluid dynamics in emerging atomtronic devices [79].

Acknowledgments. We thank A. Smerzi and A. Muñoz Mateo for valuable discussions. This work was supported by QuantERA project NAQUAS (EPSRC EP/R043434/1), EPSRC project EP/R005192/1, Fondazione Cassa di Risparmio di Firenze project QuSim2D 2016.0770, European Research Council grant agreement no. 307032 QuFerm2D and no. 637738 PoLiChroM, and European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 705269.

S.1 Numerical Methods

We study the superfluid transport of the molecular BEC through the thin barrier using two different models. Specifically, the $T=0$ dynamics is modelled by the Gross-Pitaevskii equation. At finite temperatures we instead make use of a collisionless kinetic model, in which the condensate dynamics is self-consistently coupled to a dynamical thermal cloud described by a Boltzmann equation.

S.1.1 Experiment Overview and System Parameters

Superfluids of $N\simeq 10^{5}$ atom pairs of 6Li are produced by cooling a balanced mixture of the two lowest spin states $|F=1/2,m_{F}=\pm 1/2\rangle$ to $T/T_{c}\sim 0.3(1)$ [39, 24, 25]. Interactions between fermions are parametrized by $1/(k_{F}a)$, where $k_{F}=\sqrt{2mE_{F}}/\hbar$ is the Fermi wave-vector ($m$ is 6Li atomic mass and $E_{F}$ the Fermi energy), and $a$ is the interatomic tunable $s$-wave scattering length. The focus of this work is on the regime of superfluidity of the molecular BEC, and we restrict our modelling to the case of $1/(k_{F}a)\simeq 4.6$ . To realize an atomic Josephson junction the fermionic superfluid is separated into two weakly-coupled reservoirs by focusing onto the atomic cloud a Gaussian-shaped repulsive sheet of light, yielding a trapping potential

[TABLE]

where $\omega_{x,y,z}$ are the trapping frequency along $x$, $y$ and $z$-directions, $M=2m$ is the molecular mass, $V_{0}$ is the height of the Gaussian barrier and $w\approx 2.0\pm 0.2\,\mu$m is the barrier $1/e^{2}$ width, which is just four times wider than the superfluid coherence length $\xi$. The experimental trap frequencies are $\omega_{x}\simeq 2\pi\times 15$ Hz, $\omega_{y}\simeq 2\pi\times 187$Hz, $\omega_{z}\simeq 2\pi\times 148$Hz (cigar-shaped trap), and $V_{0}$ is varied in the regime of $0.6\,\mu\lesssim V_{0}\lesssim 1.2\,\mu$ where $\mu$ denotes the chemical potential of the system.

S.1.2 Gross-Pitaevskii Equation ($T=0$)

At $T=0$ we model the system by the molecular BEC wavefunction $\psi$ obeying the time-dependent Gross-Pitaevski equation (GPE):

[TABLE]

where $g=4\pi\hbar^{2}a_{M}/M$ is the interaction strength, and $a_{\mathrm{M}}=0.6\,a\simeq 7\times 10^{-3}\,l_{\mathrm{x}}$ is the molecular scattering length. The equilibrium state is found by substituting $\psi(\vec{r},t)=\psi_{0}(\vec{r})\exp{(-i\mu t/\hbar)}$ which gives the time-independent GPE:

[TABLE]

with $\mu$ the system chemical potential. The equilibrium state is obtained numerically via imaginary time propagation in the presence of an additional linear potential, $-\epsilon x$, along the $x$-direction which sets up the desired initial population imbalance, $z_{0}$, between the two wells (i.e. initial chemical potential difference). As there is initially a larger population in the right well, the initial flow is induced along the negative $x$-direction. The BEC dynamics instead is initiated by the instantaneous linear potential removal at $t=0$. Eq. S.2 is studied in dimensionless form, with length scaled to the harmonic oscillator length along the $x$-direction, $l_{\mathrm{x}}=\sqrt{\hbar/m\omega_{\mathrm{x}}}\simeq 7.5\,\mu$m. In our numerical simulations we use grid sizes $\left[-24,24\right]l_{x},\left[-4,4\right]l_{x},\left[-4,4\right]l_{x}$ along the $x$, $y$ and $z$-directions, and $1024\times 128\times 128$ grid points respectively. Throughout this work, the barrier width is set to $w\simeq 2\mu$m $\approx 4\xi$, with $\xi=1/\sqrt{8\pi a_{M}n_{\mathrm{max}}}\approx 0.5\,\mu$m$\,\simeq 0.067\,l_{\mathrm{x}}$. To account for experimental numbers, Fig. 1(a) of the main paper considered molecule numbers in the range $(6-12)\times 10^{4}$. Throughout this Supplemental Material instead, we fix the condensate number at $N_{\mathrm{BEC}}=6\times 10^{4}$. In analyzing our results, we express the barrier height $V_{0}$ in units of the system chemical potential $\mu$. The numerically-extracted equilibrium $\mu$ is well approximated by the analytical formula in the Thomas-Fermi approximation:

[TABLE]

with $\overline{\omega}=(\omega_{x}\omega_{y}\omega_{z})^{1/3}$ and $\overline{l}=\sqrt{\hbar/M\overline{\omega}}$ the geometric mean of harmonic oscillator lengths [40]. For typical parameters ($T=0$, $N_{\mathrm{BEC}}=6\times 10^{4}$ molecules) $\mu\simeq 114\hbar\omega_{x}$ and the barrier height numerically explored is in the range $[0.6,1.22]\,\mu$.

S.1.3 The Two-mode model

In this section we discuss the validity of the commonly-used two-mode model for the specific geometry and parameter regime of our system.

The two-mode model cannot capture the observed dynamical regimes, since – with increasing population imbalance – it always predicts a transition from the Josephson oscillation regime to the Macroscopic Quantum Self-Trapping (MQST) regime, in which the population imbalance oscillates around a non-zero value.

However, as already clearly demonstrated in the main paper, for the experimental parameters probed here – namely a narrow barrier of width $w/\xi=4$ and a relatively low barrier height $V_{0}/\mu$ in the range [0.6, 1.22] – we have instead numerically found a transition from the Josephson regime (for small initial population imbalance) to a distinct, dissipative, regime associated with defect generation. Typical dynamics of the population imbalance in these two regimes are shown in the two insets to Fig. 1(a) of the main paper. Such numerical results are confirmed by their excellent agreement with the experimental data over the entire range of barrier height explored (see Fig. 1(a)).

Even though the two-mode model fails to capture the nature of the dynamical regime to which the system transitions with increasing initial population imbalance, it is nonetheless interesting to examine the extent to which such a model can predict either the location of the transition region where Josephson plasma oscillations no longer occur, or the actual frequency of the Josephson oscillations within the Josephson regime.

To address such issues, we start with a brief overview of the two-mode model ([10, 13]).

Assuming the two superfluids are well-localized in each well, the system wavefunction can be written as a linear superposition of the left and right condensate wavefunctions as

[TABLE]

where $\psi_{L}(t)=\sqrt{N_{L}}e^{i\phi_{L}}$ and $\psi_{R}(t)=\sqrt{N_{R}}e^{i\phi_{R}}$. Here

[TABLE]

with $i$, $j$ denote left, right, and $N_{L(R)}$ or $\phi_{L(R)}$ are the left (right) condensate number and phase. This simple model is called the ‘two-mode model’ as in this approximation only the ground state and the first excited state are populated. If we define the population imbalance

[TABLE]

(where $N=N_{L}+N_{R}$), and the phase difference $\Delta\phi=\phi_{L}-\phi_{R}$, the system energy in the two-mode model takes the form

[TABLE]

where $K$ is the tunneling energy and $U$ the onsite interaction energy. For $z_{0}\ll 1$, this can be written as

[TABLE]

The condition for obtaining MQST corresponds to having an initial energy $E_{0}=E(z_{0},\Delta\phi=0)>E_{cr}$, where $E_{cr}=E(z=0,\Delta\phi=\pi)$. Thus the critical imbalance for obtaining MQST is:

[TABLE]

This condition corresponds to populating the first excited state of the system, which is the antisymmetric state. Thus the tunneling energy $K$ can be extracted from the difference between the antisymmetic-state energy, $E_{a}$, and the symmetric-state energy, $E_{s}$, via

[TABLE]

Hereafter we use $E_{J}$, defined as $E_{J}=2K$. The onsite interaction energy $U$ can be calculated either from

the linear two-mode model [10], via

[TABLE]

or from the nonlinear two-mode model [41] as:

[TABLE]

To evaluate the two-mode model prediction for the location of the crossover between Josephson and MQST, we calculate $z_{cr}$ both from the linear and from the nonlinear two-mode model for each $V_{0}/\mu$, and plot it against $V_{0}/\mu$ in the experimentally relevant range $[0.6,1.22]$ in Fig.  S1. We also plot on this figure the corresponding extracted values of $z_{cr}$ from the GPE simulations [noting that in the latter case, the transition is to a different, dissipative, regime, rathen than to MQST]. We find that both linear and nonlinear two-mode model clearly significantly underestimate the numerically predicted critical value of $z_{cr}$ for all $V_{0}/\mu\leq 1.1$, yielding a correct value for the transition only at the highest probed point $V_{0}/\mu=1.22$ (see Fig. S1).

Next, we examine the extent to which the oscillation frequency of the Josephson oscillation is correctly predicted by such two-mode model.

In the regime of the Josephson plasma oscillation, where $z_{0}<z_{cr}$, and for $E_{J}\ll E_{c}$ the two-mode model prediction for the oscillation frequency is:

[TABLE]

Plotting $\omega_{J}$ against $V_{0}/\mu$ in the considered range (Fig. S2), we find the predictions of both linear and nonlinear two-mode model to be incorrect by at least a factor of 2 (for $V_{0}/\mu=1.22$), with the disagreement rapidly increasing for lower values of $V_{0}/\mu$.

The above analysis thus proves that – for the parameter regimess of the experiment under study – the two-mode model is out of its validity range, due to the considered values of the ratio $V_{0}/\mu$ and the narrow width of the junction.

S.1.4 The Collisionless “ZNG” Kinetic Model ($T>0$)

At finite temperature, the bosonic quantum gas is partially condensed and we must consider the presence of the thermal cloud. The GPE is therefore generalized to account for the thermal cloud mean field potential, $2gn_{\mathrm{th}}$, so that Eq. S.2 becomes [37]:

[TABLE]

The equilibrium $\psi_{0}$ solves the time-independent generalized GPE:

[TABLE]

where $n^{0}_{\mathrm{th}}$ is the equilibrium thermal cloud density, while $\mu(T)$ is the temperature-dependent system chemical potential counting for the thermal cloud equilibrium mean field potential. This has been used for extracting $V_{0}/\mu(T)$ ($x$-axis of phase diagrams in Fig. 1 of the main paper). To account for thermal cloud dynamics, we solve this equation self-consistently with a collisionles Boltzmann equation for the thermal molecule phase-space distribution, $f$, obeying:

[TABLE]

where $V_{eff}=V_{\mathrm{trap}}+2g[|\psi|^{2}+n_{th}]$ is the generalized mean-field potential felt by the thermal molecules, and the thermal cloud density is defined by

[TABLE]

The initial finite-temperature equilibrium distribution is obtained iteratively for a fixed total atom number, as described in Ref. [37, 42]. Our model corresponds to the collisionless limit of the “Zaremba-Nikuni-Griffin” (ZNG) kinetic theory which has been successfully used to model collective modes, vortex dynamics and evaporative cooling [65, 37, 3, 66, 67].

In choosing parameters for our finite temperature simulations, we ensure that the BEC number is equal to the corresponding $T=0$ number, fixed here to $6\times 10^{4}$ particles. In order to capture the entire thermal cloud – which resides primarily outside the BEC region – our finite temperature simulations use an extended grid length $\left[-48,48\right]l_{x},\left[-8,8\right]l_{x},\left[-8,8\right]l_{x}$ along the $x$, $y$ and $z$ directions respectively, and $2048\times 256\times 256$ grid points for the thermal cloud.

S.2 Flow Dynamics: Superfluid Current, Phase Slips, Backflow and Vortex Ring Energy

S.2.1 Extracting the superfluid current through the barrier

Figs. 1 of the main paper show numerical results for the maximum superfluid current and its temporal profile. Here we show how these results have been obtained. There are two ways to calculate the superfluid current: the first one is from the time derivative of the population imbalance $I=\dot{z}_{\rm BEC}N_{BEC}/2$ and the second one from the transverse integral of the probability current density

[TABLE]

where $R^{y(z)}_{\mathrm{TF}}$ is the Thomas-Fermi radius along the $y(z)$-direction. Here $j_{x}$ is the $x$-component of the density current of probability defined as:

[TABLE]

The latter can also be written as:

[TABLE]

where $v_{x}$ is the component of the superfluid velocity along the $x$-direction. By defining the density-weighted superfluid velocity as:

[TABLE]

we can write

[TABLE]

where $\rho_{x}=\int|\psi|^{2}(x=0,y,z)dydz$. The numerical reconstruction of Eq. (S.23) is shown in Fig. S3.

These two ways of calculating $I$ are equivalent as shown in Fig. S3 (lower plot) for the case of $V_{0}/\mu\simeq 0.8$ and $z_{0}=0.13$. We note here that the maximum of $-\left\langle v_{x}\right\rangle$ is shifted with respect to the current maximum due to the varying density (Fig. S3).

To further support the statement made in the main text that the (experimentally relevant) small thermal cloud fraction does not noticeably influence the VR nucleation process and early dynamics, we also show in Fig. S4 the superfluid current calculated from the time-derivative of the population imbalance (as in Fig. S3(c)), but for a larger population imbalance $z_{0}=0.25$, corresponding precisely to the parameter regime of Fig. 2 of the main paper. The superfluid current is seen to be unaffected by temperature, even in the case when multiple rings are generated.

S.2.2 Second order term in the Josephson current-phase relation

Here we explain the role of the second order harmonic to the current-phase relation plotted by the green shaded area in Fig. 1(b) of the main text.

The overall trend of the maximum current as a function of $V_{0}/\mu$ is quantitatively captured by extending the analytic model developed in Ref. [20] for homogeneous Bose gases and rectangular barriers, to the harmonically trapped case with a Gaussian barrier, relevant for our study. While more details will be given elsewhere, in the following we briefly summarize how such model extension has been obtained. First, we recall that based on a perturbative approach valid in the limit $V_{0}\gg\mu$, Meier and Zwerger [20] derived analytic expressions for the critical current, up to second order in the tunneling hamiltonian, for a homogeneous, $T=0$ weakly interacting BEC. In particular, they obtained explicit predictions for the first and second order contributions to the superfluid current, respectively denoted by $I_{c}$ and $J_{1}$, yielding a current-phase relation of the kind [20]:

[TABLE]

where $\bar{g}=J_{1}/I_{c}$. Notably, these terms solely depend upon the bulk condensate density and the boson tunneling amplitude, which in turn can be recast in terms of the bulk chemical potential and the single-particle transmission coefficient across the barrier [20]. As such, within this framework, predictions for the maximum current supported by a generic junction can be obtained from the knowledge of the bulk properties of the superfluid, and by evaluating the single-particle transmission coefficient associated with the specific barrier under consideration.

To proceed with obtaining the theoretical curves, we modified the analytic results given in Ref. [20] in the following ways. First, we derived an analytic expression for the tunneling amplitude $t_{cc}(V_{0},\mu)$ (using the notation of [20]) across our Gaussian barrier, by approximating the Gaussian profile with a symmetric Eckart potential of the kind $V_{Ec}=V_{0}/{\rm cosh}(x/d)^{2}$ with $d=0.6w$.

Second, we employed local density approximation to recast the first and second order currents within an integral form, in order to account for the inhomogeneous density distribution featured by our trapped samples. From the knowledge of $I_{c}$ and $J_{1}$, we then obtained the maximum current enabled by our junction, up to first and second order, respectively. These are shown in Fig. S5. In the former case, the critical current is simply given by $I_{c}$.

The yellow shaded area in Fig. S5 is delimited by the trend of $I_{c}$ predicted by our analytic model, assuming a $\pm 5\%$ relative variation of the peak chemical potential $\mu_{0}$ of the trapped Bose gas at $T=0$, around the nominal value of $\mu_{0}$ based on the measured molecule number and trapping frequencies.

While $I_{c}$ provides an excellent description of $I_{\mathrm{max}}$ derived from the experimental data and numerical simulations for $V_{0}/\mu\geq 0.9$, it systematically underestimates the maximum current found for lower barrier heights. This mismatch is expected in light of the fact that second order contributions become increasingly important for progressively lower $V_{0}/\mu$ values [20]. In order to account for second order corrections to $I_{\mathrm{max}}$, we exploited our model prediction for $J_{1}$ in connection with the analytic results obtained in Ref. [43] for a generic current-phase relation with first and second harmonics (Eq. (S.24)). For any value of $\bar{g}=J_{1}/I_{c}$, it can be shown [43] that the maximum current will read:

[TABLE]

The green shaded area in Fig. 1(b) of the main paper shows the trend of $I_{\mathrm{max}}$ based on Eq. (S.25) and the value of $\bar{g}$ derived from our analytic model, assuming a $\pm 5\%$ variation of the peak chemical potential, as for the first order case. By comparing corresponding first (yellow area) and second order (green area) contributions in Fig. S5, one can notice how inclusion of second harmonics generally increases $I_{\mathrm{max}}$, and enables to excellently reproduce both experimental and GPE results, down to barrier heights as low as $V_{0}\sim 0.6\mu$.

S.2.3 Population imbalance decay and phase slippage

Here we provide more evidence for interpreting the relation between population imbalance decay (dissipation), vortex ring (VR) nucleation [44, 45, 46], and phase slippage.

S.2.3.1 VR effects on the population imbalance profile

Fig. 1(a) of the main text plots the phase diagram in terms of a critical population imbalance $z_{\rm cr}$. This value was identified as the value of $z_{0}$ at which $z_{BEC}(t)$ firstly exhibits decay of the population imbalance to zero. Fig. S6 shows the evolution of $z_{BEC}(t)$ for $V_{0}/\mu\simeq 0.8$ and population imbalances (a) $z_{0}=0.13$ and (b) $z_{0}=0.25$.

As discussed in the main text, such dips in the otherwise rapid (and monotonic) decay of $z_{BEC}(t)$ is due to the generation of VRs, which are nucleated outside of the Thomas-Fermi surface and subsequently enter such surface – with the corresponding times indicated by vertical dashed line in Fig. S6. Specifically for $z_{0}=z_{cr}$ there is only one VR generated. For $z_{0}=0.25$, the $z_{BEC}(t)$ achieves its zero value later in time with respect to $z_{0}=0.13$. This is not just because its value is higher but also because a VR causes a backflow (see next section) everytime it is generated, i.e. larger number of VRs give a cumulative effect slowing down the population imbalance decay. Fig. S6(b) shows that a higher population imbalance leads to a larger number of vortex rings generated (an effect already shown in Fig. 3(b) of the main paper). Specifically, seven VRs are generated in the case $z_{0}=0.25$ (case discussed in Fig. 2 of main paper). The total number of VRs ($N_{\mathrm{VR}}$) propagating in the left well, which is plotted in the main paper Fig. 3(b) is found by looking at the 3D density plots. However different methods can be used as comparison to check for consistency: $N_{\mathrm{VR}}$ can be also found by counting the number of -$\left\langle v_{x}\right\rangle$ time evolution maxima. Another way to stimate $N_{\mathrm{VR}}$ is by the expression [60]:

[TABLE]

where the $\langle\cdots\rangle$ indicate time-averaged values, $(\mu_{L}-\mu_{R})$ is the chemical potential difference and $dn/dt$ is the rate of VRs crossing a certain path. The chemical potential in the left (right) well is estimated by using Eq. (S.4), upon replacing NBEC with $N_{L}=(1-z_{\mathrm{BEC}})/2$, ($N_{R}=(1+z_{\mathrm{BEC}})/2$).

Note that the first VR is generated earlier in time with higher population imbalance, as shown in Fig. S7: This is because the larger $z_{0}$, i.e. larger chemical potential difference, leads to higher initial superfluid acceleration as followed from Josephson-Anderson expression $m\dot{v}=\triangledown\mu$, i.e. the critical velocity is reached earlier in time. As mentioned in the main paper, the first generated vortex ring travels slower (with velocity extracted by a linear fit from $|x_{VR}|\geqslant 0.55\simeq 2w$) and has a lifetime which increases with increasing $z_{0}$ (Figs. 3(a) and 3(b) inset of main paper, respectively). As a result of this, the first VR for higher $z_{0}$ propagates further into the left reservoir (which has smaller condensate density) (Fig. S7(a)), maintaining a constant radius during this propagation for longer time than the smaller $z_{0}$ cases (Fig. S7(b)).

Moreover the first VR radius for $z_{0}=0.37$ during its propagation (flat area Fig. S7(b)) is larger than the one at $z_{0}=0.13$ which is consistent with its smaller propagation velocity value (shown in Fig. 3(a) main paper).

S.2.3.2 Phase slippage

In order to get more insight on the link between phase slippage and VR nucleation in our inhomogeneous BEC we show in Fig.  S8 the time evolution of the relative phase $\Delta\phi_{x}$ along the $x$-direction for $V_{0}/\mu\simeq 0.8$ and $z_{0}=0.25$ at the location of the VR when first generated. Specifically in Fig.  S8 this is calculated for $z=0$ and $y\approx 0.5l_{x}$, a value consistent with the location of the VR as seen by the density minima of the 2D and 3D density plots of Fig.  S8(a)-(b) respectively. Fig.  S8(c) shows clearly that $\Delta\phi_{x}$ jumps locally by $\sim 2\pi$ at a time $t=8.7$ms.

More details of this process are shown in Fig. S9 which plots the renormalized 2D density and phase profiles at different times during the initial VR dynamics, and specifically from the time it is nucleated, until it enters the left reservoir. The initial phase jump of around 2$\pi$ at the vortex core position is indicated by an ‘x’ at the 2D phase profile (see S9). This picture is consistent with the phase slippage concept for superfluid helium. [60].

S.2.4 Backflow

Here we provide more details about the interpretation of the drop in superfluid velocity reported in the main text in relation to Fig. 2(a).

As also noted in the main paper, we start by decomposing the superfluid velocity $\mathbf{v}$ into a superfluid potential flow $\mathbf{v}_{f}$ and in the flow generated by the superfluid singular vorticity distribution $\mathbf{v}_{\omega}$, via

[TABLE]

Given the complexity of the system studied, as a first approximation we neglect the density gradient effects on the superfluid velocity. In addition, as the role of vortex images with respect to the BEC boundaries is still partially unresolved [47, 48, 49], we do not consider the velocity field generated by the images of the VRs with respect to the boundaries of the condensate.

In order to calculate the superfluid potential flow $\mathbf{v}_{f}$ at each time $t$ in the center of the trap $O(0,0,0)$, we subtract the velocity field $\mathbf{v}_{\omega}(O,t)$ generated by the reconstructed VRs from the total superfluid velocity $\mathbf{v}(O,t)$, i.e.

[TABLE]

where $\mathbf{v}(O,t)$ is obtained via the phase gradient of the numerically computed wavefunction $\psi$ and $\mathbf{v}_{\omega}$ is calculated via the Biot-Savart integral [50]:

[TABLE]

where $\kappa=h/m$ is the quantum of circulation, $N$ is the number of VRs present in the BEC, $\mathcal{C}_{i}(t)$ is the closed curve corresponding to the $i$-th VR reconstructed via an algorithm based on the pseudo-vorticity vector [51, 52], $\mathbf{s}(\zeta,t)$ is the position of the VR line-element corresponding to arclength $\zeta$ and $\mathbf{s}^{\prime}(\zeta,t)$ its unit tangent vector. In Fig. S10 we show the temporal evolution of the axial ($x$) components of $\mathbf{v}_{f}$, $\mathbf{v}_{\omega}$ and $\mathbf{v}$ in the centre of the trap $O$ for the parameters of Fig. 2 of the main paper ($V_{0}/\mu\simeq 0.8$ and $z_{0}=0.25$).

We observe that in the event of a vortex ring generation, the superfluid flow $v_{fx}$ slows down and even reverses its sign for the first two rings generated: the nucleation of vortex rings leads to a reduction in the main flow, thus slowing down the population imbalance dynamics.

In fact, the corresponding reduction in $I$ is clearly visible in both Figs. S3 and S4, and it is particularly pronounced for the only (Fig. S3) or the first (Fig. S4) generated VR, for which it even reverses sign (i.e. $-I<0$).

S.3 Vortex Ring Energy calculation

The total energy of the system can be decomposed as following:

[TABLE]

where

[TABLE]

is the interaction energy density,

[TABLE]

is the potential energy, and

[TABLE]

is the kinetic energy density where $\displaystyle e_{q}(\vec{r},t)$ is the quantum pressure term. The quadratic term involving the velocity of the flow can be split into incompressible $\displaystyle e_{k}^{i}(\vec{r},t)$ and compressible $\displaystyle e_{k}^{c}(\vec{r},t)$ contributions as follows:

[TABLE]

where $\vec{\nabla}\cdot(\sqrt{\rho}\vec{v})^{i}=0$ and $\vec{\nabla}$x$(\sqrt{\rho}\vec{v})^{c}=0$. The fields $(\sqrt{\rho}\vec{v})^{i}$ and $(\sqrt{\rho}\vec{v})^{c}$ are calculated employing the Helmholtz decomposition [63, 53, 54, 55].

In Fig. 2(d) of the main paper we have shown the temporal evolution of the integrals of $\displaystyle e_{k}^{i}(\vec{r},t)$ and $\displaystyle e_{k}^{c}(\vec{r},t)$ (called $E_{k}^{i}$ and $E_{k}^{c}$ respectively). In Fig. S11 a) instead we show also the integral of $\displaystyle e_{q}(\vec{r},t)$ (called $E_{q}$) which is larger than $E_{k}^{i}$ and $E_{k}^{c}$ and it is slightly affected by the VR nucleation. As discussed in the main paper, when each VR enters the TF surface (vertical grey shadow areas in Fig. S11, shown also in Fig. 2 of the main paper), we see an increase of the compressible energy. This step-like increase is shown better in Fig. S11 (b) where only $E_{k}^{c}$ is shown until the second VR is nucleated. The step-like increase starts when the first VR enters the Thomas-Fermi surface and stops (reaches a ’plateau’) when the VR is at the position corresponding to maximum transverse Thomas-Fermi spatial extent, which for this case occurs around t$\sim$ 9.3 ms (shown by vertical dashed line in figure S11). Moreover, the ratio of the increase of $E_{k}^{c}$ to the compressible energy before the first VR nucleation takes the values $1.09\,,1.36\,,1.57$ for the nucleation of the $1^{\rm st}$,$2^{\rm nd}$, $3^{\rm rd}$ VRs, respectively.

In Fig. 3(a) of the main paper, we have instead used the incompressible kinetic energy of the generated VRs (shown as a function of $z_{0}$ ) to further characterize their motion. The incompressible kinetic energy of the VR $E_{k,VR}^{i}$, which depends on $v_{\omega}$ only, is obtained by the following procedure. We first integrate the incompressible kinetic energy density per unit mass $\displaystyle e_{k}^{i}$ on the volume $\mathcal{R}$ encompassing the vortex ring. The region $\mathcal{R}$ is defined as follows

[TABLE]

where $x_{V\!R}=-l_{x}$ and $\Delta x=4R_{VR}$, where $R_{VR}$ is the VR radius when $x_{V\!R}=-l_{x}$. (For improved clarity, a visualization of the region $\mathcal{R}$ is given in Fig. S12).

The values of $x_{\mathrm{VR}}$ and $\Delta x$ are chosen stemming from the following considerations. First, in a homogeneous and unbounded BEC, a cylindrical region, coaxial with respect to a VR, of height $8R_{VR}$ (the VR center being placed in the center at $4R_{VR}$) and radius $5R_{VR}$ contains more than $90\%$ of the incompressible kinetic energy of the VR. Second, the flow velocity is negligible in $\mathcal{R}$, i.e. $\mathbf{v}_{f}\approx 0$ and hence $\mathbf{v}\approx\mathbf{v}_{\omega}$. Therefore, the incompressible kinetic energy present in $\mathcal{R}$ only originates from the VR. As a consequence, $\displaystyle E_{k,VR}^{i}=\int_{\mathcal{R}}e_{k}^{i}d\mathbf{x}$ is a reasonable estimate of the VR incompressible kinetic energy. Moreover, the VR is still close to its nucleation region and therefore $E_{k,VR}^{i}$ can, at least qualitatively, be considered as proxy for the VR initial incompressible kinetic energy. Investigating the dependence of $E_{k,VR}^{i}$ on the initial population imbalance $z_{0}$, shown in Fig. 3(a) of the main paper, we find that $E_{k,VR}^{i}$ is an increasing function of $z_{0}$.

To assess the impact of confinement and inhomogeneity, we also evaluate the corresponding energies of the first nucleated VRs if they had been immersed in an unbounded and homogeneous BEC employing the analytical formula derived by Roberts and Grant [56]:

[TABLE]

where $\kappa$ is the quantum of circulation and both the density $\rho$ and the healing length $\xi$ have been evaluated on the $x$-axis for $x\sim x_{VR}$ (allowing $\rho$ to regain its bulk value not perturbed by the presence of the VR). The result of this calculation is reported in Fig. 3(a) of the main manuscript (pink triangles) showing that confinement and inhomogeneity play a fundmamental role in the energetics of the system studied [57, 58].

S.4 Dissipation of kinetic energy

In the main manuscript, to quantify the irreversible dissipation of Josephson oscillations, we introduce the acoustic dissipation $\epsilon_{c}$ and the incompressible dissipation $\epsilon_{i}$. The definition of these quantities is as follows:

[TABLE]

where $E_{k,[0,1]}^{tot}$ is the total kinetic energy flowing through the junction until the nucleation of the first VR, $E_{k,VR}^{i}$ is the incompressible kinetic energy of the first VR (as illustrated in the previous section) and $\Delta E_{k}^{c}$ is the increase of compressible energy in the system observed in correspondence of the nucleation of the first VR ([Fig. 2(d)] in the main paper).

S.5 Kelvin waves and barrier removal

Fig. 4(a)-(b) of the main paper discussed how the barrier removal process – implemented experimentally prior to time-of-flight observation – affects the VR dynamics. Specifically, we first let the system evolve for 13 ms, and subsequently we remove the barrier linearly over a period of 40 ms, and observe the motion of the fourth generated VR (i.e. VR generated around $t\sim 17.5$ ms in Fig. 2 of main paper). As commented in the paper, this process significantly extends the lifetime of this VR.

More details of this effect are given in Fig. S13, which compares the evolution of the x-position of the 4th VR in the case of barrier on (blue points in Fig. S13(a)) and barrier removal (red points), showing also a direct comparison of appropriate 3D density isosurfaces which reveal the vortex rings in both cases at times $t\approx 19.5$ ms (Fig. S13(b)-(c)) and $t\approx 23.4$ ms (Fig. S13(d)-(e)). From the slope of Fig. S13(a), we deduce that the VR velocity in the case of barrier removal is smaller than the one with barrier on. This is consistent with a larger propagating VR radius in the case of barrier removal, an effect visible by comparing the size of the VR highlighted inside the white box in the 3D density plots of Fig. S13(b)-(c). Specifically at $t=19.5$ms, the image shows two VRs (respectively the third and fourth VRs generated), while a few ms later, at 23.4ms, we see clearly that only one VR remains in the case of barrier on (Fig. S13(d)), whereas three VRs are still visible in the corresponding case of gradual barrier removal (Fig. S13(e)). Note that, in the case of barrier removal, the total number of VRs generated is less than the corresponding case with barrier on, because the decrease in the barrier height leads to an increase in the density inside the barrier, thus increasing the local speed of sound, which in turn decreases the superfluid velocity to below the speed of sound $c=\sqrt{gn/M}$ – a process which prohibits further VR generation.

Another crucial point to note in relating our findings to experimental observables is that in the experiment [25] they observe vortices propagating in the system even for initial barrier height $V_{0}>\mu$. In our simulations with fixed barrier height $V_{0}>\mu$, the VRs generated inside the barrier shrink fast without entering in the bulk, whereas the experimentally relevant barrier removal enables their detection in the bulk. This is shown in Fig. S14 for $V_{0}/\mu=1.2$ for both cases of barrier on and barrier removal.

Specifically, we observe that once the barrier height, while decreasing its value, reaches some characteristic value (here $\sim 0.9\mu$) the VR is able to escape the barrier region and propagate in the left reservoir (red points in Fig. S14(a)), unlike the corresponding case of constant barrier height $V_{0}>\mu$ (blue points). As visible in the density plots at t=20.7 ms (Fig. S14(c)-(d)), the VR in the case of barrier on (Fig. S14(c)) has a much smaller radius that the corresponding one when the barrier is removed (Fig. S14(d)), i.e. a much smaller energy. For this reason at the subsequent t=24.5 ms, the VR with barrier on (Fig. S14(e)) has already shrunk while the one generated under gradual barrier removal (Fig. S14(f)) propagates inside the superfluid. In fact the removal can in this case (t=24.5 ms) facilitate the simultaneous observation of two VRs, in stark contrast to the barrier on case which reveals none.

S.5.1 Kelvin Waves

The main paper discussed the role of Kelvin Wave (KW) excitation on the VR dynamics, already visible in Figs. 4(a)-(b) of main manuscript. Here we provide further details on this characterization. Due to the anisotropy in the transverse direction ($\omega_{y}\neq\omega_{z}$) the VR shape is elliptic when it is nucleated. During its propagation the VR shape oscillates by inverting its elliptical semiaxis, representing an $m=2$ KW excitation of the circular shape. In fact, the VR 2D profile is best fit by the function $(y/a)^{2}+(z/b)^{2}=1$. If we define $R_{VR}=(a+b)/2$, the deformation of the elliptic VR from its ideal circular shape with radius $R_{VR}$ is $(a-R_{VR})$ along the y-direction and $(b-R_{VR})$ along the z-direction. Fig. S15 shows the time evolution of these quantities in the case of the 4th VR undergoing the barrier removal process over the relevant post-generation temporal window $t\in[17.5,\,45]$ ms. After such time the VR is destroyed by interaction with the condensate boundary, an effect already studied in different contexts in [47, 48, 75, 76].

In the limit of a VR radius much larger than the core size, the period of oscillations of KW is given by the dispersion relation $\displaystyle\omega(k)\sim\kappa k^{2}/(4\pi)[\ln(2/(k\xi)-0.5772]$ [71], with $k$ the wavenumber, $\xi$ the vortex core and $\kappa=h/M$ the quantum of circulation. The KW wavenumber $k$ is found from the wavelength $\lambda=2\pi/k$, which satisfies the relation $\lambda m=2\pi R$ for $m=2$ and $R$ the VR radius. Estimating the vortex core inside the bulk as being comparable to the molecular BEC coherence length ($\xi\simeq 0.5\mu$m) and approximating $R$ as the mean value of the 4th VR radius in the time interval [24, 34] ms, the above dispersion relation predicts a KW period $\tau\simeq$ 3.3 ms,which is found to be in excellent agreement with a sinusoidal fit to our numerically extracted values of $(a-R_{VR})$ in the range [24, 34] ms (black line in inset of Fig. S15) which yields $\tau=(3.19\pm 0.03)$ ms; we note that this agreement is excellent, even though the VR radius is only five times the vortex core.

S.6 T$>$0 Vortex Ring Dynamics

The main paper has shown that, in the range $0\leq T\leq 0.4T_{\rm c}$ considered (for which $0.4V_{0}<k_{B}T<0.8V_{0}$) finite temperatures have practically no effect on the VR generation process, and the subsequent early VR dynamics, provided the condensate number is fixed to the corresponding $T=0$ value. In this case, the main effect is the addition of an extra potential to the BEC (Eq. (S.15)-(S.16)), as previously observed in the context of vortex dynamics and reconnections [65, 66, 67]. This is because such timescales are typically much shorter than corresponding ones for dissipation due to the relative motion between the BEC and the thermal cloud . At the same time, temperature has a notable cumulative effect on the overall post-generation VR dynamics: specifically, it reduces the VR lifetime (shown in Fig. 3(b) (inset) of the main paper for the first generated VR and different $z_{0}$) and breaks its motional symmetry (main paper Fig. 4(c)). This effect was already shown to become more pronounced in the cases of barrier removal, because in such cases the VR lifetimes are considerably longer. Here we provide more details on the effect of the thermal cloud on the VR dynamics.

For easier visualization, Fig. S16 shows an example of the axial condensate (black) and thermal cloud density (red) for $V_{0}/\mu\simeq 0.8$, $z_{0}=0.25$ and T$\simeq 0.4T_{c}$. In this figure, density is plotted along the x-direction (for $y=z=0$) both (a) at equilibrium, and (b) at a later time ($t=10.6$ ms) when the first VR has already entered the left well (visible for $x\sim-1.4l_{x}$). As expected, the thermal cloud has local maxima at both the barrier position, and at the edges of the condensate where the BEC density has local minima. This is because of the repulsive interaction between the thermal cloud and the condensate. Moreover when a VR is present, its core is filled by thermal molecules, as seen clearly in Fig. S16(b) (around $x\sim-1.4l_{x}$). This effect has already been reported in [65, 66, 67].

Even though the considered temperature range (with $k_{B}T<0.8V_{0}$) has a small effect on the VR generation process, nonetheless it does exert a ‘drift force’ causing the VR to go off-center while propagating along the negative $x$-direction. A clear visualization of this effect is shown in Fig. S17, showing the time evolution of the fourth VR center along the $y$ and $z$ directions for a characteristic single numerical realisation when the barrier is gradually removed (corresponding to Fig. 4(c) of the main paper) . Corresponding 3D densities are shown in Fig. S18, clearly contrasting the $T=0$ to the $T>0$ case, providing an alternative visualization to that of Fig. 4(b)-(c) of the main paper. Note that while the lifetime is well predicted within our kinetic model, the precise details of the VR trajectory – i.e. exactly how it goes off center and approaches the boundary – are sensitive to the numerical realisation.