Superflow decay in a toroidal Bose gas: The effect of quantum and thermal fluctuations

We theoretically investigate the stochastic decay of persistent currents in a toroidal ultracold atomic superfluid caused by a perturbing barrier. Specifically, we perform detailed three-dimensional simulations to model the experiment of Kumar et al. in [Phys. Rev. A 95 021602 (2017)], which observed a strong temperature dependence in the timescale of superflow decay in an ultracold Bose gas. Our ab initio numerical approach exploits a classical-field framework that includes thermal fluctuations due to interactions between the superfluid and a thermal cloud, as well as the intrinsic quantum fluctuations of the Bose gas. In the low-temperature regime our simulations provide a quantitative description of the experimental decay timescales, improving on previous numerical and analytical approaches. At higher temperatures, our simulations give decay timescales that range over the same orders of magnitude observed in the experiment, however, there are some quantitative discrepancies that are not captured by any of the mechanisms we explore. Our results suggest a need for further experimental and theoretical studies into superflow stability.


Introduction
Ultracold atomic Bose gases are versatile, highly configurable systems, in part due to their isolation from environmental effects, precise controllability with magnetic, optical, and rf fields, and the accessible imaging of many atomic observables [1,2]. This makes these systems ideal platforms for experimentally investigating superfluidity and many-body quantum phenomena [3][4][5][6]. In particular, atomic Bose-Einstein condensates (BECs) confined to multiplyconnected geometries such as a toroid are exceptionally well-suited for investigating persistent currents of superfluid flow [2,7,8], which may provide insights into the nature of supercurrents in superconducting materials. The first experimental demonstrations of persistent flow in a toroidal BEC were performed more than a decade ago [9,10]. Since then, experiments with superfluid toroidal BECs have investigated the creation and stability of persistent currents [11][12][13][14][15], atomic-gas analogs of quantum phenomena in electronic devices [16,17], quantum field dynamics in cosmic inflation [18], and compact atom interferometry in optical waveguides [19]. There have also been many recent theoretical works that have investigated superfluidity and persistent currents in one-dimensional systems [20][21][22][23][24][25], protocols for atomic-gas superfluid circuits [26,27], superflow in dipolar supersolids [28], and mechanisms for superflow decay, both within mean-field theory [29][30][31] and beyond [32,33]. Superfluid studies in toroidal atomic gases are particularly relevant to the emerging field of atomtronics, which broadly aims to develop circuit-based atomic gas devices [34,35]. Toroidal superfluids could be used to realize the matter-wave equivalent of a superconducting quantuminterference device (SQUID) [36]. However, a robust, well-functioning atomtronic SQUID requires precision control over the superflow current at the single-quantum level -as indeed does almost any atomtronic device based on superfluidity. Constructing theoretical models capable of quantitatively describing superfluid experiments is therefore essential for the future development of increasingly sophisticated atomtronic devices [37].
Here, we focus on one critical aspect of superflow: the lifetime of persistent current states in the presence of a repulsive barrier. There have been several experiments where a weak-link perturbing barrier was used in toroidal BECs to create persistent current states and alter their stability [10-13, 38, 39]. However, the experimentally-measured critical velocity of superflow significantly differs to the predictions of mean-field theory [10,39]. Indeed, despite numerous theoretical investigations into the underlying mechanisms of superflow decay in the presence of a perturbing barrier [26,[29][30][31][32][33]40], a detailed understanding of the nature and origin of superflow instability remains lacking.
Recently, an experiment by Kumar et al. studied the temperature dependence of superflow decay in a toroidal atomic superfluid [15]. Their experiments found that the rate of superflow decay was strongly dependent on temperature, in qualitative agreement with theoretical predictions [32]. To date, there have been two attempts to theoretically account for the observed decay rates in this experiment. The first was contained within the experimental report itself, where the authors used an analytic model of a solitonic vortex to estimate the energy barrier that couples the circulating and non-circulating states. From this, they found that an Arrhenius-type equation that describes thermally activated phase slips (TAPS) is inconsistent with the experimental data. Furthermore, the rate of quantum tunnelling through the energy barrier was found to be negligible for the experimental parameters studied [15]. The second theoretical analysis was performed by Kunimi et al. [31], who performed a more general and detailed TAPS calculation by estimating the energy barrier within mean-field theory. The computed decay rates were negligible compared to the experiment, thus ruling out TAPS as a quantitative description of superflow decay in the experiment of Kumar et al.
These prior approaches rely on idealised models of superflow decay based on specific decay mechanisms, and are thus unable to capture the highly non-equilibrium, three-dimensional, finite-temperature dynamics present in the experiment of Kumar et al.. To achieve a quantitative description of such a system, first-principles numerical modelling of the atomic gas is required. This is the approach we take in this work, where we perform detailed ab initio simulations of the experiment of Kumar et al. [15], within the theoretical framework of classical field (c-field) methodology [41]. In particular, our model goes beyond mean-field theory and includes both the inherent fluctuations of the multimode three-dimensional quantum state, and finite-temperature interactions with an incoherent thermal reservoir. These features are essential in order to describe the strong temperature-dependence of superflow decay observed in the experiment. Notably our model does not contain any fitted parameters, with all simulation parameters determined ab initio from the experimental atom numbers and temperatures.
Our simulations are able to capture both qualitative and quantitative features of the  Figure 1: Schematic of the key steps in the experimental procedure of Ref. [15]. (i) The atomic cloud is first prepared in a toroidal trap at temperature T . (ii) The condensate is then prepared in the l = 1 circulation state by stirring a barrier around the condensate. (iii) To induce decay of the superflow, a barrier with strength weaker than the chemical potential is raised over a period of 70ms, held constant for time t hold , and then lowered over 70ms. (iv) A measurement of the circulation is then made by releasing the atoms from the toroidal trap and subsequently observing their interference with an auxiliary disk of atoms in the center (pictured in red). experiment, with quantum and thermal fluctuations leading to a stochastic decay of the superflow. We calculate the timescale of the decay and compare this to experimental values.
The computed decay timescales range over the same orders of magnitude as the experiment and we see quantitative agreement for the lowest temperature studied. However, at higher temperatures the simulations require a larger perturbing barrier height than the experiment to achieve the experimentally-observed decay timescales. This discrepancy is largest for the highest temperature studied in the experiment, where the validity of our model is well established. This suggests that there is some aspect of the experiment that is not captured in our model. Although we have explored some possibilities for the noted discrepancies in this work, the precise origin of this effect remains unclear.

Details of the experiment
Here we briefly describe the experimental procedure of Ref. [15]. The key steps are summarised schematically in Fig. 1. In the experiment, ultracold 23 Na atoms were confined in a toroidal optical trapping potential with a mean radius of r 0 = 22.4µm. Their experiment considered Bose gases prepared at four different temperatures: T = 30, 40, 85 and 195nK. For the lowest two temperatures, vertical trapping was provided by a blue-detuned beam. For the higher temperatures, the vertical confinement was provided by a red-detuned beam, with atoms residing in the region of greatest light intensity.
To prepare the atomic superfluid in the first quantized circulation state (l = 1), a weak-link barrier slightly stronger than the chemical potential was raised adiabatically, stirred around the condensate as depicted in Fig. 1(ii), and then adiabatically lowered. In total the stirring procedure took ∼1s and prepared the desired circulation state with a fidelity of roughly 96%. The barrier itself was generated by rapidly scanning a Gaussian beam across the radial extent of the condensate, the time-average of which is approximately constant over the condensate density.
To induce decay of the superflow from the l = 1 circulation state to a non-circulating state (l = 0), the experiment introduced a stationary perturbing barrier with peak height weaker than the BEC's chemical potential. This barrier height was raised linearly over a period of 70ms, held at its maximum height V b for some variable time t hold between 0.2s and 4.6s, and then lowered linearly over another 70ms period. The raising of the barrier is not adiabatic and does lead to some small oscillations in the angular momentum that are eventually damped by the thermal cloud [32]. In order to keep the total time of the experiment constant at ∼7s, there was a variable time delay between the stirring stage and raising the barrier.
In the experiment, the final circulation state was measured by interfering the atomic cloud with a reference disk of atoms held at the central 'hole' of the toroidal trap. The measurement works as follows: the atomic cloud is released from the toroidal trap, and interferes with the reference disk as it expands during its time of flight. The resulting interference pattern is then measured, allowing the phase profile and thus the winding number to be unambiguously determined [16,42]. Further details on this measurement method are described in an earlier work by the same experimental group [39]. For a given set of parameters {T, V b , t hold }, the measurement was repeated 16−18 times to calculate the mean winding number l . The decay timescale τ was then computed by fitting l as a function of t hold to an exponential model, for each temperature T and barrier strength V b .
The resolvable range of superflow lifetimes is limited by the sample size. Given the finite number of measurements and the range of t hold values considered in the experiment, there was only a finite range of τ that could be distinguished from infinitely fast decay (τ = 0) or no decay (τ = ∞). The largest value of τ distinguishable from τ = ∞ corresponds to the case where, for the largest value of t hold , only one of the 18 measurements registers a decay event. Similarly, the smallest value of τ distinguishable from τ = 0 corresponds to the case where, for the smallest value of t hold , all but one of the 18 measurements registers a decay event. This sets limits on the possible values of l , giving the range of τ values measurable by the experiment as 70ms τ 80s.

Theoretical model
Decay of superflow due to the presence of a perturbing barrier is an inherently out-ofequilibrium scenario, for which there are limited theoretical tools capable of capturing both quantum and thermal effects in three-dimensional multimode systems. Our model of the experiment reported in Ref. [15] is formulated within the c-field theoretic framework, which is inherently non-perturbative and therefore well suited to studying a range of out-of-equilibrium phenomena, at both zero and finite temperature [41]. In this section we describe both our c-field model and our numerical simulation procedure.

Classical field methodology
The essential idea of c-field methods is that the macroscopically occupied modes of a degenerate Bose gas can be well described by an equation of motion for a classically-valued field ψ. Formally, c-field theories are constructed by dividing the full quantum field theory into a low-energy band C, which contains all modes of high occupation, and a high-energy band I which contains the remaining sparsely-occupied modes. This leads to a decomposition of the field operator as:ψ whereψ C andψ I are field operators for the C and I regions, respectively. An energy cutoff cut defines the division of the field theory into C and I regions and a projector P, defined in a convenient single-particle basis, ensures the two regions remain separated dynamically, i.e. P{ψ} =ψ C . Classical field theories treatψ C as a classical field ψ, and thus neglect the discrete nature of the atoms within the C region. In contrast, the I region is treated as a static thermal reservoir. The dynamics of the c-field ψ can be determined via a phase-space correspondence that maps the equations of motion forψ C to equations of motion for ψ [41]. Formally, within the phase-space framework, ψ is a stochastic sample of the C region's approximate phase-space distribution, with expectations of physical quantities given by ensemble averages of moments of ψ. However, an individual sample of ψ can often be loosely interpreted as the outcome of a single experimental run where, for example, the density of the Bose gas is |ψ| 2 [43]. For further details regarding c-field methodology, and examples of applications to non-equilibrium phenomena in Bose gases, see Ref. [41] and references therein. Classical field methods have successfully modelled ultracold Bose gases at both zero and finite temperature [41,44]. For systems near zero temperature, the occupation of the I region is negligible and thus the dominant beyond-mean-field effect is often inherent quantum fluctuations of the Bose gas. Within this regime, zero-temperature truncated Wigner (TW) [45][46][47] is the dominant c-field approach. For the closed-system dynamics typical of many BEC experiments, ψ is governed by a Gross-Pitaevskii equation (GPE) with initial conditions sampled from the Wigner distribution of the initial state. For many T = 0 non-equilibrium phenomena it is sufficient to treat this initial condition as a multimode coherent state that is sampled by seeding the initial mean-field condensate wavefunction with on average half an atom of vacuum noise per mode [48]. Zero temperature TW with this initial condition has successfully modelled BEC dynamics in regimes where nonclassical particle correlations become important [49][50][51][52][53][54][55][56][57][58]. For finite-temperature studies, the relevant c-field theory is the stochastic projected Gross-Pitaevskii equation (SPGPE), which describes interactions between degenerate modes of the quantum field with a static thermal reservoir [59]. The SPGPE and its sub-theories have been used extensively to study Bose gases both in and out of equilibrium, such as in Refs. [60][61][62] and Refs. [41,[63][64][65][66][67][68][69][70][71][72][73][74][75][76][77][78], respectfully. Notably, the SPGPE has been able to quantitatively describe experimental results, such as in Refs. [60,[79][80][81][82]. The SPGPE is typically applied to systems with a large thermal fraction, usually at temperatures ranging from T ∼T c /2 (where T c is the critical temperature of condensation) to just over T ∼T c .
The experiment of Kumar et al. that we model in this work investigated superflow decay in the presence of a relatively small thermal cloud; specifically, the experiment studied atomic superfluids at temperatures between T = 30nK (∼0.05T c ) and T = 195nK (∼0.4T c ) [15]. Unfortunately this is outside the typical regime of validity for both zero-temperature TW, which assumes a negligible thermal cloud, and the SPGPE, which assumes a larger proportion of thermal atoms. This is a low-temperature regime where neither quantum or thermal fluctuations truly dominate over one another. To address this regime, our model combines zero-temperature TW and SPGPE theory to include both quantum and thermal fluctuations, following the approach outlined in the Supplemental Materials of Ref. [58]. In our model, initial states are first sampled from the grand canonical ensemble of a Bose gas at thermal equilibrium using the SPGPE. Quantum fluctuations are then included by adding half a quantum of vacuum noise to each mode of the sample, as is done when sampling a coherent state in TW [48]. These initial states are then evolved using the SPGPE, with parameters describing reservoir interactions estimated from the atom number and temperature reported in the experiment. This model reduces to zero-temperature TW at very low temperatures, and to SPGPE at higher temperatures, both of which are expected to be quantitative models in their regimes of validity.
Additionally, a quantitative description of Ref. [15] requires the three-dimensional form of the SPGPE. This is because the vertical confinement in the experiment is not sufficiently large that all excitations in the z dimension are suppressed. Specifically, the trapping parameters of the experiment do not satisfy the condition ω z µ, which is required for an effective two-dimensional model to be a quantitatively correct description of the Bose gas.

SPGPE Theory
The simple-growth stochastic projected Gross-Pitaevskii equation can be written as: where µ is the chemical potential of the reservoir, dξ γ (x, t) is a complex Gaussian noise of mean zero and correlation and is the Gross-Pitaevskii mean-field operator in cylindrical coordinates. Here g = 4π 2 a s /m is the atom-atom interaction strength where a s ≈ 52a 0 for 23 Na, which was the atomic species used in the experiment. The term U t (r, z) corresponds to the toroidal trapping potential, which is approximately harmonic: where r 0 = 22.46µm is the mean radius of the ring in the experiment of Kumar et al. The form of the SPGPE we use in this work, Eq. (2), neglects number-conserving scattering interactions between the C and I regions, which are often referred to as the 'scattering' or 'energy-damping' terms [83]. These terms are commonly neglected in studies with SPGPE under the assumption that they do not significantly contribute to system dynamics near equilibrium [41]. We have confirmed that these terms have little quantitative effect on the results of this work (see Appendix D). The projector P restricts the c-field ψ to the low-energy subspace C defined by cut . This energy cutoff is typically chosen such that the highest-energy single-particle modes contained within the C region have an occupation of roughly n cut ∼1−10. We choose the energy cutoff by inverting the Bose-Einstein distribution for a fixed n cut [84]: In all our calculations, we fix n cut = 1, which means that modes with average occupation of n ≈ 1 and greater are included within the C region (see Fig. 11 in Appendix B). The dimensionless damping strength γ in Eq. (2) gives the rate of reservoir interactions, and characterises the speed at which the c-field approaches thermal equilibrium with the thermal reservoir. It can be a priori determined from the chemical potential of the reservoir µ, the reservoir temperature T , and the energy cutoff cut [63]: where is the thermal de Broglie wavelength, and Φ[z, x, a] is the Lerch transcendent. If γ = 0 and the projector is neglected, then Eq. (2) becomes the Gross-Pitaevskii equation (GPE).
We numerically implement Eq. (2) in a basis of approximate single-particle modes of the toroidal trap, extending the approach outlined for the projected GPE in Ref. [85]. A notable benefit of this approach is that the energy cutoff can be consistently defined in this single-particle basis, which diagonalises the Hamiltonian at high energies. This ensures that the occupation condition required for SPGPE validity is met even if there are regions where the spatial density vanishes. Furthermore, only a relatively small number of single-particle modes are required in order to represent the C region, compared to the much larger number of points required to represent the c-field on a three-dimensional cartesian grid. This enables the detailed finite-temperature 3D simulations performed in this work which, given the long decay timescales observed in the experiment, are impracticable using grid-based methods. Details of our simulation implementation are given in Appendix A. A description of the simulation parameters, and a justification of the parameter values chosen, is found in Appendix B. Unless otherwise stated, SPGPE simulation results reported in this work are averages over an ensemble of 96 trajectories. All reported observables are accompanied by a 95% confidence interval, estimated as twice the standard error of the ensemble averages.

The winding number
It is well known that the circulation C around some closed loop D in a superfluid is quantized in integer multiples of 2π /m [86]: where is the superfluid velocity field in cylindrical coordinates and the winding number l is an integer. The quantization of circulation is due to the relationship between velocity and the phase of the superfluid order parameter, as the phase must change by an integer multiple of 2π after 'winding' around the loop D.
For circulation around the annulus of a toroidal superfluid, the winding number is a topological quantity. This is because the torus is a multiply-connected topological space, which allows winding of the superfluid phase without the existence of vortices. In Fig. 2, we show exemplary phase profiles for the zero circulation state (l = 0), and the l = 1 circulation state with 2π phase winding.
In our numerical calculations, we compute l in each stochastic trajectory of the SPGPE and take the ensemble average to get the average winding number l . For convenience of l=0 l=1 Figure 2: A visualisation of exemplary phase profiles in the z = 0 plane in a toroidal geometry, for quantized circulation states with winding number l = 0 (left) and l = 1 (right). For the l = 1 circulation state, the phase increases by 2π as it winds around the central hole. The phase is only well defined in the region where there is significant superfluid density, and thus the center of the torus is omitted in this image.
calculation, we choose a loop around the trap minimum, which gives the following expression for the winding number (in cylindrical coordinates): where E[ * ] denotes an ensemble average over stochastic trajectories. In any given trajectory, we calculate the angular velocity field v θ by dividing the particle current by the density:

Initial state generation
The initial conditions for our simulations of Eq. (2) are stochastic samples of the C region's quantum state at thermal equilibrium in the grand canonical ensemble. It is common for finite-temperature studies in the c-field framework to neglect quantum fluctuations in the initial state, on the assumption that they are dominated by thermal fluctuations in the hightemperature regime where SPGPE is typically applied. However, this assumption is not appropriate for the regime studied in Kumar et al., where neither quantum or thermal fluctuations dominate over the other. To include both quantum and thermal effects, we add half an atom of vacuum noise per mode to initial samples of the grand canonical ensemble. This approach is similar to that used in Ref. [87] and the Supplemental Material of Ref. [58], and is akin to sampling the Wigner distribution of an incoherent mixture of coherent states that reproduce the statistics of the grand canonical ensemble. Specifically, where φ therm is the complex amplitude of a multimode coherent state |φ therm , sampled such thatρ = E[|φ therm φ therm |] is the density matrix for a thermal equilibrium state in the grand canonical ensemble 1 , and η is a complex Gaussian noise satisfying the correlation Each sample φ therm is given by evolving the simple-growth SPGPE to equilibrium. Physically, the SPGPE describes the exchange of particles and energy between the c-field ψ and a static thermal reservoir of chemical potential µ and temperature T . Consequently, it eventually evolves any initial state to thermal equilibrium in the grand canonical ensemble [41]. The value of the number-damping strength γ does not have any impact on the equilibrium properties in SPGPE theory, and thus may be chosen to give rapid convergence to equilibrium. In our simulations φ therm is sampled by evolving Eq. (2) with γ 0 = 1.0 for 100 trapping periods. Crucially, this sampling method provides the state of an interacting thermal Bose gas at equilibrium within the c-field approximation. The combination of the noise η and the stochastic noise in the thermal state sampling of φ therm gives an initial simulation state that includes both quantum and thermal fluctuations in the grand canonical ensemble. Projecting the noise η via the projector P{ * } ensures that the c-field remains within the low-energy subspace C. In practice, the projector is implicitly included when adding the noise in the single-particle basis (analogous to the projection of the noise term in the SPGPE, which is described in Appendix A).
A consequence of including quantum fluctuations in our initial state is that there will be formal corrections in the calculation of observables from the c-field ψ, due to the noncommutativity of the quantum field operatorsψ C with their conjugateψ † C [41]. However, due to the large atom numbers studied in this work (on the order of 10 5 ), these corrections are small and can thus be neglected.

Phase imprinting
In the experiment, an l = 1 circulation state is prepared by 'stirring' the perturbing barrier around the superfluid (see the Supplemental Material of Ref. [15]). Assuming this procedure perfectly prepares the metastable l = 1 circulation state, we may model this by instantaneously imprinting a 2π phase winding on our initial state: with l = 1. Modes with energy above the cutoff cut are incoherent, and so are unaffected by this transformation, and thus the thermal reservoir is treated as non-rotating.

Perturbing barrier and experimental sequence
In the experiment, the perturbing barrier is created by dithering a Gaussian beam in the radial direction such that its time average is [31]: 1 Explicitly, Here φ therm (r) is the position-space representation of |φ therm .
as a function of simulation time. In each simulation, the initial state is first generated (shaded blue) by evolving the SPGPE to equilibrium for t < t 0 . At t = t 0 , the initial state sample is seeded with quantum fluctuations and imprinted with a 2π phase winding. After some time t 1 − t 0 , the barrier is linearly ramped up over a period of time t ramp until it reaches its maximum height (χ = 1) at t = t 2 . The barrier is held at that maximum height for time t hold and then linearly ramped down over the period t 3 < t < t 4 .
where w = 6µm is the 1/e 2 half-width of the Gaussian beam and l d = 21.8µm is the width of the dither. The bracketed term ensures that the barrier vanishes at the edge of the torus. In our simulations, we choose the barrier to be centered at θ 0 = − π 2 . The time-dependent element 0 ≤ χ(t) ≤ 1 describes the raising and lowering ('ramping up/down') of the barrier. Importantly, the perturbing barrier does not affect the thermal reservoir interactions described by the SPGPE, as the peak barrier strength is weak compared to the energy cutoff, i.e. satisfies V b 2 cut /3 [63]. Indeed, all the barrier strengths considered in this work are weaker than the chemical potential of the atomic cloud.
The simulation protocol closely follows the experimental sequence and is shown in Figure 3. After the initial state is prepared as described in Sec. 3.4, the barrier is ramped up in strength over some time t ramp , held at its maximum value (χ = 1) for t hold , and then ramped down over Figure 4: Single trajectory of the SPGPE for a simulation of the experimental sequence at T = 85nK with the barrier held at its maximum strength V b = 0.82µ for t hold = 600ms. Slices of the density (peak density normalised to 1) and phase in the x − y plane are shown in (a) and (b), respectively, at different times. The winding number as a function of time is shown in (c), with the initial state generation period t < 0 (shaded blue) and barrier ramping periods (shaded red) shown. At time t = 1ms (i), the superfluid is in the metastable l = 1 circulation state, with a 2π winding of the phase profile. The barrier is raised in strength over the period 100ms < t < 170ms. At t = 170ms (ii), the barrier reaches its maximum strength and the density in the region of the barrier is visibly depleted. The superflow decays stochastically from l = 1 to l = 0 during the period where the barrier is held at its maximum 170ms < t < 770ms. After the decay event occurs (iii), the density remains depleted in the region of the barrier, however the 2π phase winding vanishes. The barrier is then ramped down over the period 770ms < t < 840ms, at the end of which (iv) the density in the region of the barrier is restored.
As noted by Mathey et al., zero-temperature mean-field theory (the GPE) predicts there is either no decay over the lifetime of the experiment or rapid decay (in less than 10ms), depending on whether the value of V b is above some critical value or not. This suggests the dynamics of the superflow is driven by Hamiltonian dynamics well captured by mean-field theory, for very small and very large barrier heights. This behaviour is to be contrasted with the superflow decay observed both in the experiment of Kumar et al. [15] and the finite- temperature simulations conducted by Mathey et al. [32], where stochastic decay events were observed at both short and long time scales, specifically across the entire time period where the barrier was held at its maximum. We observe qualitatively similar stochastic decay events in our simulations. As an example, the results from a single trajectory of the SPGPE are shown in Fig. 4, in which a decay event occurs roughly 110ms after the barrier reaches its maximum height. In Fig. 5 we present the results of SPGPE simulations for the four different temperatures T = 30, 40, 85, 195nK studied in the experiment 2 . For reference, the critical temperature is T c = 370nK for T = 40, 85nK and T c = 470nK for T = 30, 195nK, as given in the Supplemental Material of Ref. [15]. For each of the four temperatures, simulations are run for several different maximum barrier heights within a range of no more than ∼0.05µ. In any given trajectory, the transition of the winding number from l = 1 to l = 0 occurs stochastically, with the average value decaying slowly over a range of timescales ∼100ms−10s. Across all temperatures, Fig. 5 shows an extreme sensitivity to small changes in the barrier height. This is in qualitative agreement with the experimental observations, which found that tuning the barrier height over a range of ∼0.1µ changed the superflow decay rate by orders of magnitude [15].
In Fig. 6, we show that the decay of the average winding number for T = 85nK is similar to the experimental results of Kumar et al., albeit for different values of barrier strength V b . This quantitative discrepancy between the predictions of our model and the results of the experiment is analysed further in the following Section.
In both the experiment of Kumar et al. and the simulations by Mathey et al., the decay of the superflow while the barrier is at maximum strength appear well-fitted by decaying exponential trends. In Fig. 5 we can see that the results of SPGPE simulations for the lower temperatures studied do not seem to be well-suited to an exponential fit. While an exponential decay is not necessarily expected in the presence of nonlinearities, this behaviour may also be an artefact of ensemble averaging over a relatively small number of stochastic trajectories (96), or it might be suggestive of approximations made in the derivation of the SPGPE breaking down at low temperatures and long simulation times. The latter is discussed in more detail in Sec. 5.1. However, fitting these non-exponential trends to a simple exponential is still useful in estimating the timescale of the decay. In Appendix C we show that the trends are better captured by a two-timescale fit. However, it is not clear whether these 'timescales' are physically motivated, and thus it is difficult to compare them to the experimental results.

Blue-detuned trap
Red-detuned trap Figure 7: Comparison of the superflow lifetime τ between simulations (closed circles) and experiment (open circles) as a function of normalised barrier height V b /µ, for the four different temperatures studied in the experiment. Error bars on the experimental data give a 95% confidence interval. The grey shaded areas denote the timescales that sit outside the estimated range of values detectable by the experiment. The simulation and experimental data show strong quantitative agreement for T = 30nK, however there is an increasing discrepancy at higher temperatures. Specifically, the simulations suggest a larger barrier height is required to model a given decay time, and that this 'offset' increases with temperature.
Furthermore, the dominant timescales calculated within this approach are not significantly different to the single-timescale model, and thus there is no apparent advantage in using a two-timescale model in this study.

Quantitative comparison to experiment
As described in the section above, the timescale of superflow decay τ can be quantified by fitting the average winding number to the trend where c 0 and τ are positive-valued fitting parameters 3 . The fit is over the window of time t 2 ≤ t ≤ t 3 (see Fig 3), which is the period where the barrier is held at its maximum height.
In the simulations performed in this work, the barrier is held at its maximum for either t hold = 1s or t hold = 2.5s, depending on the magnitude of the barrier height. This ensures that the longest decay timescales τ accessible in the numerical data are the same order as the longest decay timescales observed in the experiment. The comparison of the computed timescale τ to that measured in the experiment is shown in Fig. 7, for each of the four temperatures. We can see that the timescale changes by four orders of magnitude over a small range of normalised barrier heights (V b /µ), with roughly exponential trends (linear trends in log-linear space). In particular, the simulation results for T = 30nK closely agree with the experiment. For higher temperatures, there is greater discrepancy between the experiment and simulation results, although the trends remain qualitatively similar. For instance, there is a discrepancy in the sensitivity of τ to changes in the barrier height (i.e. the slope of the trends in Fig. 7); this discrepancy is largest for the intermediate temperature T = 40nK and T = 85nK. More significantly, the trends in Fig. 7 suggest that our model requires a larger barrier height to reproduce superflow decay with a given timescale τ . This 'offset' in the barrier height is largest for the highest temperature, T = 195nK.
To quantify these disparities, we fit the trends in Fig. 7 to the following exponential function: The dimensionless fit parameter α captures the sensitivity of the decay timescale to changes in the barrier height. As shown in Fig. 8, the values of α for the simulation results are larger than the experimental results for all temperatures. Furthermore, the trends are qualitatively different; the experimentally-estimated values of α increase monotonically with temperature, which is not reflected in the simulation results. The quantitative discrepancy is greatest for T = 40nK, with simulations overestimating α by an order of magnitude. Agreement is strongest for the highest temperature studied T = 195nK, with the simulation value lying within 100% of the experimental value. This may be in part due to the exponential function being a poor fit to the simulation data for the lower temperatures (c.f. Fig. 5).
To quantify the discrepancy in the magnitude of the barrier height at which superflow decay begins to occur (i.e. the 'offset' seen in Fig. 7), we introduce the quantity κ t , which is the normalised barrier height for which the model predicts a one second superflow lifetime. This quantity is estimated by inverting Eq. (16) to get κ t = V b /µ for τ = 1s. In Fig. 9 we can see that the discrepancy between the value obtained from the simulations and experimental data ∆κ t increases monotonically with temperature, with strong agreement for the lowest temperature T = 30nK.

Discussion
At each temperature studied, save perhaps for the lowest temperature T = 30nK, we have observed significant discrepancy between the predictions of our model and the experimental results of Ref. [15]. For the two intermediate temperatures studied, T = 40, 85nK, the observed discrepancies could feasibly be due to the inadequacy of our model in the lower temperature limit (see Sec. 5.1.1). However, for the lowest and highest temperatures studied, we expect our model to give quantitative agreement with the experimental results, as the validity of zero-temperature TW and SPGPE theory has been well-demonstrated at low and high temperatures, respectively. This suggests that the large discrepancy observed between the results of our simulations and the experimental data for the largest temperature, T = 195nK, is likely due to additional experimental features not captured by our model.
In this section we will detail possible sources for the discrepancy, including assumptions and limitations of our theoretical model as well as technical effects within the experiment.

Validity in the low-temperature limit
As SPGPE theory is formulated in the high-temperature regime (∼T c /2−1.1T c ) where there is a significant thermal fraction of the gas, it is not a priori valid in the low-temperature regime of the experiment. Specifically, within the derivation of the SPGPE, reservoir interactions are expanded in powers of (µ − L)/k B T , and truncated at low order, the validity of which is well satisfied by the requirement k B T µ. The truncation may be valid for much lower temperatures if the system is not far from particle equilibrium with the reservoir. Since the experiment does not operate in the high-temperature regime (for all temperatures except T = 195nK, k B T < µ) it is possible this truncation discards important reservoir interactions for the two intermediate temperatures T = 40, 85nK studied in this work. In particular, omission of these higher-order reservoir interactions may be the source of the discrepancy in the sensitivity of the decay timescale on barrier height, as shown in Fig. 8. In addition, it may also account for the non-exponential nature of the trends for the lowest temperatures (T = 30, 40nK) in Fig. 5. Extensions of the SPGPE which include the higher-order terms will likely be challenging to implement numerically, however may be a useful avenue of future work.

Truncated Wigner Approximation
The SPGPE is derived by mapping a high-temperature master equation for the C region's quantum state to a partial differential equation (PDE) for the state's Wigner function. In general, this PDE suffers from the same computational intractability as the original master equation. However, by making the truncated Wigner approximation, which neglects thirdand higher-order derivatives in this PDE, the resulting equation of motion for the Wigner function takes the form of a classical Fokker-Planck equation that can be efficiently simulated via the SPGPE. The validity of our SPGPE simulations therefore depends in part upon the validity of the truncated Wigner approximation.
Although the truncated Wigner approximation is an uncontrolled approximation 4 , its validity has been verified in the classical field regime where the particle occupation per mode is high [41,89,90]. This condition is easily fulfilled in our work by the inclusion of a projector and our choice of energy cutoff. Nevertheless, the truncated Wigner approximation is only valid for a finite simulation time, since the (unquantifiable) error in the truncation will compound with time [49]. The simulations performed in this work have been over a long timescale, relative to typical cold-atom experiments, and thus it is feasible that the truncated Wigner approximation may begin to appreciably breakdown for the longest simulations (t hold =2.5s). This effect may have influenced the long-time behaviour of the T = 30nK trends shown in Fig. 5, which after about t≈2s deviate significantly from the behaviour noted for all other temperatures. However, we do not expect the truncation error to significantly contribute to the discrepancies noted at higher temperatures, since the thermal decoherence rate is likely much faster than the rate at which the truncated terms become significant [91].

Static thermal reservoir
Our SPGPE model assumes that the thermal reservoir (high-energy modes above cut ) is static and unaffected by the dynamics of the C region. The typical justification for this is that the sparsely-populated high-energy modes equilibrate rapidly compared to the macroscopicallyoccupied modes of the C region. While there are many systems for which the dynamics of the thermal cloud are an essential aspect of the physics (for example, collective modes of the condensate and thermal cloud [92]), there is currently no formulated extension of SPGPE theory that includes the dynamics of the I region. To study these effects, an alternative theoretical framework must be used, such as the coupled condensate-thermal theory of Zaremba, Nikuni, and Griffin [93]. However, it is incorrect to say the SPGPE does not treat the thermal cloud dynamically. The majority of thermal atoms are contained within the degenerate modes of the C region; it is only the high-energy sparsely-populated thermal modes that make up the thermal reservoir. In this sense, the majority of the thermal cloud is treated on the same footing as the condensate in SPGPE theory. Given the small thermal fractions considered in this work, most of which are contained within the C region, it is unlikely that neglecting non-equilibrium dynamics of the I region significantly affects the simulated decay timescales.
Despite these small thermal fractions, the atoms in the I region cannot simply be neglected. This is the approach adopted when using the projected GPE (PGPE), a c-field theory that only describes the portion of the thermal cloud within the C region [94,95]. A key attraction to this approach is its computational simplicity relative to the simple-growth SPGPE; its lack of dissipative and dynamical noise terms make it less challenging to numerically integrate than the SPGPE. In situations where the thermal cloud is small and the system is near equilibrium, the finite-temperature PGPE has been able to quantitatively describe experiments both in and out of equilibrium [96,97]. Within the context of this work, the PGPE captures the qualitative nature of the stochastic decay observed in the experiment of Kumar et al. [15]. However, quantitatively it provides a poorer description of the experiment than our SPGPE model, due to the lower rate of dissipation relative to the simple-growth SPGPE (see Fig. 18 in Appendix D). This, alongside the strong temperature-dependence of the experimentallyobserved superflow decay, suggests that it is important to retain dissipation due to scattering processes with atoms in the I region. Indeed, increasing this dissipation via increasing γ could potentially improve the agreement between our model and the experiment, although there is no a priori justification for increasing γ in this way.

Parameter estimation
In this work, the SPGPE parameters that describe reservoir interactions are estimated selfconsistently using the experimentally-measured atom numbers and temperatures. Nevertheless, to account for possible uncharacterised errors in these measurements, we investigated the sensitivity of our simulation results to changes in several parameters, including barrier width and temperature. Changes to these parameters do not strongly change our results, suggesting that the imprecise experimental calibration of these parameters cannot account for the discrepancy between our simulations and experimental data. We also varied the energy cutoff, which is not an experimental parameter, and confirmed that our results were not strongly cutoff dependent. Appendix B explains how we estimated our parameter values and fully details the effect of parameter variations on our simulation results. σ N /N = 11.1% σ N /N = 22.3% Figure 10: The effect of shot-to-shot fluctuations in the total atom number, for a barrier height of V b = 0.78µ and a temperature of T = 40nK. When the chemical potential is kept constant from shot-to-shot, µ =μ, there are fluctuations in the atom number on the order of 0.1%, due to sampling the initial state from the grand canonical ensemble. Larger shot-to-shot atom number fluctuations are modelled by randomly sampling the chemical potential µ from a Gaussian distribution with meanμ. As the fluctuations increase in magnitude, the average winding number begins to decay earlier in time. However, the overall timescale for decay is not significantly affected.

Shot-to-shot number fluctuations
Shot-to-shot fluctuations in the total atom number is one example of a common technical effect in experiments with ultracold atomic gases, the magnitude of which can certainly be temperature dependent. In fact, fluctuations in atom number were proposed as a possible cause of the discrepancy between experimental results in a related experiment (see the Supplemental Material of Ref. [39]). The simulations performed in this work include some fluctuations in the atom number, as the initial thermal state is sampled from the grand canonical ensemble. However, these fluctuations are no more than 3% for T = 195nK, and less than 1% for T = 30nK. We can increase the magnitude of shot-to-shot atom number fluctuations in our simulations by sampling the chemical potential µ from a Gaussian distribution with meanμ and width chosen to give a certain variance in the total atom number. As an example, we study the inclusion of atom number fluctuations for a fixed barrier height V b /μ = 0.78 and T = 40nK, in Fig. 10. We find that the average winding number is not significantly affected for shot-to-shot atom number fluctuations of magnitude 5−20%, which is the size observed in typical ultracold atom experiments, suggesting that such fluctuations cannot account for the noted discrepancies.

Barrier calibration
As shown in Fig. 9, the discrepancy between the experimental data and our simulations in κ t (the estimated value of κ = V b /µ that gives τ = 1s) increases with the temperature of the Bose gas. This could be caused by a temperature-dependent systematic effect in the barrier calibration, which is certainly possible as the experimental calibration of the normalised barrier strength V b /µ did not differentiate between condensate and thermal atoms.
The magnitude of this error may be estimated by computing the contribution of the thermal atoms to the chemical potential µ within a semiclassical approximation. This calculation is included in the Supplemental Material of Ref. [15], where the systematic shift in the barrier calibration was found to be approximately 3% for T = 85nK and 8% for T = 195nK. Not only is this smaller than the statistical error in the experimental calibration, but it is far smaller than the discrepancy in the range of barrier heights for which decay occurs between the experimental data and the simulation results, which is roughly 35% and 50% for T = 85nK and T = 195nK, respectively (see Fig. 9). Therefore, an unaccounted shift in the barrier calibration due to the presence of a thermal cloud cannot account for the discrepancies between the experimental results and the predictions of our model. Atom losses during the experiment could also have caused a discrepancy between the normalised barrier strengths V b /µ reported in the experiment and the values used in our simulations. The experimental report estimates that the ∼25s lifetime of the experiment reduces the chemical potential by ≈ 10%, resulting in a slight increase in the effective barrier strength [15]. However, similar to the above, the size of this effect is not sufficient to account for the magnitude of the observed discrepancy. Assuming this lifetime is vacuum limited (as were previous experiments from the same laboratory [10]), the effect of atom loss would be a common shift for all temperatures, and thus cannot account for the temperature dependence of the discrepancies observed between our simulations and the experiment.

Optical trap imperfections
In general, technical noise on the trap lasers causes heating of an optically trapped superfluid, which may lead to enhanced dissipation of superflow. Furthermore, in this experiment the superflow decay rate is particularly sensitive to small changes in barrier height. This could amplify the effect of slight violations of the trap's azimuthal symmetry in the experiment [98]. This breaks the symmetry of the ideal ring considered in the simulations, allowing more pathways for vortex escape, thus enhancing the rate of superflow decay. This may have contributed to the observation of superflow decay in the experiment at lower barrier heights than predicted by our model.
The experiments performed at higher temperatures used different optical trapping beams to the experiments performed at lower temperatures. As described in Sec. 2, the low temperature T = 30, 40nK BECs were confined in a blue-detuned dipole trap, whereas the higher temperature T = 85, 195nK BECs were confined in a red-detuned dipole trap. In red-detuned traps, atoms are located where the light intensity is maximal, and are thus intrinsically more sensitive to imperfections and noise in the trapping laser fields. The red-detuned beams used in the experiment suffered from etaloning in the vacuum cell, which resulted in fringes in the trapping potential [98,99]. The resulting effect on the measurements performed in reddetuned traps was not modelled in our higher-temperature simulations, and could account for the stronger temperature dependence of superflow decay observed in the experiment.

Conclusions
In this work, we have performed detailed three-dimensional classical-field simulations to model the experiment of Kumar et al. [15]. Our model, which combines zero-temperature TW methodology and SPGPE theory, describes the role of quantum and thermal fluctuations in the spontaneous decay of persistent currents of a superfluid BEC trapped in a toroidal geometry. We have demonstrated that our model is able to capture the essential non-equilibrium dynamics that lead to superflow decay in the presence of a perturbing barrier, in good qualitative agreement with the experimental results of Kumar et al. and a previous theoretical analysis of a related experiment [32]. Specifically, the predictions of our model are in quantitative agreement with the experiment at low temperature, and provide a qualitative description of the experiment at higher temperatures. Furthermore, our simulations predict the same range of decay timescales as observed in the experiment, across all temperatures studied. Notably, this is achieved with simulation parameters estimated solely from the experimental temperature and atom number.
For the lowest temperature Bose gas studied in the experiment, the decay timescales predicted by our model quantitatively agree with the experimental data. For the other temperatures studied, however, we have found discrepancies between the quantitative predictions of our model and the experimental data, which become most significant at the highest temperature. As discussed in Section 5.2, this is likely not solely due to limitations of the theoretical model, as there were several technical effects in the experiment that may have led to enhanced superflow decay.
In general, obtaining quantitative agreement with experimentally measured decay rates is very difficult to achieve due to the many sources of dissipation in a real superfluid. In particular, modelling the experiment of Ref. [15] was a challenging task, as a three-dimensional model was required that included the effects of both quantum and thermal fluctuations. Moreover, we have pushed the SPGPE model beyond its safe regime of applicability for studying the dynamics of a Bose gas in the low-temperature regime of the experiment, and noted a range of technical effects that could lead to enhanced dissipation of the superflow. Despite this, we have observed some level of quantitative agreement without fitted parameters for low temperatures, providing further evidence on the value of a c-field description of highly non-equilibrium dynamics in Bose gases. Nevertheless, our work suggests a need for further theory beyond SPGPE, alongside further experimental characterisation of superflow decay at higher temperatures.

A Numerical methods
Within classical-field methodology, the numerical implementation of a well-defined energy cutoff is crucial in order to make quantitative physical predictions [41]. This can be achieved by projecting the equation of motion for the classical field onto a basis where the many-body Hamiltonian is approximately diagonal. At the high energies where the cutoff is typically imposed, this is satisfied by the eigenbasis of the single-particle Hamiltonian H 0 . This corresponds to the decomposition of the classical field as: where Φ αnΓ (r, θ, z) are the single-particle modes of H 0 .
In the context of toroidal confinement, the eigenstates of the potential in Eq. (5) are not analytically known. However, recent work by Prikhodko et al. [85] has demonstrated that an approximate single-particle basis may be used: where ϕ α are the normalised Hermite-Gauss functions (in physical units): using the physicists' Hermite polynomials H n (x). As described in Ref. [85], this basis is approximately orthonormal provided it is truncated such that the highest-energy mode vanishes as r → 0. In this limit, this basis diagonalises the single-particle Hamiltonian with the energy spectrum This allows us to define the low-energy region C by only including modes with energies below the cutoff: Numerically, the projector is implemented by setting the occupations of single-particle modes with energies above cut to zero.
Casting the projected GPE (Eq. (2) with γ = 0) onto this single-particle basis is detailed explicitly in Ref. [85], so we will not describe it here. Implementing the number-damping reservoir interaction terms in this basis is a straightforward extension of their method. For the deterministic terms only the constant prefactors need to be adjusted. The number-damping noise term can be directly sampled on the single-particle basis via: where dW αnΓ is a complex Weiner noise satisfying Similarly, the quantum fluctuations seeded in the initial states (Section 3.4) can also be directly sampled on the single-particle basis: where η αnΓ is a complex Gaussian noise satisfying We numerically integrate the SPGPE using the open source XMDS2 software package [100], exploiting an adaptive fourth-fifth order Runge-Kutta algorithm. The use of a high-order adaptive Runge-Kutta algorithm is appropriate due to the additive nature of the noise term in the simple-growth SPGPE, as first noted in Appendix B of Ref. [63]. In all the simulations reported in this work, the relative error tolerance of the algorithm is set at 10 −5 .
Transformations between the single-particle basis and spatial grids (where the nonlinear |ψ| 2 ψ term is diagonal) are implemented using in-built Hermite-Gauss and Fourier transforms in XMDS2. Due to the presence of the 1/ √ r in the single-particle basis -Eq. (18) -Hermite-Gauss quadrature methods are inexact for computing spatial integrals with a finite number of grid points. To minimise this as a source of numerical error, we include an additional 16 points on the spatial quadrature grids.

B Fixing simulation parameters
The estimated simulation parameters for each of the four temperatures studied in this work are given in Table 1.

B.1 Chemical potential
In our simulations we choose the value of the chemical potential µ such that the initial thermal states have an average atom number as close as possible to the experimentally reported value, for each temperature. Although this can be determined through a semiclassical calculation (see Appendix A of Ref. [64]), we find it convenient to simply vary µ until we achieve the T (nK) k B T / ω r ω z /2π (Hz) µ/ ω r cut / ω r γ n r × n θ × n z  Table 1: Fixed parameters for the simulations performed in this work. Here n r , n θ , n z are the number of single-particle modes in the r, θ, z dimensions, respectively. desired atom number. Specifically, we first estimate the chemical potential by assuming a purely Thomas-Fermi density: We then calculate the atom number at thermal equilibrium for a range of µ around this initial estimate, and fit the resulting trend to determine the value of µ that will give an atom number closest to the experimental value. Once the chemical potential has been set, the energy cutoff cut and number-damping strength γ can be estimated using Eqs. (6) and (7) in the main text.

B.2 Energy cutoff
As described in Sec. 3.2, the energy cutoff in Eq. (6) is chosen to give an average occupation of n ≈ 1 for the single-particle modes near the energy cutoff, as is typically done in SPGPE analyses. We check this is satisfied, for simulations of the experimental sequence at T = 40nK and T = 195nK, by computing the occupation of each single-particle mode |ψ αnΓ | 2 and plotting it as a function of its energy. This is shown in Fig. 11, where it is clear that the modes near the cutoff have an average occupation on the order of n ≈ 1. As a more comprehensive check of the energy cutoff, in Fig. 12 we assess the quantitative impact of changing the cutoff value slightly around its estimated value, for the temperature T = 85nK and a barrier height of V b = 0.77µ. Although the precise value of l changes slightly as cut is varied, the decay timescale is not significantly affected. Specifically, increasing the energy cutoff by 20.0% leads to τ reducing to roughly 45% of its value at the estimated cutoff. Given that changing the barrier height by ∼0.1µ changes τ by several orders of magnitude, this variation with cut is acceptably small.

B.3 Grid size
The size of the single-particle grid used in the simulations is set by the energy cutoff. Although it is straightforward to estimate the number of radial and axial modes required to satisfy the condition set by Eq. points for the angular grid. The highest energy angular mode allowed below the energy cutoff can be calculated by assuming all energy is in the angular direction, which gives: rounded to the nearest integer. Naïvely, this would suggest that 2|n| max grid points should be used with n ∈ [−|n| max , |n| max ]. However, the angular component of the single-particle basis is a plane wave and is therefore subject to Nyquist aliasing. Formally, aliasing of modes within the C region can be avoided by using n θ ≥ 4|n| max angular grid points. However, this is typically a large number (n θ ∼ 300−400), which results in restrictive computational requirements for the long-timescale three-dimensional simulations needed for our investigation.
In practice, we choose the grid size such that only a small number of angular modes are subject to aliasing, and only those modes that have a relatively small occupation. We find it is sufficient to use n θ = 150 in all our simulations. This is confirmed in Fig. 13, which shows that the aliased modes for simulations at T = 40, 195nK have an average occupation on the order of 1, and never greater than 10. This is significantly smaller than the occupation of modes of small n, which each contain 10 3 − 10 5 atoms. We have also checked for a full simulation of the experimental procedure at T = 40nK that doubling the number of angular modes does not quantitatively change the results of the simulation.

B.4 Sensitivity of results to barrier width and temperature
We investigate the sensitivity of our results to simulation parameters that are taken directly from the experiment, specifically temperature T and barrier width w. As shown in Fig. 14 predict that stochastic superflow decay no longer occurs within t hold = 1s. This may be due to the barrier width becoming so large that the dynamics of the superfluid around the barrier maximum are suppressed. However, given that the width of the atomic density depletion due to the perturbing barrier can be measured with a precision of order 0.1µm, it is unlikely that our estimate of w deviates beyond 10% of its true value.
To investigate the sensitivity of the results on temperature, we consider the experimental data at the hottest temperature, and run simulations across a range of barrier heights at a temperature of T = 225nK which is the reported experimental value plus its standard deviation (T + σ T ). The comparison of these results to the results for T = 195nK is shown in Fig. 15. As one may expect, increasing the temperature reduces the decay timescale at lower barrier heights, improving agreement with the experimental data slightly at those barrier heights. However, for larger barrier heights, we find that increasing the temperature by this amount has very little effect, and thus there is little to no improvement in agreement for the larger barrier heights simulated. Therefore, the precise choice of the simulated temperature within the experimentally quoted confidence interval can be dismissed as a possible origin for the discrepancy observed with the experimental data.

C Multi-timescale fits of decay
We investigate the effect on our results of fitting the average winding number as a function of time with a multi-parameter exponential model, as opposed to the single timescale model Eq. (15) used in the main body of the work. In Fig. 16  Exp. Figure 15: Decay timescale as a function of normalised barrier height, for the experimental results for T = 195nK, and simulation results for T = 195nK (red diamonds) and T = 225nK (blue squares). The temperature T = 225nK is the quoted experimental value plus its standard deviation, referring to the upper limit of a 68% confidence interval on the temperature of the experiment. There is a slight quantitative difference between the two simulation trends, with the higher temperature associated with a gentler slope. However, this shift is not sufficient to explain the discrepancy with the experimental data. Error bars shown here represent a 95% confidence interval, with simulation data obtained from an ensemble average over 48 SPGPE trajectories.
temperatures with models of the form where {c i , τ i } are free parameters of the fit. In general, a two-timescale fit of the form: appears to be sufficient to describe the trends for all temperatures. This is with the exception of some T = 30nK trends, which clearly have some non-exponential qualities that are not captured by any of these fitting models. Taking the two parameters τ 1 , τ 2 to represent slow and fast timescales, respectively (we enforce τ 1 > τ 2 ), we compare the results of our simulation to the experiment in Fig. 17. This approach does not lead to clear trends, and is challenging to interpret in a meaningful way. In fact, the confidence intervals on many of the two-timescale fitted data points are very large, further demonstrating that a two-timescale fit does not deepen the analysis of the results significantly over the use of a one-timescale fit. Further, there is not a clear physical explanation for why there would be two timescales that would govern the superflow decay mechanism, and thus the addition of additional timescales is entirely ad hoc. Finally, the values for the two-timescale fits do not significantly differ from the one-timescale fits, with values spanning the same range of magnitudes as the single-parameter timescales in Fig. 7.  Overall, this analysis shows that there is no advantage in using multiple timescale fits over single timescale fits for comparing our simulations to the experiment.

D Contribution of energy-damping terms
In our study, we have neglected number-conserving reservoir interaction terms in the SPGPE often referred to as the scattering or energy-damping terms. These terms are numerically challenging to implement and are often neglected in SPGPE theory under the justification that they are expected to be dominated by the non-number-conserving γ process [41]. Below we quantitatively confirm that their neglect is justified for the simulations performed in this work. Exp. Figure 17: Slow-decay timescale τ 1 (red triangles) and fast-decay timescale τ 2 (blue inverted triangles) as a function of normalised barrier height κ = V b /µ. Error bars give 95% confidence interval. Data points with standard deviation in τ greater than 10 2 s are excluded from this figure.
When energy-damping terms are included, the SPGPE is i dΨ = P (L − µ)Ψdt + iγ(µ − L)Ψdt + i dξ γ (x, t) The energy-damping terms (33) consist of a deterministic evolution term and a noise term, both of which are non-local. The deterministic term describes particle scattering via an effective potential: which is a convolution between the particle current, and the epsilon function The energy-damping noise is real-valued, multiplicative, and non-local in space:  Figure 18: Final value of the average angular momentum per particle (which is equivalent to the final winding number) for the temperature T = 195nK and a barrier hold time of t hold = 2.5s. There are several data sets corresponding to T = 0 GPE simulations (black squares), and various subtheories of the SPGPE. Note that M is given in units of l 2 0 , where l 0 = /mω r ≈ 1.3µm is the radial lengthscale set by the trap. Each data set is fitted to a sigmoidal function: L z / N = [exp ((V b /µ − α)/β) + 1] −1 , where α, β are the fitting parameters. Error bars give a 95% confidence interval in the data. Notably, the absence of energy damping (γ = 0, M = 0) does not deviate significantly from the full SPGPE result (γ = 0, M = 0).