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

The evolution of a turbulent tangle of quantum vortices in 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 are presented and discussed. The particle Lagrangian dynamics is also studied. At large time scales, the velocity spectrum of particles is reminiscent of a classical Lagrangian turbulent behavior. At time-scales faster than the turnover time associated to the mean inter-vortex distance, the particle motion is dominated by oscillations due to Magnus effect. For light particles a non-classical 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.


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 towards small scales in a cascade process. Eventually, it reaches the smallest scales of the system, where dissipation acts efficiently. In 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 exotics 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, that 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 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][6][7][8]), unlike classical (prefect) fluids.
When energy is injected in a low-temperature superfluid at scales much larger than the mean inter-vortex distance , a classical Kolmogorov regime is expected. Such a behavior has been observed numerically [9][10][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 become relevant. The system keeps transferring energy towards small scales but through different non-classical 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 non-linear wave interactions. This energy cascade has been successfully described in the framework of weak-wave turbulence theory [15,16]. The resulting theoretical prediction have been observed numerically in vortex-filament and Gross-Pitaevskii numerical simulations [17][18][19].
Flow visualization is certainly a fundamental issue in every fluid dynamics experiment. Among the techniques which 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 been also 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][25][26]. Similar deviations from classical behaviors have been recently 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. Although it has been seen that particles unveils the dynamics of quantum vortices, it is yet not clear how much they affect the dynamics of quantum turbulent flows.
Several theoretical efforts have been made in the last 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 to study different specific problems. The interaction between one particle and one vortex has been addressed [28,29], as well the back-reaction 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 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][34][35][36][37]. Although the GP model is formally derived for dilute Bose-Einstein condensates, it is considered as 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. In this framework, the process of trapping of large active inertial particles by straight vortex lines has been thoroughly analyzed [34], and the interplay between many trapped particles and Kelvin waves has been also investigated [36].
In this 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 particles of different masses and having sizes smaller than and comparable to the mean inter-vortex distance. We also study the different regimes of the turbulent evolution from the Lagrangian point of view. The paper is organized as follows. In Section 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 Section III, we present our main results. In particular, in Section III A we address whether or not the presence of particles affects the scales of the flow at which Kolmogorov turbulence takes place. Section III B is devoted to the study of the particle dynamics inside the vortex tangle, to their trapping by vortices and to their dynamics at scales larger and smaller than the inter-vortex distance. Particle velocity and acceleration statistics are then presented in Section III C. Finally, Section IV contains our conclusions.

A. Gross-Pitaevskii equation coupled with particles
We describe a superfluid of volume V at low-temperature by using the complex field ψ, 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 q i and its classical momentum p i . Particles are modeled as strong localized potentials V p , centered at the position q j . The potentials deplete the superfluid in a sphere of radius a p , which defines the size of the particle. All the particles considered have the same size, as well as the same mass M p . The Hamiltonian of the system is given by 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πa s 2 /m. The chemical potential is denoted by µ. The particle interaction potential V ij rep is responsible for short-range repulsion between particles, so that they behave as hard spheres and do not overlap. The equations of motion that govern the superfluid field and the particle positions are obtained varying the Hamiltonian (1): 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)(3).
In the case where particles are absent, the chemical potential µ fixes the value of the ground state of the system ψ ∞ = ρ ∞ /m = µ/g. Large wavelength perturbations around this state are sound waves that propagate with the speed of sound c = gρ ∞ /m 2 , while they become dispersive at length scales smaller than the healing length ξ = 2 /2gρ ∞ . The GP model describes a superfluid with zero viscosity. Using the Madelung transformation ψ(x) = ρ(x)/m e i m φ(x) the GP equation (2) is mapped into the continuity and Bernoulli equations of a superfluid of density ρ and velocity v s = ∇φ. A superfluid flow is potential, but the phase is not defined at the nodal lines of ψ(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 ξ, and the contour integral of the superfluid velocity around a single vortex filament is the Feynman-Onsager quantum of circulation κ = h/m = 2π √ 2cξ. 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]: where ( √ ρv s ) I = P I [ √ ρv s ] and ( √ ρv s ) C = v s − ( √ ρv s ) I , the operator P I [·] being the projector onto the space of divergence-free fields. The other energies of the superfluid are the internal energy E int = (2V ) −1 (gρ − µ) 2 dx, where the energy of the ground-state is subtracted, and the interaction energy with the particles E GP P = V −1 Np i V p (|x − q i |)ρ dx, so that the total energy is given by E tot = E GP kin +E int +E GP P . 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-3) in a cubic periodic box of side L = 341ξ with N c = 512 3 collocation points by using a standard pseudo-spectral method. We use a 4 th order Runge-Kutta scheme for the time-stepping and the standard 2/3 rule for the dealiasing. In numerics, we fix c = 1 and ψ ∞ = 1.
In order 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 · 2π/L and k = 2 · 2π/L basic ABC flows: and the parameters A = 0.5196, B = 0.5774, 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 follows the ABC vortex lines. The initial mean inter-vortex distance is (t = 0) ∼ 25ξ. 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 form the vortex configuration.
The ground state for the particles consists in a number of particles (we use N p = 200 and N p = 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 the equation (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 to the particles ground state. An example of an initial field containing particles is displayed in Fig.1.d.
The potential used to model each particle is a smoothed hat-function where the parameters ζ and ∆ are set to model the particle. Their values are listed in Table I. The mass displaced by a particle is measured as where ψ p is the steady state with just one particle. 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. Since the particle boundaries are not sharp, we define the effective particle radius as a p = (3M 0 /4πρ ∞ ) 1 3 for given values of the numerical parameters ζ and ∆ a . Particles attract each other by a short range fluid mediated interaction [35,39]. In order to avoid an overlap between them we use the repulsive potential V ij rep = γ(2a p /|q i − q j |) 12 . The pre-factor γ is adjusted numerically in order to set minimum of V ij rep plus the fluid mediated attractive potential at the inter-particle distance 2a p [35]. 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 can not 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 the evaluation of the force term (3) acting on particles requires to know the value of the fields at inter-mesh points.
When the number of particles in the simulation is not large, the force can be computed with spectral accuracy using a Fourier interpolation. Such method has been used in [34][35][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. Therefore, we use a fourth-order B-spline interpolation method, which has been shown to be highly accurate with a reduced computational cost [40]. Nevertheless, some issues arising from the interpolation are discussed in Appendix.

III. PARTICLES IMMERSED IN A TANGLE OF SUPERFLUID VORTICES
Superfluid turbulence in the context of the GP model has been largely studied [9,11,38,41,42]. 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 thermal bath providing an effective dissipation acting vortices. As a consequence, vortices shrink and eventually disappear through mutual friction effects following Vinen's decay law [19,43]. In particular, it has been shown that the decrease of the incompressible kinetic energy behaves in a similar manner to decaying classical turbulence [9]. In order to make connection with decaying classical Kolmogorov turbulence, the incompressible energy dissipation or dissipation rate is usually defined in the context of GP turbulence as 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 energy spectrum satisfies the Kolmogorov prediction where C is the Kolmogorov constant which value has been found to be close to 1 in GP turbulence [11,38,42]. 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-3). In Fig.1.a, b and c the three stages of the evolution (respectively initial condition, turbulent vortex tangle and residual filaments in a bath of sound) are visualized in the case of 200 neutrally buoyant particles of radius 4ξ. Trapped particles by vortices are displayed in green, whereas free ones 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 Section 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 Bernoulli pressure. 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 to be used as probes.
In the next section, 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 units of the large-eddy-turnover time defined as T L = L/2v rms , where v rms = 2E I kin (t = 0)/3 is the root-mean-square velocity associated to the initial vortex tangle and L/2 is its characteristic length scale. We compare the case where 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 towards 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 effect is more evident if the number of large particles is increased, and could be related to an increment of the filling fraction Φ, namely the fraction of the total volume occupied by the particles. In fact, for N p = 200 particles of radius a p = 4ξ the filling fraction is Φ = 0.1%, for N p = 80 particles of radius a p = 10ξ it is Φ = 0.8% and for N p = 200 particles of radius a p = 10ξ we have Φ = 2.1%. The kinetic and repulsion energies of the particles, as well as the particle-vortex interaction E GP P are negligible compared to 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 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 where many large particles are moving in the system. In particular for this case, it is about 90% of 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 quasi-steady state and we perform the temporal average of certain physical quantities in order to improve statistical convergence.
Two other important quantity that is not affected much by the interplay between tangle and particles is the mean inter-vortex distance , whose time evolution is reported in Fig.3.a. The mean inter-vortex distance is then estimated  as = 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 2D vort (k): The spectra of momentum densities are the angle average of the norm in Fourier space of the momentum density J = ρ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 min ∼ 14ξ of the inter-vortex 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ξ), 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 E I (k) averaged around the maximum of dissipation is unaltered by particles in the system. The figure Fig.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 inter-vortex 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 with the Kolmogorov prediction E I (k) = C 2/3 max 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 particle perturbs the vortex filaments and excite 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 (a) (b) the tangle. At the times when turbulence is developed, the details of the vortex configurations are completely different (see lower row of Fig.4). Nevertheless, the statistical properties of the system in this regime remain unchanged. We stress that the inter-vortex distance in quantum turbulence experiments lies typically in the range 10 − 100 µm, which is equal 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 the believe 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 section 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), 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 where 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 wether 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 Γ = C v s · dx 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 + ∆ x ) around each particle, where ∆ 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 inter-vortex distance the majority of particles is in this situation. In the first stages of the evolution of the flow, the number of trapped particles increases rapidly, until it gets stationary always at times much smaller than one T L . The time needed to reach a stationary state slightly depends 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 trap p /N p is between 80% and 90% for small particles (2a p < ), while in average the totality of particles of size 2a p ∼ 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 trap p /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 homogenous 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 that the one observed here [31].
The circulation around each superfluid vortex filament is equal to a single quantum of circulation κ. 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 ∼ are always attached to at least two different filaments. Sometimes even more vortices pass simultaneously through the same particle, as it 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 ∆t trap spent by the particles inside 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 ∼ (∆t trap ) −α , with α ∼ 1.67. The PDF certainly vanishes much slower than an exponential decay at large ∆t trap , which would typically result form 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 peculiar of the most turbulent regime, since they persist also at the late times of the simulations (see dotted blue line in Fig.5.d). Therefore, many particles spend a time at least of the order of the simulation time (∼ 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) are short. Since 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 where B is a constant of order unity andv p (ω) is the Fourier transform of the Lagrangian particle velocity v p (t) [47,48]. Such scaling is valid in the inertial range 2π/T L ω 2π/τ η , where τ η is the Kolmogorov time scale. In our case, we build an analogous of the Kolmogorov time scale under the assumptions that the dissipation rate max is the only important physical parameter in the classical turbulence regime and that the Kolmogorov turbulent cascade ends at the inter-vortex distance min . Therefore, we define the smallest time scale of the classical turbulence regime as τ = 2 min / max 1/3 , and we expect classical turbulent phenomenology to hold for times τ t T L . In Fig.6 the measurement of the frequency spectrum of the particle velocity |v p (ω)| 2 = | q(t)e −iωt dt| 2 during the turbulent regime is shown for different species of particles, compensated with the classical scaling max ω −2 . Note that the average which defines the spectrum is meant over different realizations. In numerics we average over all the particles trajectories during the turbulent regime. At frequencies ω < τ /2π the spectra approach a plateau of value one, 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 ∼ 5τ . For comparison, we also present the velocity spectrum of a particle of size a p = 4ξ and mass M = 1, computed in a temporal window at much later times, when Kolmogorov turbulence has decayed and only few vortices are left. As expected, no Kolmogorov scaling is observed at small time scales.
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 inter-vortex distance. Unlike classical turbulence (see for instance [47]), there is still a non trivial scaling at time scales shorter than τ . Such 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 F Magnus = 3 2 ρ ∞ a p Γ × (q − v s ), where the circulation vector Γ is oriented along the vortex filament and the superfluid velocity v s contains the contributions of the mean flow and the vortex motion [36,49]. Magnus effect induces a precession of the particle about the filament with the characteristic angular velocity large particles with ap = 10ξ. 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 Ωp. Dashed red line is the scaling due to vortex reconnection or Kelvin waves ∝ |ω| −1 . Dashed golden line is the spectrum evaluated at late times in the simulation (6TL < τ < 7TL).
where the particle effective mass M eff p = M p + 1 2 M 0 = (M + 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 − 100Hz. If the Magnus force is the main force acting on a trapped particle, the Newton equation M eff pq = F Magnus implies the following expression for the frequency spectrum of the particle velocity: Independently of the external superfluid velocity, the expression (10) predicts that the spectrum |v p (ω)| 2 must be peaked around the natural frequency of trapped particles ω = Ω 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 Ω 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π/Ω p = 0.7T L for the particles of size a p = 10ξ 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.9.a, could be also related to this fact. At small time scales, a different scaling of the velocity spectrum is observed for the light particles, now in agreement with |v p (ω)| 2 ∝ |ω| −1 . This behaviour is consistent with the fact that at scales smaller than the inter-vortex distance, the typical velocities of a superfluid turbulent tangle are supposed to scale as v fast (t) ∝ κ/|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 |v p (ω)| 2 ∼v 2 fast (ω) ∝ κ|ω| −1 . Another effect that could contributes 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 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][25][26]. Such scaling has been related to the singular velocity field of quantized vortices [50,51]. At low temperatures, as Stokes drag is negligible, particles should not move with the superfluid flow and such scaling can be understood as a consequence of quantum vortex reconnections sampled by trapped particles [7,24]. Furthermore, in reference [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 inter-vortex distance. In this last section 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 ∇φ. 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 [50,52,53]. The PDF of ∇φ 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 is filtered with a Gaussian convolution in order to smooth out the noisy oscillations at frequencies ω < ω noise = 50 (2π/τ ) (see Appendix A). 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 were 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 the Discussion section.
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 reference [25], we have computed the velocity PDFs of the velocity fluctuations filtered at a given frequency ω c . The frequency was varied from values lower to larger than 2π/τ . 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, that eventually lead to a v −3 scaling. However, when using Fourier interpolation, that is an exact evaluation (up to spectral convergence of the pseudo-spectral 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 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| > σ |a| , as displayed in Fig.8.b. The core of the PDF in this case is a χ 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 [54]. In the inset of Fig.8.b, we compare our data with such distribution. For the lightest and smallest particle, the small accelerations appears to be more probable than in the classical case. Note that, as pointed out in [54], 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. Finally, in Fig.9, we show the two-point correlator of the particle acceleration, defined as: In classical Lagrangian turbulence, the decorrelation time t a (such that ρ a (t a ) = 0) is related to the Kolmogorov time scale t a = 2τ η [55]. This is not the case in quantum turbulence. Figure 9.a displays the autocorrelation ρ a (t) for all the simulations. It is apparent that the acceleration decorrelates much faster than τ , the equivalent of the Kolmogorov time scale 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 Ω p in Eq. (9). If time is normalized by Ω p (9), then t a Ω 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, like reconnection events between vortex filaments and Kelvin waves 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 presence of finite size active particles. We considered different families of particles having sizes smaller than and of the order of the mean inter-vortex 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 inter-vortex distance. Although particles are active and get captured by vortices generating Kelvin waves, there is not a significant impact at scales larger than the inter-vortex 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 particle 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 the one of Lagrangian tracers in classical turbulence, while at short time scales 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 particles 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 explaining why power-law tails are absent in our simulations. Firstly, since the simulation of each particle has an important numerical cost, the number of particles is restricted only to a couple of hundreds. 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 and turbulence could be eventually driven by particles. Although interesting, this limit is out of the scope of this work. Secondly, the GP-model is compressible and particle 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 non-linearities. 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 vanishing density at the vortex cores, in contrast to other models like the vortex filament method. As a consequence, the divergence of the superfluid velocity along the vortex lines can not be experienced by the particles. Conversely, at finite temperatures the superfluid and the normal component can be locked thanks to mutual friction, and particles might be able to sample the 1/r superfluid 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.
Appendix A: 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 f GP i = − (V p * ∇ρ) [q i ] (3) at the particle position q i using a B-spline interpolation method [40] at each time step. Such method is precise and computationally cheap, but it turns out to present some issues that we have to take care of. In order to check the reliability of the method, we re-run 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 uses the information of the full three dimensional field, that 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 the one 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 non-physical fast oscillations, but at the frequencies ω < ω noise = 50(2π/τ ) 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 ω c < ω < ω noise , where ω c is a variable infrared cut-off frequency. Then we compute the PDF of the filtered velocity for different values of ω 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 we observe power-law tails for the fast oscillations 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.