Dynamical Equilibration Across a Quenched Phase Transition in a Trapped Quantum Gas

The formation of an equilibrium quantum state from an uncorrelated thermal one through the dynamical crossing of a phase transition is a central question of non-equilibrium many-body physics. During such crossing, the system breaks its symmetry by establishing numerous uncorrelated regions separated by spontaneously-generated defects, whose emergence obeys a universal scaling law with the quench duration. Much less is known about the ensuing re-equilibrating or"coarse-graining"stage, which is governed by the evolution and interactions of such defects under system-specific and external constraints. In this work we perform a detailed numerical characterization of the entire non-equilibrium process, addressing subtle issues in condensate growth dynamics and demonstrating the quench-induced decoupling of number and coherence growth during the re-equilibration process. Our unique visualizations not only reproduce experimental measurements in the relevant regimes, but also provide valuable information in currently experimentally-inaccessible regimes.

A central emerging question is the extent to which the defect formation process described by the Kibble-Zurek mechanism can be decoupled from the dissipative evolution expected to occur when defects co-exist within a finite region, which is crucial to the modelling of any experiment. Previous work on purely 1D evaporativelycooled gases has indicated the spontaneous emergence of dark solitons which subsequently decay in time, with the coherence length of the final post-quench thermal state being set by the average distance between defects [39]. Decoupling these two dynamical features is not a trivial task, as one does not a priori know exactly when to count the defects, and, in fact, at the time at which one is supposed to count the defects, the system is so disordered, that it is not even clear whether any counting would be possible, even at a post-processing stage. Nonetheless, numerous experiments with ultracold trapped atoms have provided convincing evidence of the emergence of the predicted power in the scaling of defect number with respect to quench rate, raising the question of why this scaling should survive over such prolonged dynamical periods. The Cambridge group [19] insightfully avoided such issues by looking instead at correlation functions, which enabled them to extract critical exponents, without the need for direct defect counting.
In this work, we offer a unified analysis of quenched condensate growth dynamics in a finite elongated 3D inhomogeneous system, incorporating into a single study the dynamical evolution, for different induced quench rates, from an equilibrium thermal state above the Bose-Einstein condensation transition temperature to a near-equilibrated, low-temperature phase-coherent Bose-Einstein condensate (BEC). By demonstrating the natural emergence of symmetry-breaking in our simulations, we perform a detailed parallel analysis of both the spontaneous emergence and complex nonlinear dynamics of defects, and the related evolution of coher- Starting from a purely thermal state, with a given atom number, we linearly quench temperature to lower values and chemical potential to positive values over a ramp duration τR, to mimic experimental conditions. b Dynamical response of an equilibrium thermal gas subjected to different cooling quench rates (τR = 1440, 144, 18 ms, from top to bottom), demonstrating the equilibration route towards a finite temperature phase-coherent BEC. The characteristic regions depicted here refer to density isosurfaces of the highestpopulated mode, chosen such that n(t)/n(t → ∞) = 0.1% (yellow), or 3% (green), where n(t → ∞) describes the final equilibrated condensate density for N = 6.6 × 10 6 atoms. Different rows correspond to different durations of the constant applied external cooling ramps, from the very slow, quasi-adiabatic (top) to the very fast, nearly-instantaneous ones (bottom), with the intermediate case representing typical quenches used in experimental studies of phase transitions. For rather slow ramps, most condensate formation happens during the external cooling. For shorter quench duration, the condensate appears around the end of the ramp, and there is a small number of spontaneously-generated defects. In a fast quench, most condensate growth dynamics occurs after the end of the ramp, with the system at the end of the ramp being in a highly non-equilibrium state exhibiting large occupation of a handful of modes, and consisting of a dense random defect tangle. Such a tangle unravels in time into a phase-fluctuating condensate, or quasi-condensate, with numerous well-formed defects, whose presence perturbs the phase and opposes the formation of long-range coherence; the interacting defects gradually evolve into few long-living defects characteristic of the system geometry (quantized vortices here) which are experimentally observable after expansion, with the system asymptotically expelling all defects and acquiring phase coherence across its length.

ence.
Addressing subtle open questions in the condensate formation process, we show that coherence and number growth dynamics are in general decoupled, due to competing growth mechanisms following a quench, except for cases of adiabatically-slow growth which exhibit broadly-similar timescales. The unique insight provided by our numerical visualizations, which extend beyond experimentally-accessible time intervals, provides a natural framework for addressing the still unresolved interplay between Kibble-Zurek defect generation and coarse-graining dynamics. In this regard, our simula-tions demonstrate that the anticipated power-law defect scaling with quench rate (characteristic of Kibble-Zurek) is not significantly affected over a prolonged evolution period, despite the fact that both the number and type of defects changes considerably during this period. Specifically, the initial randomly-oriented and -shaped spontaneously-generated defects gradually relax to the canonical excitations in the given system geometry (here solitonic vortices), with their number exhibiting a rapid decrease at very early stages, followed by a more gradual decay process as the size of the condensate grows further.
These findings are consistent with late-time experimental measurements performed within our group [17,21,53], which further validate our analysis.
In order for our predictions to be directly comparable to experimental measurements in the appropriate regimes, we perform full 3D numerical simulations of the stochastic (projected) Gross-Pitaevskii equation [15,29,45,46,[49][50][51][52] in the elongated configuration of our recent experiments [17,21,53]. Such a regime is particularly interesting as it facilitates unique access to the role of coherence in a finite-size, inhomogeneous and anisotropic system, for which the phase transition process becomes position-dependent, with the condensate first emerging in the central trap regions, where the phasespace density is maximised. Examples of simulations of the BEC growth for different quench rates are given in Fig. 1, clearly showing the dependence of defect generation and subsequent dynamical equilibration on quench rate.

A. Quench Protocol
Our experiment is conducted in a cigar-shaped trap with ω x /2π = 13 Hz and ω ⊥ /2π = 131.4 Hz, where 23 Na atoms in the |F, m F = |1, −1 state are evaporatively cooled across the BEC phase transition at different rates. To faithfully mimic both changing temperature and total atom number observed in the experiments, our simulations are based on linear quenches (Fig. 1a) over a finite quench duration τ R , both in temperature (T = 790 nK → 210 nK) and chemical potential (µ = −22 ω ⊥ → µ = +22 ω ⊥ ). Consistently, the atom number goes down from N = 22 × 10 6 to N = 6.6 × 10 6 . Our parameters have been chosen such that t = 0 corresponds to the time when the system crosses the ideal gas critical temperature (µ = 0).
For each quench rate, the results are analysed over 3 to 7 individual realisations, which are sufficient for understanding the underlying physics. Our analysis is based on the characterization of the emerging condensate, defined in our simulations as the mode with the largest eigenvalue of the one-body density matrix, based on the usual Penrose-Onsager definition [49,54]. We also note that at early evolution times there are a number of approximately-equally largely populated modes, before one becomes randomly dynamically favoured by the system. Details of our experimental configuration have been reported elsewhere [17,21], while the stochastic numerical method and data analysis scheme used to model such non-equilibrium dynamics are summarised in Methods. Fig. 1b shows that for a very slow quench (top), the system grows gradually to its final equilibrium state. In doing so it remains close to the corresponding equilibrium phase-coherent BEC during its evolution, except in a narrow time window within the region of crit-ical fluctuations. Such slow evolution corresponds to the evaporatively-cooled condensate growth scenario, in which both condensate atom number and correlation function grow gradually, with a small phase-coherent condensate present shortly after the phase transition. In such cases, the system equilibration happens in real time, on the ramp, i.e., during evaporative cooling. As the ramp speed is increased (middle and bottom images), we observe more clearly a rather violent falling out of equilibrium of the dynamical system when crossing the phase transition, resulting in the spontaneous production of multiple defects: the underlying physical picture here is that coherence only forms in local patches of constant -but random -phase, and the defects separate such regions of different phases: Faster quenches (bottom) lead to a faster overall initial growth in the atom number of the quantum-degenerate state, and a larger number of spontaneously-generated defects than slower quenches (middle), with most of the condensate (and also coherence) growth occuring after the removal of the external ramp.
As the system grows, such defects are stretched out while also undergoing complicated nonlinear dynamics in the inhomogeneous background, where they interact with other defects. As a result of this, the number of defects decreases gradually (with the rate of decrease depending critically on quench rate), allowing phase coherence to spread to the entire system. For the set of experimental parameters considered here, the final equilibrated state -obtained asymptotically for the case of near-instantaneous ramps -consists of a fully phasecoherent defect-free finite-temperature BEC with condensate fraction N 0 /N ≈ 0.75. (Characteristic singlerealisation growth sequences are shown in Supplementary  Fig. 1-2 and Movies 1-2).
In the following section we analyse in detail such distinct stages in the dynamical equilibration, offering further insight to this complex non-equilibrium process, while also demonstrating very good agreement between theory and experiment, for the observables and regimes where experimental measurements have been undertaken.

B. Dynamical Crossing of the Phase Transition
The transition region dynamics is naturally incorporated within our numerical approach. As the system is cooled at a finite rate, with its atom number decreasing (as in experiments), it cannot instantaneously acquire coherence across its entire spatial extent, and so a dynamical spontaneous symmetry-breaking occurs over a relatively short temporal range; this allows the system to become infilled by a densely-packed network of randomlylocated defects shown in purple around t ≈ 30 ms. The emergence of such defects signals a stark deviation from the corresponding equilibrium system profiles. This process is demonstrated in Fig. 2 for the specific example of Above the transition temperature, we see low-amplitude random density fluctuations (yellow) within our entire simulation grid, with the density condensing in real space to reflect the underlying elongated trap geometry. During such condensation process, the system breaks its symmetry spontaneously, leading to the appearance of defects (purple) within a growing high-density region (green), signalling the emergence of a highly-fluctuating condensate. Such defects undergo their own evolution dynamics, violently perturbed by interactions with other defects (Supplementary Fig. 3 and Movie 3).
a quench time τ R = 84 ms.
Consistent with experiments, our simulations clearly show the temporal narrowing of the density distribution during the phase transition, resulting in the gradual emergence of a higher-density condensate region (green region, t ≥ 31 ms) which gradually grows towards the trap edges. The origin of this is associated with the lowering of the system temperature, the incorporated atom loss and the presence of the harmonic confinement, all of which lead to a centre-peaked position-dependent increase in local phase-space density. However, looking at the finer level during such evolution, we find the defects also interact, coalesce and decay. Such processes, which are not included in the Kibble-Zurek mechanism are nonetheless crucial to the growth of the phase-coherent regions: In the early stages, while numerous defects are present, the system can be classed as a quasi-condensate [23,55] in the sense that different regions of coherent density exhibit no common phase between them, such that the observed coherence length is considerably smaller than the system size, consistent with the large population of a number of modes.
Our advanced numerical visualisation, which encapsulates the phase transition physics, reveals the complexity of attempting to extract, whether numerically or in actual experiments, the early physics of defect formation. This is of particular relevance for addressing the interplay between the simplified Kibble-Zurek model and the "coarse-graining" dynamics governing dynamical defect evolution and decay. Our analysis shows that one cannot decouple the two processes -at least not in the context of inhomogeneous systems.
Following the crossing of the phase transition and the spontaneous formation of defects, the emerging condensate enters into a rather distinct dynamical regime, in which the physics is dominated by a complicated nonlinear combination of (i) defect stretching due to the growing condensate size, (ii) defect propagation in an inhomogeneous environment, (iii) occasional but rather violent defect interactions, including vortex reconnections, bouncing and "ejection" from the condensate [56], and (iv) additional forced relaxation in cases of slow cooling. The combination of such processes leads to the gradual dynamical equilibration of the system, which evolves from a defect-filled condensate to a state at T T C defined by the final quench parameters.

C. Condensate Growth Dynamics
Early condensate growth experiments [25,33,34] and their numerical modeling [26,28,30] revealed a number evolution N 0 (t) well described by an S-shaped curve. However, slower quenches [31] and elongated geometries [57] were observed to feature a pronounced region of critical fluctuations, leading to a "time delay", or "onset time" for condensate growth and a slower (near-linear) initial growth rate [31]; although the presence of such features is well-documented and broadly attributed to the initial emergence of the quasicondensate, there has been little quantitative discussion of this issue, which becomes particularly relevant for dynamically-driven quenches.
To address this point, we firstly compare our numerical and experimental condensate growth curves for dif-

FIG. 3.
Condensate Growth Dynamics. a-d Comparison of dynamical simulations (filled symbols, solid lines) and experimental data (open black circles, dashed lines) for different cooling rates, labelled here by the rate of change of temperature with time (dT /dt). Each subplot shows a complete averaged experimental sequence of condensate atom number N0 for a particular quench rate, with the depicted averaged numerical results corresponding to quench rates either side of the experimental values. In the simulations we assume a fixed initial (T = 790nK > TC , N = 22 × 10 6 atoms) and final equilibrium (T = 210nK < TC , N = 6.6 × 10 6 atoms) states. We also fit the numerical data with an S-shaped function (see Eq. (8) in Methods) to extract the growth timescale τG. e Dependence of τG on the ramp duration τR in the simulations, with error bars representing the 95% confidence bounds of the fits. f Condensation onset time, t bec , as a function of τR, demonstrating that slower quenches feature a delayed onset. Error bars here represent the difference in time for the condensate atom number to reach 3% (lower bound) or 7% of the total final atom number.
ferent characteristic temperature quench rates dT /dt, as shown in Fig. 3. Consistent with earlier numerical studies of condensate growth (including Kibble-Zurek dynamics [15,45]), the constant γ appearing in our theoretical model (see Eq. (2) in Methods) is treated as a free parameter, choosing its value such that it reproduces the experimental growth rate for dT /dt = 5.6 µK/s (Fig. 3b). In doing so, our numerical results are shown to accurately reproduce experimental growth curves for the entire range of quenches probed.
The precise determination of the phase-transition crossing and associated critical temperature is a rather challenging problem, both numerically and experimentally, particularly in the presence of inhomogeneous confinement. Experiments typically identify the critical re-gion as the time at which a clearly detectable condensate emerges. Following such a protocol, in our numerical simulations we identify a corresponding condensation "onset time", t bec during the early condensate growth stages, as the time at which the condensate atom number, N 0 ,defined as the largest eigenvalue of the one body density matrix -reaches 5% of the total final atom number for our chosen final equilibrium parameters; such definition is broadly consistent with our experimental measurements. Shifting our numerical growth curves by the "onset" time t bec enables a direct comparison to experimental growth curves, with the good agreement shown in Fig. 3a-d.
All our numerical condensate growth curves are well fitted by an S-shaped curve with a single free parameter, corresponding to the quenched growth timescale, τ G (Fits shown in Supplementary Fig. 4). The dependence of τ G and condensate onset time, t bec on ramp duration is shown in Fig. 3e,f. For τ R 300 ms, we find a practically constant growth timescale τ G = 15.8±0.7 ms, which is consistent with the finding of overlapping condensate growth curves once these are plotted against shifted time (t − t bec ). Note that the condensate onset time t bec is a linear-like monotonically-increasing function of τ R . This behaviour is qualitatively robust to changes in the exact definition of the condensate number/fraction chosen to mark such transition, as demonstrated by the small error bars in Fig. 3f.

D. Defect Number Dynamics and Visualization
Next, we discuss the number (Fig. 4) and nature ( Fig. 5) of the emerging defects restricting our analysis to times t ≥ t bec , with a qualitative analysis at earlier times prohibited by the densely-packed defect configurations present at early post-quench evolution times.
Defect Generation and Evolution: The number of defects already present at t bec varies significantly with ramp duration, being much higher for faster ramps, as such ramps create a more non-equilibrium initial configuration for the system. Despite the difference in absolute numbers, all defect-number curves exhibit a similar dynamical behaviour, as evident in Fig. 4a. During the first ≈ 30 ms, such curves exhibit a rather rapid initial decrease in their numbers associated with defect interactions. This is followed by a period of slower decrease with just a few defects present in the system. Such defects interact only occasionally, and mostly in pairs, because the typical distance between them during their motion is relatively-large. However, during such interactions they can violently emit a significant amount of energy, and occasionally/eventually get converted to a sound wave, or ejected by the evolving condensate near the edges of the system. Consistent with previously-reported experimental findings [53,56], for which the defect number could only be reliably measured after about a hundred ms of intrap evolution, we find the decay in the number of defects to be significantly slower once only two defects are left in  the system in any given individual run; in that case, their dynamics becomes largely decoupled, dominated by free propagation in the inhomogeneous condensate. This latter evolution stage is consistent with exponential decay, on a timescale of one hundred to few hundred ms. Relation to Kibble-Zurek Mechanism: We next investigate (Fig. 4b) the mean number of defects (vortices) as a function of ramp duration (top axis) or, equivalently, inverse temperature evolution (bottom axis). According to the idealised Kibble-Zurek scenario for defect generation, one would expect a power-law decay with a particular exponent. However, such Kibble-Zurek predictions are specific for a relatively early evolution time within the critical region, and are based on a scenario of frozen dynamics which does not account for the possibility of defect interactions, whereas experiments typically count defects after a significant in-trap evolution time and subsequent expansion imaging. Our simulations have already shown (Fig. 4a) that the defect number significantly changes during the in-trap evolution. Nonetheless, the power-law scaling of the defect decay (i.e. the slope of the decaying region in Fig. 4b) appears not to be very sensitive to the exact post-quench timing of the measurement, thus demonstrating the consistency of our numerically-identified defect dynamics with the experimental findings [21]. Specifically, a detailed analysis of defect number vs. quench rate for two different evolution times (t − t bec ) ∼ 50 ms (red diamonds) and 200 ms (blue squares), reveals a power-law decay consistent with both the predicted range and the experimental measurements (conducted at 250 ms). Our simulations also recover the experimentally-observed plateau region for fast quenches [21], a behaviour which we find to set in already at early post-quench times. At early evolution times, such behaviour is attributed to a combination of the maximum defect counting resolution of a tangled configuration in a restricted (inhomogeneous) volume, and the quench rate exceeding internal system timescales, implying that the assumption of a well-defined local temperature starts to break down. The occurrence of a plateau for fast quenches has been also discussed in recent work in the context of (2 + 1)-dimensional holographic superfluids [58].
Defect Visualization and Nature: Typical experimental measurements made after an evolution time of about 250 ms and based on integration over different (radial or axial) directions after time-of-flight (TOF) expansion are shown in Fig. 5 (left) for cases corresponding to different defect numbers. Such experimental images are in direct qualitative agreement with our numerically-generated (in situ) images at the correspondingly-long evolution times (Fig. 5, right); the latter numerically-generated results additionally enable direct visualization of the condensate phase (rightmost column), which is crucial for probing the nature of the emerging defects. From Fig. 2 it is already evident that the initially-generated defects during the critical region are highly-excited, have random shapes, sizes and orientations. An interesting question is how such random defects evolve into the experimentallyobserved solitonic vortices [59]; these are vortices with squashed 2π phase winding [60] and vortex line lying in a transverse radial plane, which corresponds to the low- est energy configuration in a highly-elongated superfluid [59,[61][62][63]. Numerical analysis of the phase evolution and corresponding density plots, appear to indicate a gradual evolution from the random distribution of different excited vortical states [64], into excited defects which are preferentially stretched along the transverse directions. Such defects gradually relax to solitonic vortices as the system equilibrates. (Supplementary Fig. 5 and Movie 4). While it is impossible to define a precise time at which defects acquire a solitonic-vortex nature, due to the gradual nature of this process which also depends on quench rate, all our findings appear broadly consistent with a solitonic-vortex nature ∼100 ms after t bec . It is worth mentioning that a comparable timescale was found in [65] for the decay of phase-imprinted dark solitons into solitonic vortices in superfluid Fermi gases in a very similar elongated geometry.

E. Coherence and Equilibration Dynamics
Having demonstrated solid agreement with experimental observations in appropriate regimes, and the ability to further interpret those through our simulations, we now use our numerical scheme to provide a deeper insight into the complicated nonlinear dynamical evolution and equilibration of quenched systems, covering also regimes where no experimental measurements are currently available.
In our quenched dynamical evolution of an initially equilibrium thermal gas, we have seen the quenched system falling out of equilibrium around the critical region, and identified a subsequent time, t bec , associated with the onset of condensation, in a manner which enables direct comparison to experimental measurements. Here we investigate the re-equilibration dynamics of such a sys-tem to a final state dictated uniquely by our final quench parameters (µ final , T final ). We show that this relaxation process depends on τ R in a nontrivial way: in particular, while the condensate-number growth dynamics depends solely on the growth timescale, τ G (which is itself a function of τ R , see Fig. 3e), the coherence growth, is additionally sensitive to details of the defect-filled state of the system following the quench. This points directly to the link between relaxation of quench-induced defects on the one hand and coherence growth and final system equilibration on the other.
From the condensate-number growth fits, we have identified two distinct dynamical regimes: for slow enough quenches (τ R 300 ms), the growth timescale is a monotonically-increasing function of the quench duration, whereas faster quenches (τ R 300 ms) all exhibit a similar number growth timescale (Fig. 3e). Nonetheless such rapid quenches lead to a notable increase in the number of spontaneously-generated defects, whose subsequent ("coarse-graining") dynamics is crucial for the evolution of coherence.
To study the growth of coherence, we follow the procedure of the Cambridge group [19] by numerically shifting the wavefunction by a fixed amount and autocorrelating this with the unshifted copy of itself. This method provides an estimate of the coherence length of the system, l coh . Due to geometrical considerations, we focus here on the axial coherence length, obtained by transversal integration (see Methods). In all cases we find that the integrated coherence length only starts increasing noticeably about 30 ms after t bec , consistent with the end of the previously identified rapid defect decay stage (Fig. 4a). The amount of vorticity present in the system sets a maximum limit to the dynamical system coherence length. This is to be expected, and has already been noted, for example, in 1D [39,66] and 2D [67]. Importantly, however, we see that the coherence length of the slowest quenches grows much more rapidly and saturates at higher values. This observation points to the importance of a system expelling practically all of its defects before it can acquire a coherence length comparable to the system size.
Starting from the same initial thermal condition, Fig. 6a shows how the dynamical correlation functions (solid lines) evolve in cases of fast (left), intermediate (middle) and slow (right) ramps and compared to the corresponding equilibrium functions (dashed lines). Our analysis here focuses on the re-equilibration dynamics, and so all dynamical data presented here correspond to times t > 0 (with µ(t) > 0) after the system has crossed the ideal-gas transition temperature, with times scaled to the quenched growth timescale τ G in order to suppress corresponding differences in condensate number growth dynamics. In the context of our adopted definition for the correlation function [19], this implies that the correlation function approaches a diagonal straight line as the coherence length approaches/exceeds the system size, which is always the case for the equilibrium systems considered here, due to their large atom numbers: this is decoupled from the fact that the equilibrium coherence length is increasing in absolute terms as the condensate size grows. Looking at times 0 < t < t bec we see that although the corresponding equilibrium correlation functions already exhibit near-perfect coherence across the probed central region, the corresponding dynamical behaviour at (t − t bec )/τ G ≈ −1.2 deviates noticeably: in the case of the faster ramp (left column), the dynamical correlation function has only mildly increased from its incoherent thermal state initial value, contrary to a much stronger increase for the slower (quasi-adiabatic) ramp, which nonetheless also exhibits a phase coherence length significantly smaller than the quantum-degenerate system size. Fig. 6b shows the dynamical coherence length as a function of the rescaled timescale (t − t bec )/τ G . For all ramps the coherence length increases with time, but it does so in a much slower manner for the fast quench. As a result, significant coherence has already built up for the slower ramp already at our identified condensation onset time, t bec , highlighting that such a system approaches the equilibrium state rapidly in the early stages after crossing the phase-transition, even before significant number growth takes place; this is in stark contrast to the faster ramp, which still exhibits hardly any coherence. Subsequent images highlight the emerging difference in the dynamical phase correlation function, g int 1 (d x ) [defined in Eq. (7)], in a most pronounced way, with the slow ramp becoming phase-coherent already for (t − t bec )/τ G 1.5, when the system is still only at the initial part of its number growth curve, as opposed to the faster ramp whose coherence length is smaller than the system size even at the much later time (t − t bec )/τ G = 12.7.
To quantify such re-equilibration dynamics further after the system has fallen out of equilibrium upon entering the critical phase-transition region, we introduce here the scaled deviation, δl coh , of the dynamical correlation length l dyn coh (t) from the corresponding equilibrium one, l equil coh (µ(t), T (t)), defined by δl coh (t) = l equil coh (µ(t), T (t)) − l dyn coh (t) l equil coh (µ(t), T (t)) . (1) Despite similar condensate growth rates (Fig. 6c) the corresponding coherence growth dynamics shown in Fig. 6d exhibit starkly distinct features, remaining strongly-dependent on the quench rate: in such relative timescales, slower ramps lead to much more rapid equilibration than the faster ramps; the latter are slowed down by the detrimental role of the defects persisting within the system. Contrary to this, slow ramps which perturb the system less lead to the emergence of a nearlydefect-free and therefore phase-coherent condensate already around t ≈ t bec . The inset in Fig. 6c highlights the rapid emergence of a single macroscopically-occupied mode for the fastest ramps (600 and 1440 ms), consistent with the rapid monotonic decrease of δl coh (t), indicating the rapid crossover to a phase-coherent condensate. The vertical error bars arising solely from our numerical averaging are significantly larger in the case of fast quenches, which is a measure of the deviation between different trajectories for the same ramp. This is easy to understand, since the more vortices there are in the system, the more likely their configuration is to be significantly different from shot to shot.

II. Discussion
Motivated by recent experiments with dilute ultracold atomic gases, we have investigated numerically the dynamics of an equilibrium thermal gas quenched over a finite timescale to deep in the phase-coherent condensate regime. Monitoring the entire evolution, we have presented an insightful graphical representation of the critical region dynamics, during which the dynamical system falls out of equilibrium with the corresponding-parameter equilibrium system, through the dynamical symmetrybreaking spontaneous emergence of defects. The emphasis of our analysis has been on the less-studied reequilibration dynamics, addressing the interplay between defect emergence and dynamical evolution, and growth of coherence.
Depending on the quench duration we have identified different emerging dynamical regimes: for fast quenches, we observed a saturation in the number of detectable defects, associated with detrimental defect interactions and the inherent difficulty in counting randomly-oriented defects within a very tight volume. Rescaling the condensate number growth for all quench rates by the characteristic growth timescale for each ramp, we demonstrated the decoupling of coherence and number growth dynamics arising from the detrimental effect of defect emergence and propagation on system coherence. Although the overall growth timescale might be longer for systems un- (1 − |dx|/L), for all times t t bec , with the correlation function for slower ramps approaching the corresponding equilibrium ones faster than for fast ramps. b Growth of the coherence length in time: systems with more than 2 vortices on average exhibit coherence smaller than the system size, indicating that for the fastest ramps the system is still in the phase-fluctuating regime. c Condensate growth curves collapsing onto a single curve when time from the condensate onset is scaled to the intrinsic system growth timescale τG. Inset shows corresponding ratio, rP O , of most-populated to next-most-populated eigenmode of the one-body density matrix. d Corresponding scaled deviation, δl coh (t), of dynamical to equilibrium coherence length [based on Eq. (1)]. At t = t bec , the slower ramps (600 ms, 1440 ms) are closer to equilibrium than the faster ones. Slower ramps, scaled to their respective τG, re-equilibrate faster to a phase-coherent BEC (for which δl coh = 0), with faster ramps only reaching full system coherence asymptotically. Statistical errors are shown by coloured bands in (a) and vertical error bars in (b) and (d). Horizontal error bars in (b)-(d) arise from the combination of the uncertainty in identifying t bec and the fitting error for τG. dergoing slower external cooling quenches, in such cases the dynamical system quickly re-approaches the corresponding equilibrium configuration, as the slow evolution enables it to mimic local equilibration during its growth, falling out of equilibrium only in a relatively small time window upon entering the region of critical fluctuations. In cases where the quench induces numerous defects, we have observed, as expected, enhanced shot-to-shot variability.
Our numerical analysis has also shed more light into the dynamical crossover from quasicondensation to true Bose-Einstein condensation, studied here in the context of an elongated geometry; such a geometry is known to lead to an enhanced decoupling of characteristic temperatures for the onset of density and phase fluctuations [55,68], with such decoupling effectively translating into different growth rates for the phase-coherent and density-coherent parts of the system, and so different growth rates for coherence and quasicondensation. Quasicondensation here refers to a defect-filled phase-fluctuating state with a coherence length smaller than the quantum-degenerate system size exhibiting suppressed density fluctuations and spanning many largelypopulated modes; this is directly contrasted to "true" condensation, which refers to a phase-coherent condensate with little, or no, defects and a single emerging macroscopically-occupied mode, a definition which holds even if the system is growing in time. Even though the quasicondensate stage in the critical region pre-empts all phase-coherent condensate growth due to the decoupling of the two corresponding characteristic temperatures, the dynamical quasicondensate regime is enhanced both in parameter space and in the temporal domain by the prolonged survival of the defects, thus being critically dependent on the quench rate responsible for their initial spontaneous generation. For slow quenches and large enough atom numbers considered here, the system may already be largely phase-coherent, i.e., a "true" condensate, as soon as it has grown to a size which makes it experimentally detectable. We note however that such findings are sensitive to the details of the system, and specifically here to both the system geometry and the chosen final state of the system, with reduced dimensionality enhancing such a distinction by enhancing the role of phase fluctuations [55,57,69].
In conclusion, we have presented a systematic numerical study of the re-equilibration dynamics of a quenched elongated 3D ultracold gas, demonstrating the decoupling of number and coherence growth in quenched quantum gases.
Such decoupling is a consequence of the emergence and dynamics of the spontaneouslygenerated defects as the system crosses the phase transition. The detailed numerical visualization of the defectfilled quenched phase-transition dynamics allowed access to a broad temporal range of dynamics not typically experimentally-accessible, demonstrating the dominant role of defect interactions and decay in the early stages of condensate formation. Moreover, our findings have been shown to be fully consistent with experimental measurements in the appropriate limits (and for the relevant quantities) where those exist.
Our work is expected to be of relevance to a broad range of future investigations with quantum-degenerate systems, and could also have technological implications for dynamical control and state-engineering of a quantum system. Given that our numerical scheme has demonstrated good qualitative description also of the muchharder-to-model approach to the phase transition, we believe that our method could in the future offer further insight into delicate features of non-equilibrium condensate dynamics, including a critical assessment and extension of the inhomogeneous Kibble-Zurek phenomenon.

III. Methods
Experiments. We produce ultracold samples of sodium atoms in the internal state |F, m F = |1, −1 in a cigarshaped harmonic magnetic trap with trap frequencies ω x /2π = 13 Hz and ω ⊥ /2π = 131.4 Hz. The thermal gas is cooled down via forced evaporative cooling and pure BECs of typically 10 7 atoms are produced. The part of the evaporation ramp in the vicinity of the transition is performed at different rates, from 50 kHz/s to 2 MHz/s. The quench ramp is followed by a variable wait time, during which a radio frequency shield is kept on to prevent from heating. After that, the atoms are released from the trap and are observed in two possible ways: Either we take simultaneous absorption images of the full atomic distribution along the radial and the axial directions [21,59] or we extract, uniformly, a small amount of atoms from the trapped sample and image it after a short time of flight [53,56]. The protocol is such that images are taken, and vortices are counted, after a fixed overall time interval from the BEC transition point, which is clearly identifiable for each quench ramp. We are also able to precisely identify the frequency ν of the RF field at T c , as well as to control the temperature variation in time, (∂T /∂t), via the speed of the evaporation ramp (∂ν/∂t) [21]. The defects that we observe at the time of imaging are quantized vortex lines which are seen as dark stripes when looking at the BEC from a radial direction after time-of-flight. The natural size of the defects in the trapped BEC, at the end of the cooling ramp, is of the order of the in-situ healing length ξ, which is of the order of 200 nm. After a long TOF, the defect size becomes larger than our imaging resolution of 3 µm. The presence of a levitating magnetic field gradient makes it possible to achieve long TOF preventing the BEC from falling. The measured vortex number is averaged over many experimental realizations in order to get good statistical samples for each experimental condition.
Numerical Model. Out study is performed by means of the (simple growth) stochastic projected Gross-Pitaevskii equation [49] (see also related model without projector [29,[50][51][52]), already demonstrated as a useful tool for the quenched crossing of the BEC phase transition [15,45,46]. In brief we simulate the low-lying highly-occupied modes of the system, denoted by the classical field (or c-field [49]) Ψ C (r, t) through the dynamical equation where P C is the projection operator truncating the modes above the c-field regime (so above an appropriately-identified energy cutoff), γ is a constant determining the condensate growth timescale, µ(t) is the time-dependent chemical potential, and . Fluctuations are included through the complex white noise, dW γ , defined by dW * (r, t)dW (r , t) = [2γk B T / ]δ C (r − r ) where δ C (r − r ) = n∈C ψ * n (r)ψ n (r ) in the chosen orthogonal basis set, {ψ n (r)}. The constant γ can be analytically approximated (at least for near-equilibrium cases) as [50,70] γ ≈ 2 m 2π 2 where g = 4π 2 a/m is the interaction strength, corresponding to the s-wave scattering length a. Although the above formula indicates the leading functional dependence of γ, the factor of "few" conceals within it the fact that this is an effective scaling, and so one should only rely on this for order-of-magnitude estimates. Consistent with earlier related analysis [15], here we treat γ as a fitting parameter. Comparing to experimental condensate growth data for dT /dt = 5.6 µK/s (Fig. 3b), we identify an "optimal" value γ = 0.005, which is about 10 times the analytically-predicted value for the initial temperature. To mimic the experimental cooling process, the chemical potential µ(t) and temperature T (t) are quenched linearly within a ramp duration τ R . The initial and final values of (µ, T ) are set as (−22 ω ⊥ , 790 nK ) and (22 ω ⊥ , 210 nK ), corresponding to a change of total equilibrium atom number (when also including above cut-off atoms under the usual assumption that they are static) from 22 × 10 6 to 6.6 × 10 6 with a 75% final condensate fraction. In our simulations, which start from a highly incoherent equilibrium state well above T c , the time t = 0 is chosen as the time when µ(t) = 0, since most interesting dynamics occurs after this time, implying that in any given simulation, the initial (equilibrium) configuration is at a time t = −τ R /2.
We solve the SPGPE with 4th-order Runge-Kutta in a plane-wave basis using a grid size L x = 54a ho,x along the x and L y = L z = 6a ho,x along the transverse directions, where a ho,x = /mω x ≈ 5.8µm is the characteristic harmonic oscillator length in the long direction (x-axis); we use a temporal discretization dt = 10 −3 /ω x and an energy cutoff fixed at 2.5 times the value of the final chemical potential (22 ω ⊥ ) in a grid consisting of N x = 1170 and N y = N z = 130 points. Simulations are run on Newcastle University's High-Performance-Computing cluster, Topsy, using 20 to 24 nodes. A single dynamical run takes between (120 − 300) CPU hours with an additional ≈ 40 CPU hours for the Penrose-Onsager diagonalization of the selected snapshots. We estimate the total amount of presented simulations to have taken over 10, 000 CPU hours.
Identification of the Condensate. The one-body density matrix is defined as where N sample is the sample number of the subensemble average and δt is set as ∆t/N sample with an appropriately short time-interval ∆t (so that the system dynamics is not masked). Such short-time averaging mimics the ensemble averaging based on the ergodicity hypothesis [49]. The notation ... used in this and the next subsections denotes the short-time subensemble average. In our simulations, N sample = 101, and the tests of probing ∆t provide a value of around 8 ms for our simulations, which is smaller than the characteristic time scale of the harmonic trap τ ho = 1/ω x ≈ 12.2 ms. The condensate, or Penrose-Onsager (PO) mode [54] at a given time t is identified as the eigenmode of the one-body density matrix ρ(r, r ; t) with the largest eigenvalue. To assess the degree of fragmentation of the condensate, in the sense of competition between different highly-occupied modes, we evaluate the ratio, r P O , of the largest to the second largest eigenvalues of ρ(r, r ; t).
Correlation Function Analysis. We follow the procedure of the Cambridge quenched-dynamics experiment [19], which measured the correlation function by interefering a displaced copy of the system with itself. Specifically we define the function where L is a chosen length. Following Ref. [19], the correlation length l coh is extracted from g H 1 by fitting it with (1−|d x |/L) exp(−|d x |/l coh ), in which the triangularshape function in the bracket arises from the integration. Here we probe the spatial coherence of our c-field wavefunction Ψ C within a region [−L/2, L/2] with L ∼ 54 µm, and numerically evaluate this through the phasephase correlation function by g 1 (d x , y, z; t) = where the operation H[f ] is defined as and φ(r, t) is the argument of Ψ C (r, t). To obtain transversal averaging we use the integrated version, where the prime integration is performed over the yellow region (i.e. density isosurface at values of 0.1% of the final c-field density.) To take finite-size effects into account we have also introduced into the above definition a density-dependent weighting function w P O (d x , y, z; t) which assigns higher weighting to the large-density regions. This is defined here through w P O (d x , y, z; t) = H[n P O (r)/n P O,peak ]/N , where N ensures the normalization condition˜ dydzw P O (d x , y, z; t) = 1 is satisfied at any given d x . (The second-order correlation function shown in Supplementary Fig. 5 and Movie 4 is defined as the onsite correlation function g 2 (r) ≡ g 2 (r, r, r, r) = |Ψ C (r)| 4 /[ |Ψ C (r)| 2 ] 2 .) Dynamical Timescales. Consistent with typical experimental measurements, in which T c is identified as the time of emergence of an observable condensate, we define here the "onset" or "delay" time for condensate growth, t bec , as the moment that the number of atoms in the condensate, N 0 , reaches 5% of the final total particle number (including in our considerations the particle number above the c-field region, which is assumed to be static). We have a posteriori verified this to provide (when used with our value of γ) an excellent description of condensate growth across all experimentally-probed regimes, and have also checked that the main findings presented in this paper are insensitive to the details of such definition. In addition to t bec , we also define the condensate growth timescale, τ G , which is extracted by fitting the condensate growth curve over the entire temporal range t with where N 0,i (N 0,f ) denote the initial (final) PO condensate atom numbers,t is the moment that N 0 (t ) reaches the mid-value (N 0,i + N 0,f )/2, and τ G is the single fitting parameter. We note here two things: firstly, that the mid-value is unique to all numerical growth curves, since we have a unique set of experimentally-relevant initial and final parameters in our simulations; moreover, we have checked that the extracted τ G values are largely insensitive to whether the fit is performed over the entire temporal range t, or whether it is constrained to values t ≥ t bec , suggesting the independence of the two timescales t bec and τ G .
Defect Identification. In our work we identify the location of vortices by the regions of high velocity, v(r) = ( /m)Im[Ψ * C (r)∇Ψ C (r)]/|Ψ C (r)| 2 , characterizing the region around the vortex core. By scanning the whole local maximum of the velocity field within the yellow region, we identify the positions of the vortex cores.
Statistical Analysis and Error Bars. For each numerically-simulated quench rate, we have analysed between N = 3 and 7 independent-noise realisations. The statistical uncertainty in the vortex number, N v , was estimated as In counting the vortex number, N v , we also estimated the possible systematic errors introduced by the use of subjective criteria in the identification of single vortices in situations where a vortex is at the boundary of the condensate or two vortex lines are very close to each other. However, we have checked that the corresponding uncertainty is significantly smaller than the statistical error defined above. Our procedure for assigning errors to the determination of the characteristic timescales is as follows: -the condensate onset time, t bec was defined as the time at which the condensate (Penrose-Onsager) atom number reaches 5% of the total final atom number. Error bars in our determination of t bec arise from shifting the (heuristic) value of 5% between 3% and 7%, values which are still consistent with the experimental growth curves reported in Fig. 3.
-the errors in the quenched growth timescales, τ G , arise from the 95% confidence bounds of the fit to our numerical growth curves with Eq. (8): the quality of the fits can be seen in Supplementary Fig. 4. Those two errors are treated as independent in the determination of temporal error bars for the scaled time (t − t bec )/τ G discussed in Fig. 6. Competing financial interests The authors declare no competing financial interests.

Materials & Correspondence
All correspondence and material requests should be addressed to N.P.P. In this Supplementary material we provide more detailed examples of the evolution presented in this manuscript.
• Supplementary Fig. 1 shows side-by-side the evolution of characteristic density isosurfaces of the highestpopulated (Penrose-Onsager) mode for a single run based on the same dynamical noise sequence and 3 different quench rates, characterised by their duration τ R . The more tangled evolution at early times, and the longer survival of defects in the fastest quench is visible.
• Supplementary Fig. 2 shows characteristic examples of evolution of the Penrose-Onsager mode for a given quench rate (τ R = 144 ms), demonstrating clearly the stochastic nature of the defect generation, and the stark differences between different numerical runs -a detailed analysis of which sets the error bars in Fig. 4 of the main paper.
• Supplementary Fig. 3 shows the evolution of the Penrose-Onsager mode in a focused spatial region and small time-intervals for a single-run corresponding to a fast quench (τ R = 18ms), clearly demonstrating the various types of processes that dominate the early physics of the system as it crosses the phase transition after which the initially-tangled-up generated defects interact and relax.
• Supplementary Fig. 4 demonstrates our fitting procedure for extracting the condensate quenched growth timescale τ G .
• Supplementary Fig. 5 shows different projections of the condensate wavefunction, phase and local secondorder correlation function, depicting clearly the defect relaxation to solitonic vortices at late evolution times.
We also provide here real-time evolution movies corresponding to the data shown in Supplementary Figures 1, 2, 3 and 5 (with more snapshots). Dependence of Condensate Growth and Equilibration on Quench Rate: Shown are single-realisation density profile evolutions for 3 typical ramps (from left to right: τR = 84, 144, 600 ms), with all subplots depicting density isosurfaces of the highest-occupied (Penrose-Onsager) modes at values 3% (green regions) and 0.1% (yellow regions). Each row corresponds to a different absolute time t from the moment the system is quenched, also depicting the corresponding shifted time (t − t bec ) which depends on τR (and so varies from column to column). All cases demonstrate an initial narrowing of the spatial distribution as the system approaches the critical region. Faster quenches are shown to have a more violent random defect tangle generation at t t bec thus exhibiting more defects at t bec [boxed plots] when the condensates have similar sizes. At late times, faster quenches exhibit more vortices, with the slowest quench considered here never acquiring a vortex within the high density region. Runs shown here are based on exactly the same dynamical noise sequence and are for illustrative purposes. Our estimated values for t bec for the three ramps are respectively (from left to right): t bec (τR = 84 ms) = 48 +4 −3 ms, t bec (τR = 144 ms) = 66 +6  (v) t=30 ms (t-t bec =−3 ms) (vi) t=33 ms (t-t bec =0 ms) (vii) t=36 ms (t-t bec =3 ms) (viii) t=42 ms (t-t bec =9 ms) (ix) t=45 ms (t-t bec =12 ms) (x) t=51 ms (t-t bec =18 ms) (xi) t=54 ms (t-t bec =21 ms) (xii) t=57 ms (t-t bec =26 ms) (xiii) t=60 ms (t-t bec =27 ms) (xiv) t=63 ms (t-t bec =30 ms) (xv) t=66 ms (t-t bec =33 ms) (xvi) t=69 ms (t-t bec =36 ms) (xvii) t=72 ms (t-t bec =39 ms) (xviii) t=75 ms (t-t bec =42 ms) (xix) t=84 ms (t-t bec =51 ms) (xx) t=90 ms (t-t bec =57 ms) (i) t=0 ms (t-t bec = −33 ms) (ii) t=15 ms (t-t bec = -18 ms) (iii) t=18 ms (t-t bec = -15 ms) (iv) t=24 ms (t-t bec = -9 ms) Emergence and Interaction of Vortices in a Growing (Quasi)condensate: Consecutive images depict gradual dynamical evolution of the density profiles of the highest-occupied (Penrose-Onsager) mode as the phase transition is crossed violently due to a rapid external quench (τR = 18 ms). Our focussed region, spanning only a subset of our numerical grid and of the final system size, depicts clearly both the initial defect tangle complexity (15 -18 4. Condensate Growth and Emerging Timescales: Typical condensate growth curves as a function of (absolute) time t, averaged over 3 realizations, shown from the moment that µ = 0, alongside corresponding fits (dashed cyan lines) with Eq. (8) from Methods for ramps with (from left to right) τR = 18, 144, 300, and 1440 ms. The errors in the Penrose-Onsager number are smaller than the size of the depicted points, reconfirming the consistency of the Penrose-Onsager analysis despite the different density profiles between runs with same τR. Corresponding extracted values of τG are 16 ± 2 , 14.9 ± 1.7,14.3 ± 1.0 and 78 ± 13 ms, with errors representing the 95% confidence level of the fits. For reference, our estimated values for t bec for the ramps depicted here are respectively (from left to right): t bec (τR = 18 ms) = 33 +3 −3 ms, t bec (τR = 144 ms) = 66 +6 −3 ms, t bec (τR = 300 ms) = 108 +7 −3 ms, and t bec (τR = 1440 ms) = 354 +21 −39 ms. (a) t=36 ms (t-t bec = -12 ms) (b) t=48 ms (t-t bec = 0 ms) (c) t=66 ms (t-t bec = 18 ms) (d) t=180 ms (t-t bec = 132 ms)

FIG. 5.
Crossover from Random Defects to Solitonic Vortices: Shown are different visualizations of density and phase of the quenched (τR = 84 ms) mode corresponding to the largest eigenvalue at four characteristic evolution times 36, 48, 60, and 180 ms, corresponding to t − t bec = −12, 0, 18 and 132 ms. In each panel, the top subpanels illustrate the highestmode (Penrose-Onsager) density isosurfaces from different viewing angles set by (i) z = 0 and (ii) y = 0, with axial view images shown in the square panels in (iii). Subsequent panels (iv) and (v) depict the column-integrated densities, respectively integrated over the z and y axes, and the phase profiles in the (vi) z = 0 and (vii) y = 0 planes. The two bottom-left panels show the corresponding local second order correlation function, g2(r, r, r, r), in the (viii) z = 0 and (ix) y = 0 planes. Transversally-integrated density of the Penrose-Onsager mode ρint(x) =˜dydznPO(r) is shown in (x), with the sliced density ρ sliced (x) = nPO(x, y = 0, z = 0, t) shown in (xi), where nPO the Penrose-Onsager density. This figure demonstrates clearly the initial emergence of random-type defects (which are most clearly pronounced in the density cuts), the transition of the phase profile from random to solitonic-vortex-like, and the survival of only a handful of well-defined solitonic vortices at late evolution times, propagating on top of an otherwise largely relaxed Thomas-Fermi-like condensate profile. These images also highlight the difficulty in counting vortices at early times from any chosen integration angle. We note that t bec (τR = 84 ms) = 48 +3 −3 ms. (A more coarse-grained dynamical evolution of projected densities is shown in Supplementary Movie 4.)