Grid superfluid turbulence and intermittency at very low temperature

Low-temperature grid generated turbulence is investigated by using numerical simulations of the Gross-Pitaevskii equation. The statistics of regularized velocity increments are studied. Increments of the incompressible velocity are found to be skewed for turbulent states. Results are later confronted with the (quasi) homogeneous and isotropic Taylor-Green flow, revealing the universality of the statistics. For this flow, the statistics are found to be intermittent and a Kolmogorov constant close to the one of classical fluid is found for the second order structure function.


I. INTRODUCTION
Superfluid turbulence has been largely studied in the last decades, especially thanks to the progress achieved in experimental technics.Today it is possible to create turbulent Bose-Einstein condensates (BEC) [1], visualize and track quantum vortices in BECs [2] and 4 He [3][4][5], and study Lagrangian dynamics by using tracers [6,7] in 4 He.As in classical 3D hydrodynamic turbulence [8], a turbulent Kolmogorov cascade is observed at scales in between the energy injection scale and mean inter-vortex distance [9,10].At scales smaller than the inter-vortex distance, the quantized vortices can not be considered as a continuous field and different mechanisms appear to be relevant to carry the energy to scales small enough to be dissipated by phonon emission [11,12].Numerical simulations of different models confirm this scenario [13,[15][16][17].In classical and superfluid three-dimensional turbulence, energy is usually injected at large scales by different types of forcing.In classical hydrodynamic turbulence, one of the most standard ways is by a fluid flowing through a grid [18][19][20].When increasing the mean flow, the fluid behind the grid develops a series of instabilities creating a turbulent wake.Such classical experiments have been also performed during the last decade using superfluids as 4 He [21,22] and 3 He [23].
In this work we investigate low-temperature superfluid turbulence generated by a moving grid using the Gross-Pitavskii equation (GPE).Statistics of (regularized) velocity increments are analyzed.The results are later confronted with homogeneous and (quasi) isotropic turbulent flow generated by the so-called Taylor-Green flow [24].Statistics of velocity increments are shown to be universal.An estimation of the dissipation energy rate leads to a measurement of the Kolmogorov constant close to the one observed in classical turbulence.Finally, the increments of the incompressible velocity are found to be intermittent.
The GPE describing a homogeneous BEC of volume V with (complex) wave-function ψ is given by where m is the mass of the condensed particles and g = 4πa 2 /m, with a the s-wave scattering length.
Madelung's transformation ψ(x, t) m exp [i m φ(x, t)] relates the wave-function ψ to a superfluid of density ρ(x, t) and velocity v = ∇φ, where φ is the phase of the wave-function.κ = h/m is the Onsager-Feynman quantum of velocity circulation around the ψ = 0 vortex lines.When Eq.( 1) is linearized around a constant ψ = ψ0 , the sound velocity is given by c = (g| ψ0 | 2 /m) 1/2 with dispersive effects taking place at length scales smaller than the coherence 2 , which also corresponds to the vortex core size [14].Using the Madelung transformation (see [13] for details) the energy term (per unit of volume) E = ( 2 /2mV ) |∇ψ| 2 dx can be rewritten as where with P inc [ • ] the projector onto the space of divergencefree fields.The super-index stand for incompressible (I), compressible (C) and quantum (Q) velocities.These fields are all regular at the vortex position as they are regularized by the term √ ρ [24].The velocity v I contains the contribution of vortices, whereas v C and v Q are related to waves, since they are by construction potential flows.The energy spectra are defined as , where v Λ is the Fourier transform of v Λ in the plane perpendicular to the mean flow at a given distance z from the grid.The superscript Λ stands for I,C and Q.
In the simulations presented in this work, the mean density mN/V is fixed to 1 and the physical constants in Eq.( 1) are determined by the values of ξ and c = 1.The quantum of circulation is given by 4πc ξ/ √ 2. Numerical integration of Eq.( 1) is performed by using a fully dealiased pseudo-spectral code.The domain is periodic in all directions.The perpendicular (respect to the direction of the mean flow) size of the domain is denoted by L ⊥ and the parallel one by L .The grid is modeled by a strong repulsive potential V Grid (x).The grid is characterized by the diameter of the rods a and the distance between the rods D. A sketch of the grid is shown in Appendix A. The fluid is initially at rest.The system is then advected with a velocity v 0 that is slowly increased from zero up to its final value.During this process, local dissipation is included far from the grid to reduce the sound emitted during the transient (see Appendix B for more details on numerics and methods).Different grids, Mach numbers M = v 0 /c and resolutions are studied (see Table I).Note that such a periodic configuration mimics recent  I. List of runs.L ⊥ and L are the sizes of the domain and N ⊥ and N the corresponding resolutions.D is the distance between the rods (see Appendix A for details).The diameter of the rods is a = 2ξ for all runs.Lengths are expressed in units of the healing length ξ.The Mach number is M = u0/c = 0.8 for all runs except for a1m and a2m with M = 0.6.For all grid runs ξ = .75L⊥ /N ⊥ .For the Taylor-Green run (tg) grid turbulence experiments with 4 He in a ring [21], and similar ideas have been used in 2D to study the possibility of an inverse cascade [25].We also study (quasi) homogenous and isotropic turbulence generated by the Taylor-Green flow (see Fig. 5 below).This standard flow consists of a number of vortex rings and it develops a turbulent tangle followed by a small-scale thermalization of sound waves (due to the Galerkin projection of GP equation).Vortices exchange momentum and energy with the thermalized waves mimicking mutual-friction effects.These generic properties of the GP model are also expected to be present in the grid simulations.Note that the thermalization and its associated effects also occur independently of the spectral cut-off if dispersive effects are important [26,27].
We first focus on simulations done using the grid.The grid generates a turbulent wake that is displayed in Figs.1.a-dby the isosurface of the density.At early times vortices are nucleated close to the grid (Figs.1.a-c, run c1), leading later to a turbulent wake (Figs.1.h-j,run b3).It is well known in the framework of 2D GP that vortex dipoles are nucleated behind a cylindrical obstacle for Mach numbers above a critical threshold M crit ≈ 0.4 [28].The equivalent in 3D are vortex rings that rapidly reconnect and create the complex tangle observed in Fig. 1.The process is identical to the one described for 3 He-B experiments using a grid reported in [23].
Two stages are observed during the development of turbulence in the wake of the grid.During the first stage, incompressible kinetic energy is injected by nucleation of rings.To account for this, the total vortex length is measured as , where θ( ) is the Heaviside function and ρ 2D (x) is the profile of a two-dimensional vortex given by the Padé approximation [30].Note that L(t) is only a rough estimate as small amplitude Kelvin waves and density oscillations along filaments modify this volume integral.The temporal evolution of L(t) is displayed in Fig. 2.a for all runs.The increase of the vortex length depends on the geometry of the grid and on the Mach number.From Fig. 1.a-g we observe two kind of structures: large elongated rings of length ∼ 2D and smaller ones close to the corners.When ξ D we expect that the main contribution to L comes from elongated rings.The vortex length can be thus estimated as L(t) ∼ D(L ⊥ /D) 2 t/T nucl , where T nucl is the typical time between two nucleations ((L ⊥ /D) 2 is proportional to the number of such rings).The characteristic time-scale of vortex injection for the grid is thus defined in terms of T nucl as τ Grid = T nucl Dξ/L 2 ⊥ (so that L/ξ ∼ t/τ Grid ).Vortex nucleation in a superflow around a disc has been extensively studied in the last 15 years [28,29,31].It was found that stable and unstable (nucleation) branches are connected through a primary saddle-node and a secondary pitchfork bifurcation.The critical Mach number M crit was found to depend on a/ξ and close to the critical point T nucl to scale as T nucl ∼ (M/M crit − 1) −1/2 , that corresponds to a dissipative saddle-node bifurcation.For the grid, geometry is more complex and M M crit , so the previous scaling is not expected to be valid.A precise determination of T nucl is out of the scope of the present work.However it can be empirically observed that T nucl ∼ (ξ/c)(D/ξ) 1/2 M −4 is compatible with data presented in this work.This is manifested by the relatively good collapse of L/ξ for the different runs displayed in the inset of Fig. 2.a.A precise study of the nucleation will be performed in a future work.reaches a maximum and then decreases, consistently with the decay of L(t) in Fig. 2.a.The time when the vortex length and energy fluctuations are the largest, corresponds to the time when the bulk of vortices reaches the opposite side of the box (respect to the grid).This fact has been checked with runs using the same parameters but with larger L (data not shown).By this time, the condensate (initially at the wavenumber k c = 0) has "jumped" to higher wavenumbers.This can be interpreted as the full system being entrained by the imposed flow.Indeed, a boost of velocity v G corresponds for GP to a multiplication of ψ by e i m vG•k [32].The mean (incompressible) kinetic energy of the condensate is thus given by 1  2 where k c is determined by the wavenumber with the largest number of particles.The temporal evolution of this energy is also shown in Fig. 2  ).For all other runs, the integration is not long enough to observe the decay.The same phenomenon is observed in equivalent 2D simulations (as the one in [25]) if the integration is performed for longer times (data not shown).
Before the decay starts, a turbulent state is observed.The energy spectrum computed at a distance z ≈ 360ξ from the grid is displayed in Fig. 2.c.A Kolmogorov scaling is expected to be observed for k k IV = 2π/ IV , where IV is the inter-vortex distance, usually estimated as IV = 1/ L/V .For the grid, turbulence is not homogeneous, but an effective volume can be obtained through a spatial average weighted by ∆E I (z) = k>0 E I k (z).A plot of ∆E I (z) and details on this average are included in the Appendix C. For the corresponding run and time of Fig. 2.c we obtain V eff = L 2 ⊥ × 169ξ that yields IV = 10.6ξ.This corresponds to k IV ξ ≈ .6,which is in good agreement with the end of the k −5/3 scaling observed in Fig. 2.c.Finally, in the second stage, rings shrink due to mutual friction effects [32].The estimations of L by a volume integral does not allow us to verify the Vinen's decay prediction [33].
One the most remarkable differences between classical and superfluid turbulence is the one-point velocity statistics [17,34,35].The probability distribution func-tion (PDF) of v = ∇φ presents power law tails ∼ v −3 unlike the Gaussian PDFs observed in the classical case.As in GP, the velocity field v is ill-defined at the vortex core because the phase is not defined, we instead look at statistics of the regularized fields (2).Non-Gaussian PDF have already been observed for v I in 2D GP simulations [27].The PDFs of the three components of the velocity are displayed in Fig. 3a-c for v I , v C and v Q at a fixed time and distance from the grid for run c1.The PDFs of v I present as in reference [27] non-Gaussian tails scaling as v I −b with b ∈ (2, 3).On the contrary v C and v Q exhibit almost Gaussian PDFs.These can be explained because sound waves (related to v C and v Q ) are indeed expected to thermalize at small scales and thus to develop Gaussian statistics [26].
Motivated by classical turbulence we define the longitudinal velocity increments as with r a unit vector.For the grid r is taken perpendicular to the mean flow.The velocity increments PDFs are displayed in Fig. 3d-e for I and Q.As in classical turbulence δv I / δv I 2 1/2 presents strongly non-Gaussian statistics that depend on the scale , manifesting the non self-similar behavior of turbulence.Figure 4.a shows a zoom of Fig. 3d together with the PDF of −δv I / δv I 2 1/2 in dashed lines.A negative skewness is apparent there.This asymmetry of the PDFs is an important property of Kolmogorov turbulence related to the energy cascade and the 4/5-law of turbulence [8].The bulk of the PDF of δv Q is Gaussian whereas the tails depend on the scale and separate from Gaussian statistics.No skewness is observed.The increment of δv C are totally Gaussian and scale independent once normalized by their rms value (not shown).
It is well known that strong velocity fluctuations in classical turbulence lead to the breakdown of the totally self-similar Kolmogorov phenomenology (K41).Intermittency is responsible for this breakdown and it is quantified by looking at the scaling of the velocity increments moments known as structure functions (average is over all directions of r).In the inertial range, i.e. at scales smaller than the integral scale L int and larger than the dissipative scale η, it is expected that S I p ( ) ∼ ζ I p .K41 predicts ζ p = p/3, whereas numerical and experimental results evidence a non-linear function [8].The deviation from ζ p = p/3 are known as intermittency corrections and ζ p as anomalous exponents.Note the (analytical) 4/5law fixes ζ 3 = 1.We now address this issue within the framework of GP.The statistics presented in Fig. 3 do not allow to obtain a clear scaling.In order to obtain a larger inertial range and better statistics, we make use of the Taylor-Green flow.This flow is known to develop a vortex tangle with a k −5/3 energy spectrum [13].A visualization of the Taylor-Green flow is presented in Fig. 5 at different times.We first compare the statistics of the TG velocity increments with those of the grid.We define a Kolmogorov (like) dissipative scale η using the quantum of circulation as η = L int (v I rms L int /κ) −3/4 , where the integral scale L int is estimated as D and L/2 for the grid and the Taylor-Green flow respectively.The corresponding values are ηGrid = 23ξ for run c1 and ηTG = 8.5ξ for Taylor-Green run.Note in Fig. 4.b that the statistics of both flows coincide, if velocity increments are compared at similar scales (in units of η).This is a manifestation universality in quantum turbulence like the one observed in classical turbulence.Slight discrepancies between the two configurations and small values of /η are due to the non unique way of defining a Kolmogorov length in quantum turbulence.
The PDFs of Taylor-Green flow velocity increments at different scales are displayed in Fig. 6.a-b for two different times.At t = 5 a clear scale dependence and skewness are observed.For instance, the skewness is equal to −0.13 for = 16ξ at t = 5.At later times (t = 30), the increments tend towards non-skewed PDFs, thought not Gaussian.The structure functions are presented in Fig. 6.c.K41 predicts S I 2 ( ) ≈ (11/3)C 2 ( ) 2/3 , with the energy dissipation rate and C 2 the Kolmogorov constant [19].In GP, can be estimated as = − dE I dt .By using this estimation we obtain (by fitting) C 2 = 2.6625 (see dashed line in Fig.The anomalous exponents are displayed in the inset of Fig. 6.c.The red dashed and blue dot-dashed lines represents K41 and She-Lévêque model respectively [36].GP intermittency is found to be stronger than the classical one (that is in general well represented by She-Lévêque model).Intermittency was already measured in 4 He by early experiments performed by Maurer et al. [9] and no difference with classical experiment was found.However, it has been observed that intermittency of the von Kármán flow (in classical fluids) is slightly stronger than other turbulent flows [37].As the Taylor-Green flow mimics the von Kármán flow, an enhancement of intermittency could also be expected.In addition, using HVBK-based shell models [38,39] a clear temperature dependence has been observed for the ζ p , presenting a maximum of intermittency around 0.6T λ (with T λ the temperature of the λ-point).These HVBK results do not directly apply to GP turbulence, that formally describes the low temperature limit of BECs.Indeed, in this limit dissipation need to be added by some ad-hoc mechanism to the HVBK model, unlike GP, where energy of vortices is naturally dissipated by phonon radiation.Furthermore, in the HVBK there is no notion of quantized vortices, as only a large-scale description is given.The results presented in this work directly apply to BECs at low temperature but are also expected to be relevant for superfluid Helium.Although today it is possible to create and track several vortex lines in BECs [2,40], a controlled experiment with such a dense turbulent vortex tangle is not still realizable.However, the large fluctuations of velocity fields reported in this work are expected to be an inherent property of turbulent BECs that could be observed in the future.
Understanding of intermittency in classical flows remains an open problem, in superfluids not enough information is available.The simulations presented here are not in a statistically steady-state which is the most suitable configuration for such a study.However, it has been shown that velocity statistics of grid turbulence are similar to these of Taylor-Green.Grid simulations could be thus used to investigate intermittency if the injection/dissipation is modified to obtain a stationary regime.Much longer simulations at higher resolutions are needed.The Galerkin projector P G takes a simple form in Fourier space: P G [ ψ (k)] = θ(k max − |k|) ψ (k) where θ is the Heaviside function and k max = N/3 is chosen following the standard 2/3 rule for a quadratic non-linearity [41].When de-aliasing is not performed as in (B1), conservation of momentum is not preserved by the discrete system.For the grid simulation, the system is not isotropic and it has finite momentum in one direction.The lack of momentum conservation typically leads to spurious and non controlled effects.The additional projector applied in Eq.(B1) costs one extra back and forth FFT per time step.An alternative to such technique is to use the standard de-aliasing but with k max = N/4 (corresponding to a cubic non-linearity) at the price of wasting half of the resolution.
The Taylor-Green flow is prepared as in reference [24] (with a different choice of parameters) but no symmetries are imposed during the temporal evolution.Symmetries are thus not preserved for all times.

Appendix C: Inter-vortex distance
The inter-vortex distance is usually estimated as IV = 1/ L/V , where L is the total vortex length and V the volume of the system.For the grid, turbulence is not homogenous therefore an estimation of the volume is needed.Vortices manly contribute to the incompressible kinetic energy.This quantity can be thus used to estimate the size of the turbulent bulk.We define energy fluctuations due to turbulence as a function of the distance to the grid z by ∆E I (z) = Finally, the inter-vortex distance is estimated as IV = V eff /L.

Figure 2 .
Figure 2.b displays the temporal evolution of the incompressible kinetic energy of run b1.The saturation of the energy is related to the growth of the mean flow, as shown by the temporal evolution of the k = 0 mode E I k=0 of the 3D energy spectrum (including the average over z).Note that the energy fluctuations ∆E I = E I − E I k=0

FIG. 2 .
FIG. 2. (Color online) a) Total vortex line for different grid runs.The inset shows the collapse of L/ξ vs t/τ Grid , with τ Grid = T nucl Dξ/L 2 ⊥ (see text).b) Temporal evolution of E I , the kinetic energy E I k=0 of the mean incompressible velocity field, the turbulent fluctuations ∆E I = E I − E I k=0 and the kinetic energy of the condensate (all energies averaged over z).Run b1.c) Averaged energy spectrum for run c1 at z/ξ ∈ (355, 365) and t = 5.5.

FIG. 3 .
FIG.3.(Color online) a-c) One-point velocity PDF of the three components of v I , v C and v Q (same color code for ac).d-e) PDFs of the velocity increments of δv I /δv I vrms and δv Q /δv Q vrms for different increments .Same scales as in e).All the data from run c1 taken at t = 5.5 and z = 360ξ.

FIG. 5 .
FIG. 5. (Color online) (Color online) 3D visualizations of the density field (rendered with the software VAPOR).Red isosurfaces are at low density values corresponding to vortices.Green density clouds correspond to sound waves (density fluctuations around ρ = 1).Temporal evolution (left to right) of Taylor-Green vortex.t = 0, 12.5, 33.5.
6.c).Note that C 2 = 2.6625, is very close to the value 2.0 ± 0.4 reported in classical turbulence.To look at the intermittency, the local slopes ζ p ( ) = d log S I p ( ) d log are presented in Fig.6.d.A power-law scaling of S I p ( ) corresponds to the plateau in ζ p ( ).The ζ I p are measured averaging the local slope for /ξ ∈ (15, 75).

FIG. 6 .
FIG. 6. (Color online) Taylor-Green run.a) PDFs of the velocity increments at t = 5 for different scales (colors as in (b)).Dashed lines correspond to PDFs of −δv I / (δv I ) 2 1/2 .b) Idem as a) but for t = 30.c) Structure functions (4) at t = 5.The black dashed line represent the K41 prediction S I 2 ( ) ≈ (11/3)C2( ) 2/3 with C2 = 2.6625 and = − dE I dt .The inset displays the anomalous exponent ζ I p , the red dashed line the Kolmogorov scaling ζp = p/3 and the point-dashed blue line the She-Lévêque model.d) Local slope of the structure functions.

FIG. 7 .
FIG. 7. (Color online) a) Scheme of the domain.b) Scheme of local dissipation.c) Protocol followed to increase the Mach number (arbitrary units).