Kolmogorov and Kelvin wave cascades in a generalized model for quantum turbulence
KKolmogorov and Kelvin wave cascades in a generalized model for quantum turbulence
Nicolás P. Müller and Giorgio Krstulovic
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange,Boulevard de l’Observatoire CS 34229 - F 06304 NICE Cedex 4, France
We performed numerical simulations of decaying quantum turbulence by using a generalizedGross-Pitaevskii equation, that includes a beyond mean field correction and a nonlocal interactionpotential. The nonlocal potential is chosen in order to mimic He II by introducing a roton minimumin the excitation spectrum. We observe that at large scales the statistical behavior of the flowis independent of the interaction potential, but at scales smaller than the intervortex distance aKelvin wave cascade is enhanced in the generalized model. In this range, the incompressible kineticenergy spectrum obeys the weak wave turbulence prediction for Kelvin wave cascade not only forthe scaling with wave numbers but also for the energy fluxes and the intervortex distance.
I. INTRODUCTION
One of the most fundamental phase transitions in lowtemperature physics is the Bose-Einstein condensation .It occurs when a fluid composed of bosons is cooled downbelow a critical temperature. In that state, the systemhas long-range order and can be described by a macro-scopic wave function. One of the most remarkable prop-erties of a Bose-Einstein condensate (BEC) is that itflows with no viscosity. Well before the first experimen-tal realization of a BEC by Anderson et al. , Kaptizaand Allen discovered that helium becomes superfluid be-low 2.17K . A couple of years later, London suggestedthat superfluidity is intimately linked to the phenomenonof Bose-Einstein condensation . Since then, superfluidhelium and BECs made of atomic gases have been ex-tensively studied, both theoretical and experimentally.In particular, the fluid dynamics aspect of quantum flu-ids has renewed interest due the impressive experimentalprogress of the last fifteen years. Today it is possibleto visualize and follow the dynamics of quantum vor-tices, one the most fundamental excitations of a quantumfluid .Quantum vortices are topological defects of the macro-scopic wave function describing the superfluid. They arenodal lines of the wave function and they manifest pointsand filaments in two and three dimensions respectively.To ensure the monodromy of the wave function, vorticeshave the topological constraint that the circulation (con-tour integral) of the flow around the vortex must be amultiple of the Feynman-Onsager quantum of circulation h/m , where h is the Planck constant and m is the massof the Bosons constituting the fluid . In superfluid he-lium their core size is of the order of 1Å whereas in atomicBECs is typically of the order of microns . Quantum vor-tices interact with other vortices similarly to the classicalones. They move thanks to their self-induced velocityand interact with each other by hydrodynamics laws .Unlike ideal classical vortices described by Euler equa-tions, quantum vortices can reconnect and change theirtopology despite the lack of viscosity of the fluid in whichthey are immersed .At scales much larger than the mean intervortex dis- tance (cid:96) , the quantum nature of vortices is not very impor-tant as many individual vortices contribute to the flow.One could expect then that flow is similar, in some sense,to classical one. Indeed, if energy is injected at largescales a classical Kolmogorov turbulent regime emerges.Such a regime has been observed numerically andexperimentally . In a three-dimensional turbulentflow, energy is transferred towards small scales in a cas-cade process . In a low temperature turbulent super-fluid, when energy reaches the intervortex distance, en-ergy keeps being transferred to even smaller scales whereit can be efficiently dissipated by sound emission. Themechanisms responsible for this are the vortex reconnec-tions and the wave turbulence cascade of Kelvin waves,that have its origin in the quantum nature of vortices .Describing a turbulent superfluid is not an easy task, inparticular for superfluid helium. One of the main reasonsis the gigantesque scale separation existing between thevortex core size and the typical size of experiments, cur-rently of the order centimeters or even meters . Theirtheoretical description began at the beginning of the 20thcentury by the pioneering works of Landau and Tiszawhere superfluid helium was modeled by two immisciblefluid components . In this two-fluid model, the ther-mal excitations constitute the so called normal fluid thatis modeled through the Navier-Stokes equations whereasa superfluid component is treated as an inviscid fluid.It was later realized that the thermal excitations inter-act with superfluid vortices through a scattering processthat leads to a coupling of both component by mutualfriction forces . Today the two-fluid description, knownas the Hall-Vinen-Bekarevich-Khalatnikov model is un-derstood as a coarse-grained model where scales smallerthan the intervortex distance are not considered. Thequantum nature aspects of superfluid vortices are there-fore lost. However, this model remains useful for de-scribing the large scale dynamics of finite temperaturesuperfluid helium. An alternative model was introducedby Schwarz , where vortices are described by vortex fil-aments interacting through regularized Biot-Savart inte-grals. However, the reconnection process between linesneeds to be modeled in an ad-hoc manner and by con-struction the model excludes the dynamics of a superfluid a r X i v : . [ c ond - m a t . o t h e r] S e p at scales smaller than the vortex core size. Finally, in thelimit of low temperature and weakly interacting BECs, amodel of different nature can be formally derived whichis the Gross-Pitaevskii (GP) equation, obtained from amean field theory . This model naturally contains vor-tex reconnections , sound emission and is knownto also exhibit a Kolmogorov turbulent regime at scalesmuch larger than the intervortex distance . Althoughthis model is expected to provide some qualitative de-scription of superfluid helium at low temperatures, itlacks of several physical ingredients. For instance, in GP,density excitations do not present any roton minimum asit does in superfluid helium, where interaction betweenboson are known to be much stronger than in GP .However, there have been some successful attempts toinclude such effect in the GP model. For instance, a ro-ton minimum can be easily introduced in GP by using anonlocal potential that models a long-range interactionbetween bosons . The stronger interaction of heliumcan also be included phenomenologically by introducinghigh-order terms in the GP Hamiltonian. Note that suchhigh order terms can be derived as beyond mean fieldcorrections . Some generalized version of the GP modelhas been used to study the vortex solutions andsome dynamical aspects such as vortex reconnections .Naively, for a turbulent superfluid, we can expect thatsuch generalization of the GP model might be importantat scales smaller than the intervortex distance and withless influence at scales at which Kolmogorov turbulenceis observed.In this work, we study quantum turbulent flows byperforming simulations of a generalized Gross-Pitaevskii(gGP) equation. We compare the effect of high-ordernonlinear terms and the effect of a nonlocal interactionpotential in the development and decay of turbulence atscales both, larger and smaller than the intervortex dis-tance. Remarkably, by modeling superfluid helium witha nonlocal interaction potential and including high-orderterms, the range where a Kelvin wave cascade is observedis extended and becomes manifest. Using the dissipation(or rate of transfer) of incompressible kinetic energy weare able to show that the weak wave turbulence results are valid not only to predict the scaling with wave num-ber but also with the energy flux.The manuscript is organized as follows. Section II in-troduces the gGP model and discusses its basic propertiesand solutions. It also discusses how the vortex profile ismodified in this generalized model. All useful definitionsto study turbulence are also given here. Section III givesa brief overview of the predictions of quantum turbulenceand the numerical methods used in this work. Also, itincludes the results of different simulations at moderateand high resolutions by varying the different parametersof the beyond mean field correction and the introductionof a nonlocal potential. Finally in Section IV we presentour conclusions. II. THEORETICAL DESCRIPTION OFSUPERFLUID TURBULENCE
On this section we introduce the generalized Gross-Pitaevskii model used in this work. We also discuss andreview some of the basic properties of the model as itselementary excitations and its hydrodynamic description.
A. Model
The Gross-Pitaevskii equation describes the low tem-perature dynamics of weakly interacting bosons of mass m i ¯ h ∂ψ∂t = − ¯ h m ∇ ψ − µψ + g | ψ | ψ, (1)where ψ is the condensate wave function, µ the chemi-cal potential and g = 4 π ¯ h a s /m is the coupling constantfixed by the s -wave scattering length a s that models a lo-cal interaction between bosons. Note that, the use of a lo-cal potential assumes a weak interaction between bosons,which certainly is not the case for other systems like HeII and for dipolar gases .A generalized model that is able to describe more com-plex systems can be obtained by considering a nonlocalinteraction between bosons. With proper modeling ,density excitations exhibit a roton minimum in theirspectrum as the one observed in He II . It also describeswell the behavior of dipolar condensates . In heliumand other superfluids, the interaction between bosonsis stronger and high order nonlinearities are needed forproper modeling. For instance, in helium high-orderterms are considered to mimic its equations of state andin dipolar BECs beyond mean field terms are needed todescribe the physics of recent supersolid experiments .We consider the generalized Gross-Pitaevskii (gGP)model written as i ¯ h ∂ψ∂t = − ¯ h m ∇ ψ − µ (1 + χ ) ψ (2) + g (cid:18)(cid:90) V I ( x − y ) | ψ ( y ) | d y (cid:19) ψ + gχ | ψ | γ ) n γ ψ. where γ and χ are two dimensionless parameters thatdetermine the order and amplitude of the high-orderterms. The interaction potential V I is normalized suchthat (cid:82) V I ( x )d x = 1 . The chemical potential and theinteraction coefficient of the high-order term have beenrenormalized such | ψ | = n = µ/g is the density ofparticles for the ground state of the system for all valuesof parameters. The GP equation is recovered by simplysetting V I ( x − y ) = δ ( x − y ) and χ = 0 .The gGP equation is not intended to be a first principlemodel of superfluid helium, but it has the advantage ofat least introducing in a phenomenological manner someimportant physical aspects of helium. B. Density waves
The dispersion relation of the GP model is easily ob-tained by linearizing equation (1) about the ground state.The Bogoliubov dispersion reads ω B ( k ) = c k (cid:114) ξ k , (3)where k is the wave number, c = (cid:112) gn /m is the speedof sound of the superfluid and ξ = ¯ h/ √ mgn is thehealing length at which dispersive effects become impor-tant. The healing length also fixes the vortex core size.A similar calculation leads to the Bogoliubov disper-sion relation in the case of the gGP model (2) ω ( k ) = ck (cid:115) ξ k V I ( k ) + χ ( γ + 1)1 + χ ( γ + 1) , (4)where ˆ V I ( k ) = (cid:82) e i k · r V I ( r )d r is the Fourier transformof the interaction potential normalized such that ˆ V I ( k = 0) = 1 . The inclusion of beyond mean-field terms anda nonlocal potential yields to a renormalized speed ofsound and healing length. They are given in terms of c and ξ by c = c (cid:112) χ ( γ + 1) (5) ξ = ξ (cid:112) χ ( γ + 1) . (6)Note that, in what concerns low amplitude density waves,the effect of high-order terms is a simple renormalizationof the healing length and the speed of sound. Dependingon the shape and properties of the nonlocal potential, thedynamics and steady solutions can be drastically modi-fied. Note that the product between c and ξ remainsconstant because it is related to the quantum of circula-tion κ = h/m = cξ π √ c ξ π √ .In order to be able to compare the systems with dif-ferent type of interactions, it is convenient to rewrite Eq.(2) in terms of its intrinsic length ξ and speed of sound c and the bulk density n . The gGP model then becomes ∂ t ψ = − i cξ √ χ ( γ + 1)) (cid:20) − (1 + χ ( γ + 1)) ξ ∇ ψ − (1 + χ ) ψ + χ | ψ | γ ) n γ ψ + ψn (cid:90) V I ( x − y ) | ψ ( y ) | d y (cid:21) . (7)In numerical simulations we will express lengths in unitof the healing length ξ . A natural time scale to study ex-citations is the fast turnover time τ = ξ/c . However, thissmall-scale based time is not appropriate for turbulentflows. For such flows, it is customary to use the large-eddy turnover time corresponding to the typical time ofthe largest coherent vortex structure and will be definedlater.
1. Modeling superfluid helium excitations
In this work, we aim at mimicking some properties ofsuperfluid helium II, in particular its roton minimum inthe dispersion relation. For the sake of simplicity, weuse an isotropic nonlocal interaction potential used inprevious works . With our normalization it reads ˆ V I ( k ) = (cid:34) − V (cid:18) kk rot (cid:19) + V (cid:18) kk rot (cid:19) (cid:35) exp (cid:18) − k k (cid:19) , (8)where k rot is the wave number associated with the rotonminimum and V ≥ and V ≥ are dimensionless pa-rameters to be adjusted to mimic the experimental dis-persion relation of helium II . The effects of differentfunctional forms of the nonlocal potential have been stud-ied in previous works, showing that only a phase-shift of ψ and the overall amplitude of the density depend on theprecise form of the interaction .In order to compare the dispersion relation (4) withthe experimental data , we plot the helium dispersionrelation in units of the helium healing length ξ He = 0 . Åand its turnover time τ He = ξ He /c He = 3 . × − s,where c He = 238 m/s is the speed of sound in He II. Themeasured helium dispersion relation is displayed in Fig.1 as dotted red lines.It was reported in Reneuve et al. that introducinga roton minimum in the GP dispersion relation thatmatches helium measurements leads to an unphysicalcrystallization under dynamical evolution of a vortex.We confirm such behavior in our simulations. In orderto avoid such spurious effect of the model, in reference the frequency associated to the roton minimum is set tohigher values to be able to study vortex reconnections.We have numerically observed that this crystallizationtakes place even for values of χ = 0 . and γ = 1 ofthe beyond mean field expansion, which correspond tothe correction of first order. For this reason, we chose ahigher order expansion of γ = 2 . for the simulations ofa nonlocal potential, value that was already used in theliterature to study the vortex density profile in superfluidhelium .The dispersion relation of a nonlinear wave system canbe measured numerically by computing the spatiotempo-ral spectrum of the wave field . As an example, in Fig. Figure 1. Spatiotemporal dispersion relation for simulationswith grid points with a nonlocal potential and beyondmean field corrections. Light zones correspond to excited fre-quencies. Figures (a) and (b) correspond to different ampli-tude of the perturbation A , both exhibiting a roton minimum.Experimental observations (red dotted line, see ) and theo-retical dispersion following equation (4) (blue dashed line) areshown. collocation points and with pa-rameters set to γ = 2 . , χ = 0 . , V = 4 . , V = 0 . and k rot ξ = 1 . (see details on numerics later in Sec.III), for two different amplitude values A . Dark zones in-dicate that no frequencies are excitated, while light zonescorrespond to the excitated ones with the total sum nor-malized to one. The parameters have been set in orderto qualitatively match the measured helium dispersionrelation. As expected for weak amplitude waves, the nu-merical and theoretical dispersion relations coincide. Forlarger wave amplitudes, theoretical prediction (4) and nu-merical measurements slightly differ together with an ap-parent broadening of the curve. This is a typical behaviorof nonlinear wave systems . In the following sections,all simulations with a nonlocal interaction are performedwith the previous set of parameters. C. Hydrodynamic description
The GP equation maps into an hydrodynamic descrip-tion by introducing the Madelung transformation ψ = (cid:112) ρ/m exp (cid:18) iφ √ cξ (cid:19) , (9)which allows the mapping of the wave function with thefluid mass density ρ = m | ψ | and with the fluid velocity v = ∇ φ . Replacing equation (9) into the gGP model (7) two hydrodynamic equations are obtained ∂ρ∂t + ∇ · ( ρ v ) = 0 (10) ∂φ∂t + 12 ( ∇ φ ) = − h [ ρ ] + ( cξ ) ∇ √ ρ √ ρ (11) h [ ρ ] = − c (1 + χ ) + c V I ∗ ρρ + c χ (cid:18) ρρ (cid:19) γ +1 , (12)where ∗ denotes is the convolution product and ρ = m | ψ | is the fluid mass density of the ground state.Those equations correspond to the continuity andBernoulli equations of a fluid with an enthalpy per unit ofmass h [ ρ ] , respectively. The last term of equation (11)is called the quantum pressure. Note that hydrodynamicpressure is given by p [ ρ ] = c ρρ (cid:20) V I ∗ ρ + χ γ + 1 γ + 2 ρ γ +1 ρ γ (cid:21) . (13)As expected, for large amplitude waves, the speed ofsound reads ∂p∂ρ (cid:12)(cid:12)(cid:12) ρ = c (1 + χ ( γ + 1)) = c .Although the fluid is potential, it admits vortices astopological defects of the wave function. A stationaryvortex solution of (7) is a zero of the wave functionwhere the circulation around it is quantized with values Γ = ± sκ , being the quantum of circulation κ = 2 π √ cξ and with s an integer. Because of this last condition,topological defects are also called quantum vortices.A quantum vortex has vortex core size of order ξ anddepends on the parameters of the gGP model. By re-placing the Madelung transformation (9) into the gGPequation (7) and solving in cylindrical coordinates, a dif-ferential equation for the vortex profile is directly ob-tained r dd r (cid:18) r d R d r (cid:19) + (cid:26) − s ξ r − V I ∗ R ++ χ (1 − R γ +2 ) (cid:27) Rξ = 0 (14)where R ( r ) = (cid:112) ρ ( r ) /ρ defines the density profile of thevortex line in the radial direction r .Figure 2 (a) displays the mass density of a two-dimensional vortex in the case where the nonlocal inter-action potential is included. The roton minimum intro-duces some density fluctuations around the center of thevortex which is a well-known pattern. Such a behaviorhas already been studied before, for example its interac-tion with an obstacle , the dynamics of vortex rings and in reconnection processes . Figure 2 (b) shows theradial dependence of the density profile of a vortex fordifferent parameters of the gGP model. Numerical simu-lations were performed with grid points with stan-dard numerical methods (see section III B for details).Even though all curves tend to collapse when plotted asa function of the healing length ξ , the vortex core size x/ξ y / ξ (a) r/ξ . . . . . ρ / ρ (b) χ = χ = γ = χ = γ = χ = γ = χ = γ = χ = γ = χ = γ = . . . . . Figure 2. (a) Mass density of a two-dimensional vortex witha nonlocal potential. (b) Density profile of a vortex for thegGP model with different values of the nonlinearity and fora nonlocal potential (yellow line). The vortex core size tendsto increase with the nonlinearity. slightly increases (in units of ξ ) when the nonlinearity ofthe system is increased. Note that for the present rangeof parameters, ξ /ξ varies in the range (1 , . . The rel-atively good collapse of the vortex core size thus justifiesthe choice of ξ to parametrize the gGP model while vary-ing the beyond mean field parameters. D. Energy decomposition and helicity insuperfluids
It is convenient to write the free energy per unit of mass F of a quantum fluid such that it vanishes when evalu-ated in the ground state of the system ( ψ = (cid:112) ρ /m = √ n ). For the gGP model in equation (7), it is given by F = c n V (cid:90) (cid:34) ξ |∇ ψ | + | ψ | n ( V I ∗ | ψ | ) − (1 + χ ) | ψ | + χ | ψ | γ +2) n γ +10 ( γ + 2) + n χn γ + 1 γ + 2 (cid:35) d r, (15)with V the volume of the fluid. Following stan-dard procedures applied in simulations of GP quantumturbulence , the free energy can be decomposed as F = E I kin + E C kin + E q + E int where E Ikin = V ρ (cid:82) (cid:0) [ √ ρ v ] I (cid:1) d r , E Ckin = V ρ (cid:82) (cid:0) [ √ ρ v ] C (cid:1) d r and E q = c ξ V ρ (cid:82) | ∇ √ ρ | d r ,with [ √ ρ v ] I the regularized incompressible velocity ob-tained via the Helmholtz decomposition and [ √ ρ v ] C = √ ρ v − [ √ ρ v ] I the compressible one. The internal energyper unit of volume is defined in the gGP model as E int = c V ρ (cid:90) (cid:20) ρ ( ρ − ρ ) V I ∗ ( ρ − ρ )+ (cid:18) ρρ (cid:19) γ +1 χργ + 2 − χρ + χ γ + 1 γ + 2 ρ (cid:21) d r. (16)Note that E int = 0 if ρ = ρ . The corresponding en-ergy spectra are defined in a straightforward way for thequadratic quantities . For the internal energy spectrum, it is defined as follows E int ( k ) = c V ρ (cid:90) (cid:34) ρ ( ρ − ρ ) (cid:86) − k ˆ V I ( k )( ρ − ρ ) (cid:86) k + χγ + 2 ˆ ρ − k (cid:16) ρρ (cid:17) γ +1 k (cid:86) + ( χρ γ +1 γ +2 − χρ ) (cid:86) k (cid:35) dΩ k , (17)where dΩ k is the element of surface of the shell | k | = k .Note that this particular choice of the spectrum is notunique and has been made so that the ground state ρ = ρ contributes with no internal energy to the system. Itis also worth noting that with this definition, the internalenergy spectrum may take negative values.Besides the energies, there is another quantity in quan-tum turbulence that presents a great interest in the dy-namics of quantum vortices , which is the central linehelicity per unit of volume H c = 1 V (cid:90) v ( r ) · ω ( r )d r. (18)Note that V H c /κ is the total number of helicity quanta.Formally, this quantity is ill defined for a quantum vortexas the vorticity is δ -supported on the filaments and ve-locity is not defined on the vortex core. However, in theGP formalism, this singularity can be removed by takingproper limits . We use the definition central line helicityproposed in reference as its numerical implementationis tedious but straightforward and well behaved for vor-tex tangles. III. EVOLUTION OF QUANTUM TURBULENTFLOWS
This section gives a brief overview about the predic-tions in quantum turbulence both at large and smallscales and details of the numerical methods used to runthe simulations. There is also a description of the flowvisualization in the presence of a nonlocal interaction po-tential, and the results of the flow evolution at moderateand high resolution are shown. In particular, it is studiedthe dependence of the different components of the energyand the helicity with beyond mean field parameters andwith the introduction of a nonlocal interaction potential.
A. A brief overview of cascades in quantumturbulence
Quantum turbulence is characterized by the disorderedand chaotic motion of a superfluid. Energy injected, orinitially contained, at large scales is transferred towardssmall scales in a Richardson cascade process . In thecontext of GP turbulence, the contribution of vortices tothe global energy can be studied by looking at the in-compressible kinetic energy E I kin and its associated spec-trum. As the system evolves, vortices interact transfer-ring energy between scales. Besides, the incompressiblekinetic energy is transferred to the quantum, internal andcompressible energy through vortex reconnections andsound emission . After some time, acoustic excita-tions thermalize and act as a thermal bath providing a(pseudo) dissipative mechanism, so vortices shrink untilthey vanish .Three-dimensional quantum turbulence presents twomain statistical properties. At scales much larger thanthe intervortex distance (cid:96) , but much smaller than the in-tegral scale L , the quantum character of vortices is notimportant and we can think as the system being coarse-grained. At such scales the system presents a behaviorthat resembles to classical turbulence with a direct en-ergy cascade, that is the transfer of energy from large tosmall structures. As a consequence, in this range, the in-compressible kinetic energy spectrum E I kin ( k ) follows theKolmogorov prediction E I kin ( k ) = C K (cid:15) / k − / , (19)where C K ∼ and (cid:15) is the dissipation rate of theflow, which in GP quantum turbulence is associatedwith the rate of change of incompressible kinetic en-ergy (cid:15) = − d E I kin / d t that is expressed in units of [ (cid:15) ] = Length /T ime .In classical three-dimensional inviscid flows, helicity(18) is also conserved. Associated to this invariant, asecond direct cascade is expected to be also present atlarge scales, obeying the scaling H ( k ) = C H η(cid:15) − / k − / , (20)where C H ∼ and η = − d H/ d t is the dissipation rateof helicity. This dual cascade, has been also observed inquantum turbulent flows described by the GP equation .At scales smaller than the intervortex distance, eachquantum vortex can be thought as if it were isolated.Hence, Its behavior can be described by the wave turbu-lence theory as such vortices admit hydrodynamic ex-citations known as Kelvin waves. Such waves propa-gate along vortices and interact nonlinearly among them-selves. As a result, energy is transferred towards smallscales through a process that can be described by thetheory of weak wave turbulence . An agitated debatearose concerning the prediction of the energy spectrum.Two independent groups leaded by L’vov & Nazarenko and Kozik & Svitsunov , starting from the same equa-tions and applying the same theory derived differentpredictions. Even though, today there is more numer-ical data supporting L’vov & Nazarenko prediction ,this issue is still debated . We present here theL’vov&Nazarenko prediction as, we will see later, it wasfound to be in agreement with our numerical data. Thistheoretical prediction is derived for an almost straightvortex of period L v and, as discussed in , some care is needed in order to apply the model to a turbulent vortextangle. We partially reproduce here and adapt to ourcase the considerations of . The wave turbulence L’vov& Nazarenko prediction is e KW ( k ) = C LN κ Λ (cid:15) / Ψ / k / , (21)with Λ = log( (cid:96)/ξ ) and C LN ≈ . . Here (cid:15) KW = − d e KW / d t is the mean energy flux per unit of length L v and density ρ . Note their respective dimen-sions are [ (cid:15) KW ] = Length /T ime and [ e KW ( k )] = Length /T ime . The dimensionless number Ψ is givenby Ψ = (12 πC LN ) / (cid:15) KW1 / κ / k / = C / ˜Ψ , (22)where k min is the smallest wave number of the Kelvinwaves, that can be associated with the wave number ofthe intervortex distance k (cid:96) = 2 π/(cid:96) in the case of a vortextangle . ˜Ψ is defined so that it is independent of theconstant C LN and proportional to Ψ .In order to compare this result with the incompress-ible kinetic energy, one can notice that the total energy ofKelvin waves is L v ρ (cid:82) e KW ( k )d k , where now L v is takenas the total vortex length in the system. As in a turbu-lent tangle the total vortex length is related to the meanintervortex distance by L v = V (cid:96) − , it follows that themean kinetic energy spectrum per unit of mass is givenby E KW ( k ) = e KW ( k ) (cid:96) − . The same logics, relates theenergy flux (cid:15) KW of the Kelvin wave cascade to the globalenergy flux (cid:15) of a tangle by (cid:15) KW = (cid:15)(cid:96) . It follows from(21) and the previous considerations that E KW ( k ) = C / κ Λ (cid:15) / (cid:96) − / ˜Ψ / k / . (23)Here we have made the assumption that the energy fluxis the same in the Kolmogorov range than in the Kelvinwave cascade. This strong assumption might be ques-tioned as energy could be already dissipated into soundby vortex reconnections at different scales diminishingthis value . Such extra sinks of energy are difficult toquantify and we will not take them into account. Finally,note that the theory of wave turbulence also predicts thevalue of the constant C LN , however in (23) several phe-nomenological considerations have been made and we donot expect an exact agreement. Nevertheless, the scalingwith the global energy flux should remain valid. B. Numerical methods
We perform numerical simulations of equation (7). Apseudo-spectral method is used for the spatial resolutionapplying the “ / rule” for dealiasing , and the nonlinearterm is dealiased twice following the scheme presentedin in order to also conserve momentum. Note in thecase of a nonlocal potential, this extra step has no extranumerical cost. A Runge-Kutta method of fourth orderis used for time stepping. All simulations were performedin a cubic L -periodic domain.To observe a Kolmogorov range in GP turbulence itis customary to start from an initial vortex configura-tion with a minimal acoustic contribution. The initialcondition for the wave function is obtained by a mini-mization process such that the resulting flow is as closeas possible to the targeted velocity field . In this workwe study the quantum Arnold-Bertrami-Childress (ABC)flow introduced in . It is obtained from the velocity field v ABC = v ( k )ABC + v ( k )ABC , where each ABC flow is given by v ( k )ABC = [ B cos( ky ) + C sin( kz )]ˆ x + (24) [ C cos( kz ) + A sin( kx )]ˆ y + [ A cos( kx ) + B sin( ky )]ˆ z. We set in this work ( A, B, C ) = V amp (0 . , , . / √ ,with V amp = 0 . c . Each ABC flow is a L -periodicstationary solution of the Euler equation with maxi-mal helicity, in the sense that ∇ × v ( k )ABC = k v ( k )ABC .The mean kinetic energy of the v ABC flow is E ABCkin = V ( A + B + C ) = 0 . c . Following , the wavefunction associated to this ABC flow is generated as ψ ABC = ψ ( k )ABC × ψ ( k )ABC , where each mode is constructedas the product ψ ( k )ABC = ψ x,y,zA,k × ψ y,z,xB,k × ψ z,x,yC,k with ψ x,y,zA,k = exp (cid:26) i (cid:20) A sin( kx ) cξ √ (cid:21) πyL + i (cid:20) A cos( kx ) cξ √ (cid:21) πzL (cid:27) (25)where the brackets [ ] indicate the integer closest to thevalue to ensure periodicity. This ansatz gives a goodapproximation for the phase of the initial condition. Inorder to set properly the mass density and the vortex pro-files, it is necessary to first evolve ψ ABC using the gener-alized Advected Real Ginzburg-Landau equation (imag-inary time evolution in a locally Galilean transformedsystem of reference) ∂ t ψ = − c ξ √ (cid:26) − ξ ∇ ψ − (1 + χ ) ψ + χ | ψ | γ ) ρ γ ψ + ψρ ( V ∗ | ψ | ) (cid:27) − i v ABC · ∇ ψ − ( v ABC ) √ cξ ψ. (26)This equation is dissipative and its final state containsa minimal amount of compressible modes. This stateis used as initial condition for the gGP equation. Unlessstated otherwise, we use a flow at the largest scales of thesystems by setting k = 2 π/L and k = 4 π/L throughoutthis work.The simulations performed in this work are summa-rized in Table I and regrouped in two different sets. Thefirst set of simulations (runs A1 - A8) have been per-formed at a moderate spatial resolution of N = 256 grid points to study the effects introduced by the beyondmean field interactions and a nonlocal potential. Each ofthem has a different value of χ and γ with a local poten-tial and were compared with a single simulation with a N χ γ L/ξ ˜ k , ˜ k InteractionpotentialA1 256 0 1 171 1,2 localA2 256 1 1 171 1,2 localA3 256 3 1 171 1,2 localA4 256 5 1 171 1,2 localA5 256 1 2.8 171 1,2 localA6 256 3 2.8 171 1,2 localA7 256 5 2.8 171 1,2 localA8 256 0.1 2.8 171 1,2 nonlocalB1 512 0 1 341 1,2 localB2 512 0.1 2.8 171 1,2 nonlocalB3 512 0.1 2.8 341 1,2 nonlocalB4 512 0.1 2.8 341 2,3 nonlocalB5 512 0.1 2.8 341 3,4 nonlocalB6 1024 0.1 2.8 683 1,2 nonlocalTable I. Table with the parameters of the different simula-tions. N is the linear spatial resolution, χ and γ are the am-plitude and order of the beyond mean field interactions, L/ξ is the scale separation between the domain size L and thehealing length ξ , ˜ k = k L/ π and ˜ k = k L/ π are the twowave numbers where the energy is concentrated for the initialcondition, and a local or a nonlocal interaction potential isused in each of them. nonlocal interaction potential. The second set (runs B1- B6) has been performed to study the scaling of the en-ergy spectra. In these runs, we used a spatial resolutionof and grid points, different scale separationsand initial conditions. These results were also comparedwith the GP model. C. Flow visualization
The introduction of a nonlocal potential, as mentionedin Sec. II C, allows the system to reproduce the rotonminimum in the excitation spectrum (see Fig. 1). Asa consequence, the density profiles close to the quantumvortices have some fluctuations around the bulk value ρ (see Fig. 2). These oscillations have been studiedfor the profile of a two-dimensional vortex and havebeen also observed in three dimensions during vortexreconnections . In the case of a helical vortex tangle,the roton minimum induces a remarkable pattern of den-sity fluctuations around a vortex line. A visualizationof the initial condition ψ ABC for run B6 is displayed inFig. 3 (a)-(b). The red structures are isosurfaces of lowdensity values ρ = 0 . ρ and thus represent the vortexlines. The greenish rendering displays density fluctua-tion of the field above the bulk value ρ , that are onlyobserved in the case of a nonlocal potential. In Fig. 3(a) we recognize the large scale structures of the ABCflow accompanied by some density fluctuations aroundthe nodal lines. Figure 3 (b) displays a zoom of the tan-gle where such fluctuations are clearly observed. Unlikethe (local) GP model, density variations around a vor-tex line have a very specific pattern, rolling around the Figure 3. (a)-(b) Visualization of an ABC flow at t = 0 and (c)-(d) for t = 1 . τ L and t = 0 . τ L respectively, for a resolutionof grid points with a nonlocal potential. The isosurfaces of a small value the mass density shown in red correspond tothe vortex lines, and in green are the values of the density fluctuations above ρ . nodal lines in a helical manner. Such helical pattern is aconsequence of the maximal helicity initial condition pro-duced by the ABC flow. Indeed, we have also produceda Taylor-Green initial condition , that has no mean he-licity, and such helical density fluctuations are absent,although they are nevertheless developed after some vor-tex reconnections, as observed in (data not shown). Fi-nally, in Fig. 3 (c)-(d) we display visualizations of thefield for times t = 1 . τ L and t = 0 . τ L respectively.The first one corresponds to a time inside the time win-dow where averages are done, and the second one to aprevious time for a better insight of the flow. As thesystem evolves, acoustic emissions are produced and thedensity fluctuations increase. In Fig. 3 (c) we observe aturbulent tangle where a large scale structure is predom-inant. Figure 3 (d) displays a zoom where reconnectionsand Kelvin waves propagating along vortices are clearlyvisible. D. Temporal evolution of global quantities
In this section we study the behavior of the globalquantities of an ABC flow described by gGP model (7)with both local and nonlocal potentials corresponding toruns A in Table I.Figure 4 shows the time evolution of the (a) incom-pressible kinetic energy and (b) the sum of the quan-tum, internal and compressible kinetic components tothe total energy. Time is expressed in units of thelarge-eddy turnover time τ L = L /v rms with v rms = (cid:113) E I kin ( t = 0) / and L its integral length scale givenby L = 2 π/k with k the largest wave number usedto generate the initial condition. We notice that in Fig.4 (a) the values of amplitude and exponent of the be-yond mean field interaction and the inclusion of rotonminimum (Runs A1-A8) have a negligible impact on theincompressible energy of the initial condition, and theireffect is very small during the temporal evolution. Onthe other hand, as the fluid can be considered to be moreincompressible due to stronger interactions, the density t/τ L . . . . . E I k i n / E A B C k i n (a) χ = 0; γ = 1 χ = 1; γ = 1 χ = 3; γ = 1 χ = 5; γ = 1 χ = 1; γ = 2 . χ = 3; γ = 2 . χ = 5; γ = 2 . Nonlocal t/τ L . . . . . ( E C k i n + E i n t + E q ) / E A B C k i n (b) t/τ L . . . . . . (cid:15) / (cid:15) m a x (c) t/τ L V H c / κ Figure 4. Time evolution of the (a) incompressible kinetic en-ergy, (b) the sum of the internal, quantum and compressiblekinetic energy and (c) the dissipation rate of incompressibleenergy for runs A in Table I. The inset in (c) shows the evo-lution of the central line helicity. variations respect to the bulk value ρ yield larger valuesof the other energy component between initial times and t ≈ τ L as displayed in Fig. 4 (b). In particular, forthe case of a nonlocal potential the larger values developthrough the whole run. Nevertheless, for all runs duringthe first large-eddy turnover times the main contributionto energy comes from vortices. At later times, energyfrom vortices is converted into sound. As stated in Sec.III A, the decay of the incompressible energy can be usedto estimate the energy dissipation rate (cid:15) . Its temporalevolution is displayed in Fig. 4 (c). As in classical decay-ing turbulent flows, for quantum flows the Kolmogorovregime is more developed at times slightly after the max-imum of dissipation is reached. The green zone in thefigure depicts the temporal window where the system is t/τ L ‘ / ξ χ = 0; γ = 1 χ = 1; γ = 1 χ = 3; γ = 1 χ = 5; γ = 1 χ = 1; γ = 2 . χ = 3; γ = 2 . χ = 5; γ = 2 . Nonlocal
Figure 5. Time evolution of the intervortex distance of thesystem in units of the healing length. All curves correspondto the runs A in Table I. considered to be in a quasi-steady state and a temporalaverage can be performed to improve statistics. The in-set of Fig. 4 (c) displays that the decay of the centralline helicity is independent of the parameters of the gGPmodel and is consistent with the one reported in .As a turbulent flow evolves, the total vortexlength L v varies in time in a competition betweenthe vortex line stretching and the reconnection pro-cess. This quantity can be obtained from the in-compressible momentum density of the flow J I ( k ) andof a two-dimensional point-vortex J D vort ( k ) as L v =2 π (cid:80) k The Kelvin wave cascade discussed in Sec. III A, isformally derived from an incompressible model in a verysimplified theoretical setting. In the context of the GPmodel, the Kelvin wave cascade was first observed in where a setting close to the theoretical prediction wasused. In the case of turbulent tangles, there was first anindirect observation of the Kelvin wave cascade by mak-ing use of the spatiotemporal spectra . In that work,the Kelvin wave dispersion relation was glimpsed anda space-time filtering of the fields was performed yield-ing a scaling in the energy spectrum compatible to theKelvin wave cascade. Then, by using an accurate track-ing algorithm of a turbulent tangle, in reference theL’vov-Nazarenko prediction was clearly observed in thespectrum of large vortex rings extracted from the tangle.Later, in Refs. , by using high-resolution numericalsimulations of the GP model, a secondary scaling rangecompatible with Kelvin wave cascade predictions was ob-served. In this section, we focus on the scaling of theincompressible energy spectra and helicity for the casewith a nonlocal potential (set of runs B) as it seems topresent a much clearer scaling at scales smaller than theintervortex distance. We vary different parameters so the − − kξ − − − − E ( k ) (a) E I kin E C kin E q Hk rot k − / − − kξ − H ( k ) η − (cid:15) / k / (b) Run B1Run B2Run B3Run B6 k rot Figure 7. (a) Helicity and energy spectra of the differentcomponents for the simulation with grid points (RunB6). (b) Helicity spectra compensated by Eq (20) for runsB1-B3 and B6 shown in Table I. range of scales (system size, intervortex distance, healinglength) and energy fluxes take different values.The spectra for the different components constitut-ing the kinetic energy and the helicity of the simulationwith grid points are shown in figure 7 (a). Clearpower laws for the Kolmogorov and Kelvin wave rangeare observed. The quantum energy shows a maximumat the scale associated with the roton minimum, whereasits contribution is negligible at large scales. The helic-ity spectrum also displays a Kolmogorov-like behaviorat large scales, while at scales between the intervortexdistance and the roton minimum it flattens. This flatrange of the helicity spectrum appears in the range wherethe Kelvin wave cascade is dominant. Whether it existsa direct relationship between the Kelvin wave cascadeand the flattening central line helicity spectrum, is stillunclear. Figure 7 (b) displays the compensated helicityspectrum according to (20) for different runs displayingdifferent scale separations and with local and nonlocal po-tentials. The parameters of these simulations correspondto the runs B1-B3 and B6 shown in Table I. At largescales all curves collapse to a constant C H ∼ , while atsmaller scales the system with a wider scale separationdisplays that the helicity contribution is more intense.To analyze further the incompressible energy spectra,we have performed two runs varying the forcing scale sothat the dissipation rate also changes (Runs B4-B5). Werecall that in classical turbulence, the energy flux (cid:15) is1 Run v rms L κ (cid:15) (cid:96) B1 0.395 L/ L/ L/ L/ L/ L/ L , the quantum of circu-lation κ , the energy dissipation rate (cid:15) and the intervortexdistance, expressed in units where the box size is L = 2 π andthe speed of sound c = 1 . fixed by the inertial range and varies as (cid:15) ∼ v /L .Our initial condition ψ ABC keeps fixed, by construction,the value of v rms . In Table II we present the values of dif-ferent physical quantities relevant for a turbulent state.Such quantities are expressed, as customary in classicalturbulence, in units of large scale quantuties. In partic-ular, the system size is L = 2 π and the speed of soundis c = 1 . With such definitions, large scale quantitiesremain almost constant when increasing the scale sep-aration between the box size and the smallest scale ofthe system, but the quantum of circulation takes smallervalues.Figure 8 (a) shows the incompressible energy spectrasimply compensated by k − / . Two plateaux are clearlyobserved but, as expected, they do not collapse becauseenergy fluxes and the intervortex distances have not beentaken into consideration.The energy spectra shown in Fig. 8 (b) have been com-pensated by the Kolmogorov law (19) and displayed asa function of k/k , with k = 2 π/L in order to empha-size the Kolmogorov regime. Once properly normalized,all runs present a plateau at large scales that collapseto values that fluctuate around a Kolmogorov constant C K ∼ , in agreement with previous simulation of theGP model . In order to emphasize the Kelvin wavecascade, we make use of the L’vov & Nazarenko waveturbulence prediction (23). Figure 8 (c) displays the in-compressible energy spectrum compensated by the theo-retical prediction as a function of k/k (cid:96) , with k (cid:96) = 2 π/(cid:96) .The collapse of the Kelvin wave cascade is remarkable.All runs having a nonlocal potential display a plateauaround a value C / ≈ . , which recovers a value of C LN ≈ . . Such value is relatively close to the pre-dicted one C LN = 0 . , in particular by considering allthe phenomenological assumptions made in Sec. III A toadapt the theoretical prediction (21) to the case of a tur-bulent tangle in Eq. 23. Although the GP run (with localinteraction potential) displays a good Kolmogorov scalingat large scales, it does not clearly exhibit a Kelvin cas-cade range at the highest resolution used in this work forthis model ( ). Note that previous works reportinga secondary k − / range in local GP have used resolu-tion of and collocation points ( and re-spectively). The incompressible kinetic energy spectrum,compensated by the Kozik & Svistunov prediction is − − kξ − E I k i n ( k ) k / (a) Run B1Run B2Run B3 Run B4Run B5Run B6 k/k E I k i n ( k ) (cid:15) − / k / (b) − k/k ‘ − E I k i n ( k ) (cid:15) − / ‘ / κ − Λ − Ψ / k / Figure 8. Compensated incompressible kinetic energy spectraby (a) k − / scaling, (b) Kolmogorov scaling and (c) L’vov-Nazarenko scaling for Kelvin waves. displayed in Appendix A. IV. CONCLUSIONS We studied the properties of the freely decaying quan-tum turbulence of the generalized Gross-Pitaevskii (gGP)model (7), that includes beyond mean field correctionsand considers a nonlocal interaction potential betweenbosons. This model pretends to give a better descriptionof superfluid helium as it reproduces a roton minimumin the excitation spectrum.The visualization of the flow with a nonlocal potentialallowed us to observe the formation of helical structuresaround the vortices produced by density fluctuations, ex-hibiting the intrinsic property of maximal helicity of anABC flow. These structures were not observed at initial2times in a flow with no helicity like a Taylor-Green flow,but they develop as the system evolves (data not shown).However, it was seen that the behavior of the helcity isindependent of the interaction potential. At large scalesthe helicity develops a spectrum that satisfies prediction(20), while at scales between the intervortex distance andthe healing length a plateau is observed. This range isusually associated with the Kelvin wave cascade regime,but it is still not known whether the formation of thisplateau is associated with Kelvin waves or not.By studying numerically the freely decaying quantumturbulence of an ABC flow, we observed that the sta-tistical behavior of the system either at large or smallscales does not depend much on the parameters of thebeyond mean field correction in the presence a local in-teraction potential between bosons. Even the introduc-tion of a nonlocal potential does not modify significantlythe behavior of the system at large scales, exhibiting aKolmogorov-like scaling law for the incompressible ki-netic energy. However, the situation changes at smallerscales when a nonlocal potential is implemented, betweenthe intervortex distance (cid:96) and the healing length ξ , rangeassociated with the Kelvin waves cascade. Here, a secondscaling of the incompressible energy spectrum is observedeven at a moderate resolution of grid points, whilein the case of a local GP model an energy spectrum com-patible with k − / scaling law begins to be recognizablefrom resolutions of collocation points , and evenin this case the range of scales where it takes place is lessthan a decade. This enhancement in the Kelvin wavescascade may be very useful for numerical and theoreti-cal studies of wave turbulence. This evident differencewith the local GP model may be used to compare if ef-fectively this model better describes the dynamics of su-perfluid helium. However, experimental observation atscales smaller than the intervortex distance still remainsa challenge.We also studied how is the scaling of the Kelvin wavespectrum with the energy flux (cid:15) and the intervortex dis-tance by varying the integral scale of the initial flow andits healing length. We observed that the different spec-tra tend to collapse to a constant according to L’vov &Nazarenko spectrum for Kelvin waves (23). The valueof the constant observed is C LN ≈ . which is close tothe predicted one C LN ≈ . . This shows how robustis this prediction in the presence of a nonlocal interac-tion potential. This is surprising given that the theoryis constructed from a single vortex line while here it isextended to a vortex tangle, with several phenomenolog-ical assumptions. The Kozik & Svistunov spectrum forKelvin waves was also studied for these energy spectra,however, by compensating them by the theory no clearplateau is observed (see Appendix A). Furthermore, inthe range of the Kelvin wave cascade the Kozik-Svistunovcascade would take values C KS ≈ . which is not of or-der one, so it might imply that the energy spectrum ofthis system is not described by this theory.The overall results of this work show that, even though the GP model does not include the roton minimum inthe excitation spectrum, both GP and gGP models de-scribe a similar behavior at large scales. Some of themare for example the Kolmogorov-like spectrum of the in-compressible kinetic energy and the scaling law observedfor the helicity. This means that the GP equation is stilla good model for describing the macroscopic behavior ofsuperfluid helium but while looking at small scales of aturbulent tangle, some care is needed. Appendix A: Kozik-Svistunov Kelvin spectrum The original Kozik & Svistunov prediction for theKelvin wave cascade was done with same geometricalconsiderations of L’vov & Nazarenko and also expressedin units of Length /T ime . Applying the same consider-ations of Sec. III A to adapt this prediction to a turbulentthree-dimensional flow leads to the following Kelvin waveenergy spectrum E KSKW ( k ) = C KS κ / Λ (cid:15) / (cid:96) − / k / . (A1)where the constant C KS could be in principle determinedby the theory if some integrals in the associated kineticequation are convergent, but its value is still unknown.Figure 9 displays the incompressible kinetic energy spec-trum compensated by prediction (A1). All the curves − k/k ‘ − − E I k i n ( k ) (cid:15) − / ‘ / κ − / Λ − k / Run B1Run B2 Run B3Run B4 Run B5Run B6 Figure 9. Compensated incompressible kinetic energy spectraby the Kozik & Svistunov prediction for Kelvin waves. tend to collapse in the range associated with Kelvinwaves, showing a proper scaling with the energy flux (cid:15) ,the intervortex distance (cid:96) and the quantum of circulation κ . However, although the Kelvin wave range is limited,a plateau is not clearly observed if the spectra are com-pensated by A1 and even though a constant cannot bewell defined, the energy spectra collapse to a constant of C KS ≈ . , which is not of order one. ACKNOWLEDGMENTS N.P.M and G.K were supported by Agence Nationalede la Recherche through the project GIANTE ANR-18-3CE30-0020-01. Computations were carried out on theMésocentre SIGAMM hosted at the Observatoire de laCôte d’Azur and the French HPC Cluster OCCIGENthrough the GENCI allocation A0042A10385. GK is also supported by the EU Horizon 2020 Marie Curie projectHALT and the Simons Foundation Collaboration grantWave Turbulence (Award ID 651471). L. P. Pitaevskii and S. Stringari, Bose-Einstein Condensa-tion and Superfluidity , Vol. 164 (Oxford University Press,2016). M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E.Wieman, and E. A. Cornell, Science, New Series , 198(1995). P. Kapitza, Nature , 74 (1938). J. F. Allen and A. Misener, Nature , 75 (1938). F. London, Nature , 643 (1938). E. Fonda, D. P. Meichle, N. T. Ouellette, S. Hormoz, andD. P. Lathrop, Proceedings of the National Academy ofSciences , 4707 (2014). S. Serafini, M. Barbiero, M. Debortoli, S. Donadello,F. Larcher, F. Dalfovo, G. Lamporesi, and G. Ferrari,Physical Review Letters , 170402 (2015). C. F. Barenghi, L. Skrbek, and K. R. Sreenivasan, Pro-ceedings of the National Academy of Sciences , 4647(2014). K. W. Schwarz, Physical Review B , 2398 (1988). J. Koplik and H. Levine, Physical Review Letters , 1375(1993). C. Nore, M. Abid, and M. E. Brachet, Physics of Fluids ,2644 (1997). A. W. Baggaley, J. Laurie, and C. F. Barenghi, PhysicalReview Letters , 205304 (2012). V. Shukla, P. D. Mininni, G. Krstulovic, P. C. di Leoni,and M. E. Brachet, Physical Review A , 043605 (2019). J. Maurer and P. Tabeling, Europhysics Letters (EPL) ,29 (1998). J. Salort, B. Chabaud, E. Lévêque, and P.-E. Roche, EPL(Europhysics Letters) , 34006 (2012). U. Frisch, Turbulence: The Legacy of A.N. Kolmogorov ,1st ed. (Cambridge University Press, 1995). W. F. Vinen, Physical Review B , 134520 (2001). B. Rousset, P. Bonnay, P. Diribarne, A. Girard, J. M.Poncet, E. Herbert, J. Salort, C. Baudet, B. Castaing,L. Chevillard, F. Daviaud, B. Dubrulle, Y. Gagne, M. Gib-ert, B. Hébral, T. Lehner, P.-E. Roche, B. Saint-Michel,and M. Bon Mardion, Review of Scientific Instruments ,103908 (2014). R. J. Donnelly, Quantized Vortices in Helium II (Cam-bridge University Press, 1991). A. Villois, D. Proment, and G. Krstulovic, Physical ReviewFluids , 044701 (2017). G. Krstulovic, M. Brachet, and E. Tirapegui, Physical Re-view E , 026601 (2008). A. Villois, D. Proment, and G. Krstulovic,arXiv:2005.02048 [cond-mat, physics:nlin, physics:physics](2020), arXiv:2005.02048 [cond-mat, physics:nlin,physics:physics]. R. J. Donnelly and C. F. Barenghi, Journal of Physical andChemical Reference Data , 1217 (1998). Y. Pomeau and S. Rica, Physical Review Letters , 247(1993). N. G. Berloff and P. H. Roberts, Journal of Physics A:Mathematical and General , 5611 (1999). J. Reneuve, J. Salort, and L. Chevillard, Physical ReviewFluids , 114602 (2018). T. D. Lee, K. Huang, and C. N. Yang, Physical Review , 1135 (1957). N. G. Berloff, M. Brachet, and N. P. Proukakis, Proceed-ings of the National Academy of Sciences , 4675 (2014). S. Villerot, B. Castaing, and L. Chevillard, Journal of LowTemperature Physics , 1 (2012). V. S. L’vov and S. Nazarenko, JETP Letters , 428(2010). T. Lahaye, C. Menotti, L. Santos, M. Lewenstein, andT. Pfau, Reports on Progress in Physics , 126401 (2009). A. Griesmaier, Journal of Physics B: Atomic, Molecularand Optical Physics , R91 (2007). L. Santos, G. V. Shlyapnikov, and M. Lewenstein, PhysicalReview Letters , 250403 (2003). S. M. Roccuzzo, A. Gallemi, A. Recati, and S. Stringari,Physical Review Letters , 045702 (2020). S. Nazarenko, Wave Turbulence , Lecture Notes in PhysicsNo. 825 (Springer, Heidelberg, 2011). P. Clark di Leoni, P. D. Mininni, and M. E. Brachet, Phys-ical Review A , 043605 (2016). S. Zuccher and R. L. Ricca, Physical Review E ,061001(R) (2015). M. W. Scheeler, D. Kleckner, D. Proment, G. L. Kindl-mann, and W. T. M. Irvine, Proceedings of the NationalAcademy of Sciences , 15350 (2014). G. Krstulovic and M. Brachet, Physical Review B ,132506 (2011). G. Krstulovic and M. Brachet, Physical Review E ,066311 (2011). S. K. Nemirovskii, Physics Reports , 85 (2013). L. Skrbek and K. R. Sreenivasan, Physics of Fluids ,011301 (2012). A. Brissaud, Physics of Fluids , 1366 (1973). P. Clark di Leoni, P. D. Mininni, and M. E. Brachet, Phys-ical Review A , 053636 (2017). E. V. Kozik and B. V. Svistunov, Journal of Low Temper-ature Physics , 215 (2009). G. Krstulovic, Physical Review E , 055301(R) (2012). L. Boué, R. Dasgupta, J. Laurie, V. L’vov, S. Nazarenko,and I. Procaccia, Physical Review B , 064516 (2011). A. W. Baggaley and J. Laurie, Physical Review B ,014504 (2014). A. Villois, D. Proment, and G. Krstulovic, Physical ReviewE , 061103(R) (2016). V. Eltsov and V. S. L’vov, JETP Letters10.1134/S0021364020070012 (2020). E. B. Sonin, arXiv:2003.09912 [cond-mat] (2020),arXiv:2003.09912 [cond-mat]. V. B. Eltsov and V. S. L’vov, JETP Letters10.1134/S002136402010001X (2020). D. Proment and G. Krstulovic, arXiv:2005.02047[cond-mat, physics:nlin, physics:physics] (2020),arXiv:2005.02047 [cond-mat, physics:nlin, physics:physics]. D. Gottlieb and S. A. Orszag, Numerical analysis of spec-tral methods: theory and applications (SIAM, 1977). P. Clark di Leoni, P. D. Mininni, and M. E. Brachet, Phys-ical Review A92