Regimes of Steady-State Turbulence in a Quantum Fluid
Tommy Z. Fischer Ashton S. Bradley



Department of Physics, University of Otago, Dunedin, New Zealand
Dodd-Walls Centre for Photonic and Quantum Technologies

Department of Physics, University of Otago, Dunedin, New Zealand 
Dodd-Walls Centre for Photonic and Quantum Technologies


§ ABSTRACT 
We simulate the Gross-Pitaevskii equation to model the development of turbulence in a quantum fluid confined by a cuboid box potential, and forced by shaking along one axis. We observe the development of isotropic turbulence from anisotropic forcing for a broad range of forcing amplitudes, and characterise the states through their Fourier spectra, vortex distributions, and spatial correlations. For weak forcing the steady-state wave-action spectrum exhibits a k^-3.5 scaling over wavenumber k; further decomposition uncovers the same power law in both compressible kinetic energy and quantum pressure, while the bulk superfluid remains phase coherent and free from extended vortices. As the forcing energy exceeds the chemical potential, extended vortices develop in the bulk, disrupting the k^-3.5 scaling. The spectrum then transitions to a k^-7/3 regime for compressible kinetic energy only, associated with dense vortex turbulence, and phase coherence limited to the healing length. The strong forcing regime is consistent with an inverse cascade of compressible energy driven by small-scale vortex annihilation. 



Regimes of Steady-State Turbulence in a Quantum Fluid
    Ashton S. Bradley
    September 9, 2024
=====================================================





§ INTRODUCTION

Turbulence is a ubiquitous phenomenon in nature, observed in a wide range of systems from the atmosphere to the oceans <cit.>, and in the dynamics of quantum fluids <cit.>.
Yet the turbulent flows understood in classical fluids <cit.> assume a different character in the quantum realm. The pursuit of understanding quantum turbulence in highly controllable Bose-Einstein condensates (BEC) has shed light on the unique quantum properties of superfluids <cit.>, and uncovered universal aspects of fluid dynamics <cit.>. Quantum fluids exhibit superfluidity and support quantum vortices <cit.>, and compressible Bogoliubov phonons <cit.>, while constituent atoms can also possess internal spin degrees of freedom <cit.>, or interact via long-range dipolar potentials <cit.>. These additional degrees of freedom offer more channels for energy transport, and can cause dramatic departure from classical behavior <cit.>. 

A cascade of energy between scales can emerge when energy is injected into a fluid at a particular length scale and dissipated at another scale. This results in the formation of ever smaller structures down to the dissipation scale (direct cascade) <cit.>, or energy accumulating into coherent structures at the system scale (inverse cascade) <cit.>. If the interaction is local in k-space, power laws can appear in the energy distribution associated with the energy cascade, providing clear signatures of particular dynamics, which may otherwise prove opaque. The most well-known example of emergent power-law behavior in fluid turbulence is the Kolmogorov-Obhukov -5/3 law <cit.>. This states that at high Reynolds number in classical incompressible fluids, the distribution of kinetic energy in Fourier-space is given by
E(k) = C ε^2/3 k^-5/3, 
where C is a universal dimensionless constant and ε is the rate of energy flux per unit mass. Recent experiments <cit.> and numerical investigations <cit.> of turbulent BECs have reported power-laws related to predictions from weak-wave turbulence (WWT) theory <cit.>, observed in the wave-action spectrum accessible in experiments. 

Wave-turbulence involves nonlinearly interacting waves, rather than vortices or eddies in typical hydrodynamic turbulence.  In particular, ideal WWT theory predicts a direct energy cascade with wave-action spectrum power-law n(k) ∝ k^-3 <cit.>. Interestingly, most of the publications report a slightly steeper scaling law than k^-3 <cit.>, conjectured to be due to the residual role of vortices (which are neglected in WWT theory), nonperturbative nonlinear interactions between waves, or quantum pressure <cit.>. Recently it has been argued that the direct energy-cascade is marginally non-local <cit.> and that a logarithmic correction can account for the steeper power-law observed for weak forcing <cit.>. However, that analysis was applied to a four-wave kinetic equation, relevant for Bose gas dynamics when there is no condensate. The need for a deeper understanding of turbulent regimes in a quantum fluid motivates further analysis of the highly condensed regime <cit.>.

In this work, we use the Gross-Pitaevskii Equation (GPE) to simulate the development of steady turbulence in a cuboid box-trap. Energy is injected via an oscillating potential that shakes the box along one axis, and damped by selective removal of high energy atoms <cit.>. By finely varying the forcing energy from weak to strong relative to the superfluid chemical potential we identify regimes of turbulence by their power laws, vortex distributions, and two point correlations. Motivated by a recent reformulation of spectral analysis for quantum fluids <cit.>, we apply those methods to the turbulent states. All quantum-phase information missing from semi-classical velocity power spectra is retained, allowing highly resolved decomposition of wave-action spectra. 

For weak forcing, we observe k^-3.5 scaling as reported in experiments and other numerical studies. Further decomposition shows that the compressible and quantum pressure components are in lock step throughout the power-law range of k, consistent with a Bogoliubov-phonon direct cascade. For strong forcing, the system contains many vortices, and approaches a k^-7/3 power law for the compressible kinetic energy only, consistent with the WWT prediction of inverse particle cascade  scaling in the non-condensed regime <cit.>. The phase coherence length of the GPE wave function drops to the healing length in this regime associated with dense vortex tangles.  

This paper is organized as follows: In secII, we present our theoretical model of forced and damped Gross-Pitaevskii turbulence. secIII details the identification of two distinct turbulent regimes through analysis of density fluctuations, vortex distributionwes, and the components of kinetic energy. In secIV, we apply a reformulated spectral analysis to these regimes. Finally, in secV we discusses the implications of our findings and present our conclusions.



§ THEORETICAL MODEL

 
Far below the critical temperature for condensation, T ≪ T_c, Bose-Einstein condensed dilute gasses are well described by the Gross-Pitaevskii Equation <cit.>. Generalizations of the GPE to finite temperature have also been developed via several different approaches, including exact phase space methods <cit.>, number-conserving Bogoliubov methods <cit.>, the Zaremba-Nikuni-Griffin (ZNG) formulation of two-fluid theory <cit.>, and the truncated-Wigner classical field method <cit.>; see <cit.> for reviews of generalized GPE approaches to the dynamics of the dilute Bose gas.

Under conditions of weak forcing, the GPE provides a quantitative treatment of the system, directly comparable with experiments. Stronger forcing involves large energy injection and requires careful consideration. Due to the numerical scale of the system and the challenges of running large-scale simulations for very long physical times (of order 10s), in this work we take the approach of constructing a minimal GPE model, similar to other recent numerical treatments <cit.>. Our primary differences are using higher working precision which allows consideration of stronger forcing, and larger spatial sizes in units of the healing length. Ideally, we aim to observe a power-law over a decade of wavenumbers in order to rule out the possibility of falsely attributing power-law behavior to a crossover <cit.>. 

The GPE is 

    iħ∂ψ (r,t)/∂ t    = [ -ħ^2/2m∇^2 + V(r,t) + g|ψ(r,t)|^2 ] ψ(r,t),

where g=4πħ^2 a/m, ħ is the reduced Planck constant, a is the S-wave scattering length, and m is the atomic mass. For a homogeneous ground state with chemical potential μ, we define the healing length ξ via

    μ   = g n_0 =ħ^2/mξ^2

The natural units of time and length are τ=ħ/μ and ξ=ħ/√(mμ) respectively. To make our results more concrete, we report our results using parameters for ^87Rb:  m=1.44× 10^-25kg and a=5.8× 10^-9m. 

In Eq. GPEdim, the potential V(r,t), consists of a hard-walled cuboid potential, as well as a time-dependent forcing potential:

    V(r,t)    = V_static(r) + V_F(r,t)

The walls of the static trapping potential are slightly smoothed rather than being piecewise defined, giving both a better description of experimental reality, and improved numerical stability due to eliminating Gibbs phenomena <cit.>. The cuboid box trap has a large flat region with dimensions (L_x,L_y,L_z)=(40,30,20)ξ.  The trap walls are modeled by shifting a steep transition potential

    V_wall(x)    = 
            V_b sin(π x/2ℓ_b)^24  for |x|≤ℓ_b/2
    
            V_bH(x)  otherwise

to coordinates x-L_x, etc, where H(x) is the Heaviside step. Over a region of width ℓ_b=6ξ centered on x=3ξ the  potential V_wall(x) ramps from zero up to V_b = 30μ, providing sufficiently steep walls for box-confined dynamics while avoiding overly sharp changes in the potential; the shape of the potential causes the flat density to drop away steeply at x=0, so that after the shift, the homogeneous region is constrained to ± L_x. Numerically finding the ground state ψ_0(𝐫) of V_static(𝐫) using imaginary time evolution for chemical potential μ = k_B× 1nK, results in a total atom number N ≃ 7.2 × 10^5. The loss of initial anisotropy of the ground state under time evolution provides provide a useful indicator of the development of steady isotropic turbulence.

Generation of steady-state turbulence requires both energy injection and  dissipation. Energy injection is provided by an oscillating linear potential applied along the z-direction in the form

    V_F(r,t)    = U_Fz/L_zsin(ω_Ft),

where ω_F is the driving frequency. U_F is the driving energy scale, here characterising the extra energy per particle introduced at the edge of the box. To explore different turbulent regimes, we use a fixed driving frequency of ω_F = 4 Hz and vary the forcing amplitude over the interval 0.1≤ U_F/μ≤ 5. The lower limit is a very weak perturbation, while the upper limit represents strong forcing. 

Dissipation occurs in the form of high-energy atoms escaping the trapping potential <cit.>. We model this numerically by introducing a dissipative `potential' outside of the box trap. The total static complex potential is    

    V_static(r)    ≡ V_box(r) - iV_diss(r)

where 

    V_diss(r)    =  
         
            2.5μ   if  |u| > 0.5L_u + 5ξ,  u ∈{x,y,z}
    
                    0    otherwise.

The momentum required for a particle to escape is a function of the height of the trap walls <cit.>, and for our trapping potential loss occurs at dimensionless dissipation wavenumber k_Dξ = √(2V_b/μ) = √(60)≃ 7.74 [This choice also gives sufficient margin between the dissipation k-value and the Nyquist k-values for a grid of 256^3 points (Appendix <ref>). High-energy particles are thus dissipated well before their spatial grid representation becomes inaccurate. ]. This energy and spatially selective loss is the only dissipation used in our simulations, and provides the energy sink required for a steady state to form under conditions of steady forcing.

To create turbulence in the forced and damped GPE, we apply V_F to the ground state ψ_0 and evolve until t=t_f≡ 10 seconds. This time interval ensures a steady state is reached for all U_F values. We use a grid of 256^3 points and time-evolve using a pseudo-spectral Runge-Kutta method, referred to as the RK4IP <cit.>. A fixed time step Δ t=10^-3s is chosen to ensure numerical convergence all simulations. The code was written in Julia and optimized to achieve a minimal memory footprint and run on GPU hardware using the CUDA.jl library <cit.>. We achieved a runtime of 38 minutes per second of dynamics on NVIDIA A100 GPU using Float64 representation. For details of the algorithm and error analysis see Appendix <ref>.



§ DEVELOPMENT OF TURBULENCE



 §.§ Vortices and density waves

To outline the general phenomenology, we first present the dynamics that occur at a forcing amplitude of U_F=μ. The initial ground state, shown in 3d_turb(a), responds to forcing by weakly oscillating between the sides of the trap, creating weak density waves, as seen in 3d_turb(b). After sufficient energy is injected, the superfluid develops many small vortex rings parallel to the x-y plane [3d_turb(c)]. After approximately two seconds, the system transitions to an isotropic density and vortex distribution involving many vortices and significant energy, density, and phase fluctuations [3d_turb(d)]. For this particular forcing, there are no extended vortex lines in the steady state, but many small vortex rings are present.
        
To study the response to anisotropic forcing along z, we consider condensate density and vortex distributions on the y=0 plane for a range of forcing amplitudes and times, shown in heatmaps. Vortices are detected by calculating the circulation of the quantum phase around each grid plaquette in the plane. We see that weak forcing excites small amplitude fluctuations, but contains insufficient energy to generate bulk vortices, as seen in heatmaps for U_F=0.2μ. As U_F increases, vortices appear at the edges of the condensate and eventually in the bulk as small rings. They do not survive as extended line vortices in the steady state until the forcing reaches U_F≈ 1.2μ. For strong forcing, e.g. heatmaps, U_F=3μ, the steady state involves a dense and isotropic vortex distribution. The transition from weak surface vortices to dense bulk vortices in the steady state is evident in the three-dimensional vortex distribution, heatmaps(b), but harder to identify due to the proliferation of surface vortices. 
 
excite_energies(a) and (b) show the normalization and energy per particle of the system during excitation. We observe a brief delay between when forcing begins and when particles have sufficient energy to escape the trap, shown in excite_energies(a). Consequently, the energy per particle of the system rapidly increases [excite_energies(b)]. Once dissipation takes effect, the energy per particle approaches an approximate steady-state, oscillating about a mean value that increases with U_F; the system thus approaches an energetic steady state while continuing to lose particles.  


 §.§ Energy decomposition

To further investigate the role of vortices and compressible excitations, we perform a decomposition of the kinetic energy <cit.>. Using a Madelung-transform, we can represent the complex wave function (away from vortex cores) in terms of two real fields: the density, ρ = |ψ|^2, and the quantum phase, Θ.

    ψ(r,t)    = √(ρ(r,t))e^iΘ(r,t)

Using the velocity field v(r,t) = ħ/m ∇Θ (r,t), we construct the density-weighted velocity field [The ordinary velocity field  𝐯 diverges near the core of vortices, and so cannot be used for decomposition. By using 𝐮 we can apply the Helmholtz decomposition to a twice continuously differentiable vector field, as required by the decomposition theorem.], u = √(ρ)v, and perform a Helmholtz decomposition of the vector field into incompressible and compressible components

    u   = u^i + u^c,

where the incompressible component is divergence-free and associated with vortices, and the compressible component is curl-free, associated with density excitations:

    ∇·u^i    = 0,    ∇×u^c = 0.

We also define a quantum pressure vector field <cit.>

    u^q    = (ħ/m) ∇√(ρ).

The quantum pressure is not a physical velocity field, but instead arises from density gradients. The total kinetic energy can then be written as the sum of the three components

    E_kin   ≡∫ d^3𝐫ħ^2|∇ψ|^2/2m=E_kin^c+E_kin^i+E_kin^q,

where 

    E_kin   =m/2∫ d^3𝐫|𝐮^α(𝐫)|^2.


Figures <ref>(c) and (d) show E_kin^i and E_kin^c during the initial period of excitation. Both components oscillate in response to the forcing, acquiring a mean value that increases with U_F. The kinetic energy components of the steady-state display a nonlinear dependence on U_F, indicating different regimes of turbulence, shown in excite_energies(e). The relative importance of compressible and incompressible energy indicates a change from wave to vortex turbulence. Significant kinetic energy is stored in density gradients even in the absence of vortices, evident in the quantum pressure.
      
Figure <ref>(e) is consistent with the behavior shown in heatmaps: from U_F = 0.2μ, E^c_kin sharply increases as density fluctuations grow [excite_energies(e)]; for the same U_F values, E^i_kin increases gradually as vortices are present only at the condensate surface. Recent reports of wave-turbulence in BECs generally correspond to these forcing amplitudes (U_F < μ) <cit.>, which is consistent with the dominance of compressible kinetic energy in this region. 

At U_F ≃1.2μ, a distinct change is evident in excite_energies(e). E^i_kin increases rapidly, associated with the formation of bulk vortices, and E^c_kin has a much weaker response to increased forcing. Finally, another change in response occurs for U_F > 2μ, where incompressible kinetic energy and quantum pressure have a weaker response to increased forcing. In heatmaps this corresponds to the steady state vortex distribution becoming homogeneous. There may be a limit to the number of vortices that can stably exist within a BEC of a certain size and density, bounding growth of incompressible kinetic energy with forcing. 

This analysis suggests distinct turbulent regimes exist for different forcing strengths. In the following sections, we apply tools of spectral analysis to relate these regimes to previous observations and theoretical predictions, aiming to identify power law behavior associated with particular regimes in excite_energies(e).


 §.§ Spectral Analysis

    Our general aim is to identify power-law behavior in kinetic energy spectra, to gain insight into the underlying mechanism of transport and the state of the turbulent Gross-Pitaevskii field. In this section, we summarize relevant results of the spectral analysis formulation presented in Ref. <cit.>, that will allow such identification.
    
For any two complex vector fields, u and v, the spectral density of their inner product ⟨𝐮𝐯⟩(k), is defined as 

    ⟨u | v⟩≡∫_0^∞ dk ⟨uv⟩(k).

It can be shown <cit.> that in three dimensions 

    ⟨u || v⟩(k)    = k^2/2π^2∫ d^3r sin(k |𝐫|)/k|𝐫| C[u,v](r),

where C[u,v](r) is the spatial two-point correlation of u and v:

    C[u,v](r) ≡∫ d^3 R ⟨u| R - r/2 ⟩⟨R + r/2|v⟩.

This reformulation of the spectral analysis problem offers several advantages. First, carrying out the angular integral in Fourier-space analytically, removes the need for numerical binning to create spectral functions of k=|𝐤|. Second, by decoupling the position and Fourier-space resolution, it allows arbitrary resolution spectra to be computed. Finally, a corresponding system-averaged two-point correlation can easily be extracted <cit.>:

    G_u v(r)   =1/4π∫ d^3 𝐫^'δ(r-|𝐫^'|) C[𝐮, 𝐯](𝐫^')
       =∫_0^∞ dk sin(k r)/kr⟨𝐮𝐯⟩(k).

    
Depending on which fields u and v are inserted into spec_IP, a variety of spectral distributions can be generated for any given wave function. For example, choosing u = v = ψ(r), we obtain

    N    = ∫_0^∞ dk ⟨ψψ⟩(k),

where the momentum spectrum, describing the particle occupation of momentum-states is

    ⟨ψψ⟩(k)   =k^2/2π^2∫ d^3𝐫sin(k |𝐫|)/krC[ψ,ψ](𝐫),

with corresponding two-point spatial correlator of the field 

    G_ψψ(r)   =∫_0^∞ dk sin(k r)/kr⟨ψψ⟩(k).


Eqs. spec_dens and N_int are one-dimensional integrals in k-space. However, spectra in WWT literature <cit.>  
are usually defined such that a three-dimensional integral returns the total inner product. This introduces a factor of 4π k^2 in the definition of the wave-action spectrum

    n(k)    ≡1/4π k^2⟨ψψ⟩ (k),

where 

    N    = ∫ d^3𝐤  n(k) = ∫_0^∞ dk   4π k^2 n(k)

        
We emphasize that there is no assumption of isotropy in the definition. Instead, formally integrating over the Fourier-space solid angle generates the sinc kernel in psispec. Another example is found by instead using the gradient of the wave function, 𝐮=𝐯=∇ψ, giving the kinetic energy spectral-density (again introducing the spherical factor 4π k^2)

    e_kin(k) ≡1/4π k^2⟨∇ψ∇ψ⟩(k),

related to the total kinetic energy by

    E_kin   = ∫ d^3k  e_kin(k) = ∫_0^∞ dk   4π k^2 e_kin(k).


We now address the decomposition problem. To set up a decomposition of the kinetic energy that includes all quantum phase information <cit.>, we define three fields:

    (w^i,w^c, w^q)    ≡ (iu^i, iu^c, u^q) e^iΘ.

These definitions furnish a linear decomposition of the gradient ∇ψ = (m / ħ)[ w^i + w^c + w^q ]. Using w^α as input fields to spec_dens, we obtain the respective components of the kinetic energy-density:

    e_kin^α(k)    ≡1/4π k^2⟨w^αw^α⟩(k), α∈i,c,q.

The kinetic energy densities calculated in e_kin_components retain all quantum phase information [The standard decomposition in quantum-turbulence studies <cit.> involves kinetic energy density calculated semi-classically as a velocity power spectrum via Plancherel's theorem (Appendix. <ref>), which discards quantum phase information.]. If instead we ignore the phase factor, as is usually done <cit.>, we would obtain the standard velocity power spectra [The usual definition does not factor out the spherical measure 4π k^2, here written explicitly for convenience in connecting with recent literature.].

    ε_kin^α(k)    ≡1/4π k^2⟨u^αu^α⟩(k) , α∈i,c,q.

An important feature of the velocity power spectra components is that while they integrate to give the respective component of kinetic energy,  E^α_kin, they do not sum locally in k-space to give the power spectral density <cit.>

    ε_kin(k)    ≠ε^i_kin(k) + ε^c_kin(k) + ε^q_kin(k).

Conversely, the components of the energy density do add locally, albeit with additional cross-terms describing redistribution of energy between scales

    e_kin(k)    = e^i_kin(k) +    e^c_kin(k) + e^q_kin(k) 
       + e^ic_kin(k) + e^iq_kin(k) + e^cq_kin(k).


The cross-terms integrate to zero and thus do not contribute to the total kinetic energy, but are required to get the correct total energy at each k <cit.>. In our simulations e^α_kin(k) and the corresponding ε^α_kin(k) were found to differ significantly, as detailed in Appendix <ref>.

As an example to check our formulation, we plot the time evolution of n(k) in k3.5. As forcing is initiated, a peak in the spectrum is evident at the wavenumber k_L_z=2π/L_z, corresponding to the forcing lengthscale. This peak then migrates towards higher k-values, as expected for a direct cascade. Shortly after reaching the dissipation length scale, n(k) remains approximately constant in time apart from a change in amplitude due to particle loss. For this particular forcing, the power law form n(k)∼ k^-3.5 is evident in k3.5, consistent with previous works <cit.>. 



§ STEADY STATE TURBULENCE

We now investigate the properties of steady-state turbulence. We  use the spectral analysis developed in the previous section to further decompose the wave-action spectrum. As kinetic energy densities and number densities only differ by a factor of ħ^2k^2/2m, the spectral decomposition allterms can be used to define the component wave-action spectra 

    n^α(k)     = 2m/ħ^2k^2 e^α_kin(k) = m/2πħ^2k^4⟨w^αw^α⟩(k),

where α∈{i,c,q}. 
We can then use the previously established decomposition of energy-densities to identify power-law behavior in distinct parts of the wave-action spectrum, considering the compressible, incompressible, and quantum pressure components [We do not consider the cross terms further in this work as they are not very intuitive to interpret. See Ref. <cit.> for a simple example.]. 

We decompose the steady-state wave-action spectrum according to nk_i and present the densities n^α(k) in decomposed_nk. The decomposed spectra reveal further information beyond what is available from n(k). The presented steady-state spectra are averaged over one forcing period, however, time averaging only alters the result near the forcing lengthscale, a further indication that the steady-state is reached. 

We do not observe any obvious power-law behaviour for U_F≪μ (not shown). Within a narrow region of forcing, 0.15μ≲ U_F ≲ 0.3μ, the compressible density approximates k^-3. However, the approximate scaling only extends over a factor of 2 in wave number, insufficient to identify power-law behaviour <cit.>.
Starting at U_F=0.4μ, we observe a clear k^-3.5 scaling over nearly a decade of kξ, in both compressible and quantum pressure momentum distributions, shown in decomposed_nk(a). This power-law remains intact for a wide interval of forcing energies, 0.4μ≤ U_F≤ 1.2μ, and as incompressible energy remains small, the power-law is also evident in the total spectrum n(k). This is consistent with the observed vortex configurations [heatmaps(a)], where extended bulk vortices remain rare. In this regime it is notable that n^c(k)≃ n^q(k) in the power law region of wavenumbers. We return to this point in discussion.
        
For U_F ≳ 1.2μ, extended vortices develop in the bulk superfluid in the steady-state, and the system enters a transitional regime without power-law behavior in any component of kinetic energy [decomposed_nk(b)]. The total wave-action spectrum gradually steepens in the ultraviolet due to an increase in the incompressible energy-density, while the compressible and quantum pressure components uncouple, and the slope of n^c(k) decreases. The appearance of bulk vortices is associated with the breakdown of k^-3.5 scaling. Eventually, we observe a wide n^c(k) ∝ k^-7/3 power-law for strong forcing, U_F≳ 2μ, in agreement with the predicted inverse particle cascade for weak forcing in the absence of a condensate [decomposed_nk(c)]  <cit.>.  

To further understand the different regimes of steady turbulence, in decomposed_nk(d) we plot the ratio of incompressible to compressible kinetic energy, with the identified power-law regimes. There is a clear change in the energy dependence on U_F at the boundary of each power-law regime. Initially, the ratio is rapidly decreasing, due to compressible excitations paired with a lack of vortex nucleation. With additional energy injected, this is succeeded by a sharp increase in the kinetic energy ratio (as surface vortices begin to proliferate), and the appearance of a k^-3.5 power-law. The breakdown of k^-3.5 involves extended vortices entering the bulk in the steady-state, and a steep rise in incompressible energy with U_F. The system enters a regime of vortex turbulence, but shows no scaling behavior. Finally, for U_f≳ 2μ, a k^-7/3 scaling emerges where vortices are approximately homogeneous and dense, and the incompressible energy responds only weakly to increased forcing. In this regime E^i_kin/E^c_kin > 1, a regime we can regard as strong quantum vortex turbulence, in contrast to weak wave turublence where E^i_kin/E^c_kin≪ 1.

While Fourier-space reveals power-law behavior, we can also examine the corresponding correlation functions in position-space. When power-law behavior weakens, the appearance of characteristic length scales can be more clearly understood in position space. Fourier-transformation of the momentum distribution gives the system-averaged spatial two-point correlation of the Gross-Piteavskii wave function, as discussed in speca. In two_point(a), we plot the system averaged two-point correlation function g_1(r)≡ G_ψψ(r)/G_ψψ(0), calculated by transforming ⟨ψψ⟩(k)=4π k^2n(k) to position space via Gpp. The two-point correlator of the ground state (U_F=0) decays over a scale comparable to the shortest cuboid length L_z. As U_F increases, the length scale is reduced. Eventually, for strong forcing the correlations only persist over a scale of order ξ. We define the correlation length l_c as the distance where g_1(r) drops to 1/e. In two_point(b) we plot the correlation length against U_F, together with the power law regimes identified. 

Initially there is a steep drop, followed by a slower decline during the k^-3.5 regime. As vortex lines develop in steady-state, we see a second steep decline in correlation length around U_F∼1.5μ. The dense vortex regime U_F≳2μ corresponds to l_c∼ξ, and the system loses all coherence beyond the vortex core scale.   



§ DISCUSSION AND CONCLUSIONS



 §.§ Discussion

We identify three distinct regimes of steady turbulence through a combination of power laws, vortex distributions, and two point correlations.

Weak wave turbulence.— For weaker forcing, 0.4μ≲ U_F≲ 1.2μ, there is a ∝ k^-3.5 power law spectrum for compressible kinetic energy, quantum pressure, and the total spectrum, and a distinct absence of extended vortex lines in the bulk superfluid. In this weak-wave turbulence regime E^i_kin/E^c_kin≪ 1. Even when small vortex rings appear in the steady state as U_F approaches μ, the k^-3.5 scaling is robust and spans a decade of k-space. The forcing amplitude reported in Ref. <cit.>, U_F = 0.8μ, is in the center of the U_F range where we identify the scaling both instantaneously [k3.5] and in time averages [decomposed_nk]. In addition, we observe a clear forcing peak near the system scale L_z in the total wave-action spectrum at early times. The wave-action spectra decomposition, combined with absence of bulk vortices implicates a direct Bogoliubov quasiparticle cascade for the k^-3.5 power law. Indeed, Bogoliubov phonons involve a quantum pressure component associated with small scale density fluctuations when kξ≳ 1. In decomposed_nk(b), we see n^c(k) and n^q(c) are closely matched within this k range, and decouple when kξ≲ 1. 

The log-corrected spectrum analysis of WWT direct cascade is strictly applicable for very weak forcing in the absence of a condensate <cit.>, and is thus not directly applicable to our system that instead involves a large condensate with dominant three wave interactions. However, the three wave kinetic equation also supports a k^-3 scaling <cit.>. Applying the log-correction to our wave-action spectra, we find good agreement if we identify the forcing scale for phonons as k_F=1/ξ, instead of the natural forcing wavenumber of the potential 2π/L_z as detailed in Appendix <ref>. Indeed, the analysis of the forcing scale by Martirosyan et al. <cit.>, reached a similar conclusion [Their result translates to k_F=√(2)/ξ in our units. However, for this estimate of the forcing scale we see a narrower flat region in the log-corrected spectrum. For further discussion see Appendix <ref>.].

Mixed turbulence.— For 1.2μ≲ U_F≲ 2.5μ the system enters a transitional regime where the k^-3.5 power-law breaks down, and power law behavior is not strongly evident in any components. The compressible and quantum pressure components decouple in k space, and the incompressible energy increases rapidly with U_F, entering a regime where E^i_kin/E^c_kin∼ 1, and the compressible and incompressible wave-action spectra become comparable in k space. This is associated with the appearance of extended vortices in the bulk superfluid. The transition is also evident in the two-point correlation function, where the correlation length decreases rapidly as vortices enter the bulk.

Strong vortex turbulence.— For 2.5μ≲ U_F, the condensate bulk becomes saturated with vortex lines and the compressible spectrum approximates a k^-7/3 power-law over more than a decade of k; this is consistent with the predicted WWT scaling for an inverse particle-cascade solution to the four-wave kinetic equation. In Ref. <cit.>, this was observed in the total spectrum for weak forcing in the absence of BEC, while here we see it in the compressible spectrum only, for strong forcing. As the power law extends far below k_ξ=1/ξ, this would be consistent with WWT if random vortices effectively destroy the condensate, and indeed the phase coherence length drops to the healing length at U_F≳ 2μ [two_point]. As our model only contains damping for high wavenumbers, a steady inverse-cascade of particles can only exist in n^c(k) if particles are transported to high k by the other components of n(k). Physical consistency thus requires a compressible particle source near the dissipation scale k∼ k_D to generate an inverse cascade; one candidate mechanism is the steady annihilation of vortices, creating bursts of sound at small scales. In Figure <ref>(d),(e) the compressible energy driving by the forcing appears to become ineffective for U_F≳ 2μ, supporting the interpretation that compressible energy is instead injected at small scales via vortex annihilation.

In the high-energy regime, the incompressible energy density is extremely stable across both U_F values and forcing phase. In all other cases, spectra visibly oscillate in the infrared with forcing, necessitating cycle-averaging of spectra. However, the dense distribution of vortices appears to depend very weakly on the phase of the forcing. Wave-turbulence closure requires the initial phase and amplitude of the Gross-Pitaevskii field be independent and random, whereas we find that for weak forcing there is still a high degree of coherence. Conversely, once the condensate becomes saturated with vortices, we see a drastic decrease in the correlation length, consistent with the assumptions underlying WWT. Despite strong forcing, wave interactions may be sufficiently weak that WWT theory remains relevant, particularly if injected energy primarily resides in the quantum pressure and incompressible components of kinetic energy.

Calculation of spectra with arbitrary infrared k-resolution and decomposition of wave-action spectra require the reformulated spectral densities <cit.>. In Appendix <ref> we compare the decomposition via the wave-action spectrum with the standard decomposition of the the velocity power spectrum used in previous studies. We find the two are similar only in a minority of cases, typically in the inertial range for well-developed turbulence.


 §.§ Conclusions

Motivated by recent experiments in forced quantum turbulence <cit.>, we have simulated a Gross-Pitaevskii model for a wide range of forcing energies, generating and analysing steady turbulence. We find evedence for regimes of weak-wave and strong-vortex turbulence separated by a mixed regime that lacks clear power-law behavior. Distinct regimes of power-law behavior are identified by application of high resolution spectral analysis <cit.>, and further studied through vortex distributions and two-point correlations. Our results are largely consistent with other recent works <cit.>, while also enabling further decomposition of the wave-action spectrum. For weak forcing, the bulk superfluid is approximately vortex free, and exhibits a power-law of ∼ k^-3.5, consistent with weak-wave turbulence predictions for a direct cascade <cit.>. The compressible and quantum pressure components also show the same power-law scaling and are in close agreement. Log-correction analysis implicates a forcing wavenumber k_F∼ξ^-1 for the cascade. In this regime the superfluid order parameter remains coherent over the system scale, with coherent only weakly degraded by small vortex rings. In contrast, for strong forcing the steady state develops a dense tangle of extended vortices that limites phase coherence to the healing length scale. Power-law scaling ∼ k^-7/3 is observed in the compressible component only, consistent with the weak-wave turbulence prediction of an inverse cascade in non-condensate four-wave kinetics <cit.>.

Several open questions would benefit from future exploration. It would be interesting to study the rapid vortex growth regime in the transition from k^-3.5 to k^-7/3 scaling, together with the different kinetic equation approaches and their associated scaling solutions. Another open question is identifying the true source of high-k compressible energy in the k^-7/3 regime. Detailed study of vortex tangles, energy fluxes, and velocity increments <cit.> would also reveal further insights into the nature and formation of quantum turbulence. 



§ ACKNOWLEDGEMENTS

We thank Gevorg Martirosyan and Tim Copland for stimulating discussions.
TF would like to thank the University of Otago for financial support.

The authors wish to acknowledge the use of New Zealand eScience Infrastructure (NeSI) high performance computing facilities as part of this research.












47
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Grant et al.(1962)Grant, Stewart, and Moilliet]grant_turbulence_1962
  author author H. L. Grant, author R. W. Stewart, and author A. Moilliet, title title Turbulence spectra from a tidal channel, https://doi.org/10.1017/S002211206200018X journal journal Journal of Fluid Mechanics volume 12, pages 241 (year 1962)NoStop
[Henn et al.(2009)Henn, Seman, Roati, Magalhães, and Bagnato]henn_emergence_2009
  author author E. A. L. Henn, author J. A. Seman, author G. Roati, author K. M. F. Magalhães, and author V. S. Bagnato, title title Emergence of Turbulence in an Oscillating Bose-Einstein Condensate, https://doi.org/10.1103/PhysRevLett.103.045301 journal journal Physical Review Letters volume 103, pages 045301 (year 2009)NoStop
[Neely et al.(2013)Neely, Bradley, Samson, Rooney, Wright, Law, Carretero-González, Kevrekidis, Davis, and Anderson]neely_characteristics_2013
  author author T. W. Neely, author A. S. Bradley, author E. C. Samson, author S. J. Rooney, author E. M. Wright, author K. J. H. Law, author R. Carretero-González, author P. G. Kevrekidis, author M. J. Davis, and author B. P. Anderson, title title Characteristics of Two-Dimensional Quantum Turbulence in a Compressible Superfluid, https://doi.org/10.1103/PhysRevLett.111.235301 journal journal Physical Review Letters volume 111, pages 235301 (year 2013)NoStop
[Kolmogorov(1941)]kolmogorov_local_1941
  author author A. Kolmogorov, title title The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds' Numbers, https://ui.adsabs.harvard.edu/abs/1941DoSSR..30..301K journal journal Akademiia Nauk SSSR Doklady volume 30, pages 301 (year 1941)NoStop
[Navon et al.(2016)Navon, Gaunt, Smith, and Hadzibabic]navon_emergence_2016
  author author N. Navon, author A. L. Gaunt, author R. P. Smith, and author Z. Hadzibabic, title title Emergence of a turbulent cascade in a quantum gas, https://doi.org/10.1038/nature20114 journal journal Nature volume 539, pages 72 (year 2016)NoStop
[Gauthier et al.(2019)Gauthier, Reeves, Yu, Bradley, Baker, Bell, Rubinsztein-Dunlop, Davis, and Neely]gauthier_giant_2019
  author author G. Gauthier, author M. T. Reeves, author X. Yu, author A. S. Bradley, author M. A. Baker, author T. A. Bell, author H. Rubinsztein-Dunlop, author M. J. Davis, and author T. W. Neely, title title Giant vortex clusters in a two-dimensional quantum fluid, https://doi.org/10.1126/science.aat5718 journal journal Science volume 364, pages 1264 (year 2019)NoStop
[Reeves et al.(2015)Reeves, Billam, Anderson, and Bradley]reeves_identifying_2015
  author author M. T. Reeves, author T. P. Billam, author B. P. Anderson, and author A. S. Bradley, title title Identifying a Superfluid Reynolds Number via Dynamical Similarity, https://doi.org/10.1103/PhysRevLett.114.155302 journal journal Physical Review Letters volume 114, pages 155302 (year 2015)NoStop
[García-Orozco et al.(2022)García-Orozco, Madeira, Moreno-Armijos, Fritsch, Tavares, Castilho, Cidrim, Roati, and Bagnato]garcia-orozco_universal_2022
  author author A. D. García-Orozco, author L. Madeira, author M. A. Moreno-Armijos, author A. R. Fritsch, author P. E. S. Tavares, author P. C. M. Castilho, author A. Cidrim, author G. Roati, and author V. S. Bagnato, title title Universal dynamics of a turbulent superfluid Bose gas, https://doi.org/10.1103/PhysRevA.106.023314 journal journal Physical Review A volume 106, pages 023314 (year 2022)NoStop
[Dogra et al.(2023)Dogra, Martirosyan, Hilker, Glidden, Etrych, Cao, Eigen, Smith, and Hadzibabic]dogra_universal_2023
  author author L. H. Dogra, author G. Martirosyan, author T. A. Hilker, author J. A. P. Glidden, author J. Etrych, author A. Cao, author C. Eigen, author R. P. Smith, and author Z. Hadzibabic, title title Universal equation of state for wave turbulence in a quantum gas, https://doi.org/10.1038/s41586-023-06240-z journal journal Nature volume 620, pages 521 (year 2023)NoStop
[Madeira et al.(2024)Madeira, García-Orozco, Moreno-Armijos, Fritsch, and Bagnato]madeira_universal_2024
  author author L. Madeira, author A. D. García-Orozco, author M. A. Moreno-Armijos, author A. R. Fritsch, and author V. S. Bagnato, title title Universal scaling in far-from-equilibrium quantum systems: An equivalent differential approach, https://doi.org/10.1073/pnas.2404828121 journal journal Proceedings of the National Academy of Sciences volume 121, pages e2404828121 (year 2024)NoStop
[Onsager(1949)]onsager_statistical_1949
  author author L. Onsager, title title Statistical hydrodynamics, https://doi.org/10.1007/BF02780991 journal journal Il Nuovo Cimento (1943-1954) volume 6, pages 279 (year 1949)NoStop
[Bogoliubov(1947)]bogoliubov_theory_1947
  author author N. N. Bogoliubov, title title On the Theory of Superfluidity, @noop  journal journal Journal of Physics (USSR) volume 11, pages 23 (year 1947)NoStop
[Tsubota and Fujimoto(2014)]tsubota_spin_2014
  author author M. Tsubota and author K. Fujimoto, title title Spin turbulence in spinor Bose-Einstein condensates, https://doi.org/10.1088/1742-6596/497/1/012002 journal journal Journal of Physics: Conference Series volume 497, pages 012002 (year 2014)NoStop
[Chomaz et al.(2022)Chomaz, Ferrier-Barbut, Ferlaino, Laburthe-Tolra, Lev, and Pfau]chomaz_dipolar_2022
  author author L. Chomaz, author I. Ferrier-Barbut, author F. Ferlaino, author B. Laburthe-Tolra, author B. L. Lev, and author T. Pfau, title title Dipolar physics: A review of experiments with magnetic quantum gases, https://doi.org/10.1088/1361-6633/aca814 journal journal Reports on Progress in Physics volume 86, pages 026401 (year 2022)NoStop
[Navon et al.(2019)Navon, Eigen, Zhang, Lopes, Gaunt, Fujimoto, Tsubota, Smith, and Hadzibabic]navon_synthetic_2019
  author author N. Navon, author C. Eigen, author J. Zhang, author R. Lopes, author A. L. Gaunt, author K. Fujimoto, author M. Tsubota, author R. P. Smith, and author Z. Hadzibabic, title title Synthetic dissipation and cascade fluxes in a turbulent quantum gas, https://doi.org/10.1126/science.aau6103 journal journal Science volume 366, pages 382 (year 2019)NoStop
[Johnstone et al.(2019)Johnstone, Groszek, Starkey, Billington, Simula, and Helmerson]johnstone_evolution_2019
  author author S. P. Johnstone, author A. J. Groszek, author P. T. Starkey, author C. J. Billington, author T. P. Simula, and author K. Helmerson, title title Evolution of large-scale flow from turbulence in a two-dimensional superfluid, https://doi.org/10.1126/science.aat5793 journal journal Science volume 364, pages 1267 (year 2019)NoStop
[Kraichnan(1967)]kraichnan_inertial_1967
  author author R. H. Kraichnan, title title Inertial Ranges in Two-Dimensional Turbulence, https://doi.org/10.1063/1.1762301 journal journal The Physics of Fluids volume 10, pages 1417 (year 1967)NoStop
[Shukla and Nazarenko(2022)]shukla_nonequilibrium_2022
  author author V. Shukla and author S. Nazarenko, title title Nonequilibrium Bose-Einstein condensation, https://doi.org/10.1103/PhysRevA.105.033305 journal journal Physical Review A volume 105, pages 033305 (year 2022)NoStop
[Zhu et al.(2023)Zhu, Semisalov, Krstulovic, and Nazarenko]zhu_direct_2023
  author author Y. Zhu, author B. Semisalov, author G. Krstulovic, and author S. Nazarenko, title title Direct and Inverse Cascades in Turbulent Bose-Einstein Condensates, https://doi.org/10.1103/PhysRevLett.130.133001 journal journal Physical Review Letters volume 130, pages 133001 (year 2023)NoStop
[Zakharov et al.(1992)Zakharov, Lvov, and Falkovich]zakharov_kolmogorov_1992-1
  author author V. Zakharov, author V. Lvov, and author G. E. Falkovich, @noop  title Kolmogorov Spectra of Turbulence 1: Wave Turbulence (publisher Springer-Verlag, address New York, year 1992)NoStop
[Nazarenko(2011)]nazarenko_wave_2011
  author author S. Nazarenko, https://doi.org/10.1007/978-3-642-15942-8 title Wave Turbulence, series Lecture Notes in Physics, Vol. volume 825 (publisher Springer, address Berlin, Heidelberg, year 2011)NoStop
[Dyachenko et al.(1992)Dyachenko, Newell, Pushkarev, and Zakharov]dyachenko_optical_1992
  author author S. Dyachenko, author A. C. Newell, author A. Pushkarev, and author V. E. Zakharov, title title Optical Turbulence - Weak Turbulence, Condensates and Collapsing Filaments in the Nonlinear Schrodinger-Equation, http://gateway.webofknowledge.com/gateway/Gateway.cgi?GWVersion=2   SrcAuth=mekentosj   SrcApp=Papers   DestLinkType=FullRecord   DestApp=WOS   KeyUT=A1992JD72100006 journal journal Physica D volume 57, pages 96 (year 1992)NoStop
[Martirosyan et al.(2024)Martirosyan, Fujimoto, and Navon]martirosyan_equation_2024
  author author G. Martirosyan, author K. Fujimoto, and author N. Navon, https://doi.org/10.48550/arXiv.2407.08738 title An Equation of State for Turbulence in the Gross-Pitaevskii model (year 2024), https://arxiv.org/abs/2407.08738 arXiv:2407.08738 [cond-mat, physics:physics] NoStop
[Bradley et al.(2022)Bradley, Kumar, Pal, and Yu]bradley_spectral_2022
  author author A. S. Bradley, author R. K. Kumar, author S. Pal, and author X. Yu, title title Spectral analysis for compressible quantum fluids, https://doi.org/10.1103/PhysRevA.106.043322 journal journal Physical Review A volume 106, pages 043322 (year 2022)NoStop
[Dalfovo et al.(1999)Dalfovo, Giorgini, Pitaevskii, and Stringari]dalfovo_theory_1999
  author author F. Dalfovo, author S. Giorgini, author L. P. Pitaevskii, and author S. Stringari, title title Theory of Bose-Einstein condensation in trapped gases, https://doi.org/10.1103/RevModPhys.71.463 journal journal Reviews of Modern Physics volume 71, pages 463 (year 1999)NoStop
[Steel et al.(1998)Steel, Olsen, Plimak, Drummond, Tan, Collett, Walls, and Graham]steel_dynamical_1998
  author author M. J. Steel, author M. K. Olsen, author L. I. Plimak, author P. D. Drummond, author S. M. Tan, author M. J. Collett, author D. F. Walls, and author R. Graham, title title Dynamical quantum noise in trapped Bose-Einstein condensates, https://doi.org/10.1103/PhysRevA.58.4824 journal journal Physical Review A volume 58, pages 4824 (year 1998)NoStop
[Gardiner and Morgan(2007)]gardiner_number-conserving_2007
  author author S. A. Gardiner and author S. A. Morgan, title title Number-conserving approach to a minimal self-consistent treatment of condensate and noncondensate dynamics in a degenerate Bose gas, https://doi.org/10.1103/PhysRevA.75.043621 journal journal Physical Review A volume 75, pages 261 (year 2007)NoStop
[Zaremba et al.(1999)Zaremba, Nikuni, and Griffin]zaremba_dynamics_1999
  author author E. Zaremba, author T. Nikuni, and author A. Griffin, title title Dynamics of Trapped Bose Gases at Finite Temperatures, https://doi.org/10.1023/A:1021846002995 journal journal Journal of Low Temperature Physics volume 116, pages 277 (year 1999)NoStop
[Sinatra et al.(2001)Sinatra, Lobo, and Castin]sinatra_classical-field_2001
  author author A. Sinatra, author C. Lobo, and author Y. Castin, title title Classical-field method for time dependent Bose-Einstein condensed gases, https://doi.org/10.1103/PhysRevLett.87.210404 journal journal Physical Review Letters volume 87, pages 210404 (year 2001)NoStop
[Graham(1973)]graham_statistical_1973
  author author R. Graham, title title Statistical Theory of Instabilities in Stationary Nonequilibrium Systems with Applications to Lasers and Nonlinear Optics, in @noop  booktitle Springer Tracts in Modern Physics, editor edited by editor G. Hohler (publisher Springer, year 1973) p. pages 1NoStop
[Polkovnikov(2010)]polkovnikov_phase_2010
  author author A. Polkovnikov, title title Phase space representation of quantum dynamics, https://doi.org/10.1016/j.aop.2010.02.006 journal journal Annals of Physics volume 325, pages 1790 (year 2010)NoStop
[Blakie† et al.(2008)Blakie†, Bradley†, Davis, Ballagh, and Gardiner]blakie_dynamics_2008
  author author P. Blakie†, author A. Bradley†, author M. Davis, author R. Ballagh, and author C. Gardiner, title title Dynamics and statistical mechanics of ultra-cold Bose gases using c-field techniques, https://doi.org/10.1080/00018730802564254 journal journal Advances in Physics volume 57, pages 363 (year 2008)NoStop
[Proukakis and Jackson(2008)]proukakis_finite-temperature_2008
  author author N. P. Proukakis and author B. Jackson, title title Finite-temperature models of Bose–Einstein condensation, https://doi.org/10.1088/0953-4075/41/20/203002 journal journal Journal of Physics B: Atomic, Molecular and Optical Physics volume 41, pages 203002 (year 2008)NoStop
[Clauset et al.(2009)Clauset, Shalizi, and Newman]clauset_power-law_2009
  author author A. Clauset, author C. R. Shalizi, and author M. E. J. Newman, title title Power-Law Distributions in Empirical Data, https://doi.org/10.1137/070710111 journal journal SIAM Review volume 51, pages 661 (year 2009)NoStop
[Piotrowska et al.(2019)Piotrowska, Miller, and Schnetter]piotrowska_spectral_2019
  author author J. Piotrowska, author J. M. Miller, and author E. Schnetter, title title Spectral methods in the presence of discontinuities, https://doi.org/10.1016/j.jcp.2019.03.048 journal journal Journal of Computational Physics volume 390, pages 527 (year 2019)NoStop
[Note1()]Note1
  note This choice also gives sufficient margin between the dissipation k-value and the Nyquist k-values for a grid of 256^3 points (Appendix <ref>). High-energy particles are thus dissipated well before their spatial grid representation becomes inaccurate.Stop
[Caradoc-Davies et al.(2000)Caradoc-Davies, Ballagh, and Blakie]caradoc-davies_three-dimensional_2000
  author author B. M. Caradoc-Davies, author R. J. Ballagh, and author P. B. Blakie, title title Three-dimensional vortex dynamics in Bose-Einstein condensates, https://doi.org/10.1103/PhysRevA.62.011602 journal journal Physical Review A volume 62, pages 011602 (year 2000)NoStop
[Besard et al.(2019a)Besard, Churavy, Edelman, and Sutter]besard_rapid_2019
  author author T. Besard, author V. Churavy, author A. Edelman, and author B. D. Sutter, title title Rapid software prototyping for heterogeneous and distributed platforms, https://doi.org/10.1016/j.advengsoft.2019.02.002 journal journal Advances in Engineering Software volume 132, pages 29 (year 2019a)NoStop
[Besard et al.(2019b)Besard, Foket, and De Sutter]besard_effective_2019
  author author T. Besard, author C. Foket, and author B. De Sutter, title title Effective Extensible Programming: Unleashing Julia on GPUs, https://doi.org/10.1109/TPDS.2018.2872064 journal journal IEEE Transactions on Parallel and Distributed Systems volume 30, pages 827 (year 2019b)NoStop
[Nore et al.(1997)Nore, Abid, and Brachet]nore_kolmogorov_1997
  author author C. Nore, author M. Abid, and author M. E. Brachet, title title Kolmogorov Turbulence in Low-Temperature Superflows, https://doi.org/10.1103/PhysRevLett.78.3896 journal journal Physical Review Letters volume 78, pages 3896 (year 1997)NoStop
[Note2()]Note2
  note The ordinary velocity field 𝐯 diverges near the core of vortices, and so cannot be used for decomposition. By using 𝐮 we can apply the Helmholtz decomposition to a twice continuously differentiable vector field, as required by the decomposition theorem.Stop
[Note3()]Note3
  note The standard decomposition in quantum-turbulence studies <cit.> involves kinetic energy density calculated semi-classically as a velocity power spectrum via Plancherel's theorem (Appendix. <ref>), which discards quantum phase information.Stop
[Note4()]Note4
  note The usual definition does not factor out the spherical measure 4π k^2, here written explicitly for convenience in connecting with recent literature.Stop
[Reeves et al.(2014)Reeves, Billam, Anderson, and Bradley]reeves_signatures_2014
  author author M. T. Reeves, author T. P. Billam, author B. P. Anderson, and author A. S. Bradley, title title Signatures of coherent vortex structures in a disordered two-dimensional quantum fluid, https://doi.org/10.1103/PhysRevA.89.053631 journal journal Physical Review A volume 89, pages 053631 (year 2014)NoStop
[Note5()]Note5
  note We do not consider the cross terms further in this work as they are not very intuitive to interpret. See Ref. <cit.> for a simple example.Stop
[Note6()]Note6
  note Their result translates to k_F=√(2)/ξ in our units. However, for this estimate of the forcing scale we see a narrower flat region in the log-corrected spectrum. For further discussion see Appendix <ref>.Stop
[Zhao et al.(2024)Zhao, Tao, and Spielman]zhao_kolmogorov_2024
  author author M. Zhao, author J. Tao, and author I. Spielman, https://doi.org/10.48550/arXiv.2408.04715 title Kolmogorov turbulence in atomic Bose-Einstein condensates (year 2024), https://arxiv.org/abs/2408.04715 arXiv:2408.04715 NoStop






§ NUMERICAL METHOD

Evaluation of the potential term in the GPE, as well as the interaction term |ψ(r,t)|^2, is straightforward using element-wise operations. However, the kinetic term contains a Laplacian, which is less straightforward. By applying a Fourier-transform to the wavefunction, we shift to momentum-space, where the kinetic operator becomes diagonal. This method is commonly known as pseudo-spectral differentiation and makes use of the Fourier-relation

    ∇^2 ψ (r,t) ↔ -k^2 ϕ(k,t),

where the arrows denote a spatial Fourier-transform, and ϕ(k,t) is the Fourier-transform of ψ(r,t). This method is efficient and offers excellent error properties compared to alternatives; however, it does come with some limitations. One limitation is that the boundary conditions imposed on the numerical grid become necessarily periodic. Another issue is that aliasing can occur when representing certain momentum values, based on the grid's density. Several tests were performed to ensure these factors were appropriately mitigated. 

We begin by ensuring that a sufficiently fine grid is being used to represent the system and that numerical error is within an acceptable margin. We test this only for the most violent and energetic case (U_F = 5μ). The error for smaller forcing amplitudes is expected to be approximately bounded by this case. Excitation was simulated across a range of step-sizes and grid points, with no damping potential applied. In the absence of damping, normalisation of the wavefunction should be conserved, allowing any deviation to be attributed to numerical error. By monitoring the normalisation over time, we can quantify the error for each case. 

As we would expect, error_tests(a) shows the error per time-step decreasing with step-size, until reaching a floor determined by the limit of Float64 precision. Reducing step-size beyond this point increases the total error [error_tests(b)], as more steps are required per second of dynamics. For this reason a step-size of Δ t = 10^-3 was chosen for calculations. 

After selecting the step-size, we determine the necessary grid density to ensure convergence. Similar to step-size, an overly dense grid is not necessarily optimal. Although it does not increase error, it significantly raises memory-usage and computation time. The number of points along each axis of the grid, M, was chosen as powers of two to ensure optimal performance under the FFTW algorithm. The total number of grid points is then M^3. To test for convergence, we performed identical excitation of a ground state using various grid sizes. We then calculate the total kinetic energy, compressible kinetic energy, and the wave-action spectrum for each case. 
    
In converge_test, results for M = 64 show significant deviations compared to results from denser grids. This suggests that the grid-size is insufficiently dense to accurately represent the system. For the three largest grid sizes, converge_test(a) and (b) show minimal variation. However, converge_test(c) reveals a notable difference between grid sizes, as the 256^3 and 128^3 grid sizes exhibit ultraviolet features that appear absent in the 512^3 grid. This feature is attributed to the Nyquist-frequency of Fourier-transforms.


Just as temporal Fourier-transforms have a maximum resolvable frequency determined by the sampling frequency, a spatial Fourier-transform has maximum resolvable k-values based on the spatial sampling frequency for each dimension. For a grid with M^3 points and lengths L_i, we have

    [ Δ x, Δ y, Δ z ] = [ L_x/M, L_y/M, L_z/M]


    k_Nyq, i = 1/22π/Δ i,    i ∈{x, y, z}

Occupation of k-values beyond the Nyquist limit (Eq. <ref>), can result in aliasing, where the condensate is misrepresented in Fourier-space. This results in the observed artifacts. Grids with more points provide higher spatial sampling, which shifts the ultraviolet artifact to higher k-values. 
    
For our simulations, we opted to use a grid size of M = 256. This allows a sufficient gap between the dissipation k-scale and the Nyquist k-scale. It also allows for efficient computation, with each second of Float64 dynamics taking approximately 38 minutes. Taking full advantage of this, we were able to test a high number of U_F values [excite_energies(e)]. 


§ SPECTRAL ANALYSIS



 §.§ Calculation of spectral densities using Plancherel's theorem

Expressing kinetic energy in terms of our density-weighted velocity field, we can utilize Plancherel's theorem to equivalently express the integral in k-space rather than position-space

    E_kin   = m/2∫ dr  |u(r)|^2 = m/2∫ dk  |u(k)|^2,

where u(k) is the three-dimensional Fourier-transform of u(r). Integrating the right-hand side of Eq. <ref> over the solid angle in k-space, we obtain a distribution which integrates over k to give the kinetic energy.

    E_kin    = ∫^∞_0 dk  m k^2/2∫^2π_0 dϕ_k ∫^π_0 dθ_k  sin(θ) |u(k, ϕ_k, θ_k)|^2 
        = ∫^∞_0 dk  ε_kin(k)

There is then a temptation to describe ε_kin(k) as the kinetic energy spectral density. However, this is somewhat misleading, as the components of ε_kin(k) do not sum locally in k-space, limiting it's interpretation as an energy density.

    ε_kin(k) ε^i_kin(k) + ε^c_kin(k) + ε^q_kin(k)



Using ε_kin(k) also implicitly introduces a semi-classical approximation, as quantum phase information is neglected when calculating |u|^2 in Eq. <ref>. We can more appropriately describe ε_kin(k) as a velocity power spectral density. For wavefunctions represented on a cartesian grid (typically the case in numerical studies), there will also be fewer points with small |k| than large, leading to a lack of resolution in the infrared region of the resulting spectra. Conversely, the k-values for the reformulated spectral-density (Eq. <ref>) are completely decoupled from the grid representing the wavefunction, allowing for arbitrarily high resolution in k-space. While k-space resolution can be arbitrarily high, the available range is still limited by Nyquist k-values. 



 §.§ Comparison of energy densities and velocity power spectra

To demonstrate the importance of including phase information, spectra_comparison compares equivalent energy density and velocity power spectra components. Both are calculated identically, in except that ε^α_kin(k) use the density-weighted velocity fields, and e^α_kin(k) use the w^α fields which incorporate quantum phase, defined in Eq. <ref>. 

In spectra_comparison, we see ε_kin(k) and e_kin(k) can have close agreement [spectra_comparison(b)] or be extremely distinct [spectra_comparison(a)] for different cases. Generally, the two distributions are similar in the inertial range for well-developed turbulence. During initial excitation and at extreme lengthscales, there is no clear pattern relating the two spectra. This indicates that for well-developed turbulence there is nontrivial phase-behaviour in the inertial range, such that neglecting quantum phase and imposing a semi-classical approximation does not introduce significant error to energy spectral densities. An interesting question then would be whether this is true for other forms of turbulence. 
        
spectra_comparison demonstrates that inclusion of phase information is clearly important in certain cases. However, this information is typically inaccessible in experimental studies of BECs. Therefore, it would be valuable to understand the relationship between ε_kin(k) and e_kin(k), and identify conditions under which they are similar enough that phase information becomes less critical. 
      


§ COMPARISON WITH LOGARITHMIC-CORRECTION MODEL

Another explanation for the k^-3.5 power law, is that the k∝ k^-3 stationary solution is marginally non-local, requiring a log-correction of the form

    n(k)    = C_d P_0^1/3 k^-3ln^-1/3(k/k_F).

The first two terms are constants and k_F is the wavenumber associated with the forcing length scale, where energy is injected. 
 

This logarithmic-correction emerges from the introduction of an infrared cutoff at the forcing wavenumber, k_F. In Ref. <cit.>, k_F is a defined parameter. With our method of energy injection this is not the case, and k_F is slightly more ambiguous. We excite the condensate along the z-direction, so one might expect k_F = k_L_z. However, in practice we find the forcing peak appears at slightly higher k-values than k_L_z [k3.5], implying the forcing length scale is slightly smaller than L_z. Regardless, we test the models fit using a wide range of k_F values.
  
Applying the logarithmic-correction to spectra from the k^-3.5 regime, we generally find poor agreement using k_F values near the actual forcing lengthscale. However, for k_F = 1/ξ, we find good agreement in the total momentum distribution and slightly stronger agreement in the compressible & quantum pressure components [logplot]. This k_F value is clearly distinct from the actual forcing wavenumber, supporting the argument presented in Ref. <cit.>, that the physical energy-injection lengthscale is distinct from the effective injection scale for the four-wave cascade. The strength of this agreement, however, is limited by the size of the inertial range, which to create a more compelling case could be increased with a greater dissipation lengthscale and a finer numerical grid to ensure the model remains stable.

The slightly improved fit of the model in n^q(k) and n^c(k), as opposed to the total momentum distribution, is likely due to incompressible excitations irrelevant to the WWT cascade being discarded. In Ref. <cit.> there is strong agreement with no partitioning of the spectrum, possibly because fewer vortices are nucleated. 

Where we use a finite box-trap with hard-wall boundaries, Ref. <cit.> uses a homogeneous condensate with periodic boundaries. Our forcing term is an emulation of experimental forcing, whereas in Ref. <cit.>, energy is injected in k-space in a narrow band about at particular k_F value. Our model also does not feature any mechanism for dissipation at large lengthscales. These factors may all play a role in introducing greater incompressible aspects to the steady turbulent state. 
