Eulerian and Lagrangian properties of turbulent flows under multi-scale helical injection

We present numerical and theoretical results concerning the properties of turbulent flows under strong multi-scale helical injection. We perform direct numerical simulations of the Navier-Stokes equations under a random helical stirring with power-law spectrum and with different intensities of energy and helicity injections. We show that there exists three different regimes where the forward energy and helicity inertial transfers are: (i) both leading with respect to the external injections, (ii) energy transfer is leading and helicity transfer is sub-leading, and (iii) both are sub-leading and helicity is maximal at all scales. As a result, the cases (ii-iii) give flows with Kolmogorov-like inertial energy cascade and tunable helicity transfers/contents. We further explore regime (iii) in a Lagrangian domain, by studying the kinetics of point-like isotropic helicoids, particles whose dynamics is isotropic but breaks parity invariance. We investigate small-scale fractal clustering and preferential sampling of intense helical flow structures. Depending on its structural parameters, the isotropic helicoids either preferentially sample co-chiral or anti-chiral flow structures. We explain these findings in limiting cases in terms of what is known for spherical particles of different densities and inertia. Furthermore, we present theoretical and numerical results for a stochastic model where Lagrangian properties can be calculated using analytical perturbation theory. Our study shows that a suitable tuning of the stirring mechanism can strongly modify the small-scale turbulent helical properties and demonstrates that isotropic helicoids are the simplest particles able to preferentially sense helical properties in turbulence.


Introduction
Helicity is an invariant of the Navier-Stokes equations in three spatial dimensions when neglecting the effects of viscous dissipation and external forcing (Frisch 1995;Moffatt & Tsinober 1992;Chen et al. 2003a;Alexakos & Biferale 2018). It is connected to the topological structure of vortex lines, characterized in terms of twist, writhe and linking numbers (Scheeler et al. 2014;Kedia et al. 2016;Laing et al. 2015;Kerr 2015). Helicity can be introduced in a flow by a stirring mechanism that breaks mirror symmetry and its effects on the turbulent energy cascade in three spatial dimensions have been widely studied since the pioneering work of Brissaud et al. (1973) (see Borue & Orszag (1997); Pelz et al. (1985); Kerr (1987); Kholmyansky et al. (1991); Kit et al. (1987) for other major contributions). In geophysical flows, helicity plays an important role in the atmospheric Ekman layer, where there exist arguments supporting a turbulent helicity cascade in the logarithmic range of the boundary layer (Deusebio & Lindborg 2014;Koprov et al. 2005;Kurgansky 2017). Recent experimental advancements allowed to produce vortex bundles with different prescribed topology (Kleckner & Irvine 2013) and the combination of shear and helicity has been studied experimentally and numerically (Herbert et al. 2012;Qu et al. 2018). Concerning the dual energy-helicity cascade, it is widely believed that for the case of Navier-Stokes equations in three spatial dimensions forced on a limited scale-range, both energy and helicity cascade forward (Chen et al. 2003a,b;Sahoo et al. 2015). This is a dual co-directional cascade according to the classification given in Alexakos & Biferale (2018). As a result, the mirror-symmetry breaking induced by the helical stirring mechanism tends to become weaker and weaker by going to smaller and smaller scales: the energy transfer is the leading mechanism and small-scale turbulence recovers a neutral statistics with zero helicity signature in average. On the contrary, if only one homochiral sector is dynamically active, one can prove that Navier-Stokes equations admit a dual counter-directional cascade (with energy flowing backward and helicity forward) the flow has global solutions (Waleffe 1992;Biferale et al. 2012;Biferale & Titi 2013) and small-scale turbulence is strongly (maximally) helical. Similarly, there are analytical and numerical hints (Linkmann 2018) that helicity induces a non-trivial decrease in the drag coefficient of turbulent flows. In this paper we further investigate the statistical properties of the dual energy-helicity transfers by adopting a power-law multi-scale stirring mechanism, which allows us to explore three different regimes concerning the relative intensity of energy and helicity injections. In particular, we show that there exists a suitable range of forcing spectral exponents, where the energy transfer is not affected by the stirring term while helicity can be controlled, leading to a turbulent realization with tunable small-scale helicity contents. Furthermore, in a regime where both small-scale energy and helicity contents are controlled by the forcing, leading to maximal-helicity flow configurations, we study the preferential concentration of isotropic helicoids (Kelvin 1871; Gustavsson & Biferale 2016), i.e. point-like particles whose dynamics is isotropic but breaks mirror-symmetry. By using both Direct Numerical Simulations (DNS) and a stochastic model for the Eulerian advecting velocity field (Gustavsson & Mehlig 2016), we show that isotropic helicoids show highly non-trivial preferential sampling of the underlying helical flow properties depending on the particle parameters. The paper is organized as follows. In Sec. 2 we describe the Eulerian part, discussing the different regimes for different helical injection power spectra and we present numerical simulations of the different regimes. In Sec. 3 we introduce the isotropic helicoids and their Lagrangian equations. We discuss the existence of two new scales of the Stokes number, St ± , which depend on the coupling between translational and rotational degrees of freedom. Furthermore, we present results on the preferential sampling of the flow helicity for different particle parameters, including two asymptotic limits where the Stokes number St is either much smaller than St + or much larger than St − . We conclude the paper in Sec. 4.

Theoretical background
We start by considering the forced Navier-Stokes equations for the fluid velocity u and the pressure p in three spatial dimensions: where ν the kinematic viscosity and f is a parity-breaking external forcing with energy injection rate = u · f and helicity injection rate h = u · (∇ × f ) + 2Ω · f , where 2Ω = ∇ × u denotes the flow vorticity. It is useful to adopt an exact decomposition of the velocity field in positive and negative Fourier helical waves (Constantin & Majda 1988;Waleffe 1992): where h ± k are the eigenvectors of the curl operator. In term of such decomposition the total energy, E = d 3 x u 2 , and the total helicity, H = 2 d 3 x u · Ω, take the forms: We can further consider the energy content of positive and negative helical modes, E ± (k) = k |k|<k+∆k |u ± k | 2 , where ∆k = 2π/L, such that the energy and helicity spectra become Supposing that there exists a dual co-directional forward cascade of energy and helicity and that the typical time at scale r ∼ k −1 is dominated by the energy eddy turn-over time τ E (r) ∼ −1/3 r 2/3 , we have for the semi-sum and semi-difference of the spectral components (Chen et al. 2003a): where C E and C H are two constants of dimension inverse length. Hence the two energy components can be written as:

Multiscale energy and helicity injections
Let us suppose a Gaussian white-in-time helical forcing, f (x, t) = k f + k h + k , whose twopoint correlation is isotropic, and with a power-law spectrum (Sain et al. 1998;Biferale et al. 2004;Kessar et al. 2015): For the sake of numerical implementation we cut off the power-law at a maximum wavenumber of the order of the Kolmogorov scale, k max ∼ k η . Using this forcing, the energy and helicity injection rates up to the scale k < k max can be estimated as: where k 0 = 2π/L is the smallest wavenumber in the box. From Eq. (2.8) we distinguish three different regimes depending on the forcing spectrum: (I) when y > 2 both energy and helicity injections are dominated by the infrared range, (k) → const. when k → ∞, and the system behaves as if it was forced at large scales only. We obtain a dual energy-helicity cascade because both quantities are transferred by the non-linear inertial terms of the Navier-Stokes equations (2.1); (II) when 1 < y < 2 the energy injection is still dominated by the infrared range, while the helicity injection depends on the wave-number, h(k) ∼ k 2−y . In this regime we obtain an energy cascade and helicity multi-scale injection; (III) when y < 1 both energy and helicity transfer are dominated by the multiscale injection, (k) ∼ k 1−y and h(k) ∼ k 2−y . We obtain an energy-helicity multiscale injection. The three regimes are summarized in Table 1. As a result, the spectral properties (2.6) are valid only for regime (I), and we can summarise the scaling for all three different regimes as follows: where the prefactors depend on the forcing as C X E ∝ |f 0 | 2 and C X H ∝ k 0 |f 0 | 2 . From the expressions (2.9) we can evaluate the mirror-symmetry recovery ratio, ) in the three regimes as: from which it follows that regime (III) is a flow with a maximal helical content at all scales where the injection is acting. Before concluding this section it is important to stress again that the prediction leading to regime III is correct under the assumption that the typical time scale guiding the transfer is the scale-dependent generalization of the eddy-turn-over time: τ E (k) ∝ k −2/3 (k) −1/3 , which is not necessarily the only possibility. In the following section we will present a first numerical investigation of the previous arguments.

Numerical simulation
In this section we show the results of a series of DNS with resolution of 512 3 grid points to explore properties of the energy and helicity of the three fluid regimes identified in the previous section. We implement a hyper-viscosity method to extend the inertial range (Borue & Orszag 1995). In particular we set ν α ∆ α u as the viscous term, with α = 2. The external forcing f in Navier-Stokes equations (2.1) has been implemented as a Langevin process with correlation time proportional to a fraction of the Kolmogorov time. As detailed in the previous section, to obtain a fully helical flow, we project the forcing only on velocity modes with positive helicity with energy injection on all wavenumber shells, up to dissipative scales k ∈ [0.5 : 60]. Three representative values for the three regimes have been selected: y = 4, 1.5, −1. Details about the 512 3 DNS set-ups are summarized in Table 2. In Fig. 1, we present four panels with the results for (a) the energy spectrum, (b) the helicity spectrum, (c) the energy flux and (d) the helicity flux as functions of k and for the three representatives values of y. The predictions for the scalings of the spectra (2.9) and for the energy fluxes (2.8) are verified with good accuracy, except for ultraviolet effects induced by the cut-off wavenumber where we stop to act with the external forcing to avoid stability issues in the code. Overall, we conclude that by changing the spectral properties of the helically-forced Navier-Stokes equations we can achieve a flow evolution with tunable energy/helicity ratios as theorized by (2.10). In particular, in the left panel  Figure 1. Direct numerical simulations for the energy and helicity spectra (a,b) and fluxes (c,d) for the three regimes y = 4 (I), y = 1.5 (II), and y = −1 (III). Parameters are given in Table 2. In panel a, the curve for y = 1.5 ( ) has been shifted with respect to the curve for y = 4 ( ) for the sake of presentation. We also superpose the scalings predicted by the relations in Eqs. (2.9) and (2.8).
of Fig. 2 we show both the positive and negative helical spectral components E ± (k) for y = −2/3 (case III). The major contribution to the energy spectrum is given by the velocity modes with positive helicity E + (k) for all wavenumbers. As a result the Navier-Stokes flow develops a dominant positive helical dynamics at all scales. The right panel of Fig. 2 shows that |E ( k) + −E − (k)| ∼ (E + (k) + E − (k)) which implies that mirrorsymmetry recovery is broken at all scales. We remark that in this regime, the scaling behaviour of |E + − E − | and E + + E − are less steep than the Kolmogorov prediction being both dominated by the external injection as predicted by (2.9).  Table 2. Right: Total energy E + (k)+E − (k) and rescaled total helicity E + (k)−E − (k) = H(k)/k. The scaling −5/9 predicted by relation (2.9) and the Kolmogorov −5/3 power laws are also shown for comparison.

Stochastic helical flows
The fully helical flow described by the regime (III) can be considered a sort of multi-scale flow dominated by the external forcing, where the Navier-Stokes non-linear evolution is sub-leading with respect to the forcing effects at all scales. As a result, it might be important to study also surrogate dynamics given by simpler stochastic evolution without any underlying structure coming from Navier-Stokes equations. This approximation is also necessary to perform analytical estimates for the dynamics of particles in the flow as discussed later. To follow this idea, we consider a random incompressible, homogeneous, and isotropic single-scale velocity field, u = ∇ × A. Here the components of the vector potential A(x, t) are independent Gaussian random functions with zero mean, a spatial correlation function decaying on a scale of order η 0 and an exponential time-correlation function with decay rate, τ 0 (see Appendix A for more details). The velocity field is normalized such that u 2 = u 2 0 . The flow is characterized by a dimensionless Kubo number Ku = u 0 τ 0 /η 0 , the ratio between the Eulerian flow decorrelation time τ 0 and the advecting time, η 0 /u 0 . The Kubo number can be seen as a dimensionless correlation time of the flow. If Ku tends to zero a white-noise flow is approached and if Ku is large a persistent flow is obtained. The latter case is important because the particle dynamics often agrees qualitatively or even quantitatively with the dynamics in a real turbulent flow (Gustavsson et al. 2015;Gustavsson & Mehlig 2016;Gustavsson et al. 2017). The former case is important because it allows for an analytical perturbative analysis in the Kubo number (Gustavsson & Mehlig 2011, that allows to understand the particle dynamics quantitatively at small values of Ku and qualitatively at large values of Ku or in DNS. In order to control the probability distribution function of the parity-breaking structures in the flow, we adopt the exact helical decomposition of each Fourier mode given by (2.2), weighting the positive modes h + k with a factor µ allows for construction of flows where positive (µ > 1, H flow > 0) or negative (µ < 1, H flow < 0) helical structures are dominant. The resulting flow has the following exponential-like distribution of helicity (see Appendix A for details):  where K ν (x) is the modified Bessel function of the second kind and H 0 is the average dimensionless helicity (2.12) Fig. 3 shows a comparison to 256 3 DNS (fourth case in Table 2) using H 0 = 0.85 (µ ≈ 1.5) to make the shape of the distribution (2.11) similar to that of the DNS described above.
In order to compare to DNS, it is necessary to take into account that the smooth length scale of the dissipation range in DNS is larger than the Kolmogorov length by a factor proportional to √ Re λ for not too large Re λ (Calzavarini et al. 2009). In our DNS we have Re λ ∼ 100 and we therefore use η 0 ∼ 10η K for the comparison. We observe that the distributions in Fig. 3 agree well for small values of H, but slightly disagrees in the right tail. This is not surprising, we cannot expect to reproduce the exact shape of the helicity distribution in Navier-Stokes equations with a single-scale stochastic flow.

Helical turbulent flows: Lagrangian Helicoids
The helical flow described in Sec. 2 breaks chirality: in configurations dominated by positive helicity, as the flow in Fig. 3, structures where the flow velocity and vorticity aligns are dominant. Heavy, inertial spherical particles are not able to distinguish the chirality of the underlying flow, they centrifuge out of vortex structures independent of their sign of helicity. We therefore study the dynamics of so called isotropic helicoids (Kelvin 1871; Happel & Brenner 2012). These are the simplest generalisation of spherical particles, their dynamics breaks parity, but remains isotropic. One example of isotropic helicoids suggested by Kelvin (Kelvin 1871) is illustrated in Fig. 4. Twelve planar wanes are attached perpendicular to the surface of a sphere at equal distances on three great circles. All vanes either form the angle +45 • (anti-chiral helicoid) or −45 • (co-chiral helicoid) with the great circle traversed clockwise, see Fig. 4. The wanes cause a coupling between translational and rotational motion. The dynamics of isotropic helicoids was studied in stationary ABC flows in Gustavsson & Biferale (2016). It was shown that the spatial distribution of isotropic helicoids depends on the relative chirality between the particle and the underlying flow. In this section we use DNS, the stochastic model and theoretical approaches to analyse the motion of isotropic helicoids in helical turbulence. The initial response to two simple flow configurations are illustrated with arrows. In response to an applied vertical difference in velocity, both helicoids accelerate in the direction of relative velocity, while their angular accelerations depend upon the sign of C0. Similarly, in response to an applied vertical difference in vorticity, the helicoids obtain the same angular acceleration, but they are accelerated in opposite directions.

Isotropic helicoids
The dynamics of an isotropic helicoid with position x, velocity v and angular velocity ω suspended in a fluid with velocity u and vorticity 2Ω = ∇ × u is governed by the following equations (Kelvin 1871; Happel & Brenner 2012;Gustavsson & Biferale 2016 (3.1) Here dots denote time derivatives and u and Ω are evaluated at the particle position x(t).
The dynamics of isotropic helicoids couples individual vector components of translational and rotational motion, but does not mix different components. The dynamics is governed by four parameters. First, τ p is a relaxation time quantifying particle inertia. In the limit of τ p → 0 the particle approaches the dynamics of a tracer, v = u and ω = Ω. Second, a = 5I 0 /2m is a measure of the particle size defined by its mass m and moment of inertia I 0 . Third, C 0 is the helicoidality. It quantifies the strength of the coupling between translational and rotational degrees of freedoms. Finally, S is the structural number that quantifies how much the rotational inertia of the isotropic helicoid differs from that of a spherical particle. When C 0 = 0, the particle dynamics is that of an isotropic particle, and if further S = 1, the dynamics is that of a spherical particle with Stokes relaxation time τ p . When C 0 = 0, invariance of the particle dynamics under mirror reflections of the particle is broken. Depending on the relative sign between C 0 and components of Ω, the particle accelerates either along the vorticity component, or opposite to it, see Fig. 4. The only constraint on the parameters is |C 0 | < √ 27S, required for the kinetic energy of the particle to remain finite. The actual size of the particle, ∼ã, should also be less than the smooth scale of the flow (a multiple of the Kolmogorov length η) for the point-particle approximation to be valid. The governing equations (3.1) exemplify why isotropic helicoids are simpler extensions to spherical particles than spheroids: the dynamics of spheroids depends on their instantaneous direction in addition to v and ω and it couples different components of the velocity and angular velocity. Moreover, in the limit of inertialess spheroids, the particle angular velocity does not simply follow Ω, but is also affected by the strain rate of the flow (Jeffery 1922).
and Ω = Ωτ η , and dropping the primes in what follows, we can write the equations of motion for each pair of components v i and ω i in dimensionless form: (3.2) Here we have introduced the dimensionless size a =ã/η and the Stokes number St = τ p /τ η . Interpreting the two-tensor D as a matrix, it has two eigenvalues d ± and corresponding eigenvectors ξ ± given by These equations are well defined for all parameter values, but the limit of isotropic particles, C 0 → 0, has a complication. Taking the limit C 0 → 0 in Eqs. (3.3) and (3.4) we need to distinguish the two cases of S < 3/10 and S > 3/10, resulting in the eigenvalues and eigenvectors given in Table 3. The reason is that when C 0 = 0, the eigenvalues cross at S = 3/10, meaning that the translational eigenvalue switches from being the largest (d + ) when S < 3/10 to the smallest (d − ) when S > 3/10. Moreover, the translational and rotational degrees of freedom decouples when C 0 = 0 and the eigenvector (1, 0) corresponding to the translational dynamics must be ξ + for S < 3/10 and ξ − for S > 3/10 with a discontinuous jump at S = 3/10. In the same way, the eigensystem corresponding to the rotational dynamics has a discontinuity at S = 3/10. When C 0 = 0 the translational dynamics has a single scale of inertia. When S ∼ 1 this scale is St ∼ 1, see  Table 4. Overview of the four parameter families in the simulations in Fig. 5. For each family, St varies over a few decades. The dynamics of the helicoids is driven by the helicoidality C0, shape factor S, and particle size a =ã/η which define the characteristic scales St−, St+ and β eff as introduced in Eqs. (3.3) and (3.14) respectively; for details see Sec. 3.1. The notion of being co-chiral or anti-chiral is made in terms of the flow helicity which is always taken in average positive in this paper.
for clustering can be explicitly related to preferential sampling. In Gustavsson & Biferale (2016) the divergence of the velocity field along the trajectory of an isotropic helicoid was derived for small values of St (St  1987). For helicoids the structures in which particles with small values of St converge are more intricate and depend on the particle parameters, a, C 0 , S, combined with the local flow helicity as expressed by the last term on the right-hand side of Eq. (3.5).
As observed by Gustavsson & Biferale (2016), for a flow region with strong helical coherence, V ∼ cA, particles cluster where (27S − 9acC 0 /5)Tr[A 2 ] is positive. As a consequence, particles of opposite helicoidality (different signs of C 0 ) may accumulate in flow regions of opposite sign of helicity c. As a result, even if the helicoids is heavier than the surrounding flow, they may cluster in vortical regions where Tr[A 2 ] < 0, similar to light spherical particles. In order to quantify the preferential sampling of helical flow structures, we simulate the dynamics (3.2) for a number of parameters summarized in Table 4 using the flows described in Sec. 2.3. For each set of parameters, once the fluid reaches its statistically stationary state it is seeded with 2.4 × 10 6 particles. The initial velocity and angular velocity of each particle are given by the fluid velocity and half the fluid vorticity evaluated at the particle position. Fig. 5a shows the fraction of particles in rotational regions of the flow. In rotational regions of a flow, the fluid gradient matrix A has complex eigenvalues, or equivalently, the sign of the discriminant is positive (Chong et al. 1990). We observe that isotropic helicoids with the same positive chirality of the underlying flow, C 0 = 1.6 (solid blue markers), depend intricately on the Stokes number: for small values of St they behave similar to light particles that oversample rotational regions where ∆ > 0, while for larger values of St they instead behave as heavy inertial particles that oversample strain regions where ∆ < 0. In contrast, for the other considered values of C 0 , the helicoids always behave as heavy particles and  Table 2. b Mean fluid helicity H = 2 u · Ω along particle trajectories as a function of St for the DNS. The data is normalized by the helicity of the flow, H flow which is always chosen to be positive in our simulations. The parameters of the simulations are given in Table 4 and the simulation results are displayed as interconnected markers. Results for neutral particles, C0 = 0, are shown as black asterisks. Results for S = 1, a = 10 helicoids are shown as hollow orange circles (anti-chiral, C0 = −5) and filled red circles (co-chiral, C0 = 5). Results for S = 0.1, a = 30 helicoids are shown as hollow lightblue boxes (anti-chiral, C0 = −1.6) and filled blue boxes (co-chiral, C0 = 1.6). Black dashed lines show P (∆ > 0) and H for tracer particles. c, d Same as a, b, but for the stochastic model flow with Ku = 10, H0 = 0.85, and ratio between the smooth length scale η0 and the Kolmogorov scale η, η0/η = 10.
oversample strain regions to different degrees depending on the particle parameters. Fig. 5b shows the mean value of fluid helicity evaluated along particle trajectories, H(x(t)) = 2 u · Ω , as a function of St. Comparing Fig. 5a and b shows that the behaviour is quite similar: helicoids with C 0 = 1.6 oversample rotational regions and have larger helicity than the underlying flow if the Stokes number is small enough. This is consistent with these particles spending long time in rotational regions of the flow where helicity is high and mainly of a given sign due to the helical nature of the underlying flow. Particles with the other investigated parameter values on the other hand, experience a fluid helicity that is lower than that of tracer particles. This is consistent with these particles aggregating in fluid strain regions where helicity is small. We also observe a transition at intermediate Stokes numbers: for small values of St, isotropic helicoids with negative values of C 0 are more likely to sample flow regions of negative helicity, while for large values of St helicoids with positive values of C 0 on average sample more negative values of helicity. In Fig. 5c and d we show the same plots of panels (a) and (b) but for the stochastic model described in Sec. 2.4. We choose parameters that are expected to correspond well with the DNS parameters. The Kubo number Ku is taken to be large, corresponding to the persistent-flow limit where the dynamics most resembles the small scales in turbulence. As described in Sec. 2.4 we take H 0 = 0.85 to match the distribution of flow helicity to that of the DNS and we take η 0 /η = 10 to compensate for the difference between the smooth length scale of the dissipation range and the Kolmogorov length in DNS. Finally, similarly to the DNS we base the Stokes number in the stochastic model on the Lagrangian time scale of tracer particles, τ η ≡ Tr(AA T ) −1/2 flow = η 0 /( √ 5u 0 ). In previous studies similar schemes have resulted in qualitative agreement between stochastic model simulations and DNS for the dynamics of spherical particles, elongated particles, and gyrotactic microswimmers (Gustavsson et al. 2015;Gustavsson & Mehlig 2016;Gustavsson et al. 2017). As shown in Fig. 5c and d we obtain a qualitative agreement with the DNS also for the isotropic helicoids. The general trends as functions of St for the different values of C 0 agree, but the detailed values disagree in some ranges. One example is the probability to find helicoids with C 0 = 1.6 in rotational regions with ∆ > 0. Although being larger than the other parameter values, it is not larger than that of the underlying flow as for the DNS case. One possible explanation for this is that the life time of vortex regions is longer in DNS than in the stochastic model.

Limiting cases
We now explain the data shown in Fig. 5

) and (3.4) to write the dynamics (3.2) in its diagonal basis asζ
(3.8) Here we have changed basis for each component i using where the columns of the 2 × 2 matrix X consist of the eigenvectors ξ − and ξ + . Eqs. (3.7) and (3.8) are only implicitly coupled through the trajectory dependence in u − and u + . We therefore expect that the two limits St St + and St St − discussed above can be taken, to a lowest-order approximation, in one of Eqs. (3.7) and (3.8) independent from the second equation. Consider first St St + with general values of St − . For the acceleration to remain finite in Eq. (3.8) in this limit, we must have 0 ∼ u +,i − ζ +,i . In terms of the original coordinates, this condition gives the following constraint: (3.10) Using Eq.(3.10) and its time derivative to replace ω andω in Eq. (3.7), and reverting to the original coordinates, we obtain: (3.11) Thus, a single equation determines the velocity of the isotropic helicoid in the limit St St + and the angular velocity is given by Eq. (3.10). It can be noted that for the case of C 0 = 0, the constraint (3.10) becomes singular because velocity and angular velocity are uncoupled. However, Eq. (3.11) still shows the same results as those obtained by letting C 0 = 0 and St St + in the original equation (3.2). When C 0 = 0 and S > 3/10, Eq. (3.11) simplifies tov = (u − v)/ St while the rotational dynamics is overdamped, ω − Ω = 0. When C 0 = 0 and S < 3/10 on the other hand, the condition St St + = 1 implies that the translational dynamics is overdamped and Eq. (3.11) simplifies tov = (u − v)/ St +u. This equation relaxes to the overdamped limit v = u after a short initial transient on the time scale of order St 1. The dynamics (3.11) can be further simplified using one or both of the following two assumptions. First, for small enough values of St / St − , we approximateu ∼ Du/Dt anḋ Ω ∼ DΩ/Dt, where D/Dt ≡ ∂ t + (u · ∇) = d/dt + (u − v) · ∇ are advective derivatives. Second, in a helical flow Ω and u tend to be aligned and we may approximate Ω ∼ cu with some proportionality constant c. Using these approximations, Eq.(3.11) simplifies to:v (3.12) If we define an effective Stokes number St eff and an effective density parameter β eff : (3.14) then Eq. (3.12) becomes identical to that for small spherical particles (Maxey & Riley 1983) with sub-dominant terms neglected: For spherical particles the density is characterised by β: 0 < β < 1 corresponds to heavy particles more dense than the fluid, β = 1 corresponds to neutrally buoyant particles and 1 < β 3 corresponds to particles lighter than the fluid. We remark that β eff in Eq. (3.14) is not constrained to the interval 0 β 3 as the case of spherical particles, β eff may also take negative values as well as values larger than 3. Fig. 6a,b compares stochastic model simulations of the approximation (3.11) to the simulation data in Fig. 5c and d. We observe a quantitative agreement of the approximation in the expected limit St St + . Fig. 6c,d show numerical simulations of (3.11) withu andΩ replaced by Du/Dt and DΩ/Dt, which approaches the results of Eq. (3.11) when St / St − 1 as expected. Finally, Fig. 6c,d also shows the approximation (3.12) using Ω = cu with c = ( Ω 2 / u 2 ) 1/2 = √ 5/2 for the stochastic model. This approximation only gives qualitative agreement, which is expected because the fluid velocity and vorticity are not perfectly aligned in the helical flow with H 0 = 0.85. However, the qualitative agreement allows us to explain the behaviour observed for small St St + in Fig. 5 in terms of what is known from spherical particles. The values of the effective density parameter β eff in Eq. (3.14) are quoted in Table 4 for our simulation parameters. Only the case C 0 = 1.6 has β eff larger than one, corresponding to spherical particles lighter than the fluid. Such particles are expected to preferentially sample rotational flow regions when the effective Stokes number St eff = St / St − is of order unity, which is consistent with the data in Fig. 5a and c. The cases of helicoids with C 0 = −1.6 or C 0 = −5 can be viewed as heavy particles because β eff is close to zero. In these cases the helicoids preferentially sample strain regions of the fluid when St eff ∼ O(1). Finally, the case C 0 = 5 has β eff ≈ 0.5, making it heavy but not as heavy as the cases with negative values of C 0 . As expected the preferential sampling of straining regions when St eff ∼ O(1) is therefore somewhat lower, see Fig. 6a. We now consider the second limiting case, St St − with general values of St + . In this limit ζ −,i in Eq. (3.7) respond slowly to changes in the flow compared to ζ +,i . Due to the symmetries of the underlying flow, we expect the averages of u −,i and consequently of ζ −,i to vanish. We therefore replace ζ −,i by its vanishing average, ζ −,i = 0, which gives the following constraint on ω Inserting this constraint and its time derivative into Eq. (3.8), gives the following equation for the velocityv As for the first limiting case we can consider a helical flow with Ω ∼ cu to obtaiṅ St − we observe a quantitative agreement. Fig. 7c,d shows that we only obtain qualitative agreement for the approximation (3.17) based on the approximation Ω = cu. This is similar to the first limiting case above and we use this limit to explain the observed data. For helicoids with negative values of C 0 , the coupling to the flow, β eff u in Eq. (3.17), is small, see Table 4. The particle motion is thus expected to be only weakly correlated to the underlying flow structures, which is consistent with the data: the particles with negative values of C 0 have approximately the same statistical properties as the flow (black dashed lines in Fig. 7). The helicoids with positive values of C 0 on the other hand have β eff ∼ 1 and are therefore expected to have preferential sampling similar to spherical particles with effective Stokes number St + / St. This is the case for simulations of Eq. (3.17), but the magnitude of the preferential sampling in Fig. 7c,d differ somewhat from the full numerical data and from Eq. (3.16) because Eq. (3.17) in only a qualitative approximation.

Small values of Ku
For the stochastic model we can solve the dynamics analytically in the limit of small Ku in terms of the full set of model parameters St, C 0 , S, a, and H 0 using the method in Gustavsson & Mehlig (2011. The calculation is outlined in Appendix B. The resulting mean helicity becomes to second order in Ku: (3.18). We found that multiplying the Stokes number by 4 and the Kubo number by 1/9 gives qualitative agreement. Fig. 8a shows that the small Kubo results have the same trends as those of the DNS and the stochastic model in Fig. 5b and d. This shows that the trends shown in Fig. 5 are robust, they do not depend on the particular nature of the underlying flow.
Using this observation, we can use the theoretical solution of the stochastic model to get an estimate of the parameter dependence of preferential sampling of helicity for general values of the five model parameters. Two examples of this dependence are illustrated in Fig. 8b and c. Fig. 8b shows how the observed Stokes-dependent preferential sampling of helicity depend on C 0 . For not too small values of |C 0 | the preferential sampling is similar to that observed in Fig. 5b and d: co-chiral helicoids oversample helicity for small Stokes numbers, while anti-chiral helicoids oversample helicity for large Stokes numbers. Helicoids with small |C 0 | behave similar to neutral particles and undersample helicity. Similar trends are observed in a neutral flow (see Fig. 8c). In a neutral flow helical structures of opposite signs are equally likely, which imposes a symmetry under the simultaneous change of H and C 0 to −H and −C 0 . This symmetry is clearly seen in Fig. 8c: upon changing C 0 to −C 0 the average helicity changes sign. As a consequence of the symmetry neutral particles (C 0 = 0) in neutral flows may not show preferential sampling of helicity, there is a thin line of no preferential sampling at C 0 = 0 in Fig. 8c.

Clustering
The previous section shows that helicoids of different chirality may go to very different regions in the flow. This is exemplified in Fig. 9 which shows a snapshot of the positions of isotropic helicoids of different chirality in a strongly helical stochastic flow. Fig. 9 clearly shows that the helicoids of opposite chirality distribute in different flow regions. It is also visible that the spatial clustering is different in nature, the helicoids with C 0 = −1.6 seems to distribute on a set which fills a larger portion of space compared to helicoids with C 0 = 1.6. It is known that spherical particles show small-scale fractal clustering due to preferential sampling and due to the dissipative nature of the dynamics. We start investigating the degree of spatial clustering of helicoids by computing the spatial  correlation dimension, D 2 . The probability distribution P (r) to find two helicoids within a spatial distance r scales as P (r) ∼ r D2 for small r. Here the exponent D 2 is the spatial correlation dimension and it quantifies the degree of spatial clustering. Fig. 10a shows the behavior of D 2 as functions of St / St − for the DNS with parameters given by the fourth case in  (Bec 2003;Toschi & Bodenschatz 2009). This explains why the helicoids with C 0 = 1.6 have a smaller fractal dimension, close to D 2 = 1.6, than the other, effectively heavy spherical particles. It also explains why the correlation dimension for the helicoids with C 0 = −1.6 and C 0 = ±5 approximately collapse on the correlation dimension of spherical particles when plotted against St / St − . The numerical data for C 0 = 1.6 shows a second peak of clustering (minimum of Figure 10. a Spatial correlation dimension D2 as a function of St / St− for the DNS with parameters given by the fourth case in Table 2. The parameters of the simulations are given in Table 4  approximation (3.17) where for β ef f ∼ O(1) we have a sort of heavy-particle behaviour. However, we remark that the result for D 2 is measured at finite separation which might not reflect the true asymptotic scaling for r → 0. Indeed, in the stochastic model the correlation dimension around St ∼ St + show a similar scaling P (r) ∼ r D2(r) with local exponent D 2 (r) < 3 for a range of r 1, while for small enough values of r, the uniform D 2 = 3 scaling is approached (not shown). This can be explained by the fact that the deviation from the approximation (3.17) is different for the two particles, which results in a uniform distribution of particles for small enough scales. Using D 2 to quantify fractal dimension of clustering may be problematic because it often converges slowly and consequently very small scales must be resolved in the distribution P (r). An alternative quantification of fractal clustering which is easier to evaluate accurately is provided by the Lyapunov dimension (Kaplan-Yorke dimension) D L (Frederickson et al. 1983). Denoting by λ 1 λ 2 · · · λ 9 the nine Lyapunov exponents of the equation system (3.1) together withẋ = v, the Lyapunov dimension is where K is the largest integer such that Fig. 10b shows numerical evaluation of the Lyapunov dimension for the stochastic model. For D L < 3 it shows similar trends as the correlation dimension D 2 in DNS. The values of D L are somewhat higher than D 2 , consistent with the fractal attractor in turbulence being a multifractal with D L D 2 (Bec 2005). Finally, we remark that the Lyapunov dimension defined in Eq. (3.19) describes the dimension of the fractal in phase-space, and it is therefore not bounded by the spatial dimension 3 as is the case for the spatial correlation dimension. For large values of St, the Lyapunov dimension is expected to approach the dimensionality D of phase-space. This is consistent with the data in Fig. 10b (D = 9 for helicoids and D = 6 for spherical particles). The long-term average of the compressibility along particle trajectories is identical to the sum over the first three Lyapunov exponents, ∇ · v = λ 1 + λ 2 + λ 3 . Using this relation, we verify in Fig. 11 the small St approximation given by Eq. (3.5). We find that the Figure 11. Stochastic model simulations of the relative amplitude of average compressibility along particle trajectories, ∇ · v = λ1 + λ2 + λ3, and the approximation (3.5), ∇ · v (1) . Parameters as in Fig. 10. approximation (3.5) tends to approach ∇ · v as St / St − is reduced and that the two expressions agrees approximately for St / St − ∼ 0.1. We end the discussion on particle clustering by remarking that we have applied the perturbation theory developed in Appendix B to calculate the Lyapunov exponents and Lyapunov dimension for small values of Ku, similar to the expansions for the Lyapunov exponents of spherical particles (Gustavsson & Mehlig 2011. The theory relies on the additional constraint St St − in order for caustic singularities to be rare. In this limit we observe good agreement between theory and numerical simulations of the stochastic model (not shown).

Discussion and conclusions
We have presented a series of numerical and theoretical results concerning the properties of turbulent flows under strong multi-scale helical injection. We performed direct numerical simulations of the Navier-Stokes equations up to resolution 512 3 and at changing the exponent of the power-law helical injection, in the limit of white-in-time noise. We first showed that there exists three different regimes for the forward energy and helicity non-linear transfers: (i) when both transfers are directed toward small-scales and the external multi-scale injection is negligible, leading to a −5/3 spectrum for both energy and helicity; (ii) when the energy cascade is fully non-linear and helicity is dominated by the forcing; and (iii) when both cascades are dominated by the forcing at all scales. For the case of turbulence under condition (iii) and for a surrogate stochastic flow we studied the evolution of isotropic helicoids, presenting a systematic assessment of preferential sampling and fractal clustering for helicoids with different properties. In particular, we showed that a suitable tuning of the chirality of the helicoids may lead to particles that behave either as being lighter or heavier than the surrounding flow. The phenomenon is enhanced for helicoids sizes larger then the Kolmogorov scale, where the approximation of point-particle might not be fully correct anymore. This apparent problem can be resolved by constructing particles with large effectiveã while the size of the particle interacting with the fluid is smaller (Gustavsson & Biferale 2016

Appendix A. Stochastic model for helical turbulence
To construct the incompressible, homogeneous, and isotropic stochastic velocity field, u = ∇ × A, used in the article, we generate the components of the vector potential A as a Fourier sum A i (r, t) = (2π) 3/4 3(1 + µ) η 5/2 Here L = 10 is the system size (we use L = 10η 0 in our simulations), the wave vector k is summed over the components k j = 2πn j /L with n j = −20, −19, . . . , +20 and j = 1, 2, 3. For each k, the vector of Fourier coefficients, a k (t), has been expanded in terms of the eigenmodes h ± k of the curl operator, weighted by a factor µ to give a bias to positive helical modes if µ > 1, and to negative helical modes if 0 µ < 1. The coefficients a i,k (t) are complex random Gaussian numbers fulfilling the condition a * i,k = a i,−k and having the statistics a i,k (t) = 0 and a i,k1 (t 1 )a * j,k2 (t 2 ) = δ ij δ k1k2 e −|t1−t2|/τ .

(A 2)
The exponential time correlation in (A 2) is generated from an underlying Ornstein-Uhlenbeck processes a i,k (t + δt) = e −δt/τ a i,k (t) + b i,k (t) , (A 3) where δt is the time step of the simulation and b i,k (t) are independent random Gaussian numbers that are white-noise in time with statistics b i,k (t) = 0 and b i,k1 (t)b * j,k2 (t) = δ ij δ k1k2 (1 − e −2δt/τ ) .

(A 5)
From this correlation function the statistics of u and its spatial derivatives follows. To obtain the distribution P 0 (H) of helicity H = 2u · Ω for the stochastic flow in Eq. (A 1), we start from the joint distribution of u and Ω P = 1 where X = (u 1 , u 2 , u 3 , Ω 1 , Ω 2 , Ω 3 ) T and C is the corresponding covariance matrix (velocity is made dimensionless in terms of u 0 and position in terms of η 0 ) obtained from Eq. (A 5) with H 0 ≡ 2/π8(µ 2 − 1)/(3(µ 2 + 1)). After a change of coordinates Ω z = (H/2 − Ω x u x − Ω y u y )/u z and integration over Ω x and Ω y the remaining joint distribution of H, u x , u y , and u z depends only on H and the combination