Helicoidal particles in turbulent flows with multi-scale helical injection

We present numerical and theoretical results concerning the properties of turbulent flows with 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 tuneable helicity transfers/contents. We further explore regime (iii) by studying its effect on 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 their 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 degrees of inertia. Furthermore, we present theoretical and numerical results for a stochastic model where dynamical 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.

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 Pelz et al. (1985), Kerr (1987), Kit et al. (1987), Kholmyansky et al. (1991), Borue & Orszag (1997) for other 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 (Koprov et al. 2005;Deusebio & Lindborg 2014;Kurgansky 2017). Recent experimental advancements allowed the production of vortex bundles with a different prescribed topology (Kleckner & Irvine 2013) and the combination of shear and helicity has been studied experimentally and numerically (Herbert et al. 2012;Qu, Naso & Bos 2018). Concerning the dual energy-helicity cascade, it is widely believed that for the case of NSE in three spatial dimensions forced on a limited range of scales, both energy and helicity cascade forward (Chen et al. 2003a,b;Sahoo, Bonaccorso & Biferale 2015). This is a dual co-directional cascade according to the classification given in Alexakis & Biferale (2018). The mirror-symmetry breaking induced by the helical stirring mechanism tends to become weaker and weaker by going to smaller and smaller spatial scales: the energy transfer is the leading mechanism and small-scale turbulence recovers a neutral statistics with zero helicity on average. On the contrary, if only one homochiral sector is dynamically active, one can prove that NSE admit a dual counter-directional cascade (with energy flowing backward and helicity forward). For this case the flow has global solutions (Waleffe 1992;Biferale, Musacchio & Toschi 2012;Biferale & Titi 2013) and small-scale turbulence is strongly (maximally) helical. In addition, 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 energyhelicity 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 tuneable small-scale helicity content. 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 1872; 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 possess highly non-trivial preferential sampling of the underlying helical flow properties depending on the particle parameters. The paper is organized as follows. In § 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 § 3 we introduce the isotropic helicoids and their dynamical 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 § 4. L. Biferale, K. Gustavsson and R. Scatamacchia 2. Helical turbulent flows: Eulerian properties 2.1. Theoretical background We start by considering the forced NSE for the fluid velocity u and the pressure p in three spatial dimensions: where ν is 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 terms 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 |u ± k | 2 , where k = 2π/L, such that the energy and helicity spectra become (2.4a,b) 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 turnover 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: It is known that for large-scale energy and helicity injection the Kolmogorov-like scaling (2.5) is observed, implying a recovery of mirror symmetry at small scales, see for example Sahoo et al. (2015), Vallefuoco et al. (2018) for recent studies about this issue with and without rotation. In order to have strong multi-scale helicity, it is necessary to resort to a power-law injection (Forster, Nelson & Stephen 1977;Seoud & Vassilicos 2007).
whose two-point correlation is isotropic, and with a power-law spectrum (Sain, Manu & Pandit 1998;Biferale, Lanotte & Toschi 2004;Kessar et al. 2015): where d is the space dimension and D 0 defines the typical forcing intensity at the smallest wavenumber that we will always assume to be k 0 = 2π/L = 1. 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: By considering spherical symmetry, the sums in (2.9) can be easily estimated and 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. and h(k) → const. when k → ∞, and the system behaves as if it is forced at large scales only. In this case, we obtain a dual energy-helicity cascade because both quantities are transferred by the nonlinear inertial terms of the NSE (2.1); (II) when 1 < y < 2 the energy injection sum is still dominated by the infrared range, while the helicity injection depends on the ultra-violet limit, 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 multi-scale injection, (k) ∼ k 1−y and h(k) ∼ k 2−y . 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 summarize the scaling for all three different regimes as follows: where the prefactors depend on the forcing intensity (2.8). From the expressions (2.10) we can evaluate the mirror-symmetry recovery ratio, L. Biferale, K. Gustavsson and R. Scatamacchia E − (k)) 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 obtained under the assumption that the typical time scale guiding the transfer is the scale-dependent generalization of the eddy turnover time: τ E (k) ∝ k −2/3 (k) −1/3 , which is not necessarily the only possibility. In order to have a quantitative assessment of the scaling properties at high Reynolds numbers one could resort to Fourier closures based on the eddy-damped quasi-normal Markovian (EDQNM) approximation as in Briard & Gomez (2017).
In the following, we resort to direct numerical simulations and we present a first numerical investigation of the flow properties under multi-scale helical injection without any approximation.

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 (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 at all wavenumbers up to dissipative scales k ∈ [1 : 70]. Three representative values for the three regimes have been selected: y = 4, 3/2, −1. Details about the 512 3 DNS set-ups are summarized in Table 2. In figure 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.
In the insets of panels (a,b) the total energy and helicity as functions of time in the stationary regime are shown. The predictions for spectra (2.10) and energy fluxes (2.9) are verified with good accuracy, except for ultraviolet effects induced by the cutoff wavenumber where we stop to act with the external forcing to avoid stability issues in the code. Note that the power-law forcing smooths down the presence of the high-wavenumber bottleneck expected in the spectrum when using hyper-viscosity (Frisch et al. 2008). Overall, we conclude that by changing the spectral properties of the helically forced NSE we can achieve a flow evolution with tuneable energy/helicity ratios as theorized by (2.11). In particular, in figure 2(a) 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 figure 2 shows that |E + (k) − E − (k)| ∼ (E + (k) + E − (k)) which implies that mirror symmetry 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 that is in both cases dominated by the external injection as predicted by (2.10). t/ †˙t/ †Ḟ IGURE 1. Time average of the energy and helicity spectra (a,b) and fluxes (c,d) for the three regimes y = 4 (I), y = 3/2 (II) and y = −1 (III). The small discontinuity at the high wavenumbers is due to the end of the range where the forcing is applied. Inset: time evolution of the total energy (a) and total helicity (b) in the stationary regime where all averages are performed. Parameters are given in Table 2. In (a), the curve for y = 3/2 ( E ) 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 (2.10) and (2.9).

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 nonlinear evolution is sub-leading with respect to the forcing effects at all scales. In order to have an analytical control and variability of the governing flow, we study also surrogate dynamics given by simpler stochastic evolution without any   Table 2. Inset: time evolution of the helical spectral components in the stationary regime. (b) 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.10) and the Kolmogorov −5/3 power laws are also shown for comparison. Inset: time evolution of total energy and rescaled total helicity. underlying structure coming from NSE. 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 (2.12) 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, and to understand the particle dynamics quantitatively at small Ku and qualitatively at large 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 µ leads to 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):  Table 2 (black crosses) and for the stochastic model (2.13) with H 0 = 0.85 (red line). The helicity is made dimensionless using the Kolmogorov scales η and τ η .
where K ν (x) is the modified Bessel function of the second kind and H 0 is the average dimensionless helicity (2.14) Figure 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.13) 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 figure 3 agree well for small values of H, but slightly disagree in the right tail. This is not surprising, we cannot expect to reproduce the exact shape of the helicity distribution in NSE with a single-scale stochastic flow.

Helical turbulent flows: suspensions of helicoidal particles
The helical flows described in § 2 break parity invariance (chiral symmetry): in configurations dominated by positive helicity, as the flow in figure 3, structures where the flow velocity and vorticity align 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 1872; Happel & Brenner 2012). These are the simplest idealized generalization of spherical particles, their dynamics breaks parity, but remains isotropic. One example of isotropic helicoids suggested by Lord Kelvin (Kelvin 1872) is illustrated in figure 4. Twelve planar vanes 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 figure 4. The vanes 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 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 C 0 . Similarly, in response to an applied vertical difference in vorticity, the helicoids obtain the same angular acceleration, but they are accelerated in opposite directions. (b) Snapshot in the stationary state for two types of helicoids of opposing helicoidality in DNS of the helical turbulent flow given by the fourth case in Table 2. Points show particle positions in a slice of height 5η. Parameters: St ≈ St − , S = 0.1, a = 30, and C 0 = −1.6 (anti-chiral, light blue) or C 0 = 1.6 (co-chiral, blue). Inset shows a zoom to highlight that particles of different chirality accumulate in different regions.
the spatial distribution of isotropic helicoids depends on the relative chirality between the particle and the underlying flow. We anticipate that this is also what happen in helical turbulence as illustrated in figure 4 where we show that isotropic helicoids move to different flow regions depending on their helicoidality also in our DNS of the NSE (2.1). In this section we use DNS, the stochastic model and theoretical approaches to analyse the motion of isotropic helicoids in helical turbulence.

Isotropic helicoids
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,ã = √ 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 figure 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) and (3.2) 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).
Rescaling to dimensionless units t = t/τ η , x = x/η, u = uτ η /η, v = vτ η /η, ω = ωτ η , 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: 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 . (3.5) These equations are well defined for all parameter values, but in the limit of isotropic particles, C 0 → 0, there is a complication. Taking the limit C 0 → 0 in (3.4) and (3.5) 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 decouple 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 3. This has been observed in simulations of inertial 3.2. Preferential sampling of vorticity and helicity Inertial spherical particles are subjected to preferential sampling of particular flow structures as well as small-scale fractal clustering (Maxey 1987;Fessler et al. 1994;Bec 2003;Gustavsson & Mehlig 2016). In the limit of small Stokes numbers the mechanism 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 St − ): where A and V are matrices with elements A ij = ∂ j u i and V ij = ∂ j Ω i . Depending on the sign of ∇ · v trajectories of close-by particles may either converge ] > 0 (Maxey 1987). For helicoids the structures in which particles with small values of St converge are more intricate and depend in addition on the particle parameters, a, C 0 , S, combined with the local flow helicity as expressed in the last term, ∝ Tr[AV ], on the right-hand side of (3.6). 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 are 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.3) for a number of parameters summarized in Table 4 Table 4 and the simulation results are displayed as interconnected markers. Results for neutral particles, C 0 = 0, are shown as black asterisks. Results for S = 1, a = 10 helicoids are shown as hollow orange circles (anti-chiral, C 0 = −5) and filled red circles (co-chiral, C 0 = 5). Results for S = 0.1, a = 30 helicoids are shown as hollow light blue boxes (anti-chiral, C 0 = −1.6) and filled blue boxes (co-chiral, C 0 = 1.6). Black dashed lines show P(∆ > 0) and H for tracer particles.
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. Figure 4 illustrates that isotropic helicoids of opposing chirality preferentially sample different flow regions in DNS of a helical turbulent flow. is positive (Chong, Perry & Cantwell 1990). We observe that isotropic helicoids with the same helicoidality (positive) of the underlying flow, C 0 = 1.6 (filled blue boxes), 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 oversample strain regions to different degrees depending on the particle parameters. Figure 5(c,d) shows the mean value of fluid helicity evaluated along particle trajectories, H(x(t)) = 2 u · Ω , as functions of 3.3. Small-scale fractal clustering The previous section shows that helicoids of different chirality may go to very different regions in the flow. This is exemplified in figure 4 which shows that the helicoids of opposite chirality distribute in different flow regions. It is also visible in figure 4 that the spatial clustering is different in nature, close-by helicoids with C 0 = −1.6 seem to distribute on different kinds of structures 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 investigate the degree of spatial clustering of helicoids by computing the spatial correlation dimension, D 2 , which defines the probability distribution P(r) to find two helicoids within a small spatial distance r: P(r) ∼ r D 2 , r 1. (3.8)

Limiting cases
We now explain the main features of the data shown in figures 5 and 6 by analysing three limiting cases of the system parameters. The first two limiting cases, St St + and St St − , apply to both DNS and the stochastic model, while the third limit is obtained for small values of Ku and therefore only applies to the stochastic model. We use (3.4) and (3.5) to write the dynamics (3.3) in its diagonal basis aṡ (3.10) Here we have changed basis for each component i using where the columns of the 2 × 2 matrix X consist of the eigenvectors ξ − and ξ + . The dimensionless parameter groups in (3.9) and (3.10) can be viewed as ratios of the time scales, τ ± /τ η ≡ St/St ± , where τ η is the Kolmogorov time of the flow and τ ± are two particle time scales (τ + τ − ) of the isotropic helicoid (τ + = τ − = τ p for a spherical particle). Equations (3.9) and (3.10) are only implicitly coupled through the trajectory dependence in u − and u + . We therefore expect that the two limiting cases St St + and St St − can be taken, to a lowest-order approximation, in one of (3.9) and (3.10) independent from the second equation. In summary, when St St + (τ + τ η ), equation (3.10) becomes overdamped and the remaining equation (3.9) gives rise to strong preferential sampling when St/St − ∼ 1 (τ − ∼ τ η ) in analogue to the case St ∼ 1 (τ p ∼ τ η ) for inertial spherical particles. In the second limit St St − (τ − τ η ), equation (3.9) becomes underdamped and ζ − can to lowest order be approximated by its mean value. The remaining equation (3.10) gives rise to strong preferential sampling when St/St + ∼ 1 (τ + ∼ τ η ). Below we discuss the two limiting cases in more detail.

Case St St +
Consider first St St + with general values of St − . For the acceleration to remain finite in (3.10) in this limit, we must have 0 ∼ u +,i − ζ +,i . In terms of the original coordinates, this condition gives the following constraint: (3.12) Using (3.12) and its time derivative to replace ω andω in (3.9), and reverting to the original coordinates, we obtaiṅ Thus, a single equation determines the velocity of the isotropic helicoid in the limit St St + and the angular velocity is given by (3.12). It can be noted that for the case of C 0 = 0, the constraint (3.12) becomes singular because velocity and angular velocity are uncoupled. However, equation ( The dynamics (3.13) can be further simplified using one or both of the following two assumptions. First, for small enough values of St/St − , we approximateu ∼ D t u andΩ ∼ D t Ω, where D t ≡ ∂ t + (u · ∇) 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, equation (3.13) simplifies tȯ (3.14) If we define an effective Stokes number St eff and an effective density parameter β eff : then (3.14) becomes identical to the equation of motion of small spherical particles (Maxey & Riley 1983) with sub-dominant terms neglected: (3.17) For spherical particles the density is characterized 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 (3.16) 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.

Comparison to DNS
Although being a crude first-order approximation, equation (3.14) allows us to use what is known from spherical particles to explain the main features of the behaviour of helicoids in DNS for St St + in figure 5(a). The values of the effective density parameter β eff in (3.16) are quoted in Table 4 for our 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 figure 5(a). 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 with strongest effect around 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 . This is consistent with a maximal preferential sampling of straining regions around St eff ∼ O(1) that is somewhat lower than for the case C 0 = −5, but inconsistent for the case of C 0 = −1.6 where the maximal preferential sampling is of the same order, see figure 5(a). The approximation (3.14) allows us to also explain the observed clustering in figure 6(a). The helicoids with C 0 = 1.6 show stronger clustering around St/St − ∼ 1 compared to the other cases. This can be explained by the observation that the dynamics of helicoids with C 0 = 1.6 is similar to that of light spherical particles and that the other types of helicoids have dynamics similar to that of heavy spherical particles with effective Stokes numbers St/St − . Light spherical particles cluster in rotational regions of the flow and show more clustering than heavy spherical particles (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 − .

Comparison to stochastic model
Below we study in detail the validity of the approximations leading to (3.14). Since DNS is slow, it is hard to reach the steady state for helicoids with St St − and to get good statistics in this limit. We therefore use the stochastic model to study the approximations. Figure 7(a,b) shows simulation results for isotropic helicoids with St St + in the stochastic model described in § 2.4 (solid lines). We choose parameters that are expected to correspond well with the DNS parameters in figure 5. We use a large Kubo number, Ku = 10, corresponding to the persistent flow limit where the dynamics most resembles the small scales in turbulence. As described in § 2.4 we  disagree in some ranges. One example is the probability of finding helicoids with C 0 = 1.6 in rotational regions with ∆ > 0. Although being larger than the probability for other parameter values, it is not larger than the probability 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. Figure 7(a,b) also compares stochastic model simulations of the approximation (3.13) to the full simulation data of (3.3). We observe a quantitative agreement of the approximation in the expected limit St St + . Figure 7(c,d) shows numerical simulations of (3.13) withu andΩ replaced by D t u and D t Ω, which approach the results of (3.13) when St/St − 1 as expected. Finally, figure 7(c,d) also shows the approximation (3.14) using Ω = cu with c = ( Ω 2 / u 2 ) 1/2 = √ 5/20 for the stochastic model (rescaled using η 0 = 10η). We observe that the prediction using spherical particles reproduces the average sampling of helicity well for St St − (figure 7d), while it does not work as well for the probability of being in vortex regions (figure 7c). Discrepancies in this approximation are expected because the fluid velocity and vorticity are not perfectly aligned in the helical flow with H 0 = 0.85. In particular, the approximation fails for the probability of being in vortex regions for the cases with C 0 ± 1.6: for C 0 = 1.6 the helicoids do not oversample vortex regions to the degree that is predicted by (3.14) and for C 0 = −1.6 the helicoids show a larger probability than is predicted. As discussed above, this is in contrast to DNS which agrees better with the trends predicted by (3.14). Figure 8 shows the positions of isotropic helicoids in a strongly helical flow in the stochastic model. Similar to the DNS in figure 4, helicoids of opposing chirality go to different regions in the flow and close-by anti-chiral particles form structures of a different kind than co-chiral helicoids. Figure 9(a) shows stochastic model results for the correlation dimension D 2 with parameters corresponding to the DNS in figure 6. We observe qualitative agreement between figure 9(a) and the DNS.
the nine Lyapunov exponents of the equation system (3.1) and (3.2) together witḣ x = v, the Lyapunov dimension is where K is the largest integer such that K i=1 λ i > 0. Figure 9(b) 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 . 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 (3.18) 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 figure 9(b) (D = 9 for helicoids and D = 6 for spherical particles).
Evaluation of the Lyapunov exponents also allows us to validate the approximation (3.6) of the local compressibility in the stochastic model. 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 figure 9(c) the small St approximation given by (3.6). We remark that the limit used to obtain (3.6) is the same as in this section (St St + ), but requires first-order corrections in St/St + to the condition in (3.12). Moreover, we need to consider St/St − 1. Consistently with this limit, we find that the average of the approximation (3.6) tends to approach ∇ · v as St/St − is reduced and that the two expressions agree approximately for St/St − ∼ 0.1.

Case St St −
We now consider the second limiting case, St St − with general values of St + . In this limit ζ −,i in (3.9) responds 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 (3.10), gives the following equation for the velocity: (3.20) As for the first limiting case we can consider a helical flow with Ω ∼ cu to obtaiṅ

Comparison to DNS
The limiting dynamics in (3.21) allows us to explain the DNS results in figure 5(b). For helicoids with negative values of C 0 , the coupling to the flow, β eff u in (3.21), 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 figure 5b). 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 what we observe for C 0 = 1.6, the shape and magnitude of the curve around St ∼ St + is similar to that of spherical particles. However, the approximation (3.21) does not work as well for C 0 = 5: even though β eff ≈ 0.5 the preferential sampling is only slightly larger than the case C 0 = −5. Since C 0 = 5 also only show small agreement with the predictions in the first limiting case St St + , we conclude that the approximations (for example Ω ∼ cu) leading to (3.14) and (3.21) may not be so accurate for C 0 = 5 in DNS. Using the approximation (3.21) to explain the fractal clustering observed in figure 6(b), we would expect that the anti-chiral helicoids, having small coupling to the flow, should show small clustering (D 2 ≈ d) and that the co-chiral helicoids should show larger clustering of the same order as the spherical particles. The numerical data for C 0 = 1.6 indeed show a second peak of clustering (minimum of D 2 ) around St ∼ St + . 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, as seen in the inset of figure 6(b), the local slope of D 2 for C 0 = 1.6 seems to drift towards larger values as r is decreased in the range of r we can resolve. This can be explained by the fact that deviations from the approximation (3.21) depend on the flow histories experienced by the particles and are therefore different for two close-by particles, which results in a uniform distribution of particles for small enough scales. In contrast, deviations from the overdamped approximation (3.14) mainly depend on the instantaneous flow and are therefore approximately the same for the two particles. We compare stochastic model simulations of the approximation (3.20) (solid lines) to simulations of (3.1) and (3.2) in figure 10. When St St − we observe a quantitative agreement for C 0 = ±1.6 and qualitative agreement when C 0 = ±5. Figure 10 also shows that the approximation (3.21) based on Ω = cu (dash-dotted lines) does not work equally well. Equation (3.21) reproduces that co-chiral helicoids have larger preferential sampling of straining regions around St ∼ St + than anti-chiral helicoids, but the degree of preferential sampling does not come out correctly in general. We conclude by remarking that, as expected, in the underdamped limit of St St + the preferential sampling for all parameter cases in Table 4 converges to that of the flow (not shown).
As seen in figure 9(a) the correlation dimension D 2 increases monotonously towards the spatial dimension d = 3 after the peak at St ∼ St − . As a consequence, D 2 does not show any clustering around St ∼ St + . The deviations from the DNS data in figure 6(b) can be explained by the fact that the observed data are better resolved in the stochastic model: the correlation dimension around St ∼ St + shows a scaling P(r) ∼ r D 2 (r) with local exponent D 2 (r) < 3 for a range of r 1 (similar to the DNS in this range), while for small enough values of r, the uniform D 2 = 3 scaling is approached (not shown).

Small values of Ku
In the stochastic model the properties of the flow can be modified by changing the value of the Kubo number (2.12). In general, this allows exploration of the robustness of results with respect to the nature of the flow. In the limit of small Ku we can solve the dynamics analytically 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, where a =ã/η 0 (velocity and position are made dimensionless in terms of u 0 and η 0 ). We found that multiplying the Stokes number by 4 and the Kubo number by 1/9 gives qualitative agreement (the effect of the former scaling is a horizontal shift of all curves and the effect of the latter scaling is a constant prefactor of the deviation from H 0 ). Figure 11(a) shows that the small Kubo results have the same trends as those of the DNS in figure 5(b) and the stochastic model in figure 7(b). This shows that the trends shown in figure 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 figure 11(b,c). Figure 11(b) shows how the observed Stokes-dependent preferential sampling of helicity depends on C 0 . For not too small values of |C 0 | the preferential sampling is similar to that observed in figure 5(b,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 figure 11c). 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 figure 11(c): upon changing C 0 to −C 0 the average helicity changes sign. As a consequence of this 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 figure 11(c). We end the discussion on small Kubo numbers by remarking that we have applied the perturbation theory developed in appendix B to calculate the Lyapunov exponents and Lyapunov dimension of particle clustering 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 NSE 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 nonlinear 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 nonlinear and helicity is dominated by the forcing; and (iii) when both cascades are dominated by the forcing at all scales. Finally, let us note that the theoretical prediction (2.10) is qualitatively well reproduced by our DNS results as shown in figure 1. Nevertheless, we must stress that for the dominant regime (case III in Table 2) the power law is not extremely clean. For the latter case, it would in future studies be important to extend the numerical resolution, such as to reduce spurious sub-leading terms.
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 small-scale 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 comparison between the turbulent and stochastic model shows very similar degree of preferential sampling for all parameters considered. Due to the different nature of the flows, this implies that the studied preferential sampling is mainly a kinematic effect: it depends on the dynamical equations (3.1) and (3.2) rather than on the existence and evolution of flow structures. This suggests that the observed effects are robust to changes in the details or nature of the flow, with interesting applications also at low or moderate Reynolds numbers. Other observables such as large-scale clustering or higher-order moments of helicity are likely to have a stronger dependency on the flow properties.
At a first glance the numerical data observed in figures 5 and 6 have sensitive and complicated parameter dependence. Nevertheless, it is remarkable that the crude approximations (3.14) and (3.21) allow us to give a first qualitative hint about the helicoid properties in terms of the relatively simpler dynamics of spherical particles, using the effective parameters St eff (3.15), β eff (3.16) and St/St + . The only term that is odd in c or C 0 in these effective parameters, as well as in (3.6), is proportional to acC 0 . This implies that helicoids with opposite helicoidality behave more differently the larger acC 0 is. For the flows and helicoids considered in this paper we have acC 0 ∼ 5. This implies an estimated helicoid sizeã larger than the smooth scale (approximately 10η), where the point-particle approximation may not be fully correct anymore. A fully systematic analysis of finite-size effects for the helicoid properties is still lacking. This apparent problem can be resolved by constructing helicoids with large effectiveã while the size of the particle interacting with the fluid remains small, for example by attaching small but heavy satellite particles to the helicoid (Gustavsson & Biferale 2016). Another solution is to consider flows or particles with larger values of c or |C 0 | (such that acC 0 is significant for small particle sizes). On one side, the magnitude of acC 0 determines how much helicoids with opposite helicoidality are different. On the other side, preferential sampling is also strongly dependent on the value of β eff , as exemplified by comparing |C 0 | = 5 and |C 0 | = 1.6 in figure 5. By optimizing β eff it is possible to find helicoids with smaller values of a that behave as light particles, i.e. have β eff > 1. We conclude by remarking that the construction in figure 4 is just one possible way to construct isotropic helicoids and what is the most general isotropic structure which breaks mirror symmetry with a given set of parameters C 0 , S and a remains an open question.
The exponential time correlation in (A 2) is generated from an underlying Ornstein-Uhlenbeck processes: a i,k (t + δt) = e −δt/τ 0 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,k 1 (t)b * j,k 2 (t) = δ ij δ k 1 k 2 (1 − e −2δt/τ 0 ). (A 4a,b) The Gaussian cutoff for large k in (A 1) ensures a Gaussian spatial correlation function with correlation length η 0 . When L η 0 , (A 1) implies the correlation function A i (r 1 , t 1 )A j (r 2 , t 1 ) = η 2 0 u 2 0 6 e −|r 1 −r 2 | 2 /(2η 2 0 )−|t 1 −t 2 |/τ 0 . (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 (A 1), we start from the joint distribution of u and Ω: 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 ) C ij = X i X j = 1 12 obtained from (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 u 2 x + u 2 y + u 2 z . Changing to spherical coordinates in u-space and integrating them away gives the final distribution of helicity, equation (2.13): where K ν (x) is the modified Bessel function of the second kind. The average helicity of the flow is determined from the helicity bias µ as follows: