# Active and finite-size particles in decaying quantum turbulence at low temperature

Umberto Giuriato, Giorgio Krstulovic

► To cite this version:

Umberto Giuriato, Giorgio Krstulovic. Active and finite-size particles in decaying quantum turbulence at low temperature. Physical Review Fluids, 2020, 5 (5), 10.1103/PhysRevFluids.5.054608. [https://doi.org/10.1103/PhysRevFluids.5.054608]

HAL Id: hal-03054341

https://hal.science/hal-03054341v1

Submitted on 28 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Active and finite-size particles in decaying quantum turbulence at low temperature

Umberto Giuriato* and Giorgio Krstulovic

Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Nice, France

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

(Received 26 February 2020; accepted 21 April 2020; published 29 May 2020)

The evolution of a turbulent tangle of quantum vortices in the presence of finite-size active particles is studied by means of numerical simulations of the Gross-Pitaevskii equation. Particles are modeled as potentials depleting the superfluid and described with classical degrees of freedom following a Newtonian dynamics. It is shown that particles do not modify the building-up and the decay of the superfluid Kolmogorov turbulent regime. It is observed that almost the totality of particles remains trapped inside quantum vortices, although they are occasionally detached and recaptured. The statistics of this process is presented and discussed. The particle Lagrangian dynamics is also studied. At large timescales, the velocity spectrum of particles is reminiscent of a classical Lagrangian turbulent behavior. At timescales faster than the turnover time associated with the mean intervortex distance, the particle motion is dominated by oscillations due to the Magnus effect. For light particles, a nonclassical scaling of the spectrum arises. The particle velocity and acceleration probability distribution functions are then studied. The decorrelation time of the particle acceleration is found to be shorter than in classical fluids, and related to the Magnus force experienced by the trapped particles.

DOI: 10.1103/PhysRevFluids.5.054608

## I. INTRODUCTION

When a fluid is stirred, energy is injected into the system exciting structures at different scales. In particular, in three-dimensional classical flows, the energy supplied at large scales is transferred toward small scales in a cascade process. Eventually, it reaches the smallest scales of the system, where dissipation acts efficiently. In the presence of a very large separation between the injection and dissipation scale, this cascade scenario proposed by Richardson leads to a fully developed turbulent state that can be described by the Kolmogorov phenomenology [1]. Kolmogorov turbulence is expected to be universal, and it is in fact commonly observed in nature, industrial applications, and in more exotic flows such as superfluids.

A superfluid is a peculiar flow, whose origin is a consequence of quantum mechanics. At finite temperature, a superfluid is considered to be a mixture of two components: the normal fluid, which can be described by the Navier-Stokes equations, and the superfluid component with zero viscosity [2]. At very low temperatures, the normal component can be neglected and the fluid becomes completely inviscid. As a consequence, an object moving at low velocities does not experience any drag from the fluid. However, when the object exceeds a critical velocity, quantum vortices are

<footnote>*Corresponding author: umberto.giuriato@oca.eu</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
nucleated [3,4]. Quantum vortices (or superfluid vortices) are the most fundamental hydrodynamical excitations of a superfluid. They are topological defects (and nodal lines) of the macroscopic wave function describing the system, and as a consequence their circulation is quantized. In superfluid helium, the core size of quantum vortices is of the order of 1 Å. Despite the lack of viscosity, quantum vortices can reconnect and change their topology (see, for instance, [5–8]), unlike classical (prefect) fluids.

When energy is injected in a low-temperature superfluid at scales much larger than the mean intervortex distance $ℓ$, a classical Kolmogorov regime is expected. Such a behavior has been observed numerically [9–11] and experimentally [12,13]. Indeed, at such scales the quantum nature of vortices is not important and the superfluid behaves like a classical fluid. At the scales of the order of $ℓ$ and smaller, the isolated nature of quantized vortices becomes relevant. The system keeps transferring energy toward small scales but through different nonclassical mechanisms [14]. An example of such mechanisms is the turbulent Kelvin wave cascade. Kelvin waves are helical oscillations propagating along quantum vortices, and the energy can be carried toward small scales thanks to nonlinear wave interactions. This energy cascade has been successfully described in the framework of weak-wave turbulence theory [15,16]. The resulting theoretical predictions have been observed numerically in vortex-filament and Gross-Pitaevskii numerical simulations [17–19].

Flow visualization is certainly a fundamental issue in every fluid dynamics experiment. Among the techniques that have been developed to sample a fluid, particle image velocimetry (PIV) and particle tracking velocimetry (PTV) are two of the most common methods [20]. The use of particles as probes has also been adapted to the study of cryogenic flows, in particular in superfluid helium $^4$He experiments [21], where micrometer-sized hydrogen and deuterium particles have been used. For instance, hydrogen ice particles have been successfully employed to visualize isolated or reconnecting vortex lines [22], as well as the propagation of Kelvin waves [23]. Moreover, the observation of power-law tails in the probability density of the particle velocity is an important difference with respect to classical turbulent states [24–26]. Similar deviations from classical behaviors have recently been reported also for the acceleration statistics [26,27]. Particles in such experiments typically have a size that can rise up to several microns, which is many orders of magnitude larger than the size of the vortex core in superfluid helium. For instance, the solidified hydrogen particles produced in the experiments [22,23] are slightly smaller than 2.7 µm, while in [25,26] their size is between 5 and 10 µm. Although it has been seen that particles unveil the dynamics of quantum vortices, it is not yet clear how much they affect the dynamics of quantum turbulent flows.

Several theoretical efforts have been made in the past decade in order to clarify what is the dynamics of particles in a superfluid and how particles interact with quantum vortices. For example, the vortex-filament (VF) method can be coupled with the classical hydrodynamical equations of a sphere, allowing us to study different specific problems. The interaction between one particle and one vortex has been addressed [28,29], as well the backreaction of tracers in a thermal counterflow [30,31]. In the context of finite-temperature superfluids, the spatial statistics of particles have been recently addressed in simulations of the Hall-Vinen-Bekarevich-Khalatnikov (HVBK) model [32].

Finally, since the work of Winiecki and Adams [4], particles described by classical degrees of freedom have been implemented self-consistently in the framework of the Gross-Pitaevskii (GP) equation [33–37]. Although the GP model is formally derived for dilute Bose-Einstein condensates, it is considered a general tool for the study of superfluid dynamics at very low temperature. Indeed, unlike the VF method or the HVBK model, it naturally contains quantum vortices as topological defects of the order parameter. It was found analytically and confirmed numerically that the GP model can reproduce the process of trapping of large active inertial particles by straight vortex lines [34], in accordance with hydrodynamical calculations [28,29]. In this framework, the interplay between many trapped particles and Kelvin waves has also been investigated [36].

In the present work, we study the influence of particles on quantum turbulent flows at very low temperature by using the GP model coupled with classical particles. In particular, we study the evolution of a free decaying superfluid turbulent vortex tangle loaded with finite-size active
particles. We consider spherical particles of different masses and having a diameter up to 20 core sizes. Such a size is about 1000 times smaller than that of solidified particles used in superfluid helium experiments. Nevertheless, it is slightly smaller than or comparable to the mean intervortex distance in our simulations, similar to current experiments. We also study the different regimes of the turbulent evolution from the Lagrangian point of view. The paper is organized as follows. In Sec. II we describe the Gross-Pitaevskii model coupled with classical particles. We also review the standard properties of the model and give the basic definitions used later to analyze the flow. We also describe the numerical method used in this work. Then, in Sec. III, we present our main results. In particular, in Sec. III A we address whether the presence of particles affects the scales of the flow at which Kolmogorov turbulence takes place. Section III B is devoted to a study of the particle dynamics inside the vortex tangle, their trapping by vortices, and their dynamics at scales larger and smaller than the intervortex distance. Particle velocity and acceleration statistics are then presented in Sec. III C. Finally, Sec. IV contains our conclusions.

## II. MODEL FOR PARTICLES IN A LOW-TEMPERATURE SUPERFLUID

### A. Gross-Pitaevskii equation coupled with particles

We describe a superfluid of volume $V$ at low temperature by using the complex field $\psi$, which obeys the GP dynamics. We consider $N_p$ particles in the system. Each particle is characterized by the position of its center of mass $\mathbf{q}_i$ and its classical momentum $\mathbf{p}_i$. The presence of a particle of size $a_p$ generates a superfluid depletion in a spherical region of radius $a_p$. This effect is reproduced by coupling the superfluid field with a strong localized potential $V_p$, which has a fixed shape and is centered at the position $\mathbf{q}_j(t)$.

All the particles considered have the same size, as well as the same mass $M_p$. The Hamiltonian of the system is given by

$$H = \int \left( \frac{\hbar^2}{2m} |\nabla \psi|^2 + \frac{g}{2} \left( |\psi|^2 - \frac{\mu}{g} \right)^2 + \sum_{i=1}^{N_p} V_p(|\mathbf{x} - \mathbf{q}_i|) |\psi|^2 \right) d\mathbf{x} + \sum_{i=1}^{N_p} \frac{\mathbf{p}_i^2}{2M_p} + \sum_{i<j}^{N_p} V_{\text{rep}}^{ij}, \quad (1)$$

where $m$ is the mass of the bosons constituting the superfluid, and $g$ is the nonlinear coupling constant between the bosons, related to the s-wave scattering length $a_s$ so that $g = 4\pi a_s \hbar^2/m$. The chemical potential is denoted by $\mu$. The particle interaction potential $V_{\text{rep}}^{ij}$ is responsible for short-range repulsion between particles, so that they behave as hard spheres and do not overlap. A detailed discussion on the inclusion of this short-range repulsion and the effect on the particle collisions in the model (1) can be found in [33]. The equations of motion that govern the superfluid field and the particle positions are obtained varying the Hamiltonian (1):

$$i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi + (g|\psi|^2 - \mu)\psi + \sum_{i=1}^{N_p} V_p(|\mathbf{x} - \mathbf{q}_i|)\psi, \quad (2)$$

$$M_p \ddot{\mathbf{q}}_i = - \int V_p(|\mathbf{x} - \mathbf{q}_i|) \nabla |\psi|^2 d\mathbf{x} + \sum_{j \neq i}^{N_p} \frac{\partial}{\partial \mathbf{q}_j} V_{\text{rep}}^{ij}. \quad (3)$$

This model has been successfully used to study vortex nucleation [4], trapping of particles by quantum vortices [34], and the interaction between particles trapped inside quantum vortices and Kelvin waves [36]. We denote by GP the Gross-Pitaevskii model without particles, and by GP-P the full coupled system (2) and (3).

In the case in which particles are absent, the chemical potential $\mu$ fixes the value of the ground state of the system $\psi_\infty = \sqrt{\rho_\infty/m} = \sqrt{\mu/g}$. Large-wavelength perturbations around this state are sound waves that propagate with the speed of sound $c = \sqrt{g\rho_\infty/m^2}$, while they become dispersive at length scales smaller than the healing length $\xi = \sqrt{\hbar^2/2g\rho_\infty}$.
The GP model describes a superfluid with zero viscosity. Using the Madelung transformation $\psi(\mathbf{x}) = \sqrt{\rho(\mathbf{x})}/m e^{i\frac{m}{\hbar}\phi(\mathbf{x})}$, the GP equation (2) is mapped into the continuity and Bernoulli equations of a superfluid of density $\rho$ and velocity $\mathbf{v}_s = \nabla\phi$. A superfluid flow is potential, but the phase is not defined at the nodal lines of $\psi(\mathbf{x})$. Therefore, the vorticity is concentrated along these filaments, which are the topological defects usually called quantum vortices. The effective size of the quantum vortex core coincides with the healing length $\xi$, and the contour integral of the superfluid velocity around a single vortex filament is the Feynman-Onsager quantum of circulation $\kappa = h/m = 2\pi\sqrt{2c\xi}$.

Using the Madelung transformation and the Helmholtz decomposition, the kinetic term of the superfluid energy density is decomposed into incompressible, compressible, and quantum energy [9]:

$$ 
\begin{aligned}
E_{\text{kin}}^{\text{GP}} &= \frac{\hbar^2}{2mV} \int |\nabla \psi|^2 d\mathbf{x} = E_{\text{kin}}^{1} + E_{\text{kin}}^{C} + E^{Q} \\
&= \frac{1}{2V} \int \left( [ (\sqrt{\rho} \mathbf{v}_s)^I ]^2 + [ (\sqrt{\rho} \mathbf{v}_s)^C ]^2 + \frac{\kappa^2}{4\pi^2} [\nabla \sqrt{\rho}]^2 \right) d\mathbf{x},
\end{aligned} 
\quad (4) 
$$

where $(\sqrt{\rho}\mathbf{v}_s)^I = \mathcal{P}_1[\sqrt{\rho}\mathbf{v}_s]$ and $(\sqrt{\rho}\mathbf{v}_s)^C = \mathbf{v}_s - (\sqrt{\rho}\mathbf{v}_s)^I$, the operator $\mathcal{P}_1[\cdot]$ being the projector onto the space of divergence-free fields. The other energies of the superfluid are the internal energy $E_{\text{int}} = (2V)^{-1} \int g(\rho/m - \mu/g)^2 d\mathbf{x}$, where the energy of the ground state is subtracted, and the interaction energy with the particles $E_{\text{GP}} = V^{-1} \int \sum_{i}^{N_p} V_p(|\mathbf{x} - \mathbf{q}_i|) \rho d\mathbf{x}$, so that the total energy is given by $E_{\text{tot}} = E_{\text{kin}}^{\text{GP}} + E_{\text{int}} + E_{\text{GP}}^{\text{GP}}$. From these definitions follow the corresponding energy spectra defined in terms of the Fourier transform of the fields [9].

## B. Numerical methods and parameters

In the simulations presented in this work, we solve the system (2) and (3) in a cubic periodic box of side $L = 341\xi$ with $N_c = 512^3$ collocation points by using a standard pseudospectral method. We use a fourth-order Runge-Kutta scheme for the time-stepping and the standard 2/3 rule for the dealiasing. In numerics, we fix $c = 1$ and $\psi_\infty = 1$.

To produce a homogeneous and isotropic tangle of quantized vortex lines, we impose an initial Arnold-Beltrami-Childress (ABC) flow, following the procedure described in [38]. In particular, we use a superposition of $k = 1 \times 2\pi/L$ and $k = 2 \times 2\pi/L$ basic ABC flows: $\mathbf{v}_{\text{ABC}} = \mathbf{v}_{\text{ABC}}^{(1)} + \mathbf{v}_{\text{ABC}}^{(2)}$, with

$$ \mathbf{v}_{\text{ABC}}^{(k)} = [B \cos(ky) + C \sin(kz)]\hat{\mathbf{x}} + [C \cos(kz) + A \sin(kx)]\hat{\mathbf{y}} + [A \cos(kx) + B \sin(ky)]\hat{\mathbf{z}}, \quad (5) $$

and the parameters $A = 0.5196$, $B = 0.5774$, and $C = 0.6351$. The basic ABC flow is a stationary (periodic) solution of the Euler equation with maximal helicity. The resulting wave function contains a tangle whose nodal lines follow the ABC vortex lines. The initial mean intervortex distance is $\ell(t = 0) \sim 25\xi$. As the flow is prepared by minimizing the energy, most of the energy of the system is in the incompressible part of the energy and resulting from the vortex configuration.

The ground state for the particles consists in a number of particles (we use $N_p = 200$ and $80$) of the same size and mass, randomly distributed in the computational box. Particles are initially at rest. This state is prepared using the imaginary-time evolution of Eq. (2). Then, the initial condition for the simulations is obtained by multiplying the wave function associated with the ABC flow and the wave function associated with the particle ground state. An example of an initial field containing particles is displayed in Fig. 1(d).

Because of the presence of a healing layer, the particle boundary is never sharp, independently of the functional form of the potential $V_p$. The superfluid field vanishes in the region where $V_p > \mu$, and at the particle boundary the fluid density passes from zero to the bulk value $\rho_\infty$ in approximately one healing length. The potential used to model each particle is a smoothed hat-function $V_p(r) = \frac{V_0}{2}(1 - \tanh[\frac{r^2 - \zeta^2}{4\Delta_a^2}])$, where the parameters $\zeta$ and $\Delta_a$ are set to model the
<figure><img src="image_3.png" /><figcaption>FIG. 1. Visualizations of the superfluid vortex tangle. Vortices are represented as isosurfaces in red of the density field ($\rho = 0.15\rho_{\infty}$), sound is rendered in blue, trapped particles in green, and free particles in purple. The upper row is without particles, the lower row is with 200 neutrally buoyant particles of radius $a_p = 4\xi$. (a,d) The ABC initial states. (b,e) The most turbulent regime ($t = 1.3T_L$). (c,f) A late time ($t = 8.1T_L$). $T_L$ denotes the large-eddy-turnover time (see the text).</figcaption></figure>

particle. Their values are listed in Table I. In particular, $\zeta$ fixes the width of the potential and it is related to the particle size, while $\Delta_a$ controls the steepness of the smoothed hat-function. The latter needs to be adjusted in order to avoid the Gibbs effect in the Fourier transform of $V_p$. Since the particle boundaries are not sharp, the effective particle radius is defined as $a_p = (3M_0/4\pi\rho_\infty)^{1/3}$, where $M_0 = \rho_\infty L^3 (1 - \int |\psi_p|^2 d\mathbf{x} / \int |\psi_\infty|^2 d\mathbf{x})$ is the fluid mass displaced by the particle and $\psi_p$ is the steady state with just one particle. Practically, given the set of numerical parameters $\zeta$ and $\Delta_a$, the state $\psi_p$ is obtained numerically with imaginary-time evolution and the excluded mass $M_0$ is measured directly. Particles attract each other by a short-range fluid mediated interaction

<table><caption>TABLE I. Simulation parameters.</caption><thead><tr><th>Run</th><th>N<sub>p</sub></th><th>a<sub>p</sub></th><th>M</th><th>ζ</th><th>Δ<sub>a</sub></th><th>V<sub>0</sub>/μ</th><th>γ/μ</th></tr></thead><tbody><tr><td>I</td><td>0</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>II</td><td>200</td><td>4.0ξ</td><td>0.125</td><td>1.5ξ</td><td>1.2ξ</td><td>20.0</td><td>1.4 × 10<sup>-4</sup></td></tr><tr><td>III</td><td>200</td><td>4.0ξ</td><td>0.25</td><td>1.5ξ</td><td>1.2ξ</td><td>20.0</td><td>1.4 × 10<sup>-4</sup></td></tr><tr><td>IV</td><td>200</td><td>4.0ξ</td><td>1.0</td><td>1.5ξ</td><td>1.2ξ</td><td>20.0</td><td>1.4 × 10<sup>-4</sup></td></tr><tr><td>V</td><td>200</td><td>4.0ξ</td><td>2.0</td><td>1.5ξ</td><td>1.2ξ</td><td>20.0</td><td>1.4 × 10<sup>-4</sup></td></tr><tr><td>VI</td><td>80</td><td>10.0ξ</td><td>1.0</td><td>8.0ξ</td><td>2.0ξ</td><td>20.0</td><td>5.8 × 10<sup>-4</sup></td></tr><tr><td>VII</td><td>200</td><td>10.0ξ</td><td>0.125</td><td>8.0ξ</td><td>2.0ξ</td><td>20.0</td><td>5.8 × 10<sup>-4</sup></td></tr><tr><td>VIII</td><td>200</td><td>10.0ξ</td><td>0.25</td><td>8.0ξ</td><td>2.0ξ</td><td>20.0</td><td>5.8 × 10<sup>-4</sup></td></tr><tr><td>IX</td><td>200</td><td>10.0ξ</td><td>1.0</td><td>8.0ξ</td><td>2.0ξ</td><td>20.0</td><td>5.8 × 10<sup>-4</sup></td></tr></tbody></table>
[33,35], thus we use the repulsive potential $V_{rep}^{ij} = \gamma(2a_p/|\mathbf{q}_i - \mathbf{q}_j|)^{12}$ in order to avoid an overlap between them. The functional form of $V_{rep}^{ij}$ is inspired by the repulsive term of the Lennard-Jones potential, and the prefactor $\gamma$ is adjusted numerically so that the interparticle distance $2a_p$ minimizes the sum of $V_{rep}^{ij}$ with the fluid-mediated attractive potential [33,35]. We express the particle mass as $M_p = MM_0$, where $M_0$ is the mass of the superfluid displaced by the particle. Namely, heavy particles have $M > 1$ and light particles have $M < 1$. In Table I all the parameters for the particles used in the simulations presented in this work are reported. In the following, we will refer to each simulation specifying the size and the mass of the particles used.

Note that although the model (1) is a minimal model for implementing particles in the GP framework, we cannot add to the system an arbitrary number of particles. Indeed, since particles have a finite size, they occupy a volume at the expense of the superfluid field, and packing effects could become important if the filling fraction is too high. Moreover, the potential $V_p$ must be updated at each time step, which is numerically costly. Finally, note that the evaluation of the force term (3) acting on particles requires us to know the value of the fields at intermesh points. When the number of particles in the simulation is not large, the force $f_i^{GP}(\mathbf{q}_i) = -(V_p * \nabla\rho)[\mathbf{q}_i]$ (3) can be computed with spectral accuracy using a Fourier interpolation. Such a method has been used in [34–36], where the particle dynamics is extremely sensitive. In this work, the use of a Fourier interpolation for each particle is numerically unaffordable, due to the large number of particles involved and the resolutions used. Instead, we use a fourth-order B-spline interpolation method, which has been shown to be highly accurate with a reduced computational cost [39] and particularly well adapted for pseudospectral codes. Indeed, the use of a Fourier interpolation to evaluate the three-dimensional force for $N_p$ particles requires $\sim 3N_pN_c$ operations and evaluations of complex exponentials ($N_c = 512^3$ in the present work). Such a cost quickly becomes too expensive at high resolutions and/or a large number of particles. On the contrary, B-spline interpolation requires just one fast Fourier transform of a field per component, and an interpolation using only four neighboring grid points per dimension [39]. Such a scheme saves a factor $\sim N_p$ of computational cost compared to Fourier interpolation. Note that in the previous discussion, we have not taken into consideration parallelization issues, where local schemes (B-splines) are much more advantageous than global ones (Fourier transforms). Nevertheless, some issues with physical quantities at small scales arising from the B-spline interpolation are discussed in the Appendix.

### III. PARTICLES IMMERSED IN A TANGLE OF SUPERFLUID VORTICES

Superfluid turbulence in the context of the GP model has been studied extensively [9,11,38,40,41]. In general, quantum turbulence develops from an initial state with a vortex configuration where the incompressible kinetic energy is mainly contained at large scale. During the evolution, vortex lines move, interact among themselves, and reconnect, creating complex vortex tangles. Through this process, sound is produced and incompressible kinetic energy is irreversibly converted into quantum, internal, and compressible kinetic energy. Eventually, the compressible energy produced in the form of acoustic fluctuations starts to dominate, thermalizes, and acts as a thermal bath providing an effective dissipation acting on the vortices. As a consequence, vortices shrink and eventually disappear through mutual friction effects following Vinen’s decay law [19,42]. In particular, it has been shown that the decrease of the incompressible kinetic energy behaves in a similar manner to decaying classical turbulence [9]. To make a connection with decaying classical Kolmogorov turbulence, the incompressible energy dissipation or dissipation rate is usually defined in the context of GP turbulence as

$$ \epsilon = - \frac{dE_{\text{kin}}^{\text{I}}}{dt}. \qquad (6) $$

As in decaying Navier-Stokes turbulence, in GP the most turbulent stage is achieved around the time when this quantity is maximal. About this time, the classical picture holds and the incompressible
<figure><img src="image_4.png" /><figcaption>FIG. 2. (a) Time evolution of the superfluid energy components in the cases with no particles (dashed line), 200 small particles (dotted line), 200 large particles (solid line), and 80 large particles (dash-dotted line). (b) Incompressible energy dissipation rate for different numbers of particles with different sizes and different masses (solid lines). Dash-dotted horizontal lines of the corresponding colors indicate the value of the maximum of dissipation, obtained averaging over the shaded region. The dissipation is expressed in units of its maximum $\epsilon_{max}$ in the case without particles.</figcaption></figure>

energy spectrum satisfies the Kolmogorov prediction

$$E_{\text{kin}}^{I} = C\epsilon^{2/3}k^{-5/3},$$

where C is the Kolmogorov constant, the value of which has been found to be close to 1 in GP turbulence [11,38,41].

The first purpose of this work is to check whether and to what extent the presence of particles in the system modifies Kolmogorov turbulence. We add to the ABC initial condition a number of randomly distributed particles and let the system evolve under the dynamics (2) and (3). In Figs. 1(a), 1(b) and 1(c), the three stages of the evolution (initial condition, turbulent vortex tangle, and residual filaments in a bath of sound, respectively) are visualized in the case of 200 neutrally buoyant particles of radius $4\xi$. See the supplemental material [43] for movies of this simulation and others with particles of a different size. Trapped particles by vortices are displayed in green, whereas free ones are displayed in purple. The algorithm to distinguish a trapped particle from a free one is based on the circulation around it and it is discussed in Sec. III B.

In Fig. 1 we observe that the building up and decay of the turbulent tangle is not strongly modified by the presence of particles. Moreover, it can be noticed how during the first stages of the evolution of the system the majority of particles gets trapped into the vortices. At zero temperature, as there is no normal component in the flow, no drag is experienced by the particles and their motion is completely driven by the pressure gradients. As a consequence, they are attracted by quantum vortices [28,34,44]. During the turbulent regime, violent and strongly nonlinear events like reconnections dominate the vortex dynamics and the flow evolution. A fundamental question is whether and how much the hydrodynamical attraction between vortices and particles is sufficient to keep them attached to the filaments. Indeed, since quantum vortices are actually the main actors of turbulence in superfluid, if particles are really able to follow them in this regime, it is a good indication that they are suitable for use as probes.

In the following subsection, we will quantitatively study the effect of particles on quantum turbulent flows. We will first focus on the large scales of the flow, where Kolmogorov turbulence takes place. Then the particle dynamics and their statistics will be addressed.

## A. The effect of particles on Kolmogorov superfluid turbulence

We shall start our analysis by comparing the temporal evolution of global quantities. In Fig. 2(a) the time evolution of the different components of the energy is displayed. Times are expressed in
<figure><img src="image_5.png" /><figcaption>FIG. 3. (a) Time evolution of the mean intervortex distance for different numbers of particles of different sizes and different masses. (b) Incompressible energy spectrum for different numbers of particles of different sizes and different masses. Inset: Compensated incompressible energy spectrum. Solid lines refer to particles of size $a_p = 4\xi$, dashed lines refer to particles of size $a_p = 10\xi$. The dotted line is the classical scaling $\epsilon_{max}k^{-5/3}$. The spectrum is computed averaging over times in the shaded region.</figcaption></figure>

units of the large-eddy-turnover time defined as $T_L = L/2v_{rms}$, where $v_{rms} = \sqrt{2E_{kin}^I(t=0)/3}$ is the root-mean-square velocity associated with the initial vortex tangle, and $L/2$ is its characteristic length scale. We compare the case in which no particles are present in the flow to the cases having particles of different sizes and of relative mass $M=1$. The net transfer of incompressible energy toward compressible, quantum, and internal energy is qualitatively unchanged in the various cases. The only difference is a slightly lower value of the incompressible energy in the case of large particles, in favor of the internal energy of the superfluid. Such an effect is more evident if the number of large particles is increased, and could be related to an increment of the filling fraction $\Phi$, namely the fraction of the total volume occupied by the particles. In fact, for $N_p = 200$ particles of radius $a_p = 4\xi$ the filling fraction is $\Phi = 0.1\%$, for $N_p = 80$ particles of radius $a_p = 10\xi$ it is $\Phi = 0.8\%$, and for $N_p = 200$ particles of radius $a_p = 10\xi$ we have $\Phi = 2.1\%$. The kinetic and repulsion energies of the particles, as well as the particle-vortex interaction $E_p^{GP}$, are negligible compared with the other energies throughout the duration of the simulations (data not shown).

The dissipation rate of the incompressible kinetic energy is reported in Fig. 2(b) for particles of different masses and different sizes. The dissipation increases in the early stages when the energy begins to be transferred to the smaller scales, it reaches a maximum when all the scales are excited, and then it starts to decay since no forcing is sustaining the turbulence. We observe that the evolution of the dissipation is clearly not significantly modified by the presence of particles. In particular, the value of the maximum of dissipation, which is the signature of the most turbulent state reached by the tangle, is slightly lower only in the case in which many large particles are moving in the system. In particular for this case, it is about 90% of $\epsilon_{max}$, the value measured in the case with no particles. The shaded region in Fig. 2(b) represents the most turbulent time of the simulations. We consider that in this short stage the system is in a quasisteady state, and we perform the temporal average of certain physical quantities in order to improve statistical convergence.

Another important quantity that is not affected much by the interplay between tangle and particles is the mean intervortex distance $\ell$, whose time evolution is reported in Fig. 3(a). The mean intervortex distance is then estimated as $\ell = \sqrt{V/L_v}$, where $L_v$ is the total vortex length in the system. This latter is estimated using the method introduced in [9], where $L_v$ is shown to be related to the proportionality constant between the incompressible momentum density $J^I(k)$ of the flow and the spectrum of a two-dimensional point-vortex $J_{vort}^{2D}(k)$:

$$ \frac{L_v}{2\pi} = \frac{\sum_k J^I(k)}{\int J_{vort}^{2D}(k) dk} \qquad (7) $$
<figure><img src="image_6.png" /><figcaption>FIG. 4. Closeup of the superfluid vortex tangle at the early stage of the simulation (upper row: $t = 0.27T_L$) and during the turbulent regime (lower row: $t = 1.1T_L$) for the cases with no particles [left column (a),(d)], small particles [central column (b),(e): $a_p = 4\xi$], and large particles [right column (c),(f): $a_p = 10\xi$]. Vortices are represented as isosurfaces of the density field ($\rho = 0.15\rho_\infty$) and rendered in red, sound is rendered in blue, trapped particles in green, and free particles in purple.</figcaption></figure>

The spectra of momentum densities are the angle average of the norm in Fourier space of the momentum density $J = \rho v_s$, and the incompressible part is obtained projecting onto the space of divergence-free fields. We have checked the validity of this formula by using the vortex filament tracking method described in [45] at some checkpoints.

In the turbulent regime, where the dissipation gets its maximum, the total length of the entangled vortices is also larger by a factor 4 compared to the initial condition, and the distance between the filaments is minimum. The value $\ell_{\min} \sim 14\xi$ of the intervortex distance in this regime will be used as a characteristic small length scale of the Kolmogorov turbulent regime. Such length is smaller than the diameter of the largest particles considered ($2a_p = 20\xi$), but nevertheless this has no appreciable repercussions on the behavior of the observables studied. Furthermore, as shown in Fig. 3(c), the scaling of the incompressible energy spectrum $\mathcal{E}^I(k)$ averaged around the maximum of dissipation is unaltered by particles in the system. Figure 3(b) displays the incompressible energy spectrum. It is apparent that the scaling of the spectrum is always compatible with classical turbulence at scales larger than the intervortex distance, and the way in which the energy is accumulated at smaller scales is not modified by the particles. In the inset of Fig. 3(b), the spectrum is compensated by the Kolmogorov prediction $\mathcal{E}^I(k) = C\epsilon_{\max}^{2/3}k^{-5/3}$ for classical hydrodynamic turbulence. The dotted horizontal black line shows that the value of the constant $C$ in the Kolmogorov law is a number of order 1 for superfluid turbulence.

The only appreciable difference observed between the case with and without particles is that in the early stages of the evolution, the trapping of particles perturbs the vortex filaments and excites Kelvin waves. A comparison between the volume renderings can be seen in the upper row of Fig. 4. Such perturbations propagate during the evolution of the tangle. At the times when turbulence is developed, the details of the vortex configurations are completely different (see the lower row of Fig. 4). Nevertheless, the statistical properties of the system in this regime remain unchanged. We stress that the intervortex distance in quantum turbulence experiments lies typically in the range 10–100 µm, which is equal to or slightly larger than the particle size [24,25,27]. In this sense, the simulations presented here are compatible with the experimental parameters. They thus support
<figure><img src="image_7.png" /></figure>

FIG. 5. (a) Fraction of trapped particles as a function of time for different numbers of particles of different sizes and different masses. Inset: The same for longer time in the case of 200 neutrally buoyant particles of size $a_p = 4\xi$. (b) Comparison between the fraction of multiply trapped particles as a function of time for neutrally buoyant particles. (c) Volume rendering of large particles ($a_p = 10\xi$) multiply trapped by quantum vortices. Vortices are rendered in red, sound in blue, particles in green. (d) Probability density function of the continuous time spent by particles inside vortices for different species of particles. The dotted blue line corresponds to the same simulation of blue circles (particles with size $a_p$ and mass $M=1$) but averaged over the full simulation times). Inset: Absolute value of the circulation around a single particle of size $a_p = 4\xi$ and mass $M=1$ as a function of time. The PDF is computed averaging over times in the shaded region.

the belief that active particles have effectively no influence on the typical development and decay of quantum turbulence. This numerical fact helps to validate past and future experiments that use particles as probes of superfluids.

On the other hand, because of the lack of a Stokes drag in the system, particles cannot be treated as simple tracers of the superfluid velocity $v_s$. Nevertheless, if they remain trapped inside the vortices they can track the evolution of the vortex filaments, which are the structures that effectively become turbulent. With the purpose of characterizing this scenario, in the next subsection we investigate the motion of particles once they are immersed in a tangle of quantum vortices.

## B. Motion of particles in the superfluid vortex tangle

Looking at the time evolution of the vortex tangle (see Fig. 1 and movies in the supplemental material), the first thing that is apparent is how particles quickly get trapped into vortex filaments. This dynamics is expected and it has been studied in the case in which vortices move slowly [34]. It is a consequence of the pressure gradients. However, it is less obvious if such behavior remains dominant when turbulence take place and reconnections become frequent.

We study the evolution of particles and compute whether they are free or trapped by vortices. The temporal evolution of the fraction of trapped particles is displayed in Fig. 5(a) for all runs. This measurement is made by computing the circulation $\Gamma = \oint_C v_s \cdot d\mathbf{x}$ of the superfluid velocity $v_s$
along contours $C$ encircling each particle, and counting for which particles it is different from zero. Specifically, we compute the circulation along many parallel square contours of side $2(a_p + \Delta_x)$ around each particle, where $\Delta_x$ is the grid spacing. If the circulation around at least one of these contours is different from zero, the particle is considered as trapped [46]. For practical reasons, due to the parallelization of the numerical code, we consider only contours perpendicular to the $z$ axis of the computational box. As a consequence, the protocol is not able to grasp vortices that are crossing the particles exactly on a plane perpendicular to the $z$ axis. This means that our estimation of the fraction of trapped particles is effectively a lower bound. However, it should be noticed that this pathological situation is an extremely rare situation that does not change the conclusions of our analysis.

In the initial condition the particles are placed randomly in the computational box. It happens then that some of them are already positioned inside a vortex. In the case of particles with a size comparable to the intervortex distance, the majority of particles are in this situation. In the first stages of the evolution of the flow, the number of trapped particles increases rapidly until it becomes stationary always at times much smaller than one $T_L$. The time needed to reach a stationary state depends slightly on the mass of the particles, as well as the fraction of trapped particles once a steady regime is reached. The steady value of $N_p^{\text{trap}}/N_p$ is between 80% and 90% for small particles ($2a_p < \ell$), while on average the totality of particles of size $2a_p \sim \ell$ is found to be trapped by vortices, independently of the filling fraction. When the system reaches the most turbulent regime (indicated by the shaded region), the fraction of trapped particles does not undergo any appreciable changing. In the inset of Fig. 5(a), $N_p^{\text{trap}}/N_p$ is also shown for late times in the case of small particles of relative mass $M = 1$. It manifestly remains stable. This means that even when the density of vortex lines is decaying (along with the intensity of turbulence), the particles stay trapped inside vortices. Note that in this work we are dealing with homogeneous and isotropic decaying quantum turbulence at low temperature. We mention that the fraction of trapped particles measured in thermal counterflow simulated by means of the VF method is lower than the one observed here [31].

The circulation around each superfluid vortex filament is equal to a single quantum of circulation $\kappa$. As a consequence, measuring the circulation along a closed line $C$ allows us to count the number of filaments in the region delimited by the line, provided that the quanta of circulation around every filament have the same sign. This is true also if the vortices are trapping particles, because their topological nature does not change. In Fig. 5(b) we show again the fraction of trapped particles, but now separating the number of particles trapped by multiple vortices. It turns out that at least the 5–10 % of the particles with size $2a_p \sim \ell$ are always attached to at least two different filaments. Sometimes even more vortices pass simultaneously through the same particle, as can be visualized in the volume plot of Fig. 5(c).

Once a particle is trapped by a vortex, it can experience violent events, for instance during vortex reconnections. In such circumstances, such a particle could be detached and expelled from the vortex until it will eventually get trapped by another vortex of the tangle. We compute the probability density function (PDF) of the continuous time intervals $\Delta t_{\text{trap}}$ spent by the particles inside the vortices regime. The PDFs for particles of different sizes and masses are displayed in Fig. 5(d). For all the species of particles examined, the probability distribution seems to follow roughly a power-law scaling in time $\sim (\Delta t_{\text{trap}})^{-\alpha}$, with $\alpha \sim 1.67$. The PDF certainly vanishes much slower than an exponential decay at large $\Delta t_{\text{trap}}$, which would typically result from a standard escape problem over energy barriers. We checked that the intermittency of the circulation and the shape of the trapping time PDF are not characteristic of the most turbulent regime, since they persist also at the late times of the simulations [see the dotted blue line in Fig. 5(d)]. Therefore, many particles spend a time at least of the order of the simulation time ($\sim 10T_L$) inside a vortex filament, i.e., the typical escape time from the vortices is virtually infinite. This observation is exemplified in the inset of Fig. 5(d), where the evolution of the circulation around a single-small neutral particle is reported (the qualitative behavior is the same for the other particles). It is also clear that the time spent by the particles with zero circulation around them (namely free from vortices) is short. Since
<figure><img src="image_8.png" /><figcaption>FIG. 6. Frequency spectrum of the particle velocity for particles of different masses and different sizes, compensated with the prediction for the Lagrangian spectrum in classical turbulence $\propto \epsilon/\omega^2$: (a) small particles with $a_p = 4\xi$; (b) large particles with $a_p = 10\xi$. The dash-dotted gray line is the frequency spectrum of a single small particle trapped in a straight vortex slightly perturbed. Dotted lines of corresponding colors are the prediction for the particle natural frequency $\Omega_p$. The dashed red line is the scaling due to vortex reconnection or Kelvin waves $\propto |\omega|^{-1}$. The dashed golden line is the spectrum evaluated at late times in the simulation ($6T_L < \tau < 7T_L$).</figcaption></figure>

we established that particles immersed in a tangle spend most of the time inside vortex filaments, in the following we study their motion once they get trapped.

At large scales, the vortex tangle seems to behave as a classical hydrodynamic turbulent system. Therefore, the first natural question is whether the particles can trace such large-scale fluctuations. In classical turbulence, it is well known that the Lagrangian velocity spectrum scales as

$$ \langle |\hat{v}_p(\omega)|^2 \rangle = B \epsilon \omega^{-2}, \qquad (8) $$

where B is a constant of order unity and $\hat{v}_p(\omega)$ is the Fourier transform of the Lagrangian particle velocity $v_p(t)$ [47,48]. Such scaling is valid in the inertial range $2\pi/T_L \ll \omega \ll 2\pi/\tau_\eta$, where $\tau_\eta$ is the Kolmogorov timescale. In our case, we build an analog of the Kolmogorov time scale under the assumptions that the dissipation rate $\epsilon_{max}$ is the only important physical parameter in the classical turbulence regime and that the Kolmogorov turbulent cascade ends at the intervortex distance $\ell_{min}$. Therefore, we define the smallest timescale of the classical turbulence regime as $\tau_\ell = (\ell_{min}^2/\epsilon_{max})^{1/3}$, and we expect classical turbulent phenomenology to hold for times $\tau_\ell \ll t \ll T_L$. In Fig. 6, the measurement of the frequency spectrum of the particle velocity $\langle |\hat{v}_p(\omega)|^2 \rangle = \langle |\int \dot{q}(t)e^{-i\omega t} dt|^2 \rangle$ during the turbulent regime is shown for different species of particles, compensated with the classical scaling $\epsilon_{max}\omega^{-2}$. Note that the average that defines the spectrum is meant over different realizations. In numerics we average over all the particle trajectories during the turbulent regime. At frequencies $\omega < \tau_\ell/2\pi$, the spectra approach a plateau of value 1, confirming that particles sample well the flow and their behavior is described by the standard classical turbulence picture at large scales. Note that the classical temporal inertial range of our simulations is pretty small, since $T_L \sim 5\tau_\ell$. For comparison, we also present the velocity spectrum of a particle of size $a_p = 4\xi$ and mass $M = 1$, computed in a temporal window at much later times, when Kolmogorov turbulence has decayed and only a few vortices are left. Note that a $\omega^{-2}$ scaling of the Lagrangian velocity spectrum has also been observed in numerical simulations of the vortex filament model [49], although not in the Kolmogorov inertial range and not related to the energy dissipation rate nor to Kolmogorov turbulence.

As expected, in our simulations no Kolmogorov scaling is observed at small timescales. Indeed, one of the most striking features of quantum turbulence is the crossover between the classical Kolmogorov regime and the physics taking place at scales smaller than the mean intervortex distance. Unlike classical turbulence (see, for instance, [47]), there is still a nontrivial scaling at
timescales shorter than $\tau_\ell$. Such a difference is a consequence of the quantum nature of the system, here manifested by the presence of quantized vortices.

When a particle is trapped by a vortex, the superfluid flow turns around it. As a consequence, while the particle moves, it experience a Magnus force. This lift force is simply expressed as $\mathbf{F}_{\text{Magnus}} = \frac{3}{2}\rho_{\infty}a_{\text{p}}\mathbf{\Gamma} \times (\dot{\mathbf{q}} - \mathbf{v}_{\text{s}})$, where the circulation vector $\mathbf{\Gamma}$ is oriented along the vortex filament, and the superfluid velocity $\mathbf{v}_{\text{s}}$ contains the contributions of the mean flow and the vortex motion [36,50]. The Magnus effect induces a precession of the particle about the filament with the characteristic angular velocity

$$ \Omega_{\text{p}} = \frac{3}{2} \frac{\rho_{\infty} a_{\text{p}}}{M_{\text{p}}^{\text{eff}}} \Gamma, \qquad (9) $$

where the particle effective mass $M_{\text{p}}^{\text{eff}} = M_{\text{p}} + \frac{1}{2}M_0 = (M + \frac{1}{2})M_0$ takes into account the added mass effect due to the mass of the superfluid displaced by the particle $M_0$. As mentioned in [36], for current experiments with hydrogen particles in superfluid helium, this frequency is of order 10–100 Hz. If the Magnus force is the main force acting on a trapped particle, the Newton equation $M_{\text{p}}^{\text{eff}}\ddot{\mathbf{q}} = \mathbf{F}_{\text{Magnus}}$ implies the following expression for the frequency spectrum of the particle velocity:

$$ \langle |\hat{\mathbf{v}}_{\text{p}}(\omega)|^2 \rangle = \frac{\Omega_{\text{p}}^2}{\Gamma^2(\omega - \Omega_{\text{p}})^2} \langle |\mathbf{\Gamma} \times \hat{\mathbf{v}}_{\text{s}}(\omega)|^2 \rangle. \qquad (10) $$

Independently of the external superfluid velocity, the expression (11) predicts that the spectrum $\langle|\hat{\mathbf{v}}_{\text{p}}(\omega)|^2\rangle$ must be peaked around the natural frequency of trapped particles $\omega = \Omega_{\text{p}}$. Such behavior has been studied in detail in the case of particles trapped inside slightly perturbed straight vortex filaments [36]. The spectrum of this simple configuration is also reported for comparison in Fig. 6(a) for a small particle of relative unit mass. A clear bump in the frequency spectrum, corresponding to $\Omega_{\text{p}}$, is still visible when particles are immersed in a complex quantum vortex tangle. For the large particles, the presence of a peak is less evident because the natural frequency is lower, and therefore a longer sampling (in time) would be necessary to resolve it properly ($2\pi/\Omega_{\text{p}} = 0.7T_L$ for the particles of size $a_{\text{p}} = 10\xi$ and mass $M=1$). Moreover, as large particles are multiply trapped by many vortices, the resulting motion is certainly more complex than a precession with a single characteristic angular frequency of one single vortex. The broadness of the peak around the Magnus frequency for the small particles in Fig. 10(a) could also be related to this fact.

At small timescales, a different scaling of the velocity spectrum is observed for the light particles, now in agreement with $\langle|\hat{\mathbf{v}}_{\text{p}}(\omega)|^2\rangle \propto |\omega|^{-1}$. This behavior is consistent with the fact that at scales smaller than the intervortex distance, the typical velocities of a superfluid turbulent tangle are supposed to scale as $v_{\text{fast}}(t) \propto \sqrt{\kappa/|t - t_0|}$, because the circulation becomes the only relevant physical parameter, and the motion of vortices is dominated by their mutual advection and reconnections. In this scenario, if particles are sufficiently light to be able to follow the fast vortex dynamics, we can substitute $\langle|\hat{\mathbf{v}}_{\text{p}}(\omega)|^2\rangle \sim \hat{v}_{\text{fast}}^2(\omega) \propto \kappa|\omega|^{-1}$. Another effect that could contribute to the same result is the attraction of particles by the vortices, since the scaling in time of the particle-vortex distance is the same as that of vortex reconnection [34]. Note that for the heaviest particles, such fast scaling is absent since their reaction is probably too slow to be sensible to the fast fluctuations of the tangle.

### C. Particle velocity and acceleration statistics

Unlike classical turbulence, where the statistics of the one-point particle velocity $v$ is known to be Gaussian [1], experiments in superfluid helium using hydrogen and deuterium particles as tracers have reported long tails, with a $v^{-3}$ power-law scaling in their velocity distribution [24–26]. Such scaling has been related to the singular velocity field of quantized vortices [51,52]. At low temperatures, as Stokes drag is negligible, particles should not move with the superfluid flow
<figure><img src="image_9.png" /><figcaption>FIG. 7. (a) Probability density function of the single-component particle velocity, for different species of particles. The dotted golden line is the Eulerian velocity field $\nabla\phi$, corresponding to the simulation without particles at the time $1.4 T_L$. The data for the particles are averaged in time between $t = 1.2T_L$ and $1.6T_L$. Inset: Standard deviation of the particle velocity as a function of the particle mass. (b) The same as (a) but with the velocities normalized by the standard deviation $\sigma_v$. Dotted lines are Gaussian, dash-dotted line is a power-law scaling $|v_i|^{-3}$.</figcaption></figure>

and such scaling can be understood as a consequence of quantum vortex reconnections sampled by trapped particles [7,24]. Furthermore, in Ref. [25], by using particle tracking velocimetry in counterflow turbulence, it was shown that while varying the sampling scale, the velocity PDFs continuously change from Gaussian statistics to power-law tails, the crossover taking place at scales of the order of the intervortex distance. In this final subsection we present measurements of particle velocity and acceleration statistics within the GP-P model.

We start the discussion by presenting the Eulerian velocity field. Formally, the velocity of the superfluid is simply given by $\nabla\phi$. This field contains the density fluctuations, as well as the divergence of the vortex velocity flow close to its core. This divergence leads to the well-observed $v^{-3}$ scaling of velocity PDF [51,53,54]. The PDF of $\nabla\phi$ is displayed in Fig. 7. We turn now to analyze the particle velocity PDFs. We compute the velocity PDFs for all runs in the turbulent regime. Data are filtered with a Gaussian convolution in order to smooth out the noisy oscillations at frequencies $\omega < \omega_{noise} = 50 (2\pi/\tau_l)$ (see Appendix). In Fig. 7 the PDF of the single-component velocity is plotted for all the species of analyzed particles. In Fig. 7(a), velocities are expressed in terms of the speed of sound $c$, whereas in Fig. 7(b) they are normalized by their root-mean-squared values. The root-mean-squared values are displayed in the inset of Fig. 7(a) as a function of the mass for the two particle sizes. It is apparent from Fig. 7(b) that the particle statistics exhibits a Gaussian distribution. Note that Gaussian velocity statistics was also observed in thermal counterflow simulations of the vortex filament method with tracers particles [30]. The absence of power-law tails could be a consequence of weak statistical sampling of large velocity fluctuations due to the low number of particles present in the system and/or by compressible effects of the GP model. We will comment more about this in Sec. IV.

We would like to remark here that high-frequency fluctuations are strongly sensitive to numerical artifacts. In the Appendix, inspired by the experimental results of Ref. [25], we have computed the velocity PDFs of the velocity fluctuations filtered at a given frequency $\omega_c$. The frequency was varied from values lower to larger than $2\pi/\tau_l$. For one simulation we have compared two different interpolation methods to evaluate the force term in Eq. (3) needed to drive the particles. It turns out that for the fourth-order B-spline method, the velocity PDFs start to develop tails while the filtering scale is varied, eventually leading to a $v^{-3}$ scaling. However, when using Fourier interpolation, which is an exact evaluation (up to spectral convergence of the pseudospectral code) of the force term, the PDFs do not develop any tail and remain Gaussian. We have decided to keep this example with spurious numerical effects in the Appendix, as it might be useful for future numerical studies
<figure><img src="image_10.png" /><figcaption>FIG. 8. (a) Probability density functions of the single-component particle acceleration. (b) Probability density functions of the norm of the particle acceleration. The dotted line is a Gaussian, the dashed line is a $\chi_3$ distribution, and the dash-dotted line is an exponential tail $e^{-|a|/\sigma_a}$. Inset: Probability density functions of the natural logarithm of the norm of the particle acceleration. The dashed golden line is a log-normal distribution.</figcaption></figure>

and data analysis of similar problems. We have checked that the results presented in the paper are independent of the interpolation scheme.

We turn now to study the acceleration statistics. As displayed in Fig. 8(a), the PDF of the acceleration presents some deviations from a Gaussian distribution at large values. The norm of the acceleration has also an exponential tail for $|a| > \sigma_a|a|$, as displayed in Fig. 8(b). The core of the PDF in this case is a $\chi_3$ distribution, which is expected for the norm of a vector with Gaussian components. In classical Lagrangian turbulence, the norm of the particle acceleration is observed to obey a log-normal distribution [55]. In the inset of Fig. 8(b), we compare our data with such distribution. For the lightest and smallest particle, the small accelerations appear to be more probable than in the classical case. Note that, as pointed out in [55], small values of the acceleration are very sensible to experimental (numerical) errors. By contrast, the large accelerations are less probable than a log-normal distribution. This observation is compatible with classical numerical calculations in the framework of the viscous vortex filament model, in which it has been shown that, because of inertia, solid particles undergo less rapid changes of velocity than fluid particles [56].

Finally, in Fig. 9, we show the two-point correlator of the particle acceleration, defined as

$$ \rho^a(t) = \frac{\langle a_i(t_0)a_i(t_0+t) \rangle - \langle a_i(t_0) \rangle \langle a_i(t_0+t) \rangle}{\sigma_a(t_0)\sigma_a(t_0+t)}. \quad (11) $$

<figure><img src="image_11.png" /><figcaption>FIG. 9. Acceleration two-point correlator, plotted vs time normalized by the dissipation timescale $\tau_\ell$ (a), and by the Magnus natural frequency $1/\Omega_p$. (b) Markers indicate the time of acceleration decorrelation $t_a$. Inset: $t_a$ normalized by $1/\Omega_p$ as a function of the particle relative mass.</figcaption></figure>
In classical Lagrangian turbulence, the decorrelation time $t_a$ [such that $\rho^a(t_a) = 0$] is related to the Kolmogorov timescale $t_a = 2\tau_\eta$ [57]. This is not the case in quantum turbulence. Figure 9(a) displays the autocorrelation $\rho^a(t)$ for all the simulations. It is apparent that the acceleration decorrelates much faster than $\tau_\ell$, the equivalent of the Kolmogorov timescale in our system. This fact is a consequence of the myriad of physical phenomena taking place at smaller scales. As most particles are trapped by vortices, they oscillate at the Magnus frequency $\Omega_p$ in Eq. (9). If time is normalized by $\Omega_p$ (9), then $t_a\Omega_p$ becomes of order 1, at least for the heaviest particles [see Fig. 9(b) and the inset therein]. For the lightest particles, the decorrelation time is even lower, meaning that they are sensible to other mechanisms, such as reconnection events between vortex filaments and Kelvin wave excitations at even smaller scales.

## IV. DISCUSSION

In this work, we used the Gross-Pitaevskii model to study free decaying quantum turbulence at zero temperature in the presence of finite-size active particles. We considered different families of spherical particles having sizes smaller than and of the order of the mean intervortex distance. We first performed a standard analysis of the observables commonly used for studying Kolmogorov turbulence, such as the energy decomposition, the temporal evolution of mean energy, the rate of incompressible kinetic energy, and the mean intervortex distance. Although particles are active and get captured by vortices generating Kelvin waves, there is not a significant impact at scales larger than the intervortex distance, where Kolmogorov turbulence takes place. Monitoring the motion of the particles in the system, we confirmed their tendency to remain trapped into vortex filaments during the evolution of the tangle, with intermittent episodes of detachment and recapture. This behavior is independent of the vortex line density. We also found that particles can be easily captured simultaneously by several quantum vortices.

We also studied turbulence from the Lagrangian point of view. In particular, we computed the power spectra of the particle velocities. At large scales the particle dynamics is compatible with that of Lagrangian tracers in classical turbulence, while at short timescales the Magnus precession around the filaments caused by the vortex circulation is dominating the motion. Such information can be extracted consistently both in the frequency spectrum of the velocity and in the decay time of the correlation of the acceleration. Furthermore, if particles are light enough, faster frequencies are also excited. This suggests (as intuitively expected) that light particles can be more sensitive to the small-scale fluctuations of the flow.

Finally, we investigated the particle velocity statistics. The distribution of the particle velocity is Gaussian, in contrast with the power-law scaling $|v_i|^{-3}$ recently observed in superfluid helium experiments [24,25]. There are several reasons why power-law tails are absent in our simulations. First, since the simulation of each particle has an important numerical cost, the number of particles is restricted only to a couple of hundred. Due to this issue, vortex reconnections might be unlikely sampled by the sparse distribution of particles. Note also that, as particles have a finite size, increasing their number keeping the size of the system constant will increase substantially the filling fraction. In this case, turbulence could even be prevented by the presence of particles. Although interesting, this limit is beyond the scope of this work. Secondly, the GP model is compressible, and particles moving at large velocities are slowed down by vortex nucleations. This certainly reduces large velocity fluctuations, perhaps limiting the development of power-law tails. It would be interesting to address such issues in generalized GP models, including a roton minimum and high-order nonlinearities. Moreover, our simulations are by definition at zero temperature, and particles do not follow the singular superfluid velocity field because of the lack of viscosity in the system. Indeed, in the GP model the pressure gradients that drive the particle dynamics are always regular because of the vanishing density at the vortex cores, unlike other models such as the vortex filament method. As a consequence, the divergence of the superfluid velocity along the vortex lines cannot be experienced by the particles. Conversely, at finite temperature the superfluid and the normal component can be locked thanks to mutual friction. In this case, since particles
would sample the normal fluid velocity because of Stokes drag, they might be able to sample the $1/r$ flow around a quantum vortex. Finally, we observed that fast velocity fluctuations are highly sensitive to interpolation and filtering methods that could even lead to power-law tails. These tails are completely spurious, and special care is needed while analyzing numerical or experimental data.

## ACKNOWLEDGMENTS

U.G. acknowledges J. I. Polanco for fruitful discussions. The authors were supported by Agence Nationale de la Recherche through the project GIANTE ANR-18-CE30-0020-01. Computations were carried out on the Mésocentre SIGAMM hosted at the Observatoire de la Côte d'Azur and the French HPC Cluster OCCIGEN through the GENCI allocation A0042A10385.

## APPENDIX: NUMERICAL ARTIFACTS ON THE PARTICLE VELOCITY STATISTICS: COMPARISON BETWEEN B-SPLINE AND SPECTRAL INTERPOLATION METHODS

As explained in the main text, we evaluate the force $\mathbf{f}_i^{GP} = -(V_p * \nabla \rho)[\mathbf{q}_i]$ (3) at the particle position $\mathbf{q}_i$ using a B-spline interpolation method [39] at each time step. Such a method is precise and computationally cheap, but it turns out to present some issues that we have to take care of. To check the reliability of the method, we rerun a simulation using Fourier interpolation for one species of particles in the time window corresponding to the turbulent regime. Fourier interpolation is exact in the sense that it uses the information of the full three-dimensional field, which is resolved with spectral accuracy (i.e., discretization errors are at most exponentially small with the number of discretization points). The numerical cost of this method is that of one Fourier transform (per particle). In Fig. 10 the velocity and acceleration spectra computed using B-spline and Fourier interpolation methods are compared. Clearly, the B-spline interpolation introduces nonphysical fast oscillations, but at the frequencies $\omega < \omega_{noise} = 50(2\pi / \tau_{\ell})$ the behavior of the spectra is unchanged. Nevertheless, some differences in the features of particle statistics are still visible at fast timescales once the noise is filtered out.

We use a Gaussian convolution to perform a filtering of the velocity time series for each particle in the frequency window $\omega_c < \omega < \omega_{noise}$, where $\omega_c$ is a variable infrared cutoff frequency. Then we compute the PDF of the filtered velocity for different values of $\omega_c$. Such PDFs are shown in Fig. 11 comparing the simulations in which Fourier and B-spline interpolation are used for the same species of particle. Surprisingly, only in the latter case do we observe power-law tails for the fast oscillation distributions. Such PDFs are similar to the ones observed experimentally [24,25], but in the present case they are just a consequence of numerical artifacts.

<figure><img src="image_12.png" /><figcaption>FIG. 10. Velocity spectra (a) and the acceleration spectra (b) for particles of size $a_p = 10\xi$ and mass $M = 0.13$, evolved using B-spline interpolation (blue lines) and spectral Fourier interpolation (green lines). The spectra are averaged over particles and over the times $1.3T_L < t < 1.5T_L$.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 11. Probability density function of the velocity filtered in the frequency window $\omega_c < \omega < \omega_{\text{noise}}$ for different values of $\omega_c$. Data refer to particles of size $a_p = 10\xi$ and mass $M = 0.13$. The dotted line is a Gaussian distribution, and the dash-dotted line is a power-law scaling $0.002(v_i - \langle v_i \rangle)^{-3}$. The data are averaged over particles and over the times $1.3T_L < t < 1.5T_L$. Different PDFs are shifted for visualization. (a) Particle force interpolated with the B-spline method. (b) Particle force interpolated with the Fourier method.</figcaption></figure>

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

[2] R. J. Donnelly, *Quantized Vortices in Helium II* (Cambridge University Press, Cambridge, 1991), Vol. 2.

[3] T. Frisch, Y. Pomeau, and S. Rica, Transition to Dissipation in a Model of Superflow, Phys. Rev. Lett. 69, 1644 (1992).

[4] T. Winiecki and C. S. Adams, Motion of an object through a quantum fluid, Europhys. Lett. 52, 257 (2000).

[5] J. Koplik and H. Levine, Vortex Reconnection in Superfluid Helium, Phys. Rev. Lett. 71, 1375 (1993).

[6] G. P. Bewley, M. S. Paoletti, K. R. Sreenivasan, and D. P. Lathrop, Characterization of reconnecting vortices in superfluid helium, Proc. Natl. Acad. Sci. USA 105, 13707 (2008).

[7] A. Villois, D. Proment, and G. Krstulovic, Universal and nonuniversal aspects of vortex reconnections in superfluids, Phys. Rev. Fluids 2, 044701 (2017).

[8] L. Galantucci, A. W. Baggaley, N. G. Parker, and C. F. Barenghi, Crossover from interaction to driven regimes in quantum vortex reconnections, Proc. Natl. Acad. Sci. USA 116, 12204 (2019).

[9] C. Nore, M. Abid, and M. E. Brachet, Decaying Kolmogorov turbulence in a model of superflow, Phys. Fluids 9, 2644 (1997).

[10] A. W. Baggaley, J. Laurie, and C. F. Barenghi, Vortex-Density Fluctuations, Energy Spectra, and Vortical Regions in Superfluid Turbulence, Phys. Rev. Lett. 109, 205304 (2012).

[11] V. Shukla, P. D. Mininni, G. Krstulovic, P. C. di Leoni, and M. E. Brachet, Quantitative estimation of effective viscosity in quantum turbulence, Phys. Rev. A 99, 043605 (2019).

[12] J. Maurer and P. Tabeling, Local investigation of superfluid turbulence, Europhys. Lett. 43, 29 (1998).

[13] J. Salort, C. Baudet, B. Castaing, B. Chabaud, F. Daviaud, T. Didelot, P. Diribarne, B. Dubrulle, Y. Gagne, F. Gauthier, A. Girard, B. Hébral, B. Rousset, P. Thibault, and P.-E. Roche, Turbulent velocity spectra in superfluid flows, Phys. Fluids 22, 125102 (2010).

[14] W. F. Vinen, Decay of superfluid turbulence at a very low temperature: The radiation of sound from a kelvin wave on a quantized vortex, Phys. Rev. B 64, 134520 (2001).

[15] V. S. L'vov and S. Nazarenko, Weak turbulence of kelvin waves in superfluid he, Low Temp. Phys. 36, 785 (2010).

[16] J. Laurie, V. S. L'vov, S. Nazarenko, and O. Rudenko, Interaction of kelvin waves and nonlocality of energy transfer in superfluids, Phys. Rev. B 81, 104526 (2010).

[17] G. Krstulovic, Kelvin-wave cascade and dissipation in low-temperature superfluid vortices, Phys. Rev. E 86, 055301(R) (2012).
[18] A. W. Baggaley and J. Laurie, Kelvin-wave cascade in the vortex filament model, Phys. Rev. B **89**, 014504 (2014).

[19] A. Villois, D. Proment, and G. Krstulovic, Evolution of a superfluid vortex filament tangle driven by the Gross-Pitaevskii equation, Phys. Rev. E **93**, 061103(R) (2016).

[20] F. Toschi and E. Bodenschatz, Lagrangian properties of particles in turbulence, Annu. Rev. Fluid Mech. **41**, 375 (2009).

[21] W. Guo, M. La Mantia, D. P. Lathrop, and S. W. Van Sciver, Visualization of two-fluid flows of superfluid helium-4, Proc. Natl. Acad. Sci. USA **111**, 4653 (2014).

[22] G. P. Bewley, D. P. Lathrop, and K. R. Sreenivasan, Superfluid helium: Visualization of quantized vortices, Nature (London) **441**, 588 (2006).

[23] E. Fonda, D. P. Meichle, N. T. Ouellette, S. Hormoz, and D. P. Lathrop, Direct observation of kelvin waves excited by quantized vortex reconnection, Proc. Natl. Acad. Sci. USA **111**, 4707 (2014).

[24] M. S. Paoletti, M. E. Fisher, K. R. Sreenivasan, and D. P. Lathrop, Velocity Statistics Distinguish Quantum Turbulence from Classical Turbulence, Phys. Rev. Lett. **101**, 154501 (2008).

[25] M. L. Mantia and L. Skrbek, Quantum, or classical turbulence?, Europhys. Lett. **105**, 46002 (2014).

[26] M. La Mantia and L. Skrbek, Quantum turbulence visualized by particle dynamics, Phys. Rev. B **90**, 014519 (2014).

[27] M. La Mantia, D. Duda, M. Rotter, and L. Skrbek, Lagrangian accelerations of particles in superfluid turbulence, J. Fluid Mech. **717**, R9 (2013).

[28] Y. A. Sergeev and C. F. Barenghi, Particles-vortex interactions and flow visualization in 4he, J. Low Temp. Phys. **157**, 429 (2009).

[29] C. F. Barenghi, D. Kivotides, and Y. A. Sergeev, Close approach of a spherical particle and a quantised vortex in helium II, J. Low Temp. Phys. **148**, 293 (2007).

[30] Y. Mineda, M. Tsubota, Y. A. Sergeev, C. F. Barenghi, and W. F. Vinen, Velocity distributions of tracer particles in thermal counterflow in superfluid $^4$He, Phys. Rev. B **87**, 174508 (2013).

[31] E. Varga, C. F. Barenghi, Y. A. Sergeev, and L. Skrbek, Backreaction of tracer particles on vortex tangle in helium ii counterflow, J. Low Temp. Phys. **183**, 215 (2016).

[32] J. I. Polanco and G. Krstulovic, Inhomogeneous distribution of particles in coflow and counterflow quantum turbulence, Phys. Rev. Fluids **5**, 032601 (2020).

[33] V. Shukla, M. Brachet, and R. Pandit, Sticking transition in a minimal model for the collisions of active particles in quantum fluids, Phys. Rev. A **94**, 041602(R) (2016).

[34] U. Giuriato and G. Krstulovic, Interaction between active particles and quantum vortices leading to kelvin wave generation, Sci. Rep. **9**, 4839 (2019).

[35] U. Giuriato, G. Krstulovic, and D. Proment, Clustering and phase transitions in a 2d superfluid with immiscible active impurities, J. Phys. A **52**, 305501 (2019).

[36] U. Giuriato, G. Krstulovic, and S. Nazarenko, How trapped particles interact with and sample superfluid vortex excitations, Phys. Rev. Research **2**, 023149 (2020).

[37] A. Griffin, S. Nazarenko, V. Shukla, and M.-E. Brachet, The vortex-particle magnus effect, arXiv:1909.11010.

[38] P. Clark di Leoni, P. D. Mininni, and M. E. Brachet, Dual cascade and dissipation mechanisms in helical quantum turbulence, Phys. Rev. A **95**, 053636 (2017).

[39] M. A. T. van Hinsberg, J. H. M. Thije Boonkkamp, F. Toschi, and H. J. H. Clercx, On the efficiency and accuracy of interpolation methods for spectral codes, SIAM J. Sci. Comput. **34**, B479 (2012).

[40] N. Sasa, T. Kano, M. Machida, V. S. L'vov, O. Rudenko, and M. Tsubota, Energy spectra of quantum turbulence: Large-scale simulation and modeling, Phys. Rev. B **84**, 054525 (2011).

[41] G. Krstulovic, Grid superfluid turbulence and intermittency at very low temperature, Phys. Rev. E **93**, 063104 (2016).

[42] W. F. Vinen, Mutual friction in a heat current in liquid helium ii iii. theory of the mutual friction, Proc. R. Soc. London, Ser. A **242**, 493 (1957).

[43] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevFluids.5.054608 for movies of the simulations.
[44] N. G. Berloff and P. H. Roberts, Capture of an impurity by a vortex line in a Bose condensate, Phys. Rev. B **63**, 024510 (2000).

[45] A. Villois, G. Krstulovic, D. Proment, and H. Salman, A vortex filament tracking method for the gross-pitaevskii model of a superfluid, J. Phys. A **49**, 415502 (2016).

[46] The circulation measured with this method is subjected to a numerical error coming from the grid spacing. Such error is removed in postprocessing, knowing that $\Gamma$ can only be an integer multiple of $\kappa$. Furthermore, extremely high values of $\Gamma$ have been removed since they are related to the ill-defined situation in which a topological defect is placed at the boundary $C$.

[47] P. K. Yeung, Lagrangian characteristics of turbulence and scalar transport in direct numerical simulations, J. Fluid Mech. **427**, 241 (2001).

[48] H. Tennekes, J. L. Lumley, J. Lumley et al., *A First Course in Turbulence* (MIT Press, Cambridge, 1972).

[49] D. Kivotides, Y. A. Sergeev, and C. F. Barenghi, Dynamics of solid particles in a tangle of superfluid vortices at low temperatures, Phys. Fluids **20**, 055105 (2008).

[50] L. Kiknadze and Y. Mamaladze, The magnus (kutta-jukovskii) force acting on a sphere, arXiv:cond-mat/0604436.

[51] A. C. White, C. F. Barenghi, N. P. Proukakis, A. J. Youd, and D. H. Wacks, Nonclassical Velocity Statistics in a Turbulent Atomic Bose-Einstein Condensate, Phys. Rev. Lett. **104**, 075301 (2010).

[52] M. S. Paoletti and D. P. Lathrop, Quantum turbulence, Annu. Rev. Condens. Matter Phys. **2**, 213 (2011).

[53] A. W. Baggaley and C. F. Barenghi, Quantum turbulent velocity statistics and quasiclassical limit, Phys. Rev. E **84**, 067301 (2011).

[54] V. Shukla, M. Brachet, and R. Pandit, Turbulence in the two-dimensional fourier-truncated gross-pitaevskii equation, New J. Phys. **15**, 113025 (2013).

[55] N. Mordant, A. M. Crawford, and E. Bodenschatz, Three-Dimensional Structure of the Lagrangian Acceleration in Turbulent Flows, Phys. Rev. Lett. **93**, 214501 (2004).

[56] D. Kivotides, C. F. Barenghi, A. J. Mee, and Y. A. Sergeev, Interaction of Solid Particles with a Tangle of Vortex Filaments in a Viscous Fluid, Phys. Rev. Lett. **99**, 074501 (2007).

[57] P. K. Yeung and S. B. Pope, Lagrangian statistics from direct numerical simulations of isotropic turbulence, J. Fluid Mech. **207**, 531 (1989).