Ultra-quantum turbulence in a quenched homogeneous Bose gas

Using the classical field method, we study numerically the characteristics and decay of the turbulent tangle of superfluid vortices which is created in the evolution of a Bose gas from highly nonequilibrium initial conditions. By analysing the vortex line density, the energy spectrum and the velocity correlation function, we determine that the turbulence resulting from this effective thermal quench lacks the coherent structures and the Kolmogorov scaling; these properties are typical of both ordinary classical fluids and of superfluid helium when driven by grids or propellers. Instead, thermal quench turbulence has properties akin to a random flow, more similar to another turbulent regime called ultra-quantum turbulence which has been observed in superfluid helium.


I. INTRODUCTION
The formation of a coherent Bose-Einstein condensate from a thermal Bose gas is a rich topic of ongoing research [1]. Recent experiments on thermally-quenched Bose gases have observed the spontaneous formation of defects in the guise of vortices [2,3] and solitonic vortices [4], confirming the occurrence of the Kibble-Zurek mechanism [5,6] in these gases. A paradigm for this non-equilibrium phase transition is the formation of a homogeneous weakly-interacting Bose gas starting from highly non-equilibrium initial conditions [7][8][9][10][11][12][13][14][15]. The gas is modelled as a classical matter field described by the Gross-Pitaevskii equation. This field undergoes a universal self-ordering into a quasi-condensate, i.e. a coherent superfluid component, and a non-condensed, thermal component. Phase dislocations that become entrapped within the quasi-condensate during this quench give rise to an irregular tangle of quantized vortex lines that permeate the system, as seen in Fig. 1. This corresponds to a state of superfluid turbulence. Over time, this turbulent state relaxes, tending towards the vortexfree, partially-condensed equilibrium state of the finitetemperature Bose gas.
The aim of our work is to explore the nature of this state of superfluid turbulence, including the characteristics of its decay. Is it similar to turbulence in ordinary (classical) fluids? The question is natural, because tangles of quantized vortices created in superfluid helium by moving grids or propellers appear to obey classical scaling laws [16][17][18].
By examining the distribution of kinetic energy over the length scales, the velocity correlation function and the temporal decay of the vortex line density, we show that the turbulence resulting from the thermal quench of a Bose gas is very different from classical turbulence and any quasi-classical regimes of superfluid turbulence. Instead, it shares key properties with a second state of * george.stagg@ncl.ac.uk turbulence [19,20] which has also been observed in superfluid helium, called 'ultra-quantum turbulence', which has unusual non-classical properties.

II. CLASSICAL FIELD METHOD
We model a weakly-interacting, homogeneous Bose gas (including thermal excitations) by means of the classical field method [7][8][9][10][11][12][13][14][15]. The gas is parametrized by a classical matter field ψ(r, t), normalized to the total number of particles N = |ψ| 2 dV , and whose evolution follows the Gross-Pitaevskii equation (GPE) The cubic term accounts for the repulsive interactions between the particles, with g = 4πh 2 a s /m, where m is the particle mass and a s is the inter-particle s-wave scattering length (here a s > 0). The total energy of the gas is The GPE is conventionally used to model a zero temperature condensate, but it is now established that, provided the modes of the gas are highly occupied, the gas evolves as an ensemble of modes, each of which follows (to leading order) the classical trajectory described by the GPE [21][22][23]. Various phenomena have been studied within this classical field formalism, including equilibration dynamics [12,13,15,24], critical temperatures [25], correlation functions [26], vortex nucleation [27][28][29][30] and decay of vortex rings [31], as well as extensions to binary condensates [32][33][34]. A classical field ψ(r, t) representing N particles is simulated in a cubic periodic box of volume D 3 . We make the GPE dimensionless using the natural units of the homogeneous system: particle number density is expressed in terms of the average value ρ = |ψ| 2 = N/D 3 , length in terms of the healing length ξ =h/ √ mgρ, speed in arXiv:1607.03719v1 [cond-mat.quant-gas] 13 Jul 2016 terms of the speed of sound c = ρg/m, energy in terms of the chemical potential µ = ρg, and time in terms of τ =h/gρ. The time evolution is computed using a fourth-order Runge-Kutta scheme with time step ∆t = 0.01τ on a 192 3 grid with isotropic grid spacing d = 0.75ξ.
To ensure that all numerically accessible modes are highly occupied, the initial condition is the highly nonequilibrium state, where the 192 3 modes of the system are labelled by the wavevector k, the coefficients a k are uniform and the phases are distributed randomly [13]. The occupation of mode k is n k = |a k | 2 . The spatial grid discretization implies that high momenta are not described; in effect an ultraviolet cutoff is introduced, n k (t) = 0 for k > k max , where k = |k| and the maximum described wave vector amplitude is k max = √ 3π/d [37]. The condensate fraction ρ 0 /ρ of the equilibrium state (where ρ 0 is the density of the quasi-condensate, to be defined in Section III A) is uniquely determined by the number density, N/D 3 , and the energy density, H /D 3 . In practice these values are controlled by a rescaling of the initial condition ψ(r, 0) so as to modify the uniform value of n k in the initial condition, with a selection of high momenta coefficients a k set to zero in order to retain the same number density after rescaling. The process can be seen for the initial condition in Figure 2.
Our results focus on the three cases: low condensate fraction (ρ 0 /ρ = 0.22), moderate condensate fraction (ρ 0 /ρ = 0.48) and high condensate fraction (ρ 0 /ρ = 0.77). The parameters used to generate these cases are listed in Table I, along with the corresponding temperature of each state, T /T c , where T c is the critical temperature for the condensation, estimated using the empirical formula provided in Ref. [31]. . Also shown is the temperature of the gas, T /Tc, estimated using the empirical formula in Ref. [31].

III. FORMATION OF THE TURBULENT VORTEX TANGLE
A. The quasi-condensate The evolution of the system in momentum space [13,34] is shown in in Fig. 2 (upper). At t = 0 (red line), the mode occupation numbers n k are distributed uniformly with k, with values of the order of unity, up to the cutoff, as per the imposed initial condition. Over time, self-ordering of the classical field leads to a smoothing of the mode occupation distribution towards a characteristic bimodal form. At low k the modes have macroscopic occupation (n k 1), characteristic of Bose-Einstein condensation and of the presence of the quasi-condensate. At high k the modes have low occupation, and are the thermal excitations. Their distribution approximately follows the Rayleigh-Jeans equilibrium distribution for non-interacting particles, n k ∼ k −2 [36]. The transition in k-space between these two components is most evident when viewing the the integral distribution of the particles over the wavenumbers, F k = k <k n k (inset of Fig. 2); here a prominent 'shoulder' in the curve marks the transition from the quasi-condensate to the thermal modes. We define this transition point between the quasi-condensate and thermal gas as k c , and focus on the nominal value k c1 = 10 (2π/D). However, to assess the effect of the cutoff on our findings we also consider a second cutoff, k c2 = 20 (2π/D).
Due to the randomized phases in the initial condition, the system is full of phase defects, and a dense vortex tangle forms in the quasi-condensate, as seen in Fig. 1 (a). The raw wavefunction ψ is too noisy to directly visualise the superfluid vortex tangle. Following [13], this problem is overcome by defining a quasi-condensate wavefunctionψ viaâ k = a k × max{1 − k 2 /k 2 c , 0}. This procedure, which filters high-frequency modes from ψ, is analogous to spatial course-grained averaging, so that ψ contains only the long-wavelength component of the classical field ψ. The quasi-condensate density is then |ψ| 2 and its value is ρ 0 = |ψ| 2 ; the condensate fraction follows as ρ 0 /ρ, i.e. the ratio of the quasi-condensate density to the total density. The evolution of the condensate fraction is shown in Fig. 2 (lower). The condensate fraction ρ 0 /ρ is approximately zero at t = 0. It immediately undergoes a rapid growth as the quasi-condensate forms. The growth then slows and asymptotes towards its equilibrium value. For our primary choice of cutoff, k c = k c1 the condensate fraction reaches 95% of its asymptotic value within a time t = 600τ ; in comparison, the decay of the vortex tangle towards the vortex-free state typically occurs on a timescale of t = 7000τ . For the second choice of cutoff, k c = k c2 , the growth of the condensate fraction is slightly faster and the final condensate fraction is slightly larger, but the qualitative behaviour is unchanged. We have verified that while our simulations are sensitive to the initial condition parameters shown in Table I, the behaviours demonstrated in Fig. 2 are consistent over different randomized initial conditions. The vortices are visualized, such as in Fig. 1, through isosurface plots of the quasi-condensate density. The iso-surface density level varies with time according to |ψ| 2 = 0.05 |ψ| 2 . This level is chosen so that the resulting vortex structures are depicted with approximately constant cross-sectional radius throughout the decay, and that this is similar to the radius of a vortex core at zero temperature.

B. Energy spectrum
One method to characterise the turbulence of the vortex tangle formed in the quasi-condensate is by studying how the kinetic energy of the vortices is distributed over wavenumber. There are three wavenumber scales of particular interest in this problem: k D = 2π/D (associated with the length scale of the computational box), k ξ = 2π/ξ (associated with the superfluid healing length), and k (associated with the typical inter-vortex spacing, ). The typical inter-vortex spacing can be estimated as ≈ 1/ √ L, where L is the vortex line-density (length of vortex line per unit volume). These scales provide the perspective required to interpret the distribution of energy.
The kinetic energy density of the quasi-condensate is defined as E kin =h 2 2m |∇ψ| 2 . Using Parseval's theorem, the corresponding kinetic energy spectrum,Ê kin (k), can be written in terms of the angle-average of F( √ E kin ) 2 [38], where F denotes the Fourier transform, so that, The kinetic energy can be further decomposed into com- pressible and incompressible parts, E kin = E i kin + E c kin , where the compressible part is associated with sound waves and the incompressible part with vortices. This decomposition can be defined through the hydrodynamic interpretation of the superfluid as √ nv = ( √ nv) c + ( √ nv) i with ∇ · ( √ nv) i = 0 and ∇ × ( √ nv) c = 0, where n = |ψ 2 | is the particle density and v is the superfluid velocity. The incompressible kinetic energy spectrum is then denotedÊ i kin (k), defined in a similar way toÊ kin (k) in Equation (4). Figure 3 shows the incompressible energy spectrum of the vortex tangle at time t/τ = 300 (after the quasicondensate and tangle have formed but before the turbulence has undergone any significant decay). This spectrum clearly shows that our turbulence is unlike turbulence in ordinary fluids. In classical turbulence, the energy is concentrated at the smallest wavenumber k D , decreasing as k −5/3 (the celebrated Kolmogorov scaling) in the region k D k k . Figure 3 shows no such pile-up of energy near k D and no Kolmogorov scaling; on the contrary, the energy peaks at length scales just above k . The visible k −3 dependence of the spectra in the region k < k < k ξ arises from the vortex cores [38]. Note that these results are insensitive to the cutoff; we see the same behaviour for our second choice of cutoff k c = k c2 (inset of Fig. 3).

IV. RELAXATION OF THE TURBULENT VORTEX TANGLE
During the decay of the superfluid turbulence, shown in Fig. 1, the vortex tangle remains random and isotropic throughout its decay to one or more vortex lines or rings, which eventually also decay, leading to a vortexfree state. Insight into the nature of the turbulent decay is obtained by monitoring the vortex line-density L, defined as the vortex length per unit volume; the vortex length is estimated as V t /A where V t is the total volume within isosurface vortex tubes, and A the circular crosssectional area of a vortex, which is constant for a steady condensate fraction. Figure 4 shows the decay of the vortex line density over time. At very early times (t < ∼ 300τ ) the tangle is in the process of forming, while at very late times (t > ∼ 3000τ ) only one or two vortices or vortex rings remain. However, in the large intervening range of time a sizeable tangle exists. During this range the decay of the vortex line density has a power-law form, L = αt −βt with α constant and β ≈ 1 (or slightly less depending on the condensate fraction). For comparison, the decay of classical Kolmogorov turbulence is significantly quicker, L ∼ t −3/2 [39][40][41], and hence Fig. 4 is consistent with Fig. 3 and the absence of a Kolmogorov spectrum in our system.
Instead, the exponent β ≈ 1 which we observe (see Table II), is consistent with the 'ultra-quantum' turbulent regime revealed by the experiments of Walmsley and Golov [19], who created the turbulence by injecting vortex rings in a sample of superfluid helium initially at rest. In another set of experiments, Walmsley and Golov found that a longer, more intense injection stage generates 'quasi-classical' turbulence which decays as L ∼ t −3/2 . Following Volovik [42], Walmsley and Golov argued that whereas quasi-classical turbulence implies the existence of an energy cascade of eddies or vortex bundles [43] similar to ordinary turbulence, ultra-quantum turbulence lacks coherent structures and is more akin to a random flow. This interpretation was confirmed by numerical calculations of Baggaley et al. [44] who simulated the experiments, reproducing both ultra-quantum (L ∼ t −1 ) and quasi-classical (L ∼ t −3/2 ) regimes, and verified the presence of a Kolmogorov energy spectrum (Ê kin (k) ∼ k −5/3 ) only in the latter.
From the fitting parameter α, assuming that the energy dissipation rate per unit mass has the classical form dE/dt = −νω 2 , where the energy per unit mass, E, is proportional to the length of vortex line, and the vorticity is estimated as ω ≈ κL, we infer that the turbulent kinematic viscosity ν is such that 0.06 < ν/κ < 0.3, in fair agreement with numerical (ν/κ = 0.06 and 0.1 [44,45]) and experimental (ν/κ = 0.1 [46]) values for superfluid helium, although our values do not capture the temperature trend in helium.
To confirm the ultra-quantum interpretation of the turbulence created by a thermal quench we compute the velocity correlation function, f (r), defined (at fixed time t) as where v is the velocity of the quasi-condensate and the ensemble average is performed over positions r. Since the flow is isotropic, we only present results for displacements along one direction, here chosen to be the x-direction.
The function f (r) is dimensionless and normalised to f (0) = 1. If the vortex lines are essentially randomly oriented then at a distance r ≈ /2 the velocity correlation should vanish. Indeed, the inset of Fig. 5 shows that at time t/τ = 300 (when we monitor the energy spectrum, and ≈ 30ξ) the correlation function has become negligible at such distances. In fluid dynamics, a convenient measure of the distance over which velocities are correlated. is the integral length scale [47], defined (at time t) as Figure 5 shows I as a function of t during the evolution of the vortex tangle. Clearly I(0) ≈ 0 at very early times, reflecting the random nature of the initial condition. At later times, I(t) increases for all condensate fractions: as the tangle decays, fewer and fewer vortices are left, and the velocity field becomes correlated over larger regions of space. The fact that I(t) remains less than at all confirms the disorganized nature of our turbulence. We find that, as with the incompressible kinetic energy and decay of vortex line-density, the qualitative behaviour of f (r) and I is largely unchanged by the choice of cutoff k c . Filtering less modes only has the effect of introducing disorder into the quasi-condensate wavefunction, slightly reducing the quantitative values of f (r) and I.

V. CONCLUSIONS
Using classical field simulations, we have modelled the evolution of a finite-temperature homogeneous Bose gas from a nonequilibrium initial condition, through the formation of a turbulent tangle of vortex lines to the relaxation to a vortex-free state. By monitoring the vortex line density, the energy spectrum and the velocity correlation function, we have determined that the superfluid turbulence created by this thermal quench process lacks the coherent structures and the Kolmogorov cascade process which are typical of ordinary (classical turbulence), and which have also been observed in superfluid helium when driven by propellers or towed grids, or flowing at high velocity along channels. Instead, thermally quenched turbulence has properties similar to another regime of superfluid turbulence called ultra-quantum turbulence, observed in superfluid helium under certain driving conditions, which is less organized and more similar to a random flow.
While our work is based the idealized paradigm of an infinite homogeneous system, such a system can be approximated experimentally through the recent advent of quasi-homogeneous Bose-Einstein condensates formed in box-liked traps [3,49]. Nonetheless, condensates are most commonly confined in harmonic potentials, and it would be interesting in future to explore how the inhomogeneity affects the nature of the turbulence induced by the thermal quench.
Data supporting this publication is openly available under an Open Data Commons Open Database License [50].