Kibble-Zurek Dynamics in a Trapped Ultracold Bose Gas

The dynamical evolution of an inhomogeneous ultracold atomic gas quenched at different controllable rates through the Bose-Einstein condensation phase transition is studied numerically in the premise of a recent experiment in an anisotropic harmonic trap. Our findings based on the stochastic (projected) Gross-Pitaevskii equation are shown to be consistent at early times with the predictions of the homogeneous Kibble-Zurek mechanism. This is demonstrated by collapsing the early dynamical evolution of densities, spectral functions and correlation lengths for different quench rates, based on an appropriate characterization of the distance to criticality felt by the quenched system. The subsequent long-time evolution, beyond the identified dynamical critical region, is also investigated by looking at the behaviour of the density wavefront evolution and the corresponding phase ordering dynamics.

In this paper we focus on the evolution of a trapped ultracold atomic gas across the transition to a Bose-Einstein condensate. We perform a detailed numerical analysis of externally-driven spontaneous symmetry breaking and dynamical growth of an elongated, harmonically confined, three-dimensional (3D) condensate by solving the stochastic (projected) Gross-Pitaevskii equation (SPGPE) in realistic experimental parameter regimes, previously identified in our quantitative analysis of the late-time relaxation dynamics probed experimentally [19,43,44,49].
Our numerical results are interpreted in terms of the homogeneous KZ mechanism by comparing the solutions of the full 3D stochastic nonlinear equation against analytical predictions of the linearized limit of the same equation. At short times from the transition, where the system is close to criticality, we find excellent agreement with the KZ scaling laws predicted by the linearized theory, with our numerical curves for different quench timescales appropriately collapsing onto a unified curve. In particular, the growth of the condensate is delayed with respect to the critical point by a delay time proportional to the KZ timescale. Remarkably, we also find that the KZ delay persists at later times as long as the system is ramped linearly in time. Specifically, density growth is found to occur along elliptically-expanding regions of phase-space which mimic the underlying trap geometry, with the rescaled expanding wavefronts collapsing to a single (non-universal) curve for different quench rates. Schematic of the homogeneous KZ mechanism, marking the interplay between the diverging system relaxation time τ (black dashed line) and the time |t| = | /˙ | to the transition (solid blue lines). The intersection points of these two curves mark the crossover times −t and +t.

arXiv:2004.09642v1 [cond-mat.quant-gas] 20 Apr 2020
Despite the inhomogeneous nature of the harmonically trapped gas, our present work seems to indicate that the temperature quenches probed in the experiments [19,43,44,49] were such that transition effectively occurs within the remit of the 'homogeneous' KZ mechanism; the predicted modifications due to the interplay of causality and geometry [11,85,86] seem not to be needed in this case.

II. Quenched Protocol and Modeling
A. Temperature quench and KZ mechanism The gas is initially prepared in a thermal state above the critical temperature T c and then ramped across the phase transition, where a symmetry breaking occurs and an order parameter appears. During such evolution, the effective distance from the critical point can be measured by a dimensionless parameter = 1 − T /T c . Close to the critical point this parameter can be assumed to be linear in time, as where τ Q is referred to as the quench time. While the system approaches T c from above, but still far from it, the evolution is adiabatic, i.e., the gas follows its adiabatic thermal equilibrium state. Such adiabaticity fails at a characteristic time, −t, when the relaxation time becomes longer than the instantaneous timescale | /˙ | = |t| at which is ramped. The relaxation time diverges as | | −zν , where ν and z are the equilibrium (correlation length) and dynamical critical exponents, respectively. From the equation |t| |t/τ Q | −zν one obtainŝ which corresponds to a deviation from criticalitŷ In the 'cartoon' version of the homogeneous KZ mechanism (see Fig. 1), during the system evolution, started at large negative initial value of , the state of the gas freezes-out at −t and subsequently remains unchanged until a time +t, when the adiabatic evolution starts again. During that period, the correlation length ξ is frozen at the valueξ of the equilibrium correlation length at −ˆ , given byξ The above scenario (adiabatic-impulse-adiabatic approximation) is of course a simplification of the actual dynamics, as physical quantities still evolve during the time the system spends in the critical (impulse) region, as qualitatively demonstrated e.g. in Ref. [19], and characterized in detail in Ref. [87]. However, the notable importance of the simplistic KZ mechanism, which also explains its broad applicability to a range of different physical systems, is that it correctly predicts the scaling of the characteristic lengthscaleξ and the timescalet with the quench time τ Q . It is noteworthy that the two scales are related byt They both diverge in the adiabatic limit, τ Q → ∞, where they become the unique relevant scales in the KZ scaling ansatz [88][89][90]. For instance, a two-point correlation function C R (t), between two points separated by a distance R, should satisfŷ Here d is the number of dimensions, η is a universal critical exponent, and G is a non-universal scaling function. Equation (6) is expected to be accurate in the long-wavelength and low-frequency limit. The adiabaticimpulse-adiabatic approximation is consistent with the scaling hypothesis (6) but it implies a particular form of the scaling function G that does not depend on t/ξ z during the freeze-out between −t and +t.
The stochastic projected Gross-Pitaevskii equation models the dynamics of the low-lying highly-populated 'classical field' ψ, through [91] where is the Gross-Pitaevskii term including single-particle evolution and mean field potential (nonlinearity), and dW denotes complex Gaussian white noise with a correlator The projection operator P in Eq. (7) restricts the dynamics below the energy cutoff E cut , which is fixed here as 2.5µ f , where µ f is the chemical potential of the gas at the end of the quench. The trapping potential The interaction strength g = 4π 2 a s /M is set by the s-wave scattering length a s . The dimensionless parameter γ controls the rate of relaxation of the classical field modes to the equilibrium state set by the chemical potential µ and temperature T of the reservoir of atoms located in thermal bath above the cutoff. The detailed numerical simulations performed in this work are based on earlier stochastic dynamics simulated by some of us in different geometries, dimensionalities, platforms and systems [19-21, 101, 106-108] (and related work [96-98, 102, 109]).

C. Parameter choice and quench protocol
Our study is performed for the parameters corresponding to a recent experiment [19,49], performed with few×10 7 23 Na atoms in the |F, m F = |1, −1 state (with a s = 2.91 nm), trapped in an anisotropic harmonic potential, , with longitudinal and transversal trap frequencies ω x = 2π × 13 Hz, ω ⊥ = 2π × 131.4 Hz, yielding a highly elongated 3D system. In the experiment, after creating a thermal cloud above the critical temperature, evaporative cooling is used to ramp the temperature down to much below T c , where the system exhibits significant condensation. The experiment performed a detailed study of the late-time evolution of vortex defects originally generated during the symmetry-breaking phase transition [49], finding a power-law decay within the range expected by the KZ mechanism. SPGPE simulations [19] were in good agreement with observations in the late-time regime where experimental data were available; however, experimental limitations could not facilitate such quantitative analysis of the system dynamics at earlier times. Here we focus on the early-time regime of the condensate formation. Our starting point is a better numerical estimate of the distance to criticality which enables us to cast the dynamics in the standard language of the KZ mechanism and characterize it in terms oft.
In our SPGPE simulations, after initially equilibrating the system via Eq. (7) to the desired initial thermal state defined by its chemical potential µ and temperature T , we linearly vary T and µ over a timescale 2τ Q , with the ramp initiated at t = −τ Q and finished at t = τ Q , based on the imposed quench protocol (for |t| ≤ τ Q ) After the end of this linear ramp at t = τ Q , the 'input' parameters T and µ remain fixed at their final values. Initial and final values are chosen to match typical numbers observed in Ref. [19,49], giving T 0 = 500 nK, ∆T = 290 nK, and µ f = 22 ω ⊥ , corresponding to initial and final atom numbers N i = 22 × 10 6 and N f = 6.6 × 10 6 . Equation (7) is solved in a ≈ 314.1 × 34.9 2 µm 3 cuboid box with 1130 × 130 2 grids (with grid size being ≈ 0.27 µm in all directions).

D. Equilibrium phase diagram
We first characterize the precise location of the equilibrium critical point ( = 0). For this purpose, we numerically calculate the equilibrium configuration of the gas for a given set of T and µ, in order to construct the corresponding equilibrium phase diagram of condensate fraction vs. temperature. This is done in the context of the self-consistent Hartree-Fock approximation [111,112]. The total atom number N tot = N c + N I includes atoms in the c-field, N c =´dr |ψ(r)| 2 , and in the thermal bath, N I . The atoms above E cut are assumed to be in a thermal reservoir at the given T and µ, with N I =´∞ Ecut d F I , where and (r, k) = For a given N tot , one can estimate the transition temperature for the corresponding ideal Bose gas, given by T c,0 ≈ 0.94 ωN . This can only serve as a useful reference value for an actual interacting system, due to the competition between thermal fluctuations and interparticle interactions [113], and the relevance of finite size effects [114,115].
Following a standard procedure [91], we can calculate the condensate atom number, N 0 , from the classical-field wavefunction ψ by means of the Penrose-Onsager criterion, through identification of the largest eigenvalue, and the corresponding eigenfunction ψ 0 (r), of the singleparticle density matrix ρ(r, r ) = ψ * (r)ψ(r ) , where · · · denotes a short-time average over 100 samples [116]. The resulting condensate fraction N 0 /N tot is plotted in Fig. 2(a) together with experimental data of the Trento group [110]. Our equilibrium simulations, in agreement with the experiment, clearly reveal that the condensate fraction arises at T /T c,0 ∼ 0.9 rather than 1, with the corresponding critical chemical potential µ shifting to a positive value due to the finite-size and interaction effects [114,117,118]. In order to better identify the critical region, in addition to the condensate fraction we calculate three further quantities exhibiting critical behaviour in a narrow temperature region [103,104,106,115,[119][120][121], namely the correlation length, l coh , the Binder cumulant, C b , and the order parameter, m. The longitudinal correlation length can be extracted by an appropriate fit to the first-order   [110] (purple diamonds) for the equilibrium condensate fraction N0/N vs. Tc,0. The dashed orange curve shows the ideal Bose gas prediction as a reference. The vertical yellow band marks the location of the numerically-identified critical region in SPGPE simulations, Tc ∼ 445.5 ± 7.3 nK; in this range, the total particle number is (5.6 ± 0.1) × 10 6 , which corresponds to the ideal gas critical temperature Tc,0 = (488 ± 5) nK. Background colour indicates the value of the chemical potential at each temperature during a quench, with µ(t) and T (t) proceeding from the rightmost to the leftmost point. (b) Filled blue squares: longitudinal correlation length l coh , extracted as in Eq. (13), during a SPGPE simulation of a quench with τQ = 150ms. Time is measured from the equilibrium critical time tc, i.e., the centre of the yellow vertical band (the same as in panel (a)). Open black circles: same quantity calculated in SPGPE simulation for equilibrium states with input values µ(t − tc) and T (t − tc). During the quench the growth of l coh is delayed with respect to the instantaneous equilibrium: we find such delay to correspond to (t − tc) ∼ 1.3t (dotted vertical cyan line), wheret (solid vertical cyan line) is the timescale predicted by the KZ mechanism. The scaled deviation, δl coh , between dynamical and equilibrium correlation lengths, defined by Eq. (15), is shown by the purple squares and exhibits a very rapid increase in the critical region, followed by a slower decay during the re-equilibration process, which reflects the phase ordering process. The end of the ramp is denoted by the vertical black dashed line. (c) Corresponding characteristic single-trajectory evolution of the Penrose-Onsager condensate density profiles (τQ = 150ms): plotted yellow and green regions respectively map out the density isosurfaces for 0.1% and 3% of the peak value of the final post-quench equilibrium condensate density; purple filaments denote region of high velocity field, corresponding to the location of spontaneously-generated vortices.
correlation function, via [19] G (1) (d x ) =¨dy dzˆL Here L x ≈ 54.4 µm denotes a central portion of the axial extent of the inhomogeneous system over which the correlation function is evaluated (for comparison, the final equilibrium condensate spans the range ≈ [−114, 114]µm). The weighting function w PO is introduced here to reduce the contribution of low-density regions in the transverse direction. Details of our procedure to calculate l coh , C b , and m are given in Appendix A, with different extraction protocols showing excellent agreement between them. As a result, we identify the equilibrium transition temperature in the range T c ∼ 445.5 ± 7.3 nK, roughly corresponding to The corresponding critical chemical potential is µ c = (4.13 ± 0.55) ω ⊥ > 0. Based on such values we can extract, for each quench τ Q , the dynamical critical time, t c , during the quench when µ(t) and T (t) reach their corresponding critical values µ c and T c . Specifically, we find We can identify the time t c as the reference time from where to measure the distance , which enables us to cast all dynamical behaviour and relation to the KZ mechanism in terms of the shifted time (t − t c ) from the equilibrium phase transition.

E. Quenched dynamics
The quenched dynamical growth of the system can be visualized by means of a particular simulation example. Building on our earlier work which focused on the (late-time) re-equilibration dynamics of a quenched Bose gas [19], Fig. 2(b)-(c) shows the evolution of the correlation length and density profiles for the particular case of τ Q = 150 ms.
Examining the evolution of the correlation length during a quench as a function of (t − t c ), shown in Fig. 2(b), we notice that the growth of the dynamical correlation length l coh (t) starts, as expected, at a later time to that of the corresponding equilibrium correlation length l equil coh , evaluated at equilibrium with the same µ(t) and T (t). In accordance with KZ mechanism, our simulations indicate a delay proportional tot.
To complement our findings, Fig. 2(c) shows corresponding 3D single-trajectory density profile snapshots during the quenched evolution. Shortly aftert, the system remains dominated by fluctuations as shown in Fig. 2(c)(i). The first evidence of condensation onset appears around 1.3t, in the form of a localized elongated higher density condensate region containing multiple spontaneously-generated defects (purple filaments) as in Fig. 2 are dominated by the interplay between condensate growth (driven by the increasing µ(t) > µ c and decreasing T (t) < T c ) and phase-ordering through defect relaxation, which was previously shown to lead to a decoupling of number and coherence growth [19]. Fig. 2(c)(v) shows a typical long-term profile after both density and coherence have saturated to their equilibrium values, which for the particular example displayed here occurs after the end of the external driving.
The evolution of coherence during the quench can also be visualized through the 'auxiliary' variable [19] δl where l equil coh (t) is the equilibrium correlation length at time t. Early on in the quench, during the adiabatic regime, the dynamical correlation length closely follows the corresponding equilibrium one, until the system enters the critical region and δl coh (t) quickly increases from 0 to 1. The value of δl coh (t) remains ≈ 1 until (t − t c ) ∼ 1.3t, after which time it clearly starts decreasing, but at a much slower rate that its initial increase: the latter decay, previously characterized in Ref. [19] is evidence of defect relaxation and phase ordering, until reaching values δl coh (t) ∼ 0, at which late time the dynamical system has grown sufficiently to become practically indistinguishable from the corresponding equilibrium one.

III. Linearized SPGPE
In the symmetric phase before the phase transition, when µ(t) < µ c , there are small thermal fluctuations around the symmetric vacuum ψ = 0. In the noninteracting limit, µ c = 0. We expect that during the non-equilibrium linear quench these fluctuations remain small until some time after the critical point. Therefore, the out of equilibrium evolution near the critical point can be reasonably described by a linearized version of Eq. (7) where the interaction term g|ψ| 2 in Eq. (8) is neglected. Furthermore, as the initial growth occurs around the centre of the trap, we assume here for simplicity that V trap (r) ≈ 0.
With these two approximations (and omitting the projector in our analytical considerations) In this framework the small fluctuations become dynamically unstable towards exponential growth of a condensate when µ(t) crosses 0 towards positive values. This can be understood by noticing that the dissipative terms on the right hand side of Eq. (7), that are proportional to γ, include a minus gradient of a Mexican-hat-like potential, When µ(t) > 0 the symmetry is broken and the potential has instantaneous minima at a ring with the dissipation driving ψ towards this instantaneous vacuum manifold.
Before proceeding with such analytical treatment below, which will guide our subsequent numerical analysis of the full nonlinear dynamics, we make two important comments. Firstly, the presented analytical discussion implies that the critical point arises exactly at t = 0. However, the experimentally relevant equilibrium phase diagram of Fig. 2(a) has already revealed a shift in time, which we will subsequently account for by replacing t by (t − t c ). Secondly, the linearized discussion neglects the role of the nonlinearity g|ψ| 2 ψ up tot, at which point it will be argued to slow down the exponential blow-up of |ψ|. However, its effect is not completely negligible even beforet in the simulation. As shown in Sec. II D, in equilibrium there is a range of values 0 < µ < µ c where the symmetry breaking Mexican hat potential is too shallow to prevent restoring the symmetry by thermal fluctuations. In this way, the actual symmetry breaking transition is shifted from the simplistic approximation of µ = 0 to the more appropriate µ = µ c > 0.
With those 'caveats' in mind, we proceed next with our analytical predictions, initially conducted for a homogeneous system, and subsequently generalized to modes beyond k = 0.

A. Uniform field
Let us first consider a uniform field ψ(t) when Eq. (16) becomes Here dW is also assumed uniform. When µ < 0 then ψ = 0 is stable and its instantaneous relaxation time is On the other hand, when µ > 0 then ψ = 0 is unstable with a Lyapunov time that is also given by formula (20). In general τ is a time-scale on which the system can adjust to the time-dependent µ(t). The time-scale diverges at the critical point µ = 0. Near the critical point the system is too slow to adjust, no matter how long τ Q is, and its state is effectively frozen between the two crossover times, ∓t, when the reaction time of the system equals the time to the transition (see Fig. 1): Solution of this equation with respect tot yields the crossover time:t This is the KZ time-scale. Near −t the uniform ψ goes out of equilibrium with the instantaneous µ(t), hence its fluctuations do not diverge at µ = 0 -as might be suggested by Eq. (19) -but remain small in consistency with the linearized approximation. The linearization remains self-consistent until near +t when ψ begins to catch up with the varying µ(t) again and the dynamical instability begins to blow up exponentially.

B. Reciprocal space
In order to go beyond the uniform case, k = 0, we consider a (modified) Fourier transform, with an extra dynamical phase pre-factor included. In the reciprocal space the linearized SPGPE (16) becomes a Wiener-like stochastic equatioṅ whereζ(t, k) is a Gaussian white noise with a correlator For k = 0 equation (24) becomes the uniform Eq. (19). When µ > 0 then all modesψ(t, k) with 2 k 2 2M < µ are dynamically unstable. At +t, when the dynamical instability begins to blow up, all modes with k up tô are already unstable. This borderlinek is a solution of 2k2 2M = µ(t). They are amplified by the dynamical instability and dominate the power spectrum near and after +t. An inverse ofk, is the KZ correlation length. The power spectrum is dominated by modes with wave lengths longer thanξ. The power laws (22), (26), (27) are consistent with the general KZ predictions (2)-(4) involving the critical exponents z and ν. Indeed, Eq. (20) implies that τ is proportional to an inverse of the distance from the critical point, here measured by |µ|, hence zν = 1. At the critical point, µ = 0, the right hand side of Eq. (24) implies relaxation with a rate ∝ k 2 , hence z = 2. Therefore, the general KZ formulas (2)-(4) predict in agreement with Eqs. (22) and (27), respectively.

C. KZ scaling hypothesis
For large τ Q , the length-scaleξ and the time-scalet become longer than any other scales and, therefore, they become the only relevant scales in the low frequency and long wave length regime. Therefore, according to the KZ scaling hypothesis, in this regime physical observables depend on time t, distance r, and wave vector k through scaled variables t/t and r/ξ, andξk, respectively. Here we verify the hypothesis for the linearized Eq. (24).
A formal solution of stochastic Eq. (24) is An equal-time correlator of these Gaussian fluctuations follows from Eq. (25) as where the spectral function is Here t k = τ Q ( 2 k 2 /2M )/µ f is the time when the Fourier modeψ(t, k) becomes dynamically unstable. The spectral function depends on t and k through a single variable This demonstrates not only the anticipated KZ scaling in the form but an even stronger relation Here F and F are non-universal scaling functions.

D. Near +t
The spectral function (31) is monotonically increasing with u. Consequently, for any time it peaks at k = 0 and for any k it is increasing with time. The peak value begins to blow up like e (t/t) 2 near t/t ≈ 1: The blow-up enhances the peak of the spectral function in a neighborhood of k = 0 where u becomes large enough for the integral in (31) to be approximated by √ π: In its regime of validityξ 2 k 2 t/t, hence it can be further simplified to a Gaussian: This Gaussian neglects fluctuations with wave lengths shorter than ξ = 2ξ(t/t) 1/2 that have not been enhanced by the blow-up yet. The Gaussian spectral function translates to a coarse-grained equal-time correlation function As anticipated, neart its range ξ becomes the KZ correlation lengthξ. It is noteworthy that for any t/t this correlation function is proportional tot/ξ 3 ∝ξ −1 which is consistent with the general scaling hypothesis (6) given that d = 3 and, in our linearized Gaussian theory, η = 0.
Setting r = r we obtain average strength of the coarsegrained fluctuations: accurate neart or later. These are also times when (39) blows up and the linearized SPGPE begins to break down. This suggests a scaling behaviour that Further growth is halted by the interaction term in the Mexican hat potential (17) that was neglected in the linearized equation. The nonlinear interaction begins to be felt already at the inflection point of the potential: Therefore, equating |ψ(t, r)| 2 to (41) is a good indicator when the linearized approximation breaks down. Thanks to the exponential nature of the blow-up (39) the breakdown time is close tot up to logarithmic corrections. It is noteworthy that, at t ≈t, the KZ correlation length equals the healing length in the instantaneous Mexican hat potential. The healing length is a width of a vortex core, hence it is not possible to stabilize a tangle of vortex lines whose separations are less than the healing length. This justifies a posteriori our coarse-graining over wave lengths shorter than the KZ coherence length. The shorter fluctuations are not relevant for formation of stable vortex lines.

E. Beyond +t
According to the linearized theory, neart the magnitude |ψ| 2 should come close to the inflection point of the Mexican hat potential. Near the inflection the potential is approximately a linear function, hence its gradient is a constant and the magnitude |ψ| should grow linearly in time. This is a significant slow-down after the initial exponential blow-up. Nevertheless, eventually |ψ| grows enough to get close to the instantaneous equilibrium magnitude (18) at the bottom of the potential (true vacuum). This equilibrium depends on time through t/τ Q , rather than t/t characteristic for the early times before and aroundt, because it follows the linear ramp that depends on t/τ Q .
However, as the equilibrium magnitude depends on time, the equilibration cannot be perfect and |ψ| 2 must be delayed with respect to the instantaneous equilibrium (18). The delay time should be proportional to a relaxation time towards the bottom of the Mexican hat potential at the moment when the magnitude's growth slows down near its inflection point. This relaxation time is proportional to the universal KZ timescale,t. Therefore, we expect that the instantaneous equilibrium (18) should be replaced by a crude formula: where α is a non-universal constant, excepted to be ∼ O(1). This is approximately valid long aftert, when the KZ scaling hypothesis no longer applies, but there is still a delay proportional to the KZ delay timet. It is worth emphasizing that even after the nearequilibration of the magnitude, the phase of ψ should remain random with a characteristic KZ coherence length ξ. The phase is the Goldstone mode for this symmetry breaking, hence it is not subject to the aforementioned relaxation. It is only in the subsequent evolution that the phase undergoes slow phase ordering kinetics [122] that proceeds by gradual annihilation of the randomlygenerated vortex networks. In this sense the KZ coherence lengthξ is a more robust imprint of the KZ physics that survives to very late times.

F. Shift of the critical point
In the proceeding discussion the nonlinearity g|ψ| 2 ψ was neglected up tot where it was argued to slow down the exponential blow-up of |ψ|. However, its effect is not completely negligible even beforet in the simulation. As shown in Sec. II D, in equilibrium there is a range of µ > 0 up to µ c ≈ 4.13 ω ⊥ where the symmetry breaking Mexican hat potential is too shallow to prevent restoring the symmetry by thermal fluctuations. In this way the actual symmetry breaking transition is shifted to µ = µ c . In addition to the shift, the equilibrium universality class is also altered with the mean-field correlation length exponent ν = 1/2 replaced by the exact ν = 0.67. Correspondingly, given the dynamical exponent z = 2, the predictedt ∝ τ zν/(1+zν) Q should be altered fromt ∝ τ 0.50 Q tot ∝ τ 0.57 Q . In the following we assume validity of the physical picture developed within the Gaussian theory but incorporate the criticality shift from t = 0 to t c , into predictions of the linearized S(P)GPE by making a replacement t → t − t c . Regarding the scaling oft with τ Q , we note that due to the statistical uncertainties, it is not possible to discriminate between the similar power laws:t ∝ τ 0.50 Q (our Gaussian approximation) and the improved scalinĝ t ∝ τ 0.57 Q , a point further discussed in Appendix B.

G. Homogeneous assumption
For the harmonically trapped system considered in this work, the instability addressed above is considered to be occurring in the volume where µ(t) − V trap (r) > 0. Due to the anisotropy, this volume is enclosed in an ellipsoid with λ ⊥ = ω ⊥ /ω x . Such an ellipsoid defines a critical volume of the system, V c ≡ 4πa x (t)a 2 ⊥ (t)/3, and expands along its principle semi-axes with velocities  3. Types of quench classified by comparing the critical wavefront a x,⊥ andâ x,⊥ for times up tot. When the quench duration is comparable witht, it becomes quasiinstantaneous, with momenta beyond 1/ξ excited. Increasing τQ leads to a broad homogeneous quench regime (blue), which fully encompasses all experimentally-relevant quenches probed in this work, whose range is marked by the vertical arrow. Slower quenches can lead to regimes where the quench is inhomogeneous only in the transverse direction (green), or in both directions (red), but both of these would require ramp durations exceeding 10s for the current trap. Note that, since the critical wavefront starts from the trap centre at tc, the quench is always homogeneous at t−tc = 0. The black dashed lines are boundaries estimated by (47) to (50).
These velocities diverge in the centre of the trap where the instability appears first at t = 0.
An investigation of possible corrections due to the system inhomogeneity requires us to compare these velocities of the critical front (v x , v ⊥ ) with the perturbation velocityv within the critical regime given bŷ The quench is effectively homogeneous when the critical front velocities (v x , v ⊥ ) are larger thanv. As the critical front velocities diverge in the center of the trap, the quench is effectively homogeneous there. In the longitudinal (transverse) direction the quench remains homogeneous until the moment when v x =v (v ⊥ =v). The latter equations can be solved with respect to a x (a ⊥ ) to respectively definê Inside the ellipsoid with semi axisâ x andâ ⊥ , where both v x and v ⊥ are faster thanv, the system is effectively homogeneous, and we find that these two quantities are respectively larger than a x and a ⊥ .
The ratio ofâ x(⊥) /a x(⊥) (t =t) for a given τ Q provides a guidance of the homogeneity of a quench. When a x(⊥) /a x(⊥) (t =t) > 1, the quench is homogeneous in the x (transverse) direction. We can compare the instability front a x(⊥) (t) andâ x(⊥) up to t =t. The conditions for a quench to be longitudinally/transversally homogeneous thus are and Identification of different criteria for a homogeneous quench across the longitudinal and transverse directions gives rise to a rich diagram of possible behaviour, based on our quenched input parameter µ(t). The types of quenches possible for the considered trapping potential, characterized in terms of their (in)homogeneity up tot are summarized in Fig. 3. In the blue region,â x(⊥) > a x(⊥) and the quenches are effectively homogeneous. As τ Q is increased to cover the green regime, the quenches are effectively longitudinally homogeneous whenâ x > a x but transversally inhomogeneous asâ ⊥ < a ⊥ . In the red region, the quenches are effectively inhomogeneous, sincê a x,⊥ < a x,⊥ . We also note here that, shortly after the transition, the quenches are all effectively homogeneous, as the critical wavefront starts growing outwards from the trap centre at t = t c .
The experimentally relevant quench parameters investigated in our present study lie well within the homogeneous quench regime, and hence the above linearized SPGPE analysis is expected to be applicable. When the quench duration becomes comparable tot, the quenches can be regarded as quasi-instaneous quench (yellow region): in such cases, the presence of the nonlinearity allows for momenta beyond 1/ξ to be excited after the termination of the fast ramp. To see the inhomogeneous effects in a quench, one could consider much slower quenches or increase the trapping frequencies and aspect ratio, tuning the boundaries according to (47) to (50).

IV. Early-time KZ scaling and SPGPE
Having demonstrated the relevance of the homogeneous KZ mechanism for the parameter regime considered in this work, we now examine the extent to which the linearized SPGPE -supplemented with the time shift (t − t c ) -can accurately explain the results of the full nonlinear SPGPE numerical simulations.
Firstly, we consider the spectral function, defined in Section III C. The time evolution of the spectral function can be extracted from the full SPGPE simulations. Figure 4 shows the evolution of the peak value f (t, 0) as a function of t − t c . The same curves are plotted in panel (b), but rescaled according to the analytic scaling law predicted by the linearized theory, Eq. (33). We see that the curves corresponding to different quench rates collapse onto each other in the approximate range (t − t c ) 2t, thus demonstrating the validity of the KZ scaling hypothesis. Furthermore, the collapsed curves blow-up near the scaled time (t − t c )/t = 1 as predicted by Eq. (35). At later times, as the field fluctuations ap-  2t. Grey lines in right column represent the parabola u = (t − tc)/t −ξ 2 k 2 = 1, along which the spectral function is predicted to be constant according to the stronger relation of Eq. (34), with t → t − tc, within their regime of validity. The horizontal pink dashed line marks the value of 1/ξ, which corresponds to the largest wave number being excited up to t − tc =t as discussed in Sec. III B.
proach the inflection point of the Mexican hat potential, the slope of the curves decreases as an effect of the nonlinear interaction term in the SPGPE not included in the linearized theory.
The evolution of the whole spectral function f (t, k) is investigated in Fig. 5, which shows raw (left column) and scaled (right column) numerical data for different values of τ Q as a function of the shifted time (t − t c ). The raw data demonstrate a strikingly different behaviour for different quench times. Nonetheless, plotting the same data scaled according to the law (33) reveal great similarity, particularly for the three slowest quenches (bottom three panels). In other words, for long enough τ Q , the scaled spectral functions collapse to a common scaling function thus confirming the KZ scaling hypothesis. In the same panels we also attempt a test of the scaling law in its stronger form of Eq. (34). The added grey parabolas satisfy u = (t − t c )/t −ξ 2 k 2 = 1. According to Eq. (34) the spectral functions should be constant along these lines, which appears to be the case here, up to statistical fluctuations.
Nonlinear effects becomes relevant when the density is large. In Fig. 6 we plot the time evolution of the number of atoms in the classical field near the center of the trap, defined as where, upon accounting for the system anisotropy, V cen has been chosen as the ellipsoid around the center within half harmonic lengths in all directions. The results for different ramps are shown in panel (a) as a function of t − t c , while in panel (b) the same curves are plotted according to the scaling law (40) predicted by the linearized theory. Again, the curves nicely collapse onto each other in the same early-time regime (t − t c ) 2t, where the spectral function also collapses, while for larger times, the scaling is less effective.
Finally, in Fig. 7, we show that the longitudinal correlation length growth defined by Eq. (13) also collapses onto a single scaling function by applying the same KZ rescaling previously used for the spectral function. All quenches exhibit an initial growth at some delayed time after (t − t c ) = 0, with faster quenches displaying a faster initial growth. In the case of the slower quenches, l coh grows smoothly to the final value of around 430 µm. However, the three fastest quenches (72, 150 and 200 ms) reveal evident fluctuations in the value of l coh during its growth. These have been previously identified [19] as being due to the persistence/dynamics of defects (vortices) within the region |x| ≤ 27.2 µm over which this correlation function is evaluated. This is more pronounced for the very fast quench [ τ Q = 72 ms (red)], for which the Evolution of the longitudinal correlation length defined by (13) and calculated in the SPGPE simulations for different quench times: (a) raw data as a function of (t − tc); (b) same curves rescaled according to the KZ scaling predicted by the linearized theory. The light blue vertical solid and dotted lines in the inset respectively mark the positions oft and 1.3t.
cooling ramp terminates at t − t c ≈ 2t (dashed black vertical lines in Fig. 5(a)(i)-(ii)), thereafter exciting higher momentum modes. Nonetheless, within (t − t c ) 2t, and after rescaling, l coh reveals excellent collapse for all curves, as evident from Fig. 7(b).

V. Late time dynamics
Up to now, we have accounted for the early-time dynamical phase transition crossing within SPGPE, interpreting the result in the context of the homogeneous KZ mechanism and the linearized theory. In this section we examine the extent to which the late-time dynamics of the nonlinear SPGPE -based on our quench protocol of fixed initial and final states, and different quench duration τ Q -are also collapsible onto a single curve in a way dictated by the KZ mechanism.
Firstly, we examine the late-time evolution of the central particle number, defined by Eq. (51), in Fig. 8(a). While the raw density evolution curves corresponding to different quench rates differ widely [see earlier Fig. (6)(a)], when plotted as a function of (t−t c −αt)/τ Q as suggested by Eq. (42), instead of (t − t c )/t, the different curves collapse nicely for the non-universal constant α = 1.3. There is only one notable, but not unexpected, outlier: the fastest quench with τ Q = 72ms. This quench hast close to the end of the linear ramp, hence its late time evolution is largely after the end of the linear ramp. A similar behaviour is also found for the evolution of the position of the density wave fronts in Fig. 8(b), determined by tracing a near constant value of the classicalfield density |ψ(r)| 2 , arbitrarily chosen here in the range [16,20] µm −3 to ensure relatively smooth curves (more details on this are given in Appendix C). Again, all data nicely collapse onto a single curve along both directions once the ellipsoidal growth mimicking the underlying anisotropic harmonic confinement is accounted for, consistent with the arguments exposed earlier in Sec. III G, namely that the growth always occurs along an ellipsoid. The slightly different dynamical behaviour of the wave front for the 72 ms ramp at intermediate times can be understood by the fact that this particular fast ramp terminates at (t − t c − 1.3t)/τ Q ∼ 0.3, as marked by the hollow red circle.

VI. Discussion and conclusions
We performed a detailed analysis of the early-stage quenched symmetry-breaking dynamics of an elongated harmonically trapped three-dimensional ultracold atomic gas evaporatively cooled from above the Bose-Einstein condensation phase transition temperature at variable rates. Our study was conducted by means of the stochastic projected Gross-Pitaevskii equation for parameters corresponding to a recent experiment, and cast in the language of the Kibble-Zurek mechanism.
Schematically, as the quenched system approaches the critical point from above, it enters a regime where it cannot follow the adiabatic evolution of the equilibrium state, due to the quench proceeding faster than the characteristic diverging relaxation time of the corresponding equilibrium system. Adiabaticity is resumed at a certain time around +t (actually we find a short delay prefactor of ∼ O(1) compared to the standard Kibble Zureck prediction) and the overall process leads to the spontaneous emergence of defects (in this case vortices), with some of those gradually becoming embedded in the growing condensate. Although the system is still evolving during its quenched evolution within the critical region -rather than remaining frozen in the 'impulse' limit of the 'cartoon' KZ version -such evolution still exhibits scaling properties predicted by the KZ mechanism.
In order to properly characterize the scaling laws for the observables in our SPGPE, we needed to first extract the equilibrium critical temperature of the interacting system numerically. Identification of the equilibrium critical point is crucial to correctly apply the KZ model to a shifted evolution time after the time t c when the system crosses the corresponding equilibrium critical point.
Then we used the analytical predictions based on the linearized form of the stochastic Gross-Pitaevskii equation and KZ ordering considerations. Such predictions were found to be valid, allowing quantities like spectral functions, correlation lengths and density growth to collapse onto unique curves for all different quench times probed here and performed in the experiments motivating this work.
At later times, the growth of the k = 0 mode and of the (ellipsoidal) density wave front proceed on similar timescales. However, the presence of highly-excited kmodes associated with the existence of defects in the growing condensate -which are more pronounced for the fastest quenches -implies that the phase-ordering process and coherence growth depend on τ Q and system geometry/inhomogeneity in a more complicated manner. This highlights the important nature of the decoupling of density and coherence degrees of freedom [19]. Although phase ordering for homogeneous systems is an established topic with known scaling laws, the presence of inhomogeneity and anisotropy introduces finite-size effects from the early stages of the evolution, making a collapse of the late-time dynamics particularly tricky even in numerical simulations.
As all our results at early times were found to be consistent with the homogeneous KZ symmetry-breaking transition, this raises the question of whether current ultracold experiments can actually observe the inhomogeneous Kibble-Zurek phenomenon, which requires the emerging driven system evolution to proceed slower than the speed of propagation of the causality front. Firstly, we note that it is very hard to extract reliable observables at early evolution times in a highly-controlled setting in the context of current harmonically trapped ultracold atom experiments. Aside from this, for the particular experiment on whose parameter this study is based this would require quench durations of many seconds, even if tighter harmonic confinement was used. This is an interesting topic for future investigations.
The numerical identification of the critical point from SPGPE equilibrium simulations is performed by using three quantities: the Binder cumulant, the correlation length and the order parameter m.
Two closely-related definitions of the Binder cumulant appear in the literature [115,119,120,[123][124][125][126]. In the first definition, appropriate for a homogeneous system, it is defined in terms of the full classical-field ψ via [125,126] The second definition has been implemented in the context of the trapped Bose gas, and extracts a similar in- indicates the value of the thermal de Broglie wavelength λ dB whose value in the critical region is ∼ 0.55µm.
(A2) In both cases, one expects a sharp jump from the value 1, below T c , to the value 2, above. The critical value of the Binder cumulant at the transition in the thermodynamic limit is C ∞ b, critical ∼ 1.2430 [119], while for trapped Bose gases it is expected to be smaller than C ∞ b, critical and affected by finite-size effects [115]. Fig. 9(a) shows our numerical results based on both definitions, with their results convincingly overlapping with each other. The jump from 1 to 2 is clearly visible and the critical value C ∞ b, critical in found in the range T /T c,0 ∈ (0.91, 0.93), corresponding to T ∈ (438, 453) nK.
In the critical region, the correlation length is also expected to diverge as |1 − (T /T c )| −ν [115,127]. Based on our chosen extraction method for the correlation length, l coh , defined by Eq. (13), we indeed find l coh starts increasing rapidly in the above probed region, as evident from Fig. 9(b). However, the inhomogeneous finite-size nature of the system, our chosen definition of an integrated coherence length, and our numerical accuracy do not allow for the identification of a sharp critical point, and thus cannot facilitate an accurate determination of the static critical exponent ν.
Finally, we also follow Refs. [125,126] and investigate the behaviour of the order parameter m defined within our computational volume V by This quantity, plotted in Fig. 9(c), is expected to be m ∼ 0 above the phase transition and m = 1 for a pure condensate at T = 0 [125,126]. Again we see that m starts increasing within the same critical region of the Binder cumulant and the correlation length. The verification that both l coh and m start increasing within the critical regime identified by the Binder cumulant, and the fact that this also coincides with the region when the condensate fraction decreases to zero, provide strong evidence for the consistency of the identification of our critical regime.
The vertical yellow area in Fig. 9 highlights the above determined range T /T c,0 ∈ (0.91, 0.93). In our system, the corresponding critical chemical potential is µ c = (4.13 ± 0.55) ω ⊥ and the critical time t c in our quench protocol is B. Mean field vs. exact critical behaviour In Section III we introduced the linearized SPGPE approach with the Gaussian approximation. Within this approach, the equilibrium critical exponent should coincide with the mean-field exponent ν = 1/2, leading to the scaling lawt ∝ τ A natural question arises about whether the scalings/collapses presented in the main text would be significantly affected upon using instead the exact value of with τ 0 = /(γµ f ) and ξ 0 = (2M µ f / 2 ) −1/2 . The differences between the above expressions and the corresponding mean-field values oft andξ are shown in Fig. 10(a). They are not significant in the considered regime of parameters. In panel (b) of the same figure, we present the correlation length and its rescaling witht andξ . Similarly to Fig. 7, except for the fastest quench with τ Q = 72ms, the curves for different quench times τ Q collapse onto a single curve for short times. Overall, this analysis suggests that our dynamical results cannot accurately distinguish between the two values of ν. The KZ exponents they imply do not differ enough.

C. Density distributions and spectral function
In the main text, we have identified the shifted time (t − t c − 1.3t) as the time when the system exits the KZ self-similar regime where evolutions corresponding to different quench times τ Q collapse onto one another. At later times, the dominant timescale governing the system is the quench time, which determines the rate at which  (29) for the same cases. In both columns, the dashed red line is the numericallytraced position of the density wave front for the density within [16,20] µm −3 , while the pink line is spectral function f of the k = 0 mode. Panels (c) and (d) show the profiles of the density distribution and the spectral function, respectively, at the time (t − tc − 1.3t)/τQ ≈ 0.9, revealing the extent of excitation still present for the faster quenches.
the system is ramped to its low-T state (with the exception of the fastest ramp, τ Q = 72 ms, which ends while the system is still well within the self-similar regime). As a complement to our previous analysis, here we further investigate the behavior of the density and the spectral function of the gas during a quench. Fig. 11(a)-(b) compares the evolution of density wavefront and k = 0 mode alongside the evolution of the full density distribution and spectral function. Specifically, Fig. 11(a) shows how the axial system density grows as a function of time. The front of the growing density area has been traced for the lowest value of the density which allows a relatively smooth curve, and is shown by the dashed red line. The corresponding wave fronts for different τ Q have already been discussed in Fig. 8(b), where they were shown not only to collapse on top of each other, but also on top of the corresponding transversal evolution wave front, once the system geometry/anisotropy were appropriately accounted for. In Fig. 11(c)-(d) we also plot the long term evolution of the the density distribution and the spectral function at time (t − t c − 1.3t)/τ Q ≈ 0. 9 We can thus draw various conclusions already briefly commented upon in the main text: (i) From Fig. 11(a)-(b), we see that the growth of the density wavefronts (dashed red lines) overlaps almost perfectly with that of the k = 0 modes (solid pink lines). This suggests that growth on this timescale is driven by the k = 0 mode, consistent with bosonic amplification. However, (ii) although density and k = 0 mode grow in parallel, the spectrum of higher excited modes looks very different on such scaled time [ Fig. 11(b)]. For comparison, the instantaneous spectral function at time (t − t c − 1.3t)/τ Q ≈ 0.9 (when densities have mostly saturated) is plotted in Fig. 11(d). We thus see that although low momentum modes are mostly excited for slow quenches, whose late-time momentum distribution is consistent with the Bose-Einstein distribution, faster quenches generate more modes with higher k, with the highest excited modes subsequently relaxing only gradually, and on a much longer timescale. Importantly, faster ramps are still in a far-from-equilibrium state at (t − t c − 1.3t)/τ Q ≈ 0.9, even though both density wavefronts and k = 0 mode occupations are close to saturating at such a time. This offers a clear perspective of the previously inferred decoupling between density and momentum/coherence relaxation.