Helicoidal particles in turbulent flows with multi-scale helical injection

Abstract

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.

1 Introduction

Helicity is an invariant of the Navier–Stokes equations (NSE) in three spatial dimensions when neglecting the effects of viscous dissipation and external forcing (Moffatt & Tsinober 1992; Frisch 1995; Chen, Chen & Eyink 2003a ; Alexakis & 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; Kerr 2015; Laing, Ricca & Witt 2015; Kedia et al. 2016). 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 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 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 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_{\pm }$ , 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.

2 Helical turbulent flows: Eulerian properties

2.1 Theoretical background

We start by considering the forced NSE for the fluid velocity $\boldsymbol{u}$ and the pressure $p$ in three spatial dimensions:

(2.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}\unicode[STIX]{x0394}\boldsymbol{u}+\boldsymbol{f},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D708}$ is the kinematic viscosity and $\boldsymbol{f}$ is a parity-breaking external forcing with energy injection rate $\unicode[STIX]{x1D716}=\langle \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{f}\rangle$ and helicity injection rate $h=\langle \boldsymbol{u}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\times \boldsymbol{f})+2\unicode[STIX]{x1D734}\boldsymbol{\cdot }\boldsymbol{f}\rangle$ , where $2\unicode[STIX]{x1D734}=\unicode[STIX]{x1D735}\times \boldsymbol{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):

(2.2) $$\begin{eqnarray}\displaystyle \boldsymbol{u}(\boldsymbol{x},t)=\mathop{\sum }_{k}[u_{\boldsymbol{k}}^{+}(t)\boldsymbol{h}_{\boldsymbol{k}}^{+}+u_{\boldsymbol{k}}^{-}(t)\boldsymbol{h}_{\boldsymbol{k}}^{-}]\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{h}_{k}^{\pm }$ are the eigenvectors of the curl operator. In terms of such decomposition the total energy, $E=\int \text{d}^{3}x\,\boldsymbol{u}^{2}$ , and the total helicity, $H=2\int \text{d}^{3}x\,\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D734}$ , take the forms:

(2.3a,b ) $$\begin{eqnarray}\displaystyle E=\mathop{\sum }_{\boldsymbol{k}}|u_{\boldsymbol{k}}^{+}|^{2}+|u_{\boldsymbol{k}}^{-}|^{2},\quad H=\mathop{\sum }_{\boldsymbol{k}}k(|u_{\boldsymbol{k}}^{+}|^{2}-|u_{\boldsymbol{k}}^{-}|^{2}). & & \displaystyle\end{eqnarray}$$

We can further consider the energy content of positive and negative helical modes, $E^{\pm }(k)=\sum _{|\boldsymbol{k}|=k}|u_{\boldsymbol{k}}^{\pm }|^{2}$ , where $\unicode[STIX]{x0394}k=2\unicode[STIX]{x03C0}/L$ , such that the energy and helicity spectra become

(2.4a,b ) $$\begin{eqnarray}\displaystyle E(k)=E^{+}(k)+E^{-}(k),\quad H(k)=k[E^{+}(k)-E^{-}(k)]. & & \displaystyle\end{eqnarray}$$

Supposing that there exists a dual co-directional forward cascade of energy and helicity and that the typical time at scale $r\sim k^{-1}$ is dominated by the energy eddy turnover time $\unicode[STIX]{x1D70F}_{E}(r)\sim \unicode[STIX]{x1D716}^{-1/3}r^{2/3}$ , we have for the semi-sum and semi-difference of the spectral components (Chen et al. 2003a ):

(2.5a,b ) $$\begin{eqnarray}\displaystyle E^{+}(k)+E^{-}(k)\sim C_{E}\unicode[STIX]{x1D716}^{2/3}k^{-5/3},\quad E^{+}(k)-E^{-}(k)\sim C_{H}h\unicode[STIX]{x1D716}^{-1/3}k^{-8/3}, & & \displaystyle\end{eqnarray}$$

where $C_{E}$ and $C_{H}$ are two constants of dimension inverse length. Hence the two energy components can be written as:

(2.6) $$\begin{eqnarray}\displaystyle E^{\pm }(k)\sim C_{E}\unicode[STIX]{x1D716}^{2/3}k^{-5/3}\pm C_{H}h\unicode[STIX]{x1D716}^{-1/3}k^{-8/3}. & & \displaystyle\end{eqnarray}$$

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).

2.2 Multi-scale energy and helicity injections

Let us suppose a Gaussian white-in-time helical forcing,

(2.7) $$\begin{eqnarray}\displaystyle \boldsymbol{f}(\boldsymbol{x},t)=\mathop{\sum }_{\boldsymbol{k}}f_{\boldsymbol{k}}^{+}(t)\boldsymbol{h}_{\boldsymbol{k}}^{+}\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}, & & \displaystyle\end{eqnarray}$$

whose two-point correlation is isotropic, and with a power-law spectrum (Sain, Manu & Pandit 1998; Biferale, Lanotte & Toschi 2004; Kessar et al. 2015):

(2.8) $$\begin{eqnarray}\displaystyle \langle f_{\boldsymbol{k}}^{+}(t)f_{\boldsymbol{k}^{\prime }}^{+}(t^{\prime })\rangle =D_{0}k^{1-d-y}\unicode[STIX]{x1D6FF}(t-t^{\prime })\unicode[STIX]{x1D6FF}_{\boldsymbol{k},\boldsymbol{k}^{\prime }}, & & \displaystyle\end{eqnarray}$$

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\unicode[STIX]{x03C0}/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}\sim k_{\unicode[STIX]{x1D702}}$ . Using this forcing, the energy and helicity injection rates up to the scale $k<k_{max}$ can be estimated as:

(2.9a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}(k)\sim \mathop{\sum }_{|\boldsymbol{q}|<k}|\boldsymbol{q}|^{1-d-y},\quad h(k)\sim \mathop{\sum }_{|\boldsymbol{q}|<k}|\boldsymbol{q}|^{2-d-y}. & & \displaystyle\end{eqnarray}$$

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, $\unicode[STIX]{x1D716}(k)\rightarrow \text{const.}$ and $h(k)\rightarrow \text{const.}$ when $k\rightarrow \infty$ , 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)\sim 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, $\unicode[STIX]{x1D716}(k)\sim k^{1-y}$ and $h(k)\sim k^{2-y}$ . The three regimes are summarized in Table 1.

Table 1. Energy and helicity injection regimes depending on the forcing spectrum in (2.8).

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:

(2.10) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}E^{\pm }(k)\sim C_{E}^{I}\unicode[STIX]{x1D716}^{2/3}k^{-5/3}\pm C_{H}^{I}h\unicode[STIX]{x1D716}^{-1/3}k^{-8/3},\quad y>2,\\ E^{\pm }(k)\sim C_{E}^{II}\unicode[STIX]{x1D716}^{2/3}k^{-5/3}\pm C_{H}^{II}k^{2-y}\unicode[STIX]{x1D716}^{-1/3}k^{-8/3},\quad 1<y<2,\\ E^{\pm }(k)\sim C_{E}^{III}k^{2(1-y)/3}k^{-5/3}\pm C_{H}^{III}k^{2-y}k^{-(1-y)/3}k^{-8/3},\quad y<1,\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where the prefactors depend on the forcing intensity (2.8). From the expressions (2.10) we can evaluate the mirror-symmetry recovery ratio, $R(k)=|E^{+}(k)-E^{-}(k)|/(E^{+}(k)+E^{-}(k))$ in the three regimes as:

(2.11) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}R^{I}(k)\sim k^{-1},\quad y>2,\\ R^{II}(k)\sim k^{1-y},\quad 1<y<2,\\ R^{III}(k)\sim \text{const.},\quad y<1,\end{array}\right\} & & \displaystyle\end{eqnarray}$$

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: $\unicode[STIX]{x1D70F}_{E}(k)\propto k^{-2/3}\unicode[STIX]{x1D716}(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.

2.3 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 $\unicode[STIX]{x1D708}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D6E5}^{\unicode[STIX]{x1D6FC}}u$ as the viscous term, with $\unicode[STIX]{x1D6FC}=2$ . The external forcing $\boldsymbol{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\in [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^{\pm }(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)|\sim (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).

Figure 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$ (○) 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).

Table 2. Parameters of the numerical simulations: grid resolution $N^{3}$ , Kolmogorov length scale $\unicode[STIX]{x1D702}$ in simulation units (SU), grid spacing $\unicode[STIX]{x0394}x=2\unicode[STIX]{x03C0}/N$ (SU), time step $\unicode[STIX]{x0394}t$ (SU), kinematic viscosity $\unicode[STIX]{x1D708}$ (SU), Kolmogorov eddy turnover time $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716})^{1/2}$ with $\unicode[STIX]{x1D716}$ the energy dissipation rate (SU), forcing correlation time $\unicode[STIX]{x1D70F}_{S}$ (in units of $\unicode[STIX]{x0394}t$ ), forcing power-law exponent $y$ , hyper-viscosity parameter $\unicode[STIX]{x1D6FC}$ , number of helicoids per each family $N_{h}$ .

Figure 2. (a) Time average of positive and negative helical spectral components for the direct numerical simulations with parameters of the fourth parameter set in 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.

2.4 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 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, $\boldsymbol{u}=\unicode[STIX]{x1D735}\times \boldsymbol{A}$ . Here the components of the vector potential $\boldsymbol{A}(\boldsymbol{x},t)$ are independent Gaussian random functions with zero mean, a spatial correlation function decaying on a scale of order $\unicode[STIX]{x1D702}_{0}$ and an exponential time-correlation function with decay rate, $\unicode[STIX]{x1D70F}_{0}$ (see appendix A for more details). The velocity field is normalized such that $\langle \boldsymbol{u}^{2}\rangle =u_{0}^{2}$ . The flow is characterized by a dimensionless Kubo number

(2.12) $$\begin{eqnarray}\displaystyle Ku=u_{0}\unicode[STIX]{x1D70F}_{0}/\unicode[STIX]{x1D702}_{0}, & & \displaystyle\end{eqnarray}$$

the ratio between the Eulerian flow decorrelation time $\unicode[STIX]{x1D70F}_{0}$ and the advecting time, $\unicode[STIX]{x1D702}_{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, 2016), 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 $\boldsymbol{h}_{\boldsymbol{k}}^{+}$ with a factor $\unicode[STIX]{x1D707}$ leads to flows where positive ( $\unicode[STIX]{x1D707}>1$ , $\langle H\rangle _{flow}>0$ ) or negative ( $\unicode[STIX]{x1D707}<1$ , $\langle H\rangle _{flow}<0$ ) helical structures are dominant. The resulting flow has the following exponential-like distribution of helicity (see appendix A for details):

(2.13) $$\begin{eqnarray}\displaystyle P_{0}(H)=\frac{9}{\unicode[STIX]{x03C0}}\frac{\unicode[STIX]{x1D702}_{0}^{2}}{u_{0}^{4}}\frac{\displaystyle |H|\exp \left[\frac{3H_{0}}{5-H_{0}^{2}}\frac{H\unicode[STIX]{x1D702}_{0}}{u_{0}^{2}}\right]K_{1}\left[\frac{3\sqrt{5}}{5-H_{0}^{2}}\frac{|H|\unicode[STIX]{x1D702}_{0}}{u_{0}^{2}}\right]}{\sqrt{5[5-H_{0}^{2}]}}, & & \displaystyle\end{eqnarray}$$

where $K_{\unicode[STIX]{x1D708}}(x)$ is the modified Bessel function of the second kind and $H_{0}$ is the average dimensionless helicity

(2.14) $$\begin{eqnarray}\displaystyle H_{0}\equiv \frac{\unicode[STIX]{x1D702}_{0}}{u_{0}^{2}}\langle H\rangle _{flow}=\frac{8}{3}\sqrt{\frac{2}{\unicode[STIX]{x03C0}}}\frac{\unicode[STIX]{x1D707}^{2}-1}{\unicode[STIX]{x1D707}^{2}+1}. & & \displaystyle\end{eqnarray}$$

Figure 3 shows a comparison to $256^{3}$ DNS (fourth case in Table 2) using $H_{0}=0.85$ ( $\unicode[STIX]{x1D707}\approx 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 $\sqrt{Re_{\unicode[STIX]{x1D706}}}$ for not too large $Re_{\unicode[STIX]{x1D706}}$ (Calzavarini et al. 2009). In our DNS we have $Re_{\unicode[STIX]{x1D706}}\sim 100$ and we therefore use $\unicode[STIX]{x1D702}_{0}\sim 10\unicode[STIX]{x1D702}_{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.

Figure 3. Distribution of helicity of the flow $P_{0}(H)$ for DNS with parameters given by the fourth case in 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 $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ .

3 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^{\circ }$ (anti-chiral helicoid) or $-45^{\circ }$ (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 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.

Figure 4. (a) Illustration of two isotropic helicoids with opposite helicoidality as suggested by Lord Kelvin (Kelvin 1872), anti-chiral ( $C_{0}<0$ , left panel) and co-chiral ( $C_{0}>0$ , right panel). 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\unicode[STIX]{x1D702}$ . Parameters: $St\approx 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.

3.1 Isotropic helicoids

The dynamics of an isotropic helicoid with position $\boldsymbol{x}$ , velocity $\boldsymbol{v}$ and angular velocity $\unicode[STIX]{x1D74E}$ suspended in a fluid with velocity $\boldsymbol{u}$ and vorticity $2\unicode[STIX]{x1D734}=\unicode[STIX]{x1D735}\times \boldsymbol{u}$ is governed by the following equations (Kelvin 1872; Happel & Brenner 2012; Gustavsson & Biferale 2016):

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{v}}=\frac{1}{\unicode[STIX]{x1D70F}_{p}}\left[\boldsymbol{u}(\boldsymbol{x}(t),t)-\boldsymbol{v}+\frac{2\tilde{a}}{9}C_{0}(\unicode[STIX]{x1D734}(\boldsymbol{x}(t),t)-\unicode[STIX]{x1D74E})\right], & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D74E}}=\frac{1}{\unicode[STIX]{x1D70F}_{p}}\left[\frac{10}{3}S(\unicode[STIX]{x1D734}(\boldsymbol{x}(t),t)-\unicode[STIX]{x1D74E})+\frac{5}{9\tilde{a}}C_{0}(\boldsymbol{u}(\boldsymbol{x}(t),t)-\boldsymbol{v})\right]. & \displaystyle\end{eqnarray}$$

Here dots denote time derivatives and $\boldsymbol{u}$ and $\unicode[STIX]{x1D734}$ are evaluated at the particle position $\boldsymbol{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, $\unicode[STIX]{x1D70F}_{p}$ is a relaxation time quantifying particle inertia. In the limit of $\unicode[STIX]{x1D70F}_{p}\rightarrow 0$ the particle approaches the dynamics of a tracer, $\boldsymbol{v}=\boldsymbol{u}$ and $\unicode[STIX]{x1D74E}=\unicode[STIX]{x1D734}$ . Second, $\tilde{a}=\sqrt{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 $\unicode[STIX]{x1D70F}_{p}$ . When $C_{0}\neq 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 $\unicode[STIX]{x1D734}$ , the particle accelerates either along the vorticity component, or opposite to it, see figure 4. The only constraint on the parameters is $|C_{0}|<\sqrt{27S}$ , required for the kinetic energy of the particle to remain finite. The actual size of the particle, ${\sim}\tilde{a}$ , should also be less than the smooth scale of the flow (a multiple of the Kolmogorov length $\unicode[STIX]{x1D702}$ ) 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 $\boldsymbol{v}$ and $\unicode[STIX]{x1D74E}$ 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 $\unicode[STIX]{x1D734}$ , but is also affected by the strain rate of the flow (Jeffery 1922).

Rescaling to dimensionless units $t^{\prime }=t/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ , $\boldsymbol{x}^{\prime }=\boldsymbol{x}/\unicode[STIX]{x1D702}$ , $\boldsymbol{u}^{\prime }=\boldsymbol{u}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D702}$ , $\boldsymbol{v}^{\prime }=\boldsymbol{v}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D702}$ , $\unicode[STIX]{x1D74E}^{\prime }=\unicode[STIX]{x1D74E}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ , and $\unicode[STIX]{x1D734}^{\prime }=\unicode[STIX]{x1D734}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ and dropping the primes in what follows, we can write the equations of motion for each pair of components $v_{i}$ and $\unicode[STIX]{x1D714}_{i}$ in dimensionless form:

(3.3a,b ) $$\begin{eqnarray}\displaystyle \left(\begin{array}{@{}c@{}}\dot{v}_{i}\\ \dot{\unicode[STIX]{x1D714}}_{i}\end{array}\right)=\unicode[STIX]{x1D63F}\left(\begin{array}{@{}c@{}}u_{i}-v_{i}\\ \unicode[STIX]{x1D6FA}_{i}-\unicode[STIX]{x1D714}_{i}\end{array}\right),\quad \unicode[STIX]{x1D63F}=\frac{1}{St}\left(\begin{array}{@{}cc@{}}1 & \displaystyle \frac{2C_{0}a}{9}\\ \displaystyle \frac{5C_{0}}{9a} & {\displaystyle \frac{10}{3}}S\end{array}\right). & & \displaystyle\end{eqnarray}$$

Here we have introduced the dimensionless size $a=\tilde{a}/\unicode[STIX]{x1D702}$ and the Stokes number $St=\unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ . Interpreting the two-tensor $\unicode[STIX]{x1D63F}$ as a matrix, it has two eigenvalues $d_{\pm }$ and corresponding eigenvectors $\unicode[STIX]{x1D743}_{\pm }$ given by

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle d_{\pm }=\frac{1}{18St}\left(9+30S\pm \sqrt{40C_{0}^{2}+9(3-10S)^{2}}\right)\equiv \frac{St_{\pm }}{St}, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D743}_{\pm }=\frac{1}{\sqrt{(2C_{0}a)^{2}+81(St_{\pm }-1)^{2}}}\left(\begin{array}{@{}c@{}}2C_{0}a\\ 9(St_{\pm }-1)\end{array}\right). & \displaystyle\end{eqnarray}$$

These equations are well defined for all parameter values, but in the limit of isotropic particles, $C_{0}\rightarrow 0$ , there is a complication. Taking the limit $C_{0}\rightarrow 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 $\unicode[STIX]{x1D709}_{+}$ for $S<3/10$ and $\unicode[STIX]{x1D709}_{-}$ 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\sim 1$ this scale is $St\sim 1$ , see Table 3. This has been observed in simulations of inertial particles, where the most interesting dynamics occurs around values of $St$ of order unity, see for example Fessler, Kulick & Eaton (1994), Bec et al. (2007), Falkovich & Pumir (2007). The translational dynamics of isotropic helicoids on the other hand has two characteristic inertial scales $St_{-}$ and $St_{+}$ that depend on the helicoid parameters $C_{0}$ and $S$ . These scales may be well separated in the meaning that $St_{+}/St_{-}$ can take arbitrarily large values. We therefore expect that isotropic helicoids may show significantly different behaviour depending on whether the Stokes number $St$ is of the order of $St_{-}$ or $St_{+}$ . Below, we illustrate this by numerical simulations and analysis of two different limiting cases. We remark that all statistical measurements have been made after that the particle dynamics and the flow velocity reached stationarity. Moreover, all considered statistical quantities are related to clustering in sub-viscous scales where we expect weak dependence on the Reynolds number (Bec et al. 2007).

Table 3. Rescaled eigenvalues $St_{\pm }=d_{\pm }St$ and corresponding eigenvectors $\unicode[STIX]{x1D709}_{\pm }$ in (3.4) and (3.5) for isotropic particles, $C_{0}=0$ .

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\ll St_{-}$ ):

(3.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}\sim -\frac{St}{27S-C_{0}^{2}}\left(27S\text{Tr}[\unicode[STIX]{x1D63C}^{2}]-\frac{9aC_{0}}{5}\text{Tr}[\unicode[STIX]{x1D63C}\unicode[STIX]{x1D651}]\right)+o(St), & & \displaystyle\end{eqnarray}$$

Figure 5. (a,b): fraction of particles in rotational flow regions ( $\unicode[STIX]{x1D6E5}>0$ ) as functions of (a) $St/St_{-}$ and (b) $St/St_{+}$ for the DNS given by the fourth case in Table 2. (c,d): mean fluid helicity $\langle H\rangle =2\langle \boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D734}\rangle$ along particle trajectories as functions of (c) $St/St_{-}$ and (d) $St/St_{+}$ for the DNS. The data are normalized by the helicity of the flow, $\langle H\rangle _{flow}$ which is chosen to be positive in all 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, $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(\unicode[STIX]{x1D6E5}>0)$ and $\langle H\rangle$ for tracer particles.

where $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D651}$ are matrices with elements $A_{ij}=\unicode[STIX]{x2202}_{j}u_{i}$ and $V_{ij}=\unicode[STIX]{x2202}_{j}\unicode[STIX]{x1D6FA}_{i}$ . Depending on the sign of $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}$ trajectories of close-by particles may either converge ( $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}<0$ ) or diverge ( $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}>0$ ). It is expected that particles cluster in regions where $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}<0$ , i.e. where $27S\text{Tr}[\unicode[STIX]{x1D63C}^{2}]>9aC_{0}\text{Tr}[\unicode[STIX]{x1D63C}\unicode[STIX]{x1D651}]/5$ . For heavy spherical particles $C_{0}$ is zero and particles cluster in straining regions of the flow where $\text{Tr}[\unicode[STIX]{x1D63C}^{2}]>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, $\propto \text{Tr}[\unicode[STIX]{x1D63C}\unicode[STIX]{x1D651}]$ , on the right-hand side of (3.6). As observed by Gustavsson & Biferale (2016), for a flow region with strong helical coherence, $\unicode[STIX]{x1D651}\sim c\unicode[STIX]{x1D63C}$ , particles cluster where $(27S-9acC_{0}/5)\text{Tr}[\unicode[STIX]{x1D63C}^{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 $\text{Tr}[\unicode[STIX]{x1D63C}^{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 using the flows described in § 2.3. For each set of parameters, once the fluid reaches its statistically stationary state it is seeded with $2.4\times 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. Figure 5(a,b) shows the fraction of particles in rotational regions of the flow. The data are plotted against $St/St_{-}$ in figure 5(a) and $St/St_{+}$ in figure 5(b), i.e. against the inverse of the two eigenvalues (3.4) of the dynamics (when $C_{0}=0$ the data are only plotted against $St/St_{-}$ because for this case the preferential sampling cannot depend on $St_{+}$ ). In rotational regions of a flow, the fluid gradient matrix $\unicode[STIX]{x1D63C}$ has complex eigenvalues, or equivalently, the sign of the discriminant

(3.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}=\left(\frac{\text{det}\unicode[STIX]{x1D63C}}{2}\right)^{2}-\left(\frac{\text{tr}[\unicode[STIX]{x1D63C}^{2}]}{6}\right)^{3}, & & \displaystyle\end{eqnarray}$$

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 $\unicode[STIX]{x1D6E5}>0$ , while for larger values of $St$ they instead behave as heavy inertial particles that oversample strain regions where $\unicode[STIX]{x1D6E5}<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, $\langle H(\boldsymbol{x}(t))\rangle =2\langle \boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D734}\rangle$ , as functions of $St/St_{-}$ (figure 5 c) and $St/St_{+}$ (figure 5 d). Comparing figure 5(a) and (c) 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 with low degree of helicity, while for large values of $St$ helicoids with positive values of $C_{0}$ on average sample lower degree of helicity.

Table 4. Overview of the five parameter families in the simulations in figure 5. For each family, $St$ varies over a few decades. The dynamics of the helicoids is driven by the helicoidality $C_{0}$ , shape factor $S$ and particle size $a=\tilde{a}/\unicode[STIX]{x1D702}$ which define the characteristic scales $St_{-}$ , $St_{+}$ and $\unicode[STIX]{x1D6FD}_{eff}$ as introduced in (3.4) and (3.16) respectively; for details see § 3.1. The notion of being co-chiral or anti-chiral is made in terms of the flow helicity which is always taken as on average positive in this paper. The parameter $c$ entering in the definition of $\unicode[STIX]{x1D6FD}_{eff}$ is obtained as $c=(\langle \unicode[STIX]{x1D734}^{2}\rangle /\langle \boldsymbol{u}^{2}\rangle )^{1/2}=0.1$ in DNS.

3.3 Small-scale fractal clustering

Figure 6. Spatial correlation dimension $D_{2}$ as function of (a) $St/St_{-}$ and (b) $St/St_{+}$ for the DNS and parameters of figure 5. Insets show local slope $\text{d}\log (P(r))/\text{d}\log (r)$ .

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$ :

(3.8) $$\begin{eqnarray}\displaystyle P(r)\sim r^{D_{2}},\quad r\ll 1. & & \displaystyle\end{eqnarray}$$

Figure 6(a,b) shows the behaviours of $D_{2}$ as a function of (a) $St/St_{-}$ and (b) $St/St_{+}$ for DNS with the same parameters as in figure 5. Starting at small $St/St_{-}$ , the correlation dimension is close to the spatial dimension and clustering is weak. As $St/St_{-}$ is increased, the clustering increases until it reaches a maximum around $St/St_{-}\sim 1$ . Finally, as $St/St_{-}$ is further increased, the clustering becomes weaker and saturates at the spatial dimension for large values of $St/St_{-}$ . The helicoids with $C_{0}=1.6$ show stronger clustering around $St/St_{-}\sim 1$ compared to the other cases.

3.4 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\ll St_{+}$ and $St\gg 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 as

(3.9) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D701}}_{-,i} & = & \displaystyle \frac{St_{-}}{St}(u_{-,i}-\unicode[STIX]{x1D701}_{-,i}),\end{eqnarray}$$

(3.10) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D701}}_{+,i} & = & \displaystyle \frac{St_{+}}{St}(u_{+,i}-\unicode[STIX]{x1D701}_{+,i}).\end{eqnarray}$$

Here we have changed basis for each component $i$ using

(3.11a,b ) $$\begin{eqnarray}\displaystyle \left(\begin{array}{@{}c@{}}v_{i}\\ \unicode[STIX]{x1D714}_{i}\end{array}\right)=\unicode[STIX]{x1D653}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D701}_{-,i}\\ \unicode[STIX]{x1D701}_{+,i}\end{array}\right)\quad \text{and}\quad \left(\begin{array}{@{}c@{}}u_{i}\\ \unicode[STIX]{x1D6FA}_{i}\end{array}\right)=\unicode[STIX]{x1D653}\left(\begin{array}{@{}c@{}}u_{-,i}\\ u_{+,i}\end{array}\right), & & \displaystyle\end{eqnarray}$$

where the columns of the $2\times 2$ matrix $\unicode[STIX]{x1D653}$ consist of the eigenvectors $\unicode[STIX]{x1D743}_{-}$ and $\unicode[STIX]{x1D743}_{+}$ .

The dimensionless parameter groups in (3.9) and (3.10) can be viewed as ratios of the time scales, $\unicode[STIX]{x1D70F}_{\pm }/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\equiv St/St_{\pm }$ , where $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ is the Kolmogorov time of the flow and $\unicode[STIX]{x1D70F}_{\pm }$ are two particle time scales ( $\unicode[STIX]{x1D70F}_{+}\leqslant \unicode[STIX]{x1D70F}_{-}$ ) of the isotropic helicoid ( $\unicode[STIX]{x1D70F}_{+}=\unicode[STIX]{x1D70F}_{-}=\unicode[STIX]{x1D70F}_{p}$ for a spherical particle). Equations (3.9) and (3.10) are only implicitly coupled through the trajectory dependence in $\boldsymbol{u}_{-}$ and $\boldsymbol{u}_{+}$ . We therefore expect that the two limiting cases $St\ll St_{+}$ and $St\gg 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\ll St_{+}$ ( $\unicode[STIX]{x1D70F}_{+}\ll \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ ), equation (3.10) becomes overdamped and the remaining equation (3.9) gives rise to strong preferential sampling when $St/St_{-}\sim 1$ ( $\unicode[STIX]{x1D70F}_{-}\sim \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ ) in analogue to the case $St\sim 1$ ( $\unicode[STIX]{x1D70F}_{p}\sim \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ ) for inertial spherical particles. In the second limit $St\gg St_{-}$ ( $\unicode[STIX]{x1D70F}_{-}\gg \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ ), equation (3.9) becomes underdamped and $\unicode[STIX]{x1D701}_{-}$ can to lowest order be approximated by its mean value. The remaining equation (3.10) gives rise to strong preferential sampling when $St/St_{+}\sim 1$ ( $\unicode[STIX]{x1D70F}_{+}\sim \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ ). Below we discuss the two limiting cases in more detail.

3.4.1 Case $St\ll St_{+}$

Consider first $St\ll St_{+}$ with general values of $St_{-}$ . For the acceleration to remain finite in (3.10) in this limit, we must have $0\sim u_{+,i}-\unicode[STIX]{x1D701}_{+,i}$ . In terms of the original coordinates, this condition gives the following constraint:

(3.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D74E}-\unicode[STIX]{x1D734}=\frac{9}{2C_{0}a}(1-St_{-})(\boldsymbol{u}-\boldsymbol{v}). & & \displaystyle\end{eqnarray}$$

Using (3.12) and its time derivative to replace $\unicode[STIX]{x1D74E}$ and $\dot{\unicode[STIX]{x1D74E}}$ in (3.9), and reverting to the original coordinates, we obtain

(3.13) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{v}}=\frac{St_{-}}{St}(\boldsymbol{u}-\boldsymbol{v})+\frac{1}{3(St_{-}-St_{+})}\left[3(St_{-}-1)\dot{\boldsymbol{u}}-\frac{2C_{0}a}{3}\dot{\unicode[STIX]{x1D734}}\right]. & & \displaystyle\end{eqnarray}$$

Thus, a single equation determines the velocity of the isotropic helicoid in the limit $St\ll 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 (3.13) still shows the same results as those obtained by letting $C_{0}=0$ and $St\ll St_{+}$ in the original equation (3.3). When $C_{0}=0$ and $S>3/10$ , equation (3.13) simplifies to $\dot{\boldsymbol{v}}=(\boldsymbol{u}-\boldsymbol{v})/St$ while the rotational dynamics is overdamped, $\unicode[STIX]{x1D74E}-\unicode[STIX]{x1D734}=0$ . When $C_{0}=0$ and $S<3/10$ on the other hand, the condition $St\ll St_{+}=1$ implies that the translational dynamics is overdamped and (3.13) simplifies to $\dot{\boldsymbol{v}}=(\boldsymbol{u}-\boldsymbol{v})/St+\dot{\boldsymbol{u}}$ . This equation relaxes to the overdamped limit $\boldsymbol{v}=\boldsymbol{u}$ after a short initial transient on the time scale of order $St\ll 1$ .

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 approximate $\dot{\boldsymbol{u}}\sim \text{D}_{t}\boldsymbol{u}$ and $\dot{\unicode[STIX]{x1D734}}\sim \text{D}_{t}\unicode[STIX]{x1D734}$ , where $\text{D}_{t}\equiv \unicode[STIX]{x2202}_{t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})$ are advective derivatives. Second, in a helical flow $\unicode[STIX]{x1D734}$ and $\boldsymbol{u}$ tend to be aligned and we may approximate $\unicode[STIX]{x1D734}\sim c\boldsymbol{u}$ with some proportionality constant $c$ . Using these approximations, equation (3.13) simplifies to

(3.14) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{v}}=\frac{St_{-}}{St}(\boldsymbol{u}-\boldsymbol{v})+\frac{9(St_{-}-1)-2cC_{0}a}{9(St_{-}-St_{+})}\text{D}_{t}\boldsymbol{u}. & & \displaystyle\end{eqnarray}$$

If we define an effective Stokes number $St_{eff}$ and an effective density parameter $\unicode[STIX]{x1D6FD}_{eff}$ :

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle St_{eff}=\frac{St}{St_{-}}, & \displaystyle\end{eqnarray}$$

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FD}_{eff}=\frac{9(St_{-}-1)-2cC_{0}a}{9(St_{-}-St_{+})}, & \displaystyle\end{eqnarray}$$

then (3.14) becomes identical to the equation of motion of small spherical particles (Maxey & Riley 1983) with sub-dominant terms neglected:

(3.17) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{v}}=\frac{1}{St}(\boldsymbol{u}-\boldsymbol{v})+\unicode[STIX]{x1D6FD}\text{D}_{t}\boldsymbol{u}. & & \displaystyle\end{eqnarray}$$

For spherical particles the density is characterized by $\unicode[STIX]{x1D6FD}$ : $0\leqslant \unicode[STIX]{x1D6FD}<1$ corresponds to heavy particles more dense than the fluid, $\unicode[STIX]{x1D6FD}=1$ corresponds to neutrally buoyant particles and $1<\unicode[STIX]{x1D6FD}\leqslant 3$ corresponds to particles lighter than the fluid. We remark that $\unicode[STIX]{x1D6FD}_{eff}$ in (3.16) is not constrained to the interval $0\leqslant \unicode[STIX]{x1D6FD}\leqslant 3$ as the case of spherical particles, $\unicode[STIX]{x1D6FD}_{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\ll St_{+}$ in figure 5(a). The values of the effective density parameter $\unicode[STIX]{x1D6FD}_{eff}$ in (3.16) are quoted in Table 4 for our parameters. Only the case $C_{0}=1.6$ has $\unicode[STIX]{x1D6FD}_{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 $\unicode[STIX]{x1D6FD}_{eff}$ is close to zero. In these cases the helicoids preferentially sample strain regions of the fluid with strongest effect around $St_{eff}\sim O(1)$ . Finally, the case $C_{0}=5$ has $\unicode[STIX]{x1D6FD}_{eff}\approx 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}\sim 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_{-}\sim 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}=\pm 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\ll St_{-}$ and to get good statistics in this limit. We therefore use the stochastic model to study the approximations.

Figure 7. Comparison of stochastic model simulations for (a) the probability of being in a rotational region, $P(\unicode[STIX]{x1D6E5}>0)$ , and (b) the average helicity, $\langle H\rangle$ , to results in the limiting case $St\ll St_{+}$ discussed in § 3.4.1. Markers show results from simulations of the full dynamics (3.1) and (3.2) and solid lines show the approximation (3.13) evaluated using stochastic model simulations. The upper bound $St/St_{-}\sim 20$ corresponds to $St\sim 0.2St_{+}$ . (c,d) An enlargement of the region where $St\ll St_{-}\ll St_{+}$ . Solid lines and markers as in (a) and (b). Dashed lines show the approximation (3.13) with $\dot{\boldsymbol{u}}$ and $\dot{\unicode[STIX]{x1D734}}$ replaced by $\text{D}\boldsymbol{u}/\text{D}t$ and $\text{D}\unicode[STIX]{x1D734}/\text{D}t$ . Dash-dotted lines show the approximation (3.14). Helicoid parameters corresponding to Table 4 and figure 5: $S=1$ , $a=10$ , hollow orange circles (anti-chiral, $C_{0}=-5$ ) and filled red circles (co-chiral, $C_{0}=5$ ).  $S=0.1$ , $a=30$ hollow light blue boxes (anti-chiral, $C_{0}=-1.6$ ) and filled blue boxes (co-chiral, $C_{0}=1.6$ ). Black dashed lines show results for tracer particles.

Figure 7(a,b) shows simulation results for isotropic helicoids with $St\ll 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 fix $H_{0}=0.85$ to match the distribution of flow helicity to that of the DNS and we take $\unicode[STIX]{x1D702}_{0}/\unicode[STIX]{x1D702}=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, $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\equiv \langle \text{Tr}(\unicode[STIX]{x1D63C}\unicode[STIX]{x1D63C}^{\text{T}})\rangle _{flow}^{-1/2}=\unicode[STIX]{x1D702}_{0}/(\sqrt{5}u_{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). Comparing figure 7(a) and (b) to figure 5(a) and (c) 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 of finding helicoids with $C_{0}=1.6$ in rotational regions with $\unicode[STIX]{x1D6E5}>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\ll St_{+}$ . Figure 7(c,d) shows numerical simulations of (3.13) with $\dot{\boldsymbol{u}}$ and $\dot{\unicode[STIX]{x1D734}}$ replaced by $\text{D}_{t}\boldsymbol{u}$ and $\text{D}_{t}\unicode[STIX]{x1D734}$ , which approach the results of (3.13) when $St/St_{-}\ll 1$ as expected. Finally, figure 7(c,d) also shows the approximation (3.14) using $\unicode[STIX]{x1D734}=c\boldsymbol{u}$ with $c=(\langle \unicode[STIX]{x1D734}^{2}\rangle /\langle \boldsymbol{u}^{2}\rangle )^{1/2}=\sqrt{5}/20$ for the stochastic model (rescaled using $\unicode[STIX]{x1D702}_{0}=10\unicode[STIX]{x1D702}$ ). We observe that the prediction using spherical particles reproduces the average sampling of helicity well for $St\ll St_{-}$ (figure 7 d), while it does not work as well for the probability of being in vortex regions (figure 7 c). 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}\pm 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. Snapshot of positions of isotropic helicoids in stochastic model simulations for a flow with $Ku=10$ and maximal helicity $H_{0}\approx 2.1$ ( $\unicode[STIX]{x1D707}\rightarrow \infty$ ), and particle parameters $St=0.9St_{-}$ , $a=30$ , $S=0.1$ and $C_{0}=-1.6$ (a) and $C_{0}=1.6$ (b). The history of the underlying flow is identical for the two simulations. Coordinate axes are non-dimensionalized using $\unicode[STIX]{x1D702}_{0}$ .

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.

Figure 9. Stochastic model simulations of (a) correlation dimension $D_{2}$ , (b) phase-space Lyapunov dimension $D_{L}$ and (c) the relative amplitude of average compressibility along particle trajectories estimated using the sum of the first three Lyapunov exponents, $\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}\rangle =\unicode[STIX]{x1D706}_{1}+\unicode[STIX]{x1D706}_{2}+\unicode[STIX]{x1D706}_{3}$ and the approximation (3.6),  $\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}^{(1)}\rangle$ . Parameters as in figure 7 with the addition of spherical particles (black asterisks).

Using $D_{2}$ to quantify the 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, Kaplan & Yorke 1983). Denoting by $\unicode[STIX]{x1D706}_{1}\geqslant \unicode[STIX]{x1D706}_{2}\geqslant \ldots \geqslant \unicode[STIX]{x1D706}_{9}$ the nine Lyapunov exponents of the equation system (3.1) and (3.2) together with $\dot{\boldsymbol{x}}=\boldsymbol{v}$ , the Lyapunov dimension is

(3.18) $$\begin{eqnarray}\displaystyle D_{L}=K+\mathop{\sum }_{i=1}^{K}\unicode[STIX]{x1D706}_{i}/|\unicode[STIX]{x1D706}_{K+1}|, & & \displaystyle\end{eqnarray}$$

where $K$ is the largest integer such that $\sum _{i=1}^{K}\unicode[STIX]{x1D706}_{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}\geqslant 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, $\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}\rangle =\unicode[STIX]{x1D706}_{1}+\unicode[STIX]{x1D706}_{2}+\unicode[STIX]{x1D706}_{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\ll St_{+}$ ), but requires first-order corrections in $St/St_{+}$ to the condition in (3.12). Moreover, we need to consider $St/St_{-}\ll 1$ . Consistently with this limit, we find that the average of the approximation (3.6) tends to approach $\langle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}\rangle$ as $St/St_{-}$ is reduced and that the two expressions agree approximately for $St/St_{-}\sim 0.1$ .

3.4.2 Case $St\gg St_{-}$

We now consider the second limiting case, $St\gg St_{-}$ with general values of $St_{+}$ . In this limit $\unicode[STIX]{x1D701}_{-,i}$ in (3.9) responds slowly to changes in the flow compared to $\unicode[STIX]{x1D701}_{+,i}$ . Due to the symmetries of the underlying flow, we expect the averages of $u_{-,i}$ and consequently of $\unicode[STIX]{x1D701}_{-,i}$ to vanish. We therefore replace $\unicode[STIX]{x1D701}_{-,i}$ by its vanishing average, $\unicode[STIX]{x1D701}_{-,i}=0$ , which gives the following constraint on $\unicode[STIX]{x1D74E}$ :

(3.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D74E}=\frac{9(St_{+}-1)}{2aC_{0}}\boldsymbol{v}. & & \displaystyle\end{eqnarray}$$

Inserting this constraint and its time derivative into (3.10), gives the following equation for the velocity:

(3.20) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{v}}=\frac{St_{+}}{St}\left[\frac{3(3St_{+}-10S)\boldsymbol{u}+2aC_{0}\unicode[STIX]{x1D734}}{9(St_{+}-St_{-})}-\boldsymbol{v}\right]. & & \displaystyle\end{eqnarray}$$

As for the first limiting case we can consider a helical flow with $\unicode[STIX]{x1D734}\sim c\boldsymbol{u}$ to obtain

(3.21) $$\begin{eqnarray}\displaystyle \dot{\boldsymbol{v}}=\frac{St_{+}}{St}[\unicode[STIX]{x1D6FD}_{eff}\boldsymbol{u}-\boldsymbol{v}], & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}_{eff}$ is the parameter in (3.16) occurring in the limit of $St\ll St_{-}$ . Thus, in the limit $St\gg St_{-}$ the equation of motion for helicoids in a helical flow is like a Stokes drag with effective Stokes number $St/St_{+}$ and with a rescaled amplitude of the fluid velocity.

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, $\unicode[STIX]{x1D6FD}_{eff}\boldsymbol{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 5 b). The helicoids with positive values of $C_{0}$ on the other hand have $\unicode[STIX]{x1D6FD}_{eff}\sim 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\sim St_{+}$ is similar to that of spherical particles. However, the approximation (3.21) does not work as well for $C_{0}=5$ : even though $\unicode[STIX]{x1D6FD}_{eff}\approx 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\ll St_{+}$ , we conclude that the approximations (for example $\unicode[STIX]{x1D734}\sim c\boldsymbol{u}$ ) 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}\approx 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\sim 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\rightarrow 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.

Comparison to stochastic model

Figure 10 shows simulation results in the range $St\gg St_{-}$ for isotropic helicoids in the stochastic model, plotted against $St/St_{+}$ . Similar to the case $St\ll St_{+}$ , comparison between the full model data in figure 10(a) and (b) to the corresponding DNS data in figure 5(b) and (d) shows similar trends as functions of $St$ , while the details differ in some ranges. The co-chiral helicoids, having $\unicode[STIX]{x1D6FD}_{eff}\sim 1$ show preferential sampling of straining regions with $St\sim St_{+}$ , while the anti-chiral helicoids have negligible preferential sampling of straining regions.

Figure 10. Comparison of stochastic model simulations for (a) the probability of being in a rotational region, $P(\unicode[STIX]{x1D6E5}>0)$ , and (b) the average helicity, $\langle H\rangle$ , to results in the limiting case $St\gg St_{-}$ discussed in § 3.4.2. Markers show results from simulations of the full dynamics (3.1) and (3.2), same parameters as in figure 7, solid lines show the approximation (3.20) and dash-dotted lines show the refined approximation (3.21) for flows where fluid vorticity and velocity align. The lower bound $St/St_{+}\sim 0.05$ corresponds to $St\sim 5St_{-}$ .

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\gg St_{-}$ we observe a quantitative agreement for $C_{0}=\pm 1.6$ and qualitative agreement when $C_{0}=\pm 5$ . Figure 10 also shows that the approximation (3.21) based on $\unicode[STIX]{x1D734}=c\boldsymbol{u}$ (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\sim 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\gg 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\sim St_{-}$ . As a consequence, $D_{2}$ does not show any clustering around $St\sim 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\sim St_{+}$ shows a scaling $P(r)\sim r^{D_{2}(r)}$ with local exponent $D_{2}(r)<3$ for a range of $r\ll 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).

3.4.3 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, 2016). The calculation is outlined in appendix B. The resulting mean helicity becomes to second order in $Ku$ ,

(3.22) $$\begin{eqnarray}\displaystyle \langle H\rangle & = & \displaystyle H_{0}+\frac{3Ku^{2}St}{4(5C_{0}^{2}-27(1+2St)(5S+3St))(10C_{0}^{2}-27(1+St)(10S+3St))^{2}}\nonumber\\ \displaystyle & & \displaystyle \times \,\left\{-(25+6H_{0}^{2})C_{0}\overline{a}\left[50C_{0}^{4}+135C_{0}^{2}(10S(St-2)+3St(5St-1))\right.\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\,729(50S^{2}(St-1)+18St^{3}+15SSt(5St-1))\right]\nonumber\\ \displaystyle & & \displaystyle +\,30H_{0}\left[\overline{a}^{2}(90C_{0}^{4}St-486C_{0}^{2}(5SSt-3St^{3}))+50C_{0}^{4}(10S+9St)\right.\nonumber\\ \displaystyle & & \displaystyle \left.\left.-\,675C_{0}^{2}(40S^{2}+42SSt+9St^{2})+729(5S+3St)(10S+3St)^{2}\right]\right\},\end{eqnarray}$$

where $\overline{a}=\tilde{a}/\unicode[STIX]{x1D702}_{0}$ (velocity and position are made dimensionless in terms of $u_{0}$ and $\unicode[STIX]{x1D702}_{0}$ ). Figure 11(a) shows the analytical solution for the mean helicity together with data for $Ku=10$ . In order to compensate for the different magnitudes of relevant time scales in flows with $Ku=10$ and flows with small values of $Ku$ , the parameters $Ku$ and $St$ which depend on the correlation time of the flow have been rescaled in (3.22). 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 11 c). 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).

Figure 11. Evaluation of theory (3.22) for mean helicity $\langle H\rangle$ for small values of $Ku$ with rescaled parameters $Ku$ and $St$ , see text. (a) Comparison of the theory (3.22) to simulation data for $Ku=10$ . Markers correspond to the data in figure 7(b) (without division by $H_{0}$ and plotted against $St$ ). Solid lines show (3.22). The five parameter combinations of $C_{0}$ , $S$ and $a$ are plotted as five lines with the colours of the corresponding markers. (b) Heat map of the theory (3.22) for the deviation of average helicity, $\langle H\rangle$ , from that of the flow, $\langle H\rangle _{flow}=H_{0}$ , plotted against $St$ and $C_{0}$ for $S=0.1$ , $a=30$ and $H_{0}=0.85$ . Dashed lines correspond to the curves with $C_{0}=-1.6$ (light blue), $C_{0}=0$ (black) and $C_{0}=1.6$ (blue) shown in (a). (c) Same as (b) but for a neutral flow, $H_{0}=0$ .

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, 2016). The theory relies on the additional constraint $St\ll 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).

4 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), $\unicode[STIX]{x1D6FD}_{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}\sim 5$ . This implies an estimated helicoid size $\tilde{a}$ larger than the smooth scale (approximately $10\unicode[STIX]{x1D702}$ ), 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 $\tilde{a}$ 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 $\unicode[STIX]{x1D6FD}_{eff}$ , as exemplified by comparing $|C_{0}|=5$ and $|C_{0}|=1.6$ in figure 5. By optimizing $\unicode[STIX]{x1D6FD}_{eff}$ it is possible to find helicoids with smaller values of $a$ that behave as light particles, i.e. have $\unicode[STIX]{x1D6FD}_{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.