Persistent current formation in double-ring geometries

Quenching an ultracold bosonic gas in a ring across the Bose–Einstein condensation phase transition is known, and has been experimentally observed, to lead to the spontaneous emergence of persistent currents. The present work examines how these phenomena generalize to a system of two experimentally accessible explicitly two-dimensional co-planar rings with a common interface, or to the related lemniscate geometry, and demonstrates an emerging independence of winding numbers across the rings, which can exhibit flow both in the same and in opposite directions. The observed persistence of such findings in the presence of dissipative coupled evolution due to the local character of the domain formation across the phase transition and topological protection of the randomly emerging winding numbers should be within current experimental reach.


Introduction
At the critical point of a second order phase transition, the symmetry of a system is spontaneously broken, with both relaxation time and correlation length diverging [1]. In the transition region, the resulting state is not perfectly ordered, rather building a mosaic pattern of frozen coherent domains, with defects spontaneously emerging at their boundaries, thus becoming embedded within the system's growing coherent global state. 5 Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
The study of emerging defect dynamics and their subsequent relaxation has been a topic of active research over the course of many years.
This effect was first discussed by Kibble in a cosmological context, setting an upper bound on the domain regions [2]. Zurek extended this by quantifying the emerging domain size on the basis of the universality of critical slowing down [1]. Consideration of the role of the quench timescale in finite-duration and linear quenches led to the universal Kibble-Zurek mechanism, which has been observed in a variety of complex systems [3], including liquid crystals [4], liquid helium [5,6], superconducting loops [7][8][9][10], ion chains [11], Bose-Einstein condensates (BECs) [12][13][14][15][16][17][18][19][20][21], and, recently, in Rydberg lattices [22]. In ring geometries, like the original configuration considered by Zurek [1], the frozen phase of the wave function may lead to emergence of a supercurrent of integer topological charge q, i.e. a 2πq phase winding along the ring. The focus of numerous previous studies has been to show how the supercurrent charge scales with the quench time through the phase transition [8-10, 17, 23, 24], thus verifying the applicability of the Kibble-Zurek scaling law. In particular, the numerical work of reference [23] provided a detailed visualization of, and highlighted, the role of local phase formation and subsequent evolution to a stable persistent current in the context of a one-dimensional model. Although such Kibble-Zurek scaling only applies to finite duration quenches, the underlying phenomenon of spontaneous symmetry breaking and generation of defects are at the heart of any crossing through a second-order phase transition, also including the numerically simpler instantaneous quenches, which provide crucial information for the formation dynamics of coherence in macroscopic systems and on the critical universal properties of the system.
Dynamical quenches (whether instantaneous or gradual) are in fact a critical ingredient of envisaged circuits in the emerging field of atomtronics, a highly-promising interdisciplinary field at the interface between matter-wave optics and the photonics/semiconductor technologies [25,26], which has also been suggested as a potential platform for quantuminformation devices, such as qubits [27]. The primary goal of atomtronics is to use the precise control, admitted by the ultracold quantum matter, to generate highly coherent circuits of neutral atoms, which are analogous to conventional solid-state systems, but with the potential to offer much improved and/or otherwise inaccessible technological applications [28,29]. Available atomtronic circuits include, in particular, variable-resistance RLC schemes [30][31][32], Josephsonjunction SQUIDs [33][34][35][36], and diodes [37].
In this work, we explore the formation and structure of spontaneously generated supercurrents induced by an instantaneous crossing of the phase transition to the BEC state in a co-planar, side-by-side, double-ring geometry. This setup has been chosen as a cold-atom analogue of a qubit made from adjacent superconducting loops, known as the Mooij-Harmans qubit [38][39][40]. Its operation relies on the use of coherent quantum-phase slips to coherently transport vortices through Josephson links. A theoretical proposal to implement this phenomenon in two-component quantum gases has been made [41], with Rabi coupling driving the phase slips with experimentally accessible spinor two-component gases [42] constituting a potential platform for such work. Gaining control over the tunnelling of persistent currents in a double-ring geometry would closer emulate the geometry of the original theoretical scheme. Recently, several proposals for two parallel/stacked rings have shown that the winding number can tunnel between rings due to the creation of fluxons at their boundary [43,44], and at the singleparticle level persistent current tunnelling has been predicted between arrays of adjacent rings and in similar configurations [45][46][47][48].
The aim of this work is to explore the full two-dimensional (2D) formation dynamics of spontaneous supercurrents in such a side-by-side geometry, and demonstrate the resulting (perhaps somewhat counter-intuitive) independence of the two connected rings, despite their density overlap and coupled dynamics. After presenting our model, in the form of the commonly used stochastic projected Gross-Pitaevskii equation (SPGPE) in section 2, we analyze the formation of persistent currents in the double-ring trap geometry (section 3), and discuss the distribution of observable persistent-current states. We then explicitly demonstrate how the double-ring setting can be reduced to that for two independent 2D annuli (section 4): doing so, we are able to draw analogies to Zurek's arguments, and also demonstrate the extension of previous 1D SPGPE single-ring simulations to an explicitly 2D setting. Our findings for the double-ring structures are shown to be robustly insensitive to details of the geometry, by explicitly verifying the findings in the lemniscate (figure-of-eight) configurations in section 5. Identifying regimes which are optimal for experimental observation paves the way for the use of such geometries for the design of potential qubit operations.

The theoretical model
Quench dynamics in quantum gases are well modelled by the above-mentioned stochastic (projected) GPE, which provides a numerically tractable effective field theory for low-lying 'coherent', or 'classical' modes of the system. First proposed in [49] to study the experimentally-observed reversible crossing of the BEC phase transition [50], it was developed independently with the addition a projection procedure in order to separate 'coherent' and 'incoherent' constituents of the dynamics [51][52][53][54]. This approach has become the workhorse for modelling the condensate formation across different platforms [13,20,23,[55][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70][71]. These include a number of simulations explicitly performed in a ring trap geometry [23,62,66], of direct relevance to the present analysis. Specifically, 1D SPGPE simulations by Das et al [23], under controlled finite-duration quenches to condensation from noisy initial conditions, revealed very clearly the local character of the formation of phase, and its long-term dynamics leading to the spontaneous formation of persistent currents (in agreement with Kibble-Zurek). The same model was used to study spontaneous Josephson vortex formation across two linearlycoupled 1D ring traps [63]. In 2D rings, vortex decay and persistent-current formation were considered in the presence of external stirring [62], closely matching experimental data [72], while spontaneous emergence of persistent currents in 2D ring traps was discussed in reference [66], as a limiting factor in the study of controlled generation and stability of counterpropagating dark soliton pairs. As the SPGPE models the dynamics of the condensate and low-lying modes of the classical field ψ (r, t), it does not include purely quantum effects like entanglement between persistent current states, which should be kept in mind when considering future atomtronic applications.
In our current implementation, building on our earlier works [20,65,70,71], individual trajectories in the coherent sector of the dynamics are governed by the stochastic equation of motion [54], describing their coupling to the incoherent sector, wherê is the Gross-Pitaevskii operator, and m is the atomic mass of 87 Rb atoms. Here, g 2D = √ 8π 2 a s /ml z is the two-body interaction strength in the two-dimensional (2D) geometry, determined by the s-wave interaction strength a s , μ is the chemical potential, and l z = /mω z is the transverse confinement scale imposed by the harmonic-oscillator trap V(z) = mω 2 z z 2 /2. In our explicitly 2D study, we adopt tight transverse confinement ω z = 2π × 1000 Hz, in order to explicitly satisfy the quasi-2D condition ω z μ. The complex Gaussian noise is characterized by correlations η(r, t)η(r , t) = 0 and η (r, t)η(r , t ) = 2γ k B Tδ(r − r )δ(t − t ), where T is the bath temperature, γ is the growth rate, and the asterisk stands for the complex conjugate. In this work, we fix γ = 0.05, and conclude via numerical testing that our steadystate results are insensitive to the exact value of γ, a statement which we have explicitly numerically confirmed here in the range 0.001 γ 0.2, which constitute values consistent with both theoretical expectations [53,54] and previous successful comparisons to experiments [20,62,70,71]. ProjectorP implements the energy cut-off, ensuring that the occupation of the largest included mode has average occupation of order unity. The energy cut-off is adopted as cut (μ, T) = k B T log(2) + μ-derived by setting the Bose-Einstein distribution f( cut ) = 1 and solving for cut [73]-setting the spatial numerical grid for simulations with spacing Δx π/ √ 8 m cut [53]. The kinetic energy for a winding number n w around a ring of radius R is given by 2 n 2 w /2mR 2 , thus for the parameters we consider here the adopted energy cut-off limits the winding numbers to |n w | 80, far beyond the range we address in this work.
To simulate the quench-induced dynamics, we simulated equation (1) using the software package XMDS [76], starting from initial condition ψ(x, y) = 0, with a random realization of the initial noisy field η(x, y), which is devoid of any phase coherence. This is akin to an input condition with N atoms, all with the initial energy larger than the cut-off energy cut , which enter the cut-off region at a rate governed by γ until the gas reaches thermal equilibrium (actual values of N are given below). This can be seen as an instantaneous thermal quench (cooling) from T T ∞ BKT to T = 10 nK T ∞ BKT . To confirm the broad validity of our findings, we have also explicitly verified that temperature quenches from pre-formed equilibrated 2D thermal clouds at temperatures T > T ∞ BKT (rather than from an initial noisy condition) yield the same qualitative findings and practically identical long-term winding number combinations (histograms; see subsequent figure 3(b)).

Instantaneous quench in the double-ring geometry
We consider the in-plane side-by-side double-ring (dr) configuration shown in figure 1(a). This structure is also known as the two-torus. It is defined by the 2D potential where ρ(x, y) = x 2 + y 2 , and the centres of the two rings with radius R and width w are set at δ 0 being a shift that separates the two side-by-side rings. The choice R > w ξ, for healing length ξ = / √ mμ, gives the rings 2D character and thus justifies the use of T ∞ BKT as a measure of critical temperature. Further, the height of the potential, V 0 = 27.5 k B nK > μ, is fixed throughout the work. This potential is the double-ring extension of the single-ring trap used in reference [77]. However, as the two-torus is not homeomorphic to the torus, due to the differing topology, one might expect the system dynamics within each section of the double-ring to be different from the single-ring case. We have verified that approximating the single-ring potential as V ext ∝ (ρ − R) 2 , instead of the expression adopted in equation (4), does not conspicuously affect results presented below. Panel (i) in figure 1(a) shows potential (4) for R = 25 μm, w = 6 μm and δ = 0, and the blue curve in (ii) shows a cross section along y = 0. These trap parameters are chosen to fit known experimental ranges, typically R = (12-70) μm and w = (3-12) μm [17,78]. Experiments with even larger radii, reaching R = 262 μm, have been performed with time-averaged potentials; however the condensate was not coherent across the whole ring in this setup [79].
Following the instantaneous quench, the atomic density gradually grows to the equilibrium value, while the phase is relaxing to a steady-state configuration. The process is random, with each numerical run proceeding differently. Vortices, spontaneously created in the course of the growth of the condensate, eventually decay, potentially leaving a persistent current with winding numbers, n L and n R , in the left and right rings, respectively. Figure 1(b) displays these dynamics, as the density equilibrates in panels (i)-(v), and the coalescence of the phase patterns is observed in (vi)-(x). For clarity's sake, the phase plots are masked by the Heaviside function, Θ(0.9V 0 − V dr (x, y)), set at 90% of the potential's height, as shown by the grey shaded area in panel (ii) of figure 1(a). In this example, the final state evolves towards an equilibrium density with about N ∼ 2 × 10 5 atoms and with phase winding numbers (n L , n R ) = (1, 1), where we use the convention of positive winding numbers for clockwise circulation [62]. The temporal evolution of n L (n R ) is shown in blue (red) in figure 1(a)-(iii), for a single numerical run. Each simulation leads to a random observation of winding numbers in the left and right rings, and we seek to quantify the effect that the presence of one ring has on the other through the stochastic distribution of these winding numbers. A movie of the time evolution for this example is provided as supplementary material to this work.
Extracting the winding number from the phase pattern is a two step-process. First, we apply high-pass filtering to the wave function in the momentum space, settingψ(k) = 0 for k > k cutoff = π/ξ, thus removing excitations with the spatial scale smaller than a healing length, ξ, including sound waves, vortices, and thermal noise. Then, we extract the phase, φ(ρ, θ), from the Madelung representation of the wave , at a radial distance ρ = R, and count the number of jumps Δφ = 2π around the ring. An example is shown in figures 2(i) and (ii), where a clear 2π jump is visible in the azimuthal phase of both rings. At early times the phase is random, as shown by the dotted and dashed lines, for both the azimuthal and radial phase profiles. At later times (after the thermalisation process) a smooth radial profile appears, whilst azimuthally the phase shows evidence of a phase winding of 2π with a spatially varying gradient. The nonlinear gradient is evident in figure 2(i) around θ = 0 and (ii) around θ = ±π. As we will discuss, this effect is pronounced for larger winding numbers, and a more extreme example can be found in appendix A and in the supplemental videos provided.
The code has been rigorously tested by manually imprinting persistent current states up to n w = 30 and comparing to the numerically obtained count. Of 10 000 tests with varying degrees of numerical noise, none was counted incorrectly. Therefore we do not include any estimate of error in results that follow.
Experimentally, the winding number is usually measured through a variety of destructive techniques. In reference [80], a ring was populated by two hyperfine states of 87 Rb, and the interference pattern between the rotating and non-rotating states was measured. A commonly employed method is to measure the size of the central hole after a time-of-flight expansion, either directly after removal of the trap [33,77], or after transforming the ring trap into a simply connected sheet first, before turning off the trap [77,80]. Recent advances in the application of this technique have been achieved by the inclusion of a small stationary BEC disk inside the annulus and measuring the interference pattern between the disk and annulus after the expansion [17,24]. A minimally destructive technique has been employed to find the winding number through measuring the Doppler shift of standing phonon modes [81], allowing for repeated measurements in the course of one experiment. There is also a recent proposal to allow a small number of atoms tunnel into a linear waveguide adjacent to the ring to monitor the persistent current in time [82] (similar to the setup well known in optics, with a microring resonator coupled to a straight waveguide [83]).
We summarize results of the simulations by means of a 2D histogram of distributions of winding numbers for the double-ring geometry. To approximate the true winding number distribution, we have performed 5000 numerical runs, each with random initial noise. Figure 3 show that (n L , n R ) stay constant at γt 10 ms. The distribution of these winding numbers fits a bivariate normal form with no correlation between the left and right rings: this is seen through calculation of the sample Pearson correlation coefficient 6 −1 < r < 1 with r = 0 corresponding to an uncorrelated state, with our numerical data we find r ∼ 10 −16 .
Despite the lack of correlation in the formation of persistent currents, an effect of the neighbouring ring on the other one is still visible in velocity fields v(x, y) = ( /m)∇φ(x, y), seen in figure 3(ii) subplots. When the inter-ring shift in equation (5) is δ = 0, the speed of atoms across the overlap region, where the two ring traps are abutting on each other, may be approximated by the relation The validity of this relation is most easily observed in panel (ii) figure 3(d), with (n L , n R ) = (2, 2). The speed is clearly zero 6 The sample Pearson correlation coefficient is defined as where P(n L = i) is the probability of observing state n L = i,n denotes the mean probability, and the summation, with respect to index i, is performed from j = min(n L , n R ) to k = max(n L , n R ).
at the centre due to the counter-flow between the two rings; nevertheless, the total phase winding around each ring is still 4π. We conclude that when, n L = n R , atoms are quiescent at the centre due to the counter-flow. The atom fluxes flow along outside path, still maintaining the total phase winding around each ring. However, when n L = −n R , there is a co-flow across the centre, with the speed doubled in the central reservoir. This scenario is shown in panel (ii) of figure 3(c), where the velocity for the case of (n L , n R ) = (2, −2) is shown as a function of the angle around the ring. The velocity is indeed doubled across the overlap of the two rings, and correspondingly suppressed around the rest of the rings, such that values 2πn w of the total phase windings are maintained in the system.
The atom flow for all other cases, with |n L | = |n R |, leads to an exchange of atoms between the rings, determined by the magnitude of |n L − n R |. Perhaps the most interesting case is displayed in panel (ii) of figure 3(f), where the supercurrent in the left ring remains zero, despite constantly exchanging atoms with the right counterpart. This is evidenced by the phase gradients present around the left ring, while the maintained total phase accumulation is zero.
The impact of these results is apparent when considering potential experimental measurements of the current. Destructive measurements are unlikely to be possible, given the close proximity of the two rings. However, the minimally destructive measurement techniques, currently applied to single-ring geometries, will be affected by the angular dependence of the velocity around each ring, shown in figure 3(ii). This will induce a spatial dependence on the Doppler shift of phonon modes [81] and the atom flux entering an adjacent linear waveguide [82].

Comparison to instantaneous thermal quenches
The use of a random noise initial condition is a numerically less demanding approximation to computing an instantaneous thermal quench. To justify our choice for the former, we compare the distribution obtained through a dynamical quench to that of one obtained from a random initial state. Performing a thermal quench simulation requires equilibrating to a thermal cloud with temperature T = 300 nK, then instantaneously quenching to the target temperature T = 10 nK at t = 0. As before, the winding numbers are measured at γt = 30 ms. Summing over the observed probabilities with fixed n L (or n R ) from the histogram yields a marginal distribution for n R (or n L ). Figure 3(b) shows that both methods produce practically identical results, and thus we choose to equilibrate from noise for the rest of this work. Error bars represent a 95% confidence interval containing the true probability, found by fitting with a normal distribution.

Benchmarking against instantaneous quenches in a single ring
The observation of the formation of uncorrelated persistent currents in the two rings in both instantaneous quenches from noise or from a thermal initial state, despite the evident density overlap, suggests that such a distribution may in fact be explained in terms of well-known results for independent single-ring settings. In quenches from a noisy initial condition, this result might have been anticipated by extrapolation of the 1D findings of Das et al [23], who highlighted the critical importance of local phase formation and evolution, over the global evolution around the ring. However in the present case, this was by no means a priori guaranteed for two reasons: firstly, the simulations in reference [23] were limited to an explicitly 1D ring, as opposed to the 2D rings numerically simulated here, which also account for the role of radial phase fluctuations (see figure 2 and movies); secondly, and most importantly, given that the two-torus is not homeomorphic to the torus, the final results in the two cases need not be mappable onto each other. Indeed, while the overall distribution of winding numbers remains the same across the two non-homeomorphic cases, the actual azimuthal phase gradient is non-uniform in the case of the double ring, in stark contrast to the uniform phase gradient of a single ring. In the end, our explicit 2D numerical simulations confirm that the local character of phase evolution also dominates in the 2D case in the determination of the long-term winding number combinations across the two rings, while the short-term evolution is random.
To make a connection to the single-ring case, we consider a single ring, defined by potential For ease of comparison, we keep here all parameters the same as in the previous section, and observe the distribution of winding numbers after an instantaneous quench in such a ring. In particular, for μ kept fixed, the arising atom number in the single trap is found to be N ∼ 1.1 × 10 5 . Our findings are presented in figure 4, showing explicitly the potential (panel figure 4(a)), and the equilibrium density at γt = 30 ms in panel figure 4(b).
The obtained distribution of winding numbers, n w , after 5000 numerical runs (figure 4(c)) reveals a (univariate) normal distribution N w ∼ N(0, σ 2 ), with σ = 1.5338 ± 0.037 within a 95% confidence interval; this is plotted as a red curve in the histogram of figure 4(c). There are several factors that can reduce the observed value of σ, including a trap's tilt angle away from the azimuthal plane [23], a longer thermal quench time [17], and smaller chemical potentials [24].
Zurek's original work considered how the thermal quench through a phase transition would leave behind a superfluid circulation in an annulus [1]. In the course of the phase transition, the condensate forms N ≈ C/d uniformly spaced, independent regions of coherent phase around the ring, with circumference C and defect size d. Taking d ∼ w [1] gives N = 2πR/w regions. Paraoanu derived [84] that for N independent condensates with uniform phases, placed alongside one another, the maximum stable winding number is N/4. Thus, for these parameters we would expect winding numbers n w < 7, which is consistent with our finding of n w = 6 being the largest recorded winding number.
Further confidence in our methodology and presented numerical predictions is provided by explicitly comparing our numerical results to findings of a recent experiment with ultracold atoms on a single ring. Specifically, Corman et al [17] considered finite-duration quenches on a quasi-2D ring, in a direct test of the Kibble-Zurek scaling law. Simulating such a finite-duration cooling quench protocol by means of the present scheme and starting from an initial thermal state, we obtain excellent agreement (see appendix B).
To compare the results obtained in the framework of the double-ring geometry, to those considered in the single-ring geometry, we compare to the marginal distributions from figure 3. These results are practically identical to the distribution of the winding numbers for the single ring, as clearly seen in figure 4(d). We thus conclude that a robust property, inherent to each ring, is that, under quenched spontaneous condensate growth dynamics, the presence of the second ring does not alter the steady-state distributions of persistent currents in a given one, due to the prevailing importance of local phase formation and evolution even in a purely 2D setting. Despite the non-homeomorphic nature of the potentials, the close agreement between the final 'steady-state' results for the two rings and the single one is understood by the fact that the current is a topologically protected number associated with each ring (it may be thought of as having a 'ghost' vortex in the centre of the ring).

Dependence on the trap geometry
The robustness of the winding number can also be tested by replacing the two identical rings by the Bernoulli lemniscate, or figure-of-eight configuration 7 . It is defined by the in-plane potential where x lem and y lem are described parametrically for t ∈ (0, 2π) as x lem (t) = 2R cos t sin 2 t + 1 and y lem (t) = 2R cos t sin t sin 2 t + 1 . 7 This was suggested to us by J Dalibard.
The choice of an effective radius R = a/ √ 2, where a is the length from the origin to the foci, gives two rings of comparable radius to the work above.
The results corresponding to this geometry are shown in figure 5. Although the overlap area between the two ringlike structures is significantly different in this case, with flow from each ring structure in the central region directly passing through each other, we still observe that the winding numbers in each ring are completely uncorrelated. Although, in this geometry, the velocity fields show exchange of atoms between the two circulation loops, there is no transfer of persistent currents between them.
We also considered the case of a finite shift between the two coupled ring traps, i.e., δ = 0 in equation (5). In the case of co-flows (n L = −n R ), patterns of the fluid flow are similar to those in the case of δ = 0. However, for realizations exhibiting counter-flows (n L = n R ) varying δ changes the interpretation of the flow circuit. At 0 < δ < w, the shear flow between the rings can create a vortex in the low-density overlap region, or several vortices, if the winding number is large enough. In this scenario, the Kelvin-Helmholtz instability may develop at the interface [85], as seen in simulations of the merger of two stacked rings with circulation [86,87], see also [88]. An example of this for δ = w/2 is shown in appendix C.
All cases considered above uphold our conclusion that, in our geometry, the resulting winding numbers in both rings are produced by the random phase profiles spontaneously emerging in the course of the condensate growth, and the topological stability of the established states is maintained by the presence of persistent currents in each closed ring-shaped geometry.

Discussion
We have explored the spontaneous growth of persistent currents in a quenched 2D Bose gas in the co-planar side-byside double-ring geometry. The emerging persistent currents are stable and long-lived, and do not transfer between the rings. The overall value of the winding number in each ring behaves independently, which can be directly attributed to the importance of the local nature of phase establishment in the azimuthal direction, an observation that remains valid even when allowing for radial random phase variations. While the distribution of winding numbers across the two rings can therefore be directly mapped onto the well-known results for the single-ring setting, it is important to note that the azimuthal phase gradient in each ring (for a given non-zero winding number) is not constant, in direct contrast to the constant phase gradient of the single ring case. This is more clearly visible in cases of large |n L | + |n R |, examples of which are shown in appendix A and movies of the real-time evolution of figures 1 and A1. This is one manifestation of the non-homeomorphism between a single torus and a two-torus. Remarkably, however, this is not necessarily an impediment to potential future atomtronic devices, which could utilize the observed robust nature of the superfluid current configuration. The independence of the ring winding numbers and stability of the formed supercurrents may facilitate a well-controlled current transfer protocol between rings. Such deterministic transfer of winding numbers across multiple-loop atomtronic architectures is a promising direction for future research.
Varying the separation of the rings demonstrates the ability to modify the atom transfer between the rings and generate vortices between them. However, there is no sign of angular momentum transfer between the rings. The robustness of the persistent currents after the quench holds steadfast even against the change of trap geometry, such as replacement of the double ring by the lemniscate potential, where one might naïvely expect uncorrelated persistent currents in individual rings to be suppressed due to atomic motion along a 'figure-of-eight' path. We also varied the phenomenological damping parameter γ in a broad range of values, and extended the simulation time, which produced no evidence of decay, confirming the efficiency of the topological protection of the spontaneously generated winding numbers.
A natural question suggested by the present results is the transfer of the winding numbers between the rings, which will be the subject of further work. The ability to control the transfer of winding numbers could be envisaged as a prototypical atomtronic switch, providing an avenue for controllable realization of coherent quantum phase slips required for the Mooij-Harmans qubit [38,39]. A hot topic of current studies of toroidal BECs is the supercurrent decay mechanism, an understanding of which could help one to control the transfer of current states between the rings. The Gross-Pitaevskii equation and its many extensions [53][54][55] do not capture supercurrent decay, even at finite temperature, without an external barrier. However, due to roughness of the potential, the decay time for large winding numbers decay in experiments is on the order of seconds, whereas n w ∼ 1 may be stable on times on the order of minutes [80]. Theoretical studies of the decay mechanism so far have all relied on effects along the annulus, caused by a repulsive barrier, the decay being visualized as vortices crossing the barrier region radially. In reference [88] it was shown that the temperatureinduced decay does not fit the Caldeira-Leggett superconductivity model, i.e., the observed decay rate does not match simple models of quantum tunnelling and thermal activation of phase slips. Using a truncated Wigner approximation, reference [90] attributed some of the disagreements between the theory and experiment to thermal fluctuations, however exact identification of the decay mechanism remains an open question.
Data supporting this publication is openly available under an Open Data Commons Open Database License [91].  can be seen in the supplemental material provided with this work.

Appendix B. Comparison of experimental results and theoretical predictions
In 1985 Zurek considered how a thermal quench through a phase transition would leave behind a superfluid circulation in an annulus [1]. He postulated that during the phase transition the condensate forms N ≈ C/d independent regions of coherent phase around the ring, with circumference C and defect size d. In previous works, the resulting average winding number, calculated at the end of the thermal quench, has been shown to scale as where we assume that the width of the ring is close to the defect's size, d ∼ w [1].

B.1. Direct confirmation of Kibble-Zurek scaling for finite-duration quenches
In the related work of Das et al [23] they tested the dependence of relation (B.1) on a linear temperature quench in time in an explicit 1D setting, and found excellent agreement. In this appendix we carry out finite duration quenches from a high to a low temperature, in an explicitly 2D geometry, thus assessing the 2D nature through the radial width w. In figure B1(a) it is seen that our analysis agrees with this prediction, indicating optimal experimental geometries for observing larger winding numbers. Specifically here we have chosen w = (3, 6) μm ξ and varied R = (12.5, 18.75, 25, 37.5, 50)μm using a finite duration quench with the ramp time τ r = 0.001 s from an initial thermal state at T = 250 nK (T ∼ 1.25T ∞ BKT ) to T = 10 nK (T ∼ 0.05T ∞ BKT ).

B.2. Comparison to experiment of reference [17]
Next we compare our results to the recent work by Corman et al [17], which addressed temperature quenches with different ramp rates in the single-ring geometry with parameters R = 12 μm and w = 3 μm. We take μ = 12.5 k B nK to match the atom number N = 36 000. At the end of each temperature quench, the winding number was measured by turning off the trap potential and looking at the ensuing interference pattern, produced by the interplay of the ring with internal stationary disk. We have carried out simulations of one of those experiments, for the finite duration thermal quench from T = 300 nK (T ∼ 2T ∞ BKT ) down to T = 10 nK (T ∼ 0.07T ∞ BKT ), with the ramp time τ r = 0.025 s. The experiment was repeated 36 times, with the aim to create a histogram of the resulting winding numbers. In figure B1 we compare the histograms representing the experimental findings and our numerical results (red and blue columns, respectively), the latter ones produced by 5000 simulations of the SPGPE. The figure demonstrates excellent agreement. Due to a relatively low number of experimental realizations, the respective histogram is not exactly symmetric about n w = 0, although it features n w ≈ 0, confirming the stochasticity of the distribution. The experimentally produced average absolute winding number corresponding to this dataset is |n w | = 0.6, while our simulations yield |n w | = 0.5926. The experiments did not feature winding numbers |n w | > 2, while 5 of our 5000 simulations yielded |n w | = 3. It may be that still larger values are possible in this geometry, but with a probability < 0.001. Comparing to the relation from reference [84], we expect that the maximum permitted value is |n w | = 6, although this prediction does not account for the temperature ramp rate.

Appendix C. Effect of finite ring separation distance δ on winding number histogram
Addressing effects of the finite separation between the rings (δ > 0 in equation (5)), in figure C1 we present results of 5000 simulations of the SPGPE with δ = w/2. The histogram of steady-state winding numbers does not display any significant difference from the δ = 0 case. This conclusion was verified by calculating the relative error between the two observed probability distributions, Δ(probability) = |P(n w (δ = 0)) − P(n w (δ = w/2))|/P(n w (δ = 0)), the result being that the variation is < 7%. The correlation of the winding numbers between the two separated rings is still effectively zero (r ∼ 10 −15 ). Note that, in the limit of large δ, the rings become completely independent systems, for which the above results for the single ring are directly relevant.
For intermediate values of δ (0 < δ < w) the velocity fields are remarkably similar to those presented for δ = 0. However, as stated in the main text, states with n L = n R exhibit shear flow between the rings which may create vortices in the low density overlap region [85][86][87].
At δ > w the rings are spatially separated, with little density overlap. Recent works have found that, even in this case, the angular-momentum states (truncated to |n w | 1) can couple to one another and tunnel, at the single-particle level [45][46][47][48]. However, the nonlinearity in the Gross-Pitaevskii equation couples the setting to higher-order angular-momentum states, and destroys the simple picture. By tuning the nonlinearity to be negligible through the Feshbach resonance [92], it may be possible to create a superfluid state admitting tunnelling of angular-momentum states.