# Grid superfluid turbulence and intermittency at very low temperature

Giorgio Krstulovic

<figure><img src="image_1.png" /></figure>

## To cite this version:

Giorgio Krstulovic. Grid superfluid turbulence and intermittency at very low temperature. Physical Review E, 2016, 93 (6), pp.063104. 10.1103/PhysRevE.93.063104. hal-03526708

HAL Id: hal-03526708

https://hal.science/hal-03526708

Submitted on 17 Jan 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Grid superfluid turbulence and intermittency at very low temperature

Giorgio Krstulovic

Laboratoire J.L. Lagrange, UMR7293, Université de la Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, B.P. 4229, 06304 Nice Cedex 4, France

(Received 23 November 2015; revised manuscript received 18 May 2016; published 8 June 2016)

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.

DOI: 10.1103/PhysRevE.93.063104

## I. INTRODUCTION

Superfluid turbulence has been largely studied in recent decades 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–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 intervortex distance [9,10]. At scales smaller than the intervortex distance, the quantized vortices cannot 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–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–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 past 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-Pitaevskii 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.

## II. EQUATIONS AND NUMERICAL SETUP

The GPE describing a homogeneous BEC of volume V with (complex) wave-function $\psi$ is given by

$$i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi + g|\psi|^2 \psi, \quad (1)$$

where m is the mass of the condensed particles and g = $4\pi a\hbar^2/m$, with a the s-wave scattering length. Madelung's transformation $\psi(\mathbf{x},t) = \sqrt{\frac{\rho(\mathbf{x},t)}{m}} \exp[i\frac{m}{\hbar}\phi(\mathbf{x},t)]$ relates the wave-function $\psi$ to a superfluid of density $\rho(\mathbf{x},t)$ and velocity $\mathbf{v} = \nabla\phi$, where $\phi$ is the phase of the wave function. $\kappa = h/m$ is the Onsager-Feynman quantum of velocity circulation around the $\psi = 0$ vortex lines. When Eq. (1) is linearized around a constant $\psi = \hat{\psi}_0$, the sound velocity is given by $c = (g|\hat{\psi}_0|^2/m)^{1/2}$ with dispersive effects taking place at length scales smaller than the coherence length $\xi = (\hbar^2/2m|\hat{\psi}_0|^2g)^{1/2}$, which also corresponds to the vortex core size [14]. Using the Madelung transformation (see Ref. [13] for details) the energy term (per unit of volume) $E = (\hbar^2/2mV) \int |\nabla\psi|^2 d\mathbf{x}$ can be rewritten as $E = E^1 + E^C + E^Q = (1/2V) \int (|\mathbf{v}^1|^2 + |\mathbf{v}^C|^2 + |\mathbf{v}^Q|^2) d\mathbf{x}$, where

$$\mathbf{v}^1 = \mathcal{P}_{\text{inc}}[\sqrt{\rho}\nabla\phi], \quad \mathbf{v}^C = \sqrt{\rho}\nabla\phi - \mathbf{v}^1, \quad \mathbf{v}^Q = \frac{\hbar}{m}\nabla\sqrt{\rho}, \quad (2)$$

with $\mathcal{P}_{\text{inc}}[\cdot]$ the projector onto the space of divergence-free 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 $\sqrt{\rho}$ [24]. The velocity $\mathbf{v}^1$ contains the contribution of vortices, whereas $\mathbf{v}^C$ and $\mathbf{v}^Q$ are related to waves, since they are by construction potential flows. The energy spectra are defined as $E_k^\Lambda(z) = \frac{1}{2} \int_{|\mathbf{p}|=k} |\tilde{\mathbf{v}}^\Lambda(\mathbf{p},z)|^2 d^2\mathbf{p}$, where $\tilde{\mathbf{v}}^\Lambda$ is the Fourier transform of $\mathbf{v}^\Lambda$ in the plane perpendicular to the mean flow at a given distance z from the grid. The superscript $\Lambda$ 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 $\xi$ and $c = 1$. The quantum of circulation is given by $4\pi c \xi / \sqrt{2}$. Numerical integration of Eq. (1) is performed by using a fully dealiased pseudospectral 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_\perp$ and the parallel one by $L_\parallel$. The grid is modeled by a strong repulsive potential $V_{\text{Grid}}(\mathbf{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
<table><caption>TABLE 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 = u_0/c = 0.8$ for all runs except for a1m and a2m with $M = 0.6$. For all grid runs $ξ = 0.75L_⊥/N_⊥$. For the Taylor-Green run (tg) $L_⊥ = L_‖ = L$, $N = N_⊥ = N_‖$, and $ξ = L/N$. No symmetries are enforced.</caption><thead><tr><th>Run</th><th>L<sub>⊥</sub></th><th>L<sub>‖</sub></th><th>L<sub>⊥</sub>/D</th><th>N<sub>⊥</sub></th><th>N<sub>‖</sub></th><th>Run</th><th>L<sub>⊥</sub></th><th>L<sub>‖</sub></th><th>L<sub>⊥</sub>/D</th><th>N<sub>⊥</sub></th><th>N<sub>‖</sub></th></tr></thead><tbody><tr><td>a1</td><td>170</td><td>683</td><td>9</td><td>128</td><td>512</td><td>b1</td><td>341</td><td>683</td><td>17</td><td>256</td><td>512</td></tr><tr><td>a2</td><td>170</td><td>683</td><td>5</td><td>128</td><td>512</td><td>b2</td><td>341</td><td>683</td><td>9</td><td>256</td><td>512</td></tr><tr><td>a3</td><td>170</td><td>683</td><td>7</td><td>128</td><td>512</td><td>b3</td><td>341</td><td>683</td><td>7</td><td>256</td><td>512</td></tr><tr><td>a4</td><td>170</td><td>683</td><td>3</td><td>128</td><td>512</td><td>b4</td><td>341</td><td>683</td><td>5</td><td>256</td><td>512</td></tr><tr><td>c1</td><td>683</td><td>683</td><td>7</td><td>512</td><td>512</td><td>tg</td><td>512</td><td>512</td><td>-</td><td>512</td><td>512</td></tr></tbody></table>

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 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 cutoff if dispersive effects are important [26,27].

### III. GENERATION AND DECAY OF TURBULENCE

We first focus on simulations done using the grid. The grid generates a turbulent wake that is displayed in Figs. 1(a)–1(d) by the isosurface of the density. At early times vortices are nucleated close to the grid [Figs. 1(a)–1(c), run c1], leading later to a turbulent wake [Figs. 1(h)–1(j), run b3]. See Supplemental Material [28] for a movie showing the vortex nucleation close to the grid. 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} \approx 0.4$ [29]. 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 Ref. [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 $L(t) = \int \theta(0.2 - \rho(\mathbf{x},t))d^3\mathbf{x}/\int \theta(0.2 - \rho^{2D}(\mathbf{x}))d^2\mathbf{x}$, where $\theta(\cdot)$ is the Heaviside function and $\rho^{2D}(\mathbf{x})$ is the profile of a two-dimensional vortex given by the Padé approximation [31]. 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 Figs. 1(a)–1(g)

<figure><img src="image_3.png" /><figcaption>FIG. 1. 3D visualizations of the density field (rendered with the software VAPOR). Red isosurfaces are at low density values corresponding to vortices. Blue and green density clouds correspond to sound waves (density fluctuations around $\rho = 1$). (a–g) Zoom close to the grid at early times $t = 0.1, 0.6, 1.1, 1.6, 2.4, 3.1$, and 6.4 (from left to right) for run c1. (h–j) Run b3 at $t = 21, 70, 126$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 2. (a) Total vortex line for different grid runs. The inset shows the collapse of $\mathcal{L}/\xi$ vs. $t/\tau_{\text{Grid}}$, with $\tau_{\text{Grid}} = T_{\text{nucl}} D\xi / L_{\perp}^2$ (see text). (b) Temporal evolution of $E^1$, the kinetic energy $E_{k=0}^1$ of the mean incompressible velocity field, the turbulent fluctuations $\Delta E^1 = E^1 - E_{k=0}^1$, and the kinetic energy of the condensate (all energies averaged over $z$) of run b1. (c) Averaged energy spectrum for run c1 at $z/\xi \in (355,365)$ and $t = 5.5$.</figcaption></figure>

we observe two kinds of structures: large elongated rings of length ~2D and smaller ones close to the corners. When $\xi \ll D$ we expect that the main contribution to $\mathcal{L}$ comes from elongated rings. The vortex length can be thus estimated as $\mathcal{L}(t) \sim D(L_{\perp}/D)^2 t/T_{\text{nucl}}$, where $T_{\text{nucl}}$ is the typical time between two nucleations ($(L_{\perp}/D)^2$ is proportional to the number of such rings). The characteristic timescale of vortex injection for the grid is thus defined in terms of $T_{\text{nucl}}$ as $\tau_{\text{Grid}} = T_{\text{nucl}} D\xi / L_{\perp}^2$ (so that $\mathcal{L}/\xi \sim t/\tau_{\text{Grid}}$). Vortex nucleation in a superflow around a disk has been extensively studied in the past 15 years [29,30,32]. 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_{\text{crit}}$ was found to depend on $a/\xi$ and close to the critical point $T_{\text{nucl}}$ to scale as $T_{\text{nucl}} \sim (M/M_{\text{crit}} - 1)^{-1/2}$, which corresponds to a dissipative saddle-node bifurcation. For the grid, geometry is more complex and $M \gg M_{\text{crit}}$, so the previous scaling is not expected to be valid. A precise determination of $T_{\text{nucl}}$ is out of the scope of the present work. However, it can be empirically observed that $T_{\text{nucl}} \sim (\xi/c)(D/\xi)^{1/2}M^{-4}$ is compatible with data presented in this work. This is manifested by the relatively good collapse of $\mathcal{L}/\xi$ 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.

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_{k=0}^1$ of the 3D energy spectrum (including the average over $z$). Note that the energy fluctuations $\Delta E^1 = E^1 - E_{k=0}^1$ reaches a maximum and then decreases, consistently with the decay of $\mathcal{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_{\parallel}$ (data not shown). By this time, the condensate (initially at the wavenumber $k_c = 0$) has “jumped” to higher wave numbers. 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 $\psi$ by $e^{i\frac{m}{\hbar}v_G \cdot k}$ [33]. The mean (incompressible) kinetic energy of the condensate is thus given by $\frac{1}{2}(\hbar/m)^2|k_c|^2$, where $k_c$ is determined by the wave number with the largest number of particles. The temporal evolution of this energy is also shown in Fig. 2(b). The full system moving at velocity $k_c\hbar/m$ has an effective Mach number $(v_0 - k_c\hbar/m)/c$ (e.g., ~4.6 for run b1). At this Mach number ring nucleation stops, leading to the later decay of $\mathcal{L}(t)$ (runs b1–b4). 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 Ref. [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 \approx 360\xi$ from the grid is displayed in Fig. 2(c). A Kolmogorov scaling is expected to be observed for $k \ll k_{\text{IV}} = 2\pi/\ell_{\text{IV}}$, where $\ell_{\text{IV}}$ is the intervortex distance, usually estimated as $\ell_{\text{IV}} = 1/\sqrt{\mathcal{L}/V}$. For the grid, turbulence is not homogeneous, but an effective volume can be obtained through a spatial average weighted by $\Delta E^1(z) = \sum_{k>0} E_k^1(z)$. A plot of $\Delta E^1(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_{\text{eff}} = L_{\perp}^2 \times 169\xi$ that yields $\ell_{\text{IV}} = 10.6\xi$. This corresponds to $k_{\text{IV}}\xi \approx .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 [33]. The estimations of $\mathcal{L}$ by a volume integral does not allow us to verify the Vinen’s decay prediction [34].

## IV. VELOCITY STATISTICS AND INTERMITTENCY

One the most remarkable differences between classical and superfluid turbulence is the one-point velocity statistics [17,35,36]. The probability distribution function (PDF) of $v = \nabla\phi$ presents power law tails $\sim 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 Eq. (2). Non-Gaussian PDF have already been observed for $v^1$ in 2D GP simulations [27]. The PDFs of the three components of the velocity are displayed in Figs. 3(a)–3(c) for $v^1$, $v^C$, and $v^Q$ at a fixed time and distance from the grid for run c1. The PDFs of $v^1$ present as in Ref. [27] non-Gaussian tails scaling as $v^{1-b}$ with $b \in (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

$$ \delta v_{\ell}^{\Lambda} = [\mathbf{v}^{\Lambda}(\mathbf{x} + \ell \hat{\mathbf{r}}) - \mathbf{v}^{\Lambda}(\mathbf{x})] \cdot \hat{\mathbf{r}}, \quad (3) $$

with $\hat{\mathbf{r}}$ a unit vector. For the grid $\hat{\mathbf{r}}$ is taken perpendicular to the mean flow. The velocity increments PDFs are displayed
<figure><img src="image_5.png" /><figcaption>FIG. 3. (a–c) One-point velocity PDF of the three components of $\mathbf{v}^1$, $\mathbf{v}^C$, and $\mathbf{v}^Q$ (same color code for a–c). (d, e) PDFs of the velocity increments of $\delta v^I / \delta v_{\text{rms}}^I$ and $\delta v^Q / \delta v_{\text{rms}}^Q$ for different increments $\ell$. Same scales as in (e). All the data from run c1 taken at $t = 5.5$ and $z = 360\xi$.</figcaption></figure>

in Figs. 3(d) and 3(e) for I and Q. As in classical turbulence $\delta v_{\ell}^I / \langle \delta v_{\ell}^I \rangle^{1/2}$ presents strongly non-Gaussian statistics that depend on the scale $\ell$, manifesting the non self-similar behavior of turbulence. Figure 4(a) shows a zoom of Fig. 3(d) together with the PDF of $-\delta v_{\ell}^I / \langle \delta v_{\ell}^I \rangle^{1/2}$ in dashed lines. A negative skewness is apparent there. This asymmetry of the PDFs is an important property of Kolmogorov turbulence

<figure><img src="image_6.png" /><figcaption>FIG. 4. (a) Zoom of Fig 3(a) together with the corresponding PDFs of $-\delta v_{\ell}^I / \langle \delta v_{\ell}^I \rangle^{1/2}$ in dashed lines. (b) Comparison between grid (run c1) and Taylor-Green turbulence.</figcaption></figure>

related to the energy cascade and the 4/5-law of turbulence [8]. The bulk of the PDF of $\delta v_{\ell}^Q$ is Gaussian, whereas the tails depend on the scale and separate from Gaussian statistics. No skewness is observed. The increment of $\delta v_{\ell}^C$ are totally Gaussian and scale independent once normalized by their root mean square 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,

$$S_p^I(\ell) = \langle \| (\mathbf{v}^I(\mathbf{x} + \ell\hat{\mathbf{r}}) - \mathbf{v}^I(\mathbf{x})) \| \rangle^p, \quad (4)$$

known as structure functions (average is over all directions of $\hat{\mathbf{r}}$). In the inertial range, i.e., at scales smaller than the integral scale $L_{\text{int}}$ and larger than the dissipative scale $\eta$, it is expected that $S_p^I(\ell) \sim \ell^{\zeta_p}$. K41 predicts $\zeta_p = p/3$, whereas numerical and experimental results evidence a nonlinear function [8]. The deviation from $\zeta_p = p/3$ are known as intermittency corrections and $\zeta_p$ as anomalous exponents. Note the (analytical) 4/5-law fixes $\zeta_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 $\tilde{\eta}$ using the quantum of circulation as $\tilde{\eta} = L_{\text{int}}(v_{\text{rms}}^I L_{\text{int}}/\kappa)^{-3/4}$, where the integral scale $L_{\text{int}}$ is estimated as $D$ and $L/2$ for the grid and the Taylor-Green flow, respectively. The corresponding values are $\tilde{\eta}_{\text{Grid}} = 23\xi$ for run c1 and $\tilde{\eta}_{\text{TG}} = 8.5\xi$ 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 $\tilde{\eta}$). 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 $\ell/\tilde{\eta}$ are due to the nonunique way of defining a Kolmogorov length in quantum turbulence.

The PDFs of Taylor-Green flow velocity increments at different scales are displayed in Figs. 6(a) and 6(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 $\ell = 16\xi$ at $t = 5$. At later times ($t = 30$), the increments tend towards nonskewed PDFs, though not Gaussian. The structure functions are presented in Fig. 6(c). K41 predicts $S_2^I(\ell) \approx (11/3)C_2(\epsilon\ell)^{2/3}$, with $\epsilon$ the energy dissipation rate and $C_2$ the Kolmogorov constant [19]. In GP, $\epsilon$ can be estimated as $\epsilon = -\frac{dE^I}{dt}$. By using this estimation we obtain (by fitting) $C_2 = 2.6625$ [see dashed line in Fig. 6(c)]. Note that $C_2 = 2.6625$ is very close to the value $2.0 \pm 0.4$ reported in classical turbulence.

To look at the intermittency, the local slopes $\zeta_p(\ell) = \frac{d \log S_p^I(\ell)}{d \log \ell}$ are presented in Fig. 6(d). A power-law scaling of $S_p^I(\ell)$ corresponds to the plateau in $\zeta_p(\ell)$. The $\zeta_p^I$ are measured averaging the local slope for $\ell/\xi \in (15,75)$. The anomalous exponents are displayed in the inset of Fig. 6(c). The red dashed and blue dot-dashed lines represents K41 and
<figure><img src="image_7.png" /><figcaption>FIG. 5. 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 $\rho = 1$). Temporal evolution (left to right) of Taylor-Green vortex. $t = 0, 12.5, 33.5$.</figcaption></figure>

She-Lévêque model, respectively [37]. GP intermittency is found to be stronger than the classical one (that is in general well represented by She-Lévêque model).

## V. DISCUSSION

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 [38]. 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 [39,40] a clear temperature dependence has been observed for the $\zeta_p$, presenting a maximum

of intermittency around $0.6T_\lambda$ (with $T_\lambda$ the temperature of the $\lambda$-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,41], 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 and dissipation is modified to obtain a stationary regime. Much longer simulations at higher resolutions are needed.

<figure><img src="image_8.png" /><figcaption>FIG. 6. 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 $-\delta v_t^l / ((\delta v_t^l)^2)^{1/2}$. (b) Idem as (a) but for $t = 30$. (c) Structure functions (4) at $t = 5$. The black dashed line represents the K41 prediction $S_2^l(\ell) \approx (11/3)C_2(\epsilon\ell)^{2/3}$ with $C_2 = 2.6625$ and $\epsilon = -\frac{dE}{dt}$. The inset displays the anomalous exponent $\zeta_p^1$, the red dashed line the Kolmogorov scaling $\zeta_p = p/3$ and the point-dashed blue line the She-Lévêque model. (d) Local slope of the structure functions.</figcaption></figure>

## ACKNOWLEDGMENTS

The author acknowledges useful scientific discussions with J. Bec, M. E. Brachet, V. Shukla. Computations were carried out at Mésocentre SIGAMM hosted at the Observatoire de la Côte d'Azur.

## APPENDIX A: MODEL AND PROCEDURE

We consider the Gross-Pitaevskii equation. The grid is modeled by a strong repulsive potential $V_{\text{grid}}(\mathbf{x})$ and an advection term is added in the left-hand side to impose the mean flow. The final equation reads

$$i\hbar \left( \frac{\partial \psi}{\partial t} + \mathbf{v}_0 \cdot \nabla \psi \right) = -\frac{\hbar^2}{2m} \nabla^2 \psi + g |\psi|^2 \psi - \tilde{\mu} \psi + V_{\text{grid}}(\mathbf{x}) \psi. \quad (A1)$$
<figure><img src="image_9.png" /><figcaption>FIG. 7. (a) Scheme of the domain. (b) Scheme of local dissipation. (c) Protocol followed to increase the Mach number (arbitrary units).</figcaption></figure>

When $\vec{v}_0 = 0$ and $V_{\text{grid}}(\mathbf{x}) = 0$, Eq. (A1) conserves the total energy $H = \int (\frac{\hbar^2}{2m} |\nabla \psi|^2 + \frac{g}{2} |\psi|^4) d\mathbf{x}$ and the total number of particles $\int |\psi|^2 d\mathbf{x}$. The grid potential $V_{\text{grid}}(\mathbf{x})$ is defined as follows

$$ 
\begin{aligned}
V_{\text{grid}}(x,y,z) &= V_0 V_{\text{rod1D}}(z,z_0,a_{\parallel},L_{\parallel}) V_{\text{grid2D}}(x,y), \\
V_{\text{grid2D}}(x,y) &= \min [V_{\text{grid1D}}(x) + V_{\text{grid1D}}(y), 1], \\
V_{\text{grid1D}}(x) &= \sum_{i=1}^{L_{\perp}/D} V_{\text{rod1D}} \left( x, \left(i - \frac{1}{2}\right)D, a_{\perp}, L_{\perp} \right), \\
V_{\text{rod1D}}(x,x_0,a,L) &= \frac{1}{2} + \frac{1}{2} \tanh \left[ \frac{\cos^2 \frac{\pi(x-x_0)}{L} - \cos^2 \left( \frac{\pi a}{2L} \right)}{\epsilon^2/2^2} \right].
\end{aligned}
\quad (A2) $$

The distance between the rods of the grid is given by the mesh size $D$, the diameter of each rod of the grid is given by $a_{\parallel}$ and $a_{\perp}$ in the parallel and perpendicular direction, respectively (respect to the flow). The dimensions of the box are $L_{\perp} \times L_{\parallel} \times L_{\parallel} \cdot z_0$ is the position of the grid. See Fig. 7(a) for an illustration of the grid and the box. For all simulations $a_{\parallel} = a_{\perp} = 2\xi$, $\epsilon = 0.5\xi$, $z_0 = L_{\perp}/8$, and $V_0 = 20c^2$. See Table I for the values of the other parameters and the list of the different runs.

The initial condition containing the grid and the fluid at rest is obtained by using a Newton-Raphson method that ensures a perfect and clean initial condition [29]. However, if Eq. (A1) is abruptly integrated with $v_0 \neq 0$, a lot of sound is emitted. In order to minimize the initial emission of sound a local dissipative version of Eq. (A1) is used for the short times, typically up to $t_W \sim L_{\parallel}/c$, which corresponds to the time of sound waves take to traverse the box. The dissipative GP is a mix of the real Ginzburg-Landau and GPE, namely

$$ 
\begin{aligned}
i\hbar \left( \frac{\partial \psi}{\partial t} + \mathbf{v}_0 \cdot \nabla \psi \right) &= [1 - ia(\mathbf{x})] \left( -\frac{\hbar^2}{2m} \nabla^2 \psi + g|\psi|^2 \psi - \tilde{\mu} \psi + V_{\text{grid}}(\mathbf{x}) \psi \right), \\
\end{aligned}
\quad (A3) $$

where $a(\mathbf{x})$ is zero almost everywhere but close to the faces opposite to the grid, as displayed in Fig. 7(b). Furthermore, the Mach number $M = v_0/c$ is increased by steps as schematized in Fig. 7(c). After $t_W$ only the advective GP is used.

## APPENDIX B: NUMERICAL INTEGRATION

We use a standard pseudospectral code with Runge-Kutta of order 2 for time stepping. Note that the potential in Eq. (A3) has been carefully chosen periodic. To have a fully dealiased code, the scheme proposed in Ref. [33] is used (see Appendix of that reference). It consists on applying the Galerkin projector $\mathcal{P}_G$ as follows:

$$ 
\begin{aligned}
i\hbar \left( \frac{\partial \psi}{\partial t} + \mathbf{v}_0 \cdot \nabla \psi \right) &= \mathcal{P}_G \left[ -\frac{\hbar^2}{2m} \nabla^2 \psi + g \mathcal{P}_G [|\psi|^2] \psi - \tilde{\mu} \psi + V_{\text{grid}}(\mathbf{x}) \psi \right].
\end{aligned}
\quad (B1) $$

The Galerkin projector $\mathcal{P}_G$ takes a simple form in Fourier space: $\mathcal{P}_G[\psi(\hat{\mathbf{k}})] = \theta(k_{\text{max}} - |k|)\psi(\hat{\mathbf{k}})$, where $\theta$ is the Heaviside function and $k_{\text{max}} = N/3$ is chosen following the standard 2/3 rule for a quadratic nonlinearity [42]. When dealiasing is not performed as in Eq. (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 noncontrolled 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 dealiasing but with $k_{\text{max}} = N/4$ (corresponding to a cubic nonlinearity) at the price of wasting half of the resolution.

The Taylor-Green flow is prepared as in Ref. [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: INTERVORTEX DISTANCE

The intervortex distance is usually estimated as $\ell_{IV} = 1/\sqrt{\mathcal{L}/V}$, where $\mathcal{L}$ is the total vortex length and $V$ the volume of the system. For the grid, turbulence is not homogeneous, therefore an estimation of the volume is needed. Vortices mainly 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

$$ \Delta E^I(z) = \sum_{k>0} E_k^I(z), \quad (C1) $$
where $E_k^1(z)$ is the energy spectrum computed with the wave vectors perpendicular to the mean flow for a fixed value $z$. This quantity is displayed in Fig. 8 for run c1 at $t = 5.5$. The bulk is clearly visible. The effective volume containing the bulk is defined as $V_{\text{eff}} = L_{\perp}^2 \times L_{\text{eff}}$ with $L_{\text{eff}} = (\int_0^{L_{\parallel}} z \Delta E^1(z) dz)/(\int_0^{L_{\parallel}} \Delta E^1(z) dz)$:

$$L_{\text{eff}} = \frac{\int_{0}^{L_{\parallel}} z \Delta E^{1}(z) dz}{\int_{0}^{L_{\parallel}} \Delta E^{1}(z) dz}. \qquad (C2)$$

Finally, the intervortex distance is estimated as $\ell_{\text{IV}} = \sqrt{V_{\text{eff}}/\mathcal{L}}$.

<figure><img src="image_10.png" /><figcaption>FIG. 8. Spatial dependence of the incompressible kinetic energy $\Delta E^1(z)$. Run c1 at $t = 5.5$.</figcaption></figure>

[1] E. Henn, J. Seman, G. Roati, K. Magalhães, and V. Bagnato, Phys. Rev. Lett. 103, 045301 (2009).

[2] S. Serafini, M. Barbiero, M. Debortoli, S. Donadello, F. Larcher, F. Dalfovo, G. Lamporesi, and G. Ferrari, Phys. Rev. Lett. 115, 170402 (2015).

[3] G. P. Bewley, D. P. Lathrop, and K. R. Sreenivasan, Nature 441, 588 (2006).

[4] M. La Mantia and L. Skrbek, Phys. Rev. B 90, 014519 (2014).

[5] M. La Mantia, D. Duda, M. Rotter, and L. Skrbek, J. Fluid Mech. 717, R9 (2013).

[6] J. Gao, a. Marakov, W. Guo, B. T. Pawlowski, S. W. Van Sciver, G. G. Ihas, D. N. McKinsey, and W. F. Vinen, Rev. Sci. Instrum. 86, 093904 (2015).

[7] D. E. Zmeev, F. Pakpour, P. M. Walmsley, A. I. Golov, W. Guo, D. N. McKinsey, G. G. Ihas, P. V. E. McClintock, S. N. Fisher, and W. F. Vinen, Phys. Rev. Lett. 110, 175303 (2013).

[8] U. Frisch, Turbulence: The Legacy of A. N. Kolmogorov (Cambridge University Press, Cambridge, 1995).

[9] J. Maurer and P. Tabeling, Europhys. Lett. 43, 29 (1998).

[10] J. Salort, C. Baudet, B. Castaing, B. Chabaud, F. Daviaud, T. Didelot, P. Diribarne, B. Dubrulle, Y. Gagne, F. Gauthier et al., Phys. Fluids 22, 125102 (2010).

[11] W. F. Vinen and J. J. Niemela, J. Low Temp. Phys. 128, 167 (2002).

[12] P. M. Walmsley, A. I. Golov, H. E. Hall, A. A. Levchenko, and W. F. Vinen, Phys. Rev. Lett. 99, 265302 (2007).

[13] C. Nore, M. Abid, and M. E. Brachet, Phys. Rev. Lett. 78, 3896 (1997).

[14] Using a Padé approximation for the vortex profile (see Ref. [31]), the coherence length $\xi$ corresponds to a vortex core of size $a$ defined as $\rho(r=a) = 0.27\rho_{\infty}$, where $\rho_{\infty}$ is the value of the density far away from the vortex.

[15] J. Yepez, G. Vahala, L. Vahala, and M. Soe, Phys. Rev. Lett. 103, 084501 (2009).

[16] J. Salort, B. Chabaud, E. Lévêque, and P.-E. Roche, Europhys. Lett. 97, 34006 (2012).

[17] A. W. Baggaley and C. F. Barenghi, Phys. Rev. E 84, 067301 (2011).

[18] G. Comte-Bellot and S. Corrsin, J. Fluid Mech. 25, 657 (1966).

[19] S. B. Pope, Turbulent Flows (Cambridge University Press, Cambridge, 2000).

[20] S. R. Stalp, L. Skrbek, and R. J. Donnelly, Phys. Rev. Lett. 82, 4831 (1999).

[21] D. E. Zmeev, P. M. Walmsley, A. I. Golov, P. V. E. McClintock, S. N. Fisher, and W. F. Vinen, Phys. Rev. Lett. 115, 155303 (2015).

[22] D. I. Bradley, S. N. Fisher, A. M. Guénault, R. P. Haley, M. Kumar, C. R. Lawson, R. Schanen, P. V. E. McClintock, L. Munday, G. R. Pickett, M. Poole, V. Tsepelin, and P. Williams, Phys. Rev. B 85, 224533 (2012).

[23] D. I. Bradley, D. O. Clubb, S. N. Fisher, A. M. Guénault, R. P. Haley, C. J. Matthews, G. R. Pickett, V. Tsepelin, and K. Zaki, Phys. Rev. Lett. 95, 035302 (2005).

[24] C. Nore, M. Abid, and M. E. Brachet, Phys. Fluids 9, 2644 (1997).

[25] M. T. Reeves, T. P. Billam, B. P. Anderson, and A. S. Bradley, Phys. Rev. Lett. 110, 104501 (2013).

[26] G. Krstulovic and M. Brachet, Phys. Rev. Lett. 106, 115303 (2011).

[27] V. Shukla, M. Brachet, and R. Pandit, New J. Phys. 15, 113025 (2013).

[28] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevE.93.063104 for a movie showing the vortex nucleation close to the grid.

[29] C. Huepe and M.-E. Brachet, Physica D: Nonlin. Phenom. 140, 126 (2000).

[30] T. Frisch, Y. Pomeau, and S. Rica, Phys. Rev. Lett. 69, 1644 (1992).

[31] L. M. Pismen, Vortices in Nonlinear Fields: From Liquid Crystals to Superfluids, from Nonequilibrium Patterns to Cosmic Strings, Vol. 100 (Oxford University Press, Oxford, 1999).

[32] K. Sasaki, N. Suzuki, and H. Saito, Phys. Rev. Lett. 104, 150404 (2010).

[33] G. Krstulovic and M. E. Brachet, Phys. Rev. E 83, 066311 (2011).

[34] W. F. Vinen, Proc. Roy. Soc. London A: Math. Phys. Eng. Sci. 242, 493 (1957).

[35] M. La Mantia and L. Skrbek, Europhys. Lett. 105, 46002 (2014).

[36] A. C. White, C. F. Barenghi, N. P. Proukakis, A. J. Youd, and D. H. Wacks, Phys. Rev. Lett. 104, 075301 (2010).
[37] Z.-S. She and E. Leveque, Phys. Rev. Lett. **72**, 336 (1994).

[38] J. Salort, B. Chabaud, E. Lévêque, and P.-E. Roche, J. Phys. Conf. Ser. **318**, 042014 (2011).

[39] L. Boué, V. L'vov, A. Pomyalov, and I. Procaccia, Phys. Rev. Lett. **110**, 014502 (2013).

[40] V. Shukla and R. Pandit, arXiv:1508.00448.

[41] G. Lamporesi, S. Donadello, S. Serafini, F. Dalfovo, and G. Ferrari, Nat. Phys. **9**, 656 (2013).

[42] D. Gottlieb and S. A. Orszag, *Numerical Analysis of Spectral Methods: Theory and Applications*, Vol. 26 (Siam, Philadelphia, 1977).