Stabilizing persistent currents in an atomtronic Josephson junction necklace

Arrays of Josephson junctions are at the forefront of research on quantum circuitry for quantum computing, simulation, and metrology. They provide a testing bed for exploring a variety of fundamental physical effects where macroscopic phase coherence, nonlinearities, and dissipative mechanisms compete. Here we realize finite-circulation states in an atomtronic Josephson junction necklace, consisting of a tunable array of tunneling links in a ring-shaped superfluid. We study the stability diagram of the atomic flow by tuning both the circulation and the number of junctions. We predict theoretically and demonstrate experimentally that the atomic circuit withstands higher circulations (corresponding to higher critical currents) by increasing the number of Josephson links. The increased stability contrasts with the trend of the superfluid fraction – quantified by Leggett’s criterion – which instead decreases with the number of junctions and the corresponding density depletion. Our results demonstrate atomic superfluids in mesoscopic structured ring potentials as excellent candidates for atomtronics applications, with prospects towards the observation of non-trivial macroscopic superpositions of current states.

An array of junctions in a multiply-connected geometry forms a Josephson junction necklace (JJN).In this configuration, the Josephson effect is used to control the current of persistent states, implementing robust dynamical regimes characterized by the competition between tunneling and interaction energies [15].JJNs with one or two junctions realize common quantum interference devices (SQUIDs) [17,18], which find applications in rotation sensing with superfluid gyroscopes [19,21] and magnetic field sensing with superconducting rings [17,22].Furthermore, JJNs are key elements of atomtronic circuits [23][24][25][26].Ultracold atoms in toroidal traps with a single junction or a weak link have been explored for the experimental realization of ultra-stable circulation states [27][28][29][30], including the study of various superfluid decay phenomena [31][32][33], current-phase relations [34] and hysteresis [35].These experiments have stimulated several theoretical studies that have mainly focused on the analysis of different instability phenomena in ring superfluids with various types of defects and potentials [36][37][38][39][40][41][42][43].In addition, double-junction atomtronic SQUIDs have enabled the observation of different regimes of Josephson dynamics [44], resistive flow [45] and quantum interference of currents [46].Interestingly, as conjectured by Feynman [47], further intriguing * These authors contributed equally to this work quantum coherence effects can arise -due to the stiffness of the superfluid phase -in ring systems hosting arrays of multiple junctions.However, despite advancements both in manufacturing mesoscopic nanostructured multi-linkcircuits [48][49][50][51][52] and in engineering atomic trapping potentials [25,[53][54][55], the realization of tunable JJNs with arbitrary number of junctions remains technologically and experimentally challenging, and so far elusive in both superconducting and superfluid platforms.
In this work, we investigate supercurrent states in an atomtronic JJN.We analytically predict the stabilization of persistent currents against decay by increasing the number of junctions, n.We support this surprising prediction by numerical simulations and we demonstrate it experimentally in a bosonic superfluid ring with n up to 16.Such an effect is a direct consequence of the single-valuedness of the order parameter, reflecting the macroscopic phase coherence of the superfluid state.Increasing the number of Josephson links leads to a decrease of the superfluid speed across each junction and to the corresponding increase of the global maximum (critical) current in the ring.Furthermore, the density depletion associated to an increasing n determines a decrease of the superfluid fraction according to Leggett's formulation [56,57] that, however, does not result in a decrease of the critical current.The full control of our atomtronic circuit opens exciting prospects toward the realization of non-trivial quantum superpositions of persistent currents [58][59][60][61][62].

Critical current in a multi-junction Josephson necklace.
A steady superfluid state can be described by a collective wavefunction ψ(r) = ρ(r)e iϕ (r) , with ρ(r) and ϕ(r) being the density and the phase of the superfluid, respectively [63].arXiv:2311.05523v1[cond-mat.quant-gas]9 Nov 2023 The latter is related to the superfluid speed by υ(r) = ℏ m ∇ϕ(r), where m is the atomic mass and ℏ the reduced Planck constant.To ensure a single-valued wavefunction, the integral of ∇ϕ(r) calculated around any arbitrary closed path Γ must be a multiple of 2π, where the integer (winding) number w is a topological invariant.In a multiply-connected geometry (e.g. in a toroidal superfluid), Eq. ( 1) defines a series of quantized persistentcurrent states labeled by w [64,65].Although the ground state is w = 0, metastable finite-circulation states (w 0) can be generated, as first demonstrated with liquid helium [66,67] and more recently with ultracold atomic gases [27,28,30,[68][69][70].
Let us consider, for the sake of illustration, a onedimensional (1D) JJN of radius R with n equivalent junctions modelled as narrow Gaussian potential barriers, rotating with angular velocity Ω (see [71]).In the rotating frame, the current of stationary states is given by where θ is the azimuthal angle along the ring.Each junction induces a dip in the particle density ρ(θ), shown as the blue line in Fig. 1(a) and (b) as calculated from the stationary state of the one-dimensional Gross-Pitaevskii equation (GPE (c) Maximum, υ max (solid orange line), and bulk, υ bulk (dotted black line), superfluid speed as a function of the number of junctions.Results in all panels are obtained from the stationary state of the onedimensional GPE with w = 1 and Ω = 0.
[71]).We emphasize that the barrier height is larger than the chemical potential and the barriers width is of the order of the superfluid healing length [71], with the density not vanishing inside the barrier.Due to the conservation of mass-flow [see Eq. ( 2)] a density dip implies a local increase of the superfluid speed υ(θ) [orange lines in Fig. 1(a) and (b)].Comparing the panels (a) and (b) of Fig. 1, obtained for the same value of the circulation w and for different number of junctions, n = 1 and n = 6, respectively, we observe that the maximum superfluid speed, υ max , drops by increasing n.This is a consequence of the topological invariance expressed by Eq. ( 1).This is seen by writing υ(θ) = υ bulk + υ n−peaks (θ), where υ bulk is the bulk speed, given by the minimum velocity along the ring, and υ n−peaks (θ) describes the n peaks of the superfluid speed.
Replacing this expression for υ(θ) into Eq.( 1), we find The bulk contribution in Eq. ( 3) is expected to change only slightly when adding sufficiently-narrow junctions to the JJN [see the dotted black line in Fig. 1(c)].On the contrary, the second term in Eq. ( 3) is proportional to nυ max .Therefore, for a given w, υ max must decrease roughly as 1/n in order to keep the integral in Eq. ( 3) constant.The decrease of υ max is confirmed by the results of GPE simulations reported in Fig. 1(c) [solid orange line].The reduction of the superfluid velocity at each barrier implies a decrease of the phase gain δϕ across each junction, upon increasing n.For a more quantitative study, we use Eq. ( 3) and notice that υ bulk = JR/ρ bulk + ΩR, where ρ bulk is the bulk density, given by the maximum density along the ring, and nδϕ = (mR/ℏ) dθ υ n−peaks (θ).We find the relation where J = J/J R , J R = ℏ/(mR 2 ) is the current quantum in the homogeneous (no junctions) case and w = w − Ω/J R is an effective circulation in the rotating frame.Varying Ω allows to address non-integer w and thus continuous values of the current.Furthermore, by inserting Eq. (2) into Eq.( 1), we obtain where f ( w, n) = (2π) 2 dθ/ρ(θ) −1 .We note that f ( w, n) ≤ f s , where f s ∈ [0, 1] is Leggett's superfluid fraction [56,57,[72][73][74][75].The latter expresses the phase rigidity of the system, quantified by the kinetic-energy response to a phase twist of the superfluid order parameter.Since f ( w, n) = f s for w = 0 and in the limit Ω → 0 [71], Eq. ( 5) connects the superfluid fraction to the current in the ring.It is possible to see that f s decreases with n as far as the junctions do not overlap substantially [71], therefore, according to Eq. ( 5) the current decreases as well.
On the other hand, by combining Eqs. ( 4) and ( 5), the phase across each junction reads   6), where f ( w, n) and ρ bulk are obtained with GPE calculations.Symbols are obtained for w = 1.44 (downward triangles), w = 2.15 (squares) and w = 3 (upward triangles), which correspond to the maximum value of w for n = 1, 3 and 5, respectively, for which a stable solution can be found.For larger values of w, for the given n, the system is unstable due to the nucleation of solitons.Lines are guides to the eye.In particular, the solid black line connects maxima of δϕ obtained for different w, separating the stable (blue) from the unstable (orange) region.The inset shows the superfluid phase ϕ as a function of the angle θ along the ring, for n = 1 (dotted green line) and n = 6 (solid blue line).(b) Critical current as a function of the number n of junctions.The analytic formula Eq. ( 7) (large black dots) superpose to the numerical calculation of the maximum current.Small white dots show the current J calculated for Ω = 0 and different values of w, ranging from w = 1 (lower) to w = 8 (upper).Solid and dotted lines are guides to the eye.The orange region corresponds to values of the current above Jc and are thus inaccessible in the system.Inset: phase across each junction as a function of the current (symbols) for n = 1 (green squares) and n = 6 (blue circles).The solid lines are the current-phase relations δϕ = arcsin( J/ Jc ) − 2πℓ J, with Jc and ℓ extracted from fitting.
In Fig. 2(a), we plot δϕ as a function of n, Eq. ( 6), where the quantities f ( w, n) and ρ bulk are obtained from the stationary states of the 1D GPE.Symbols refer to different values of w.We clearly see that δϕ decreases with n.This implies that the condition δϕ ≈ π/2 [15,16] -that determines the maximum (or critical) current Jc in the JJN -for increasing n is met for higher values of w.We find an explicit expression for Jc , by considering the usual current-phase relation δϕ = sin −1 ( J/ Jc ) − 2πℓ J [15,31], with ℓ an adimensional kinetic inductance.The condition J = Jc provides a critical phase δϕ c = π/2 − 2πℓ Jc .Replacing this value into Eq.( 6) and using Eq. ( 5), we find where f c and ρ c are the values of f ( w, n) and ρ bulk obtained when J = Jc .Neglecting the small inductance (nℓ ≪ 2π), we find that the critical current is mainly determined by the competition between n and f c (n).In Figure 2(b) we plot the critical current obtained from the GPE solution as a function of n.Numerical values agree with Eq. ( 7) (black dots, with the solid line being a guide to the eye), where ℓ is extracted from a fit of the numerical current-phase relation, e.g.shown in the inset for n = 1 (green squares) and n = 6 (blue circles).Furthermore, small white dots in Fig. 2(b) show the current of metastable states in the case Ω = 0, where J assumes only quantized values, see Eq. ( 5) with w = w.[41,42] for a study of the case n = 1).
Although the above discussion is restricted to a 1D geometry, the predicted effects are expected to hold also in higher dimensions, due to the general validity of Eq. (1).To confirm this expectation, we have performed 3D time-dependent Gross-Pitaevskii simulations [71].We prepare the ground state in an annular trap, impose a circulation w 0 , and observe the dynamics of the system in the presence of n junctions.Consistently with the results of Fig. 2, we observe a decrease with n of both the superfluid speed and the timeaveraged phase gain across each junction [71].The results of numerical simulations are schematically summarized as in Fig. 3(a).If the number of junctions is below a critical value n c that depends on w 0 , then vortices are emitted symmetrically from each barrier, causing phase slippage and a decay of both the current and the winding number in time [71].This vortex emission is the 3D analogue of the observed simultaneous nucleation of n solitons in 1D simulations in the unstable regime.
If n is increased above n c , then the emission of vortices is suppressed and the circulation remains constant in time.A higher stable circulation corresponds to a larger critical current.
Experimental system and persistent current states.We investigate experimentally the predicted increase of current stability in JJNs by realizing a Bose-Einstein condensate (BEC) of 6 Li molecules in an annular trap equipped with a variable number (n ≤ 16) of static planar junctions.Both the ring-shaped trap and the array of junctions are produced by the same digital micromirror device (DMD) illuminated with blue-detuned light to provide a repulsive optical potential.Using the high resolution of the DMD projection setup, we create a dark ring-shaped region in the x-y plane delimited by hard walls whose height is much larger than the chemical potential of the superfluid (given by µ/h ≃ 850 Hz in the clean ring), with R in ≃ 11.7 ± 0.2 µm and R rout ≃ 20.6 ± 0.2 µm being the inner and outer radius of the annulus.The potential is completed by a tight harmonic confinement along the vertical z direction, of trapping frequency ω z = 2π × (383 ± 2) Hz.The junctions can be modelled as Gaussian peaks of initial height V 0 ≃ (1.3 ± 0.2) µ and 1/e 2 -width σ = (1.2 ± 0.2) ξ, with ξ ≈ 0.68 µm being the healing length (see Ref. [71] for details on the barrier characterization).We initially trap ≃ 6.8 × 10 3 condensed atom pairs inside the ring with a shot-to-shot stability around 5%. Due to the finite lifetime of our molecular BEC, the pair number decreases over the course of the current decay by at most 20%, causing a decrease of the chemical potential of the superfluid.Consequently the value of V 0 /µ increases by up to ∼ 15% depending on the holding time.
We initialize the superfluid ring in a quantized circulation state with winding number w 0 ∈ {1, 2, 3, 4}.Following the procedure described in Ref. [30], different values of w 0 are obtained on-demand by shining a DMD-made azimuthal light intensity gradient onto the ring over a duration t I ≪ ℏ/µ, i.e. shorter than the characteristic density response time, ℏ/µ.In this way, we imprint a phase Φ(θ) = U 0 (θ) × t I /ℏ to the condensate wavefunction without modifying the atomic density [70], where U 0 (θ) is the spin-independent potential exerted by the light field on the atomic states that varies linearly with θ [30].After the imprinting, we wait 300 ms to let the cloud reach equilibrium, allowing the possible density excitations following the imprinting procedure to damp out [43].We then progressively ramp up the n Gaussian junctions over approximately 1 ms (corresponding to ≈ 6 ℏ/µ).

Stability phase diagram.
To measure the winding w in the ring, we exploit an interferometric probe [30,34,76]: we equip the atomic superfluid with a central disk acting as a phase reference [see panels (i) and (iv) in Fig. 3(b)] and measure the relative phase between the disk and the ring from the interference pattern arising after a short time-of-flight.The number of spiral arms in the interferogram provides access to the value of the circulation (winding number) at time t, w(t).The different panels of Fig. 3(b) display typical examples of experimental images.In panels (i) and (iv) we show the in-situ atomic density profile at t = 0.The atomic density (averaged over 10 experimental images) is characterized by a homogeneous bulk both in the azimuthal and radial directions.The n = 2 (i) and n = 4 (iv) junctions are clearly visible and are associated to local dips in the density, similarly as in Fig. 1 and Fig. 3(a).In panels (ii) and (iii) we show examples of spiral interference patterns emerging for an unstable dynamics, namely w(t) decreasing in time below w 0 (here, w 0 = 2 and n = 2): in (ii) t = 1 ms and w(t) = 2, while in (iii) t = 7 ms and w(t) = 1.In particular, panel (iii) shows the presence of a vortex identified as a localized low-density defect and marked by the orange arrow.The vortex emission signals the decrease of w by one quantum.In panels (v) and (vi) we show instead the interferograms for stable dynamics, namely w(t) = w 0 (here, w 0 = 2 and n = 4).A non-circular, polygonal interference pattern is visible both at short [(v), t = 1 ms] and at long [(vi), t = 20 ms] times due to the sharp phase gain at the junctions.
By averaging the winding number over ∼ 15 experimental realizations under the same conditions, we extract the evolu- tion of the mean circulation w(t) for various n.We study the dynamics up to 250 ms, which is sufficient to observe steady current states at long-times while still limiting particle losses.The measured w(t) is shown in Fig. 4(a) for w 0 = 2.We fit each curve with an exponential decay given by w(t) = w f + ∆w exp (−Γt).The fitting parameters w f , ∆w and Γ allow us to characterize the mean supercurrent.As w(t) is obtained from statistical averaging, the figure shows that the number of realizations w(t) that remain stable in time increases with the number of junctions.In particular, the number of stable realizations increases substantially when changing the number of junctions from n = 2 (red diamonds) to n = 4 (yellow squares).For n = 10 (blue circles), all realizations are stable: this demonstrates the experimental capability to create stable finite-circulation states in a JJN.
Figure 4(b) summarizes the results obtained for different w 0 and n, in the form of a stability phase diagram.In particular, we plot the quantity Γ = ∆w Γ/ max n (∆w Γ), where each horizontal line of the phase diagram is normalized to its maximum value for fixed w 0 .This quantity combines information on the difference between the initial and the final winding numbers, ∆w, namely how much the currents decay, and on the timescale over which this decay takes place, Γ.Values of Γ ≈ 1 (red regions) are obtained when most of the realizations w(t) rapidly decay towards values of the circulation lower than the initial w 0 .On the contrary, small values of Γ ≈ 0 (blue regions) are obtained when most of the realizations are stable over time, namely w(t) = w 0 .The phase diagram clearly shows that, on average, the system supports a higher number of stable realizations when increasing the number of junctions [71].By the choice of normalization, Γ shows a sharp transition from Γ ≈ 1 to Γ ≈ 0 when increasing n.The dashed white line in Fig. 3(b) denotes the critical winding number w c (n) and the corresponding current (right axes) as a function of n, as computed numerically from 3D GPE simulations.The numerical critical curve w c (n) is obtained for V 0 /µ = 1.8 and match the experimental phase diagram well.The need for a larger V 0 /µ in numerical simulations with respect to the one estimated in the experiment, is consistent with the finite lifetime of the sample (which implies that V 0 /µ increases during the dynamics) and the finite resolution of the DMD potential, which makes the barriers not perfectly identical [71].Anyway, we note that the only effect of a change of V 0 /µ on the critical line w c (n) is to provide a linear shift, meaning that the particular choice of V 0 /µ does not affect its trend, which well reproduce the experimental findings.
Given that Jc (n) ∼ n f c (n) from Eq. ( 7), a significant decrease of the superfluid fraction f s ≥ f c would overshadow the stabilization mechanism arising from increasing n.For this reason, in Fig. 4(c), we study the dependence of f s on n and indeed find a mildly decreasing trend, which is insufficient to disrupt the enhanced stability of currents for large n.According to a variational calculation by Leggett [56,57], the superfluid fraction f s can be bounded experimentally from the in-situ density profile [73][74][75]: where the density ρ(r, θ, z) is calculated from the ground state of the 3D GPE.The bounds in Eq. ( 8) are computed by restricting the azimuthal angle θ over a unit cell of size d = 2π/n and using the normalization dz dr r cell dθ ρ(r, θ, z) = 1 [56,57,72].In Fig. 4 we plot the upper (dashed red line) and lower (dashed blue line) bounds in Eq. (8).They are very close to each other as our system is approximately separable in the transverse spatial directions [73] and coincide in 1D, where f s = lim w=0, Ω→0 f ( w, n) [71].Increasing n enhances the size of the density dip relative to the unit cell length and thus decreases both the lower and upper limits in Eq. ( 8), see Fig. 4(c).Experimentally, for each value of n, we compute Leggett's upper bound on 10 different images of the experimental density.We compute the integral on the right-hand side of Eq. ( 8) by summing over all pixels inside an annular region with inner and outer radii r cut1 > R in and r cut2 < R out respectively.We have numerically verified that the values of the bounds do not depend on the exact size of this region.The corresponding mean values and standard deviations are shown as circles in Fig. 4(a).The deviations from f s = 1 in the clean torus (n = 0) are mainly due to noise in the experimental images, as well as the finite pixel size of our imaging sensor.Experimental results are well reproduced when taking into account the finite resolution of the imaging system (solid blue and red lines) and clearly show a decrease of f s with n.

II. DISCUSSION
Our work showcases the first experimental observation of ring supercurrents in periodic arrays of Josephson junctions.Such stable currents can be experimentally observed only for a sufficiently large number of links, as predicted by our theory modeling.In particular, our work shows that the maximum current flowing across the atomtronic circuit is due to a cooperative mechanism involving all the junctions rather than only to the properties of the single Josephson link.We expect the mechanism demonstrated in this manuscript to apply to any superfluids and superconductors as it soleley depends on the single valuedness of the wavefunction in a multiply-connected topology.
Therefore, a natural extension of our work will be to investigate whether the same effect stabilizes supercurrents in other annular systems, such as atomic Fermi superfluids [29,30] and supersolids [77].In the former case, the condensate fraction differs from unity even at T = 0 [78] and additional dissipative effects, such as Cooper pair-breaking [79,80] may compete with the stabilization mechanism.In the latter, intrinsic density modulations realize an array of self-induced Josephson junctions -as recently demonstrated in Ref. [75] for an elongated atomic system -which can be controlled by tuning the confinement parameters.
Finally, the exquisite controllability offered by our platform opens the way toward realizing exotic quantum superposition of superflow states [58][59][60][61][62] with possible implications in both atomtronic and quantum technologies.Numerical methods.We discuss here numerical methods used to obtain the results discussed in the main text.
1D GPE.The 1D simulations shown in Fig. 1 and 2 refer to static solutions of the GPE equation: where f (θ) = ρ(θ), J = 2π(w − Ω)/ 2π 0 dθ/ρ(θ) and Ω = Ω/J R .Here energies are rescaled in units of ℏ 2 /(mR 2 ), μ is the rescaled chemical potential, g is the interaction strength, Ṽ is the sum of Gaussian barriers centered at θ j = 2π j/n, and θ is the azimuthal angle along the ring.The free parameters g, σ and Ṽ0 are chosen in order to match the experimental conditions: σ/ξ = 1.2, Ṽ0 / μ0 = 1.4 and ξ/R = 0.056 (with R = 12 µm being approximately the inner radius of the experimental system), where µ 0 is the chemical potential obtained in the homogeneous case (without barriers) and for w = Ω = 0.For a given number of barriers, the solution of Eq. ( 9) is obtained by imaginary time evolution.
3D GPE.In order to better capture the experimental procedure and the dynamics of the system, in 3D we solve numerically the time-dependent GPE for static barriers, with ψ(r, t) being the condensate order parameter, M the molecule mass, g = 4πℏ 2 a/M the interaction strength, a = 1010 a 0 the s-wave scattering length and a 0 the Bohr radius.The external trapping potential is V(r) = V harm (r) + V ring (r) + V barr (r).Here, V harm (r) = M(ω 2 ⊥ r 2 + ω 2 z z 2 )/2 is an harmonic confinement with {ω ⊥ , ω z } = 2π × {2.5 , 396} Hz.The hardwall potential creating the ring confinement in the x-y plane is given by (11) with R in = 10.09µm and R out = 21.82µm being the inner and outer radius, respectively.The parameter d = 1.1 µm characterizes the stiffness of the hard walls, fixed such that the numerical density profiles match the in-situ experimental ones.We take V r larger than the chemical potential µ such that the density goes to zero at the boundary.The n barriers are modelled as identical Gaussian peaks of trapping potential (12) with constant width σ = 0.8 µm.We first find the system ground state by solving the GPE by imaginary time evolution and in the presence of n barriers.We then instantaneously imprint a current of winding w 0 by multiplying the ground state wavefunction by the phase factor exp(−i2πw 0 θ), where θ is the azimuthal angle.We finally study the system dynamics by solving the time-dependent GPE.For a particle number N = 6.8 × 10 3 (corresponding to the experimental condensate number), we obtain µ = 1.09 kHz leading to a value of the healing length ξ = 0.59 µm.Equation 10 is solved numerically by the Fourier split-step method on a Cartesian grid of {N x , N y , N z } = {256, 256, 80} points dividing a grid size of length −34.846 µm ≤ r ≤ 34.846 µm and −11.0 µm ≤ z ≤ 11.0 µm in the radial plane and axial direction, respectively.The time step is set to ∆t = 1 × 10 −5 ω −1 ⊥ .3D simulations results.We characterise the condensate dynamics by studying the winding number, w(t) calculated at at z = 0 and averaged over closed circle paths ranging from the inner to the outer radius.
Stable configuration.As discussed in the main text, for a fixed initial circulation w 0 we find the transition from unstable (decaying w(t)) to stable current (time-independent evolution of the winding) when the number of barriers n exceeds a critical value n c (w 0 ), see Fig. 5(a).In the stable configuration, 3D simulations reproduce the typical finding of the 1D case, namely that both the maximum of the superfluid speed and the phase gain at each junction decrease with n as shown in Fig. 5(b)-(c).In particular, figure 5 θ. Figure 5(c) instead shows the time-averaged phase gain across each junction, δϕ, as a function of n (symbols), together with a 1/n fit (dashed lines).Finally, in Fig. 5(d), we plot the time-averaged δϕ as a function of the time-averaged current (symbols).The dotted lines are the 1D current-phase relation obtained for the same number of junctions and rescaled to be consistent with the data.We see that the 3D results are consistent with the trend of an increasing critical current with n found in 1D.
Unstable configuration.If n < n c (w 0 ), we find that both w and the current decay in time via vortex emission.In Fig. 6 we show the numerical densities illustrating the microscopic mechanism of the vortex emission process.Vortices are emitted symmetrically from each barrier: they enter the ring from the central part, propagate along the transverse direction close to the barrier position until they enter the bulk and travel at the outer edge of the ring.Each vortex entering the bulk through the barrier causes a global decrease of the winding number by one.In particular, for the considered case of n = 4, the winding at t = 5.7 ms is equal to zero.We note that the detailed vortex emission process depends on the value of the barrier height considered.
Critical circulation.The stability phase diagram can be either characterized by the critical number of barriers n c (w 0 ) for fixed w 0 , as discussed above, or by the critical circulation w c (n).For a given n, the critical circulation is the largest value of w 0 for which we find a stable dynamics.In Fig. 7 we plot w c (n) as a function of n and for different values of V 0 /µ.Interestingly, all the curves are parallel to each other and thus   identical, but their height and size is distributed around the mean values mentioned in the main.To check whether nonidentical barriers could affect the current stability in the JJN, we study this case with 3D numerical simulations.In particular, in Fig. 8 we report the time evolution of the winding for n = 4 barriers and w 0 = 3, which is stable for identical barriers, when each barrier height and width are randomly selected from Gaussian distributions of mean values V 0 /µ 0 = 1.4 and σ = 0.8 µm and standard deviation ∆V 0 /µ 0 = 0.2 and ∆σ = 0.12 µm, respectively.The mean and standard deviation values correspond to the ones measured from the experimental characterization of the barriers.In the figure, we show the winding number as a function of time for 8 different runs (symbols), corresponding to different configurations of the barriers.We also plot the statistical mean value (solid line).These simulations reproduce qualitatively the experimental findings: for some barrier configurations, the winding number remains constant in time, for some others it de-cays, eventually also reaching negative values.Correspondingly, also the average winding number (solid line) decays in time.To summarize, having non-identical barriers is observed to reduce the stability of currents in the JJN for the same w 0 , explaining the discrepancy in V 0 /µ found between the experimental and numerical phase diagram.Finally, we have also checked numerically -by solving the collisionless Zaremba-Nikuni-Griffin model [20,81] for an experimentally estimated condensed fraction of 80% -that finite-temperature dissipation does not affect the critical winding number and it affects only slightly the decay time.
Experimental methods.We discuss here additional details of the methods employed to obtain the experimental data presented in the main text.
Characterization of the tunneling barriers.Due to the finite resolution of the DMD-projecting setup, the barriers of experimental JJNs are not identical.We characterize the properties of each barrier in the different configurations at various n by acquiring an image of the DMD-created light profile by means of a secondary camera, and calibrating the optical potential via the equation of state of a BEC in a well characterized 3D harmonic trap [78].Then, we extract the height and 1/e 2 -width by fitting the radially-averaged profile of each barrier with a Gaussian.From this set of data, we extract the mean values and standard deviation of barrier height V 0 ≃ (1.3 ± 0.2) µ and width σ = (1.2 ± 0.2) ξ.Error bars denote the standard deviation of the parameters over the set of barriers.Even though the barriers are not strictly identical, the obtained results show that it is possible to create similar barriers with fluctuations on V 0 and σ that are only a fraction of the chemical potential and healing length, respectively.
Experimental phase profile in the JJN.As already commented in the main text, the interferograms associated with stable realizations, namely in which w = w 0 , show interference fringes with a clear polygonal structure (e.g.squared for n = 4), which are a manifestation of the phase jump at each Josephson junction.Thanks to the high resolution of the imaging setup, we can extract the local relative phase between the ring and the reference central disc, ϕ, as a function of the azimuthal angle θ, as reported in Fig. 11 To enhance the fringe contrast, for each value of θ we consider a radial slice of the polar interferogram averaged over ∆θ = 0.47 rad.The image is obtained by averaging 5 similar experimental spiral patterns and unwrapping the resulting image.The lower panel is the fitted azimuthal trend of ϕ, showing a phase jump in correspondence of each barrier.The unwrapped phase here is averaged over the profiles of all experimental images for (w 0 = 2, n = 4) where w(t) = w 0 at times t > 10 ms.
polar coordinates [(a)] display a characteristic step-like shape of the fringes, closely resembling the predicted behavior of the JJN phase by the analytical model and from the numerical simulations (see the inset of Fig. 2(a) for comparison).We then quantitatively extract the value of ϕ(θ) as the phase shift in the sinusoidal fit of a slice of the polar interferogram at constant θ.As shown in Fig. 11 (b), the ϕ(θ) trend clearly deviates from the linear behavior expected in a clean ring [30], but it rather exhibits a number of jumps in correspondence of the barriers in the JJN.
Experimental stability phase diagram.In Fig. 9 we provide additional experimental data regarding the statisticallyaveraged winding number as a function of time and for different n and w 0 .In the case w 0 = 1, ⟨w⟩ is found to be constant in time up to 250 ms for any n.In particular, in Fig. 9(a) we plot the case n = 16, averaged over about 20 realizations.The inset of Fig. 9(a) shows ⟨w⟩ at time t = 250 ms and for n ranging from 1 to 16.Only in the case n = 12, we found a single experimental realization (out of 18 independent runs) with w = 0.In Fig. 9(b) and (c) we plot the cases w 0 = 3 and w 0 = 4 [the case w 0 = 2 is shown in Fig. 4(a)].
In Fig. 10 we report the analogue of the stability phase diagram of 4(a) here plotting w f /w 0 , where w f is the average circulation at long time, as obtained from a fit (see main text).This shows the experimental deterministic realization of sta- ble circulation states for w = 1 and w = 2 in a toroidal trap with up to n = 16 junctions.Superfluid fraction and the f ( w, n) function.The superfluid fraction for neutral atoms in a ring trap rotating at an angular velocity Ω can be defined as [56,57] where L is the expectation value of the angular momentum and I cl is the classical moment of inertia.In 1D, we have where ρ(θ) is normalized to one over the unit cell of azimuthal size d.Equation ( 14) is derived by noticing that the two bounds in Eq. ( 8) coincide in 1D.In our case, restricting to the unit cell as in Ref. [56,57] is not necessary and Eq. ( 14) is unchanged if we write f s = 1 (2π) 2 [ dθ ρ(θ) ] −1 with ρ(θ) normalized to one over the full circle, even in the presence of n junctions.In particular, we have f s = lim w=0, Ω→0 f ( w, n), where f ( w, n) is related to the current according to Eq. ( 5).In Fig. 12(a) we plot f s (circles) and f c [corresponding to f ( w = wc , n), dots] as a function of n.Both functions decrease with n until the barriers start to overlap.In Fig. 12(b) we plot f ( w, n) as a function of w for n = 6.
To compare numerical and experimental data in Fig. 4(c), we have taken into account the finite spatial resolution of the imaging system, characterized by a Point Spread Function (PSF) of full-witdh-half-maximum FWHM = 0.83 µm [78].To estimate the theoretical curves of Fig. 4(c), we first integrate the 3D numerical densities along the z direction, Then, we account for the finite experimental resolution by convolving the integrated numerical densities with a two-dimensional Gaussian with a FWHM matching the experimental PSF.This procedure leads to a decrease in the resolution of the density modulation, which causes the estimated superfluid fraction to increase and yields results in good agreement with experimentally extracted values [see Fig. 4(c)].

FIG. 1 .
FIG.1.Superfluid speed in a JJN.Panels (a) and (b) show the particle density ρ (blue line) and the superfluid speed υ (orange line) in a 1D ring, divided by the density (ρ 0 ) and speed (υ 0 ) in the homogeneous ring, respectively.The two panels correspond to a onedimensional JJN with n = 1 (a) and n = 6 (b) junctions, respectively.(c) Maximum, υ max (solid orange line), and bulk, υ bulk (dotted black line), superfluid speed as a function of the number of junctions.Results in all panels are obtained from the stationary state of the onedimensional GPE with w = 1 and Ω = 0.

15 FIG. 2 .
FIG.2.Superfluid phase and critical current in a JJN.(a) Phase gain across each junction as a function of n, Eq. (6), where f ( w, n) and ρ bulk are obtained with GPE calculations.Symbols are obtained for w = 1.44 (downward triangles), w = 2.15 (squares) and w = 3 (upward triangles), which correspond to the maximum value of w for n = 1, 3 and 5, respectively, for which a stable solution can be found.For larger values of w, for the given n, the system is unstable due to the nucleation of solitons.Lines are guides to the eye.In particular, the solid black line connects maxima of δϕ obtained for different w, separating the stable (blue) from the unstable (orange) region.The inset shows the superfluid phase ϕ as a function of the angle θ along the ring, for n = 1 (dotted green line) and n = 6 (solid blue line).(b) Critical current as a function of the number n of junctions.The analytic formula Eq. (7) (large black dots) superpose to the numerical calculation of the maximum current.Small white dots show the current J calculated for Ω = 0 and different values of w, ranging from w = 1 (lower) to w = 8 (upper).Solid and dotted lines are guides to the eye.The orange region corresponds to values of the current above Jc and are thus inaccessible in the system.Inset: phase across each junction as a function of the current (symbols) for n = 1 (green squares) and n = 6 (blue circles).The solid lines are the current-phase relations δϕ = arcsin( J/ Jc ) − 2πℓ J, with Jc and ℓ extracted from fitting.

Figure. 2
(b) clearly shows that Jc increases with the number of junctions.When J > Jc , the current enters the unstable regime [red regions in Fig. 2(a)-(b)], characterized by the simultaneous emission of n solitons from the barriers (see Refs.

FIG. 3 .
FIG.3.Sketch of the experiment and observables.(a) After preparing an initial persistent current state with circulation w 0 , the n junctions are ramped up (see text).The 3D density plots are isosurfaces obtained from 3D GPE numerical simulations of the experimental set-up.If n is below a critical value n c depending on w 0 , the initial current is dissipated via the nucleation of vortices (here n = 2 and vortices are highlighted by orange cycling arrows in the upper right plot).Conversely, if n ≥ n c (here n = 4), the system remains stable with w = w 0 (lower right plot).(b) Examples of single-shot experimental in-situ images and interferograms obtained for w 0 = 2 and for the same number of junctions n as in (a): n = 2 (unstable configuration), at t = 0 (i), t = 1 ms (ii) and t = 7 ms (iii); and n = 4 (stable configuration) for t = 0 ms (iv), t = 1 ms (v) and t = 20 ms (vi).In the case (iii), the circulation has decayed (w(t) < w 0 ) and the vortex emission is identified by the single spiral arm and the presence of a localized region of low density, i.e. a vortex.

10 FIG. 4 .
FIG. 4. Stability phase diagram of an atomtronic JJN.(a) Mean circulation as a function of time, for w 0 = 2 and different number of barriers, n (symbols), with averages and error bars obtained from ∼ 15 repeated measurements for each point.The dashed lines are exponential fits, w(t) = w f + ∆w exp (−Γt).(b) Effective decay rate Γ ∝ ∆w Γ (colormap), extracted from the exponential fits as in panel (a) as a function of w 0 and n.Γ quantifies the stability of an initial finite-circulation state w 0 .The dashed white line is the critical circulation w c (n) and the corresponding current (right axis) as a function of n, obtained from 3D GPE simulations.(c) Upper (dashed red line) and lower (dashed blue line) bounds to the superfluid fraction f s , Eq. (8), as a function of the number of junctions.Bounds are obtained from the ground state density of the numerical GPE.The solid lines are the bounds evaluated by including the finite resolution of the experimental imaging system.Circles are the upper bound evaluated using experimental in-situ images and averaged over 10 realizations.

12 FIG. 5 .
FIG. 5. Results of the 3D GPE numerical simulations.(a) Winding number as a function of time for fixed w 0 = 4 and at different values of n (see legend).For the parameters considered in these simulations n c (w 0 ) = 8.(b) Absolute value of the superfluid velocity extracted at the mean radius R = 15 µm for z = 0, as a function of the azimuthal angle θ.Black (red) lines as obtained for n = 4 (n = 12), an initial winding number w 0 = 2 and at time t = 5 ms.(c) Time averaged phase gain across each junctions as a function of the number of barriers n and for different values of w 0 .The dashed line is the fit function that goes as 1/n.(d) Time-averaged phase-current values extracted from the 3D time-dependent GPE simulations for n = 4 (green squares) and n = 12 (blue).The dashed lines are the corresponding phase-current curves, obtained with 1D simulations [as in the inset Fig. 2(a)].

FIG. 6 .
FIG. 6. Density profiles obtained for an unstable configuration of the system (w 0 = 4, n = 4 and barrier height V 0 /µ = 1.4) from 3D GPE numerical simulations.Each panel represents the superfluid density integrated along the z-direction.The figure clearly shows the simultaneous nucleation of n vortices, each vortex being emitted from the junction edge.

FIG. 7 .FIG. 8 .
FIG. 7. Effect of the barrier heights on the critical circulation.(a) Critical circulation w c (n) as a function of n and for different values of the barrier height V 0 /µ.(b) By substracting from w c (n) the critical value for n = 2, the various curves collapse onto a single curve and we obtain a curve independent of V 0 /µ.

FIG. 9 .
FIG. 9. Statistically-averaged winding number as a function of time and for different n (symbols).The different panels refer to different values of w 0 : w 0 = 1 (a), w 0 = 3 (b) and w 0 = 4 (c).Dashed lines represent the exponential fit of each dataset, using the same fitting function as in Fig. 4(a).The inset of panel (a) reports the averaged winding number at t = 250 ms as a function of n for w 0 = 1.

FIG. 10 .
FIG. 10.Measured ratio w f /w 0 (colormap) as a function of w 0 and n.The white dashed line is the same as in Fig. 4(b).

15 FIG. 11 .
FIG.11.Extraction of the azimuthal profile of the JJN phase from the experimental interferograms.The upper panel shows the interferogram image unwrapped into polar coordinates.A sinusoidal fit of each azimuthal slice is performed by using the function A cos(αr+ϕ).To enhance the fringe contrast, for each value of θ we consider a radial slice of the polar interferogram averaged over ∆θ = 0.47 rad.The image is obtained by averaging 5 similar experimental spiral patterns and unwrapping the resulting image.The lower panel is the fitted azimuthal trend of ϕ, showing a phase jump in correspondence of each barrier.The unwrapped phase here is averaged over the profiles of all experimental images for (w 0 = 2, n = 4) where w(t) = w 0 at times t > 10 ms.

FIG. 12 .
FIG. 12. Function f ( w, n) calculated for stationary states of the 1D GPE.Panel (a) plots f ( w, n) as a function of n for two interesting cases: f c (n) = f ( w = wc , n) (dots dots) and lim w=0,Ω→0 f ( w, n) = f s (circles), corresponding to Leggett's superfluid fraction, Eq. (14).Panel (b) shows f ( w, n) as a function of w and for n = 6 (dots).In both panels, lines are guides to the eye.