Helicoidal particles in turbulent flows with multi-scale helical injection

TL;DR

This study investigates how multi-scale helical injection influences turbulent flow regimes, particle dynamics, and flow structure, revealing tunable helicity transfer and preferential sampling by isotropic helicoids through numerical and theoretical analysis.

Contribution

It introduces a comprehensive analysis of turbulent flows with multi-scale helical injection, including new regimes and particle-flow interactions, supported by numerical simulations and theoretical models.

Findings

Three distinct inertial transfer regimes identified.

Flow regimes with Kolmogorov-like cascades and tunable helicity.

Isotropic helicoids preferentially sample co-chiral or anti-chiral structures.

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.

Figures

Tables

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

(I)	(II)	(III)
$y>2$
$ϵ\left(k\right)=const.$	$ϵ\left(k\right)=const.$	$ϵ\left(k\right)\sim k^{1-y}$
$h\left(k\right)=const.$	$h\left(k\right)\sim k^{2-y}$	$h\left(k\right)\sim k^{2-y}$

Table 2. Table 2: Parameters of the numerical simulations: grid resolution N 3 superscript 𝑁 3 N^{3} , Kolmogorov length scale η 𝜂 \eta in simulation units (SU), grid spacing Δ ​ x = 2 ​ π / N Δ 𝑥 2 𝜋 𝑁 \Delta x=2\pi/N (SU), time step Δ ​ t Δ 𝑡 \Delta t (SU), kinematic viscosity ν 𝜈 \nu (SU), Kolmogorov eddy turn-over time τ η = ( ν / ϵ ) 1 / 2 subscript 𝜏 𝜂 superscript 𝜈 italic-ϵ 1 2 \tau_{\eta}=(\nu/\epsilon)^{1/2} with ϵ italic-ϵ \epsilon the energy dissipation rate (SU), forcing correlation time τ S subscript 𝜏 𝑆 \tau_{S} (in units of Δ ​ t Δ 𝑡 \Delta t ), forcing power law exponent y 𝑦 y , hyper-viscosity parameter α 𝛼 \alpha , number of helicoids per each family N h subscript 𝑁 ℎ N_{h} .

Regime	$N^3$	$\eta$	$\Delta x$	$\Delta t$	$\nu$	${\tau }_{\eta }$	${\tau }_S$	$y$	$\alpha$	$N_h$
(I)	${512}^3$	0.008	0.012	0.0003	1.9 $\times {10}^{-7}$	$0.03$	$1$	4	2	-
(II)	${512}^3$	0.008	0.012	0.0003	1.9 $\times {10}^{-7}$	$0.03$	$1$	3/2	2	-
(III)	${512}^3$	0.008	0.012	0.0003	1.9 $\times {10}^{-7}$	$0.012$	$1$	$-1$	2	-
(III)	${256}^3$	0.016	0.024	0.0006	0.0052	$0.025$	$500$	$-2/3$	1	$2.4\times {10}^6$

Table 3. Table 3: Rescaled eigenvalues St ± = d ± ​ St subscript St plus-or-minus subscript 𝑑 plus-or-minus St \operatorname{St}_{\pm}=d_{\pm}\operatorname{St} and corresponding eigenvectors ξ ± subscript 𝜉 plus-or-minus \xi_{\pm} in Eqs. ( 16 ) and ( 17 ) for isotropic particles, C 0 = 0 subscript 𝐶 0 0 C_{0}=0 .

Condition	${\mathrm{St}}_{-}$	${\mathrm{St}}_{+}$	${\xi }_{-}$	${\xi }_{+}$
$S>3/10$	$1$	$10S/3$	$(1,0)$	$(0,1)$
	$10S/3$	$1$	$(0,-1)$	$(1,0)$

Table 4. Table 4: Overview of the five parameter families in the simulations in Fig.  5 . For each family, St St \operatorname{St} varies over a few decades. The dynamics of the helicoids is driven by the helicoidality C 0 subscript 𝐶 0 C_{0} , shape factor S 𝑆 S , and particle size a = a ~ / η 𝑎 ~ 𝑎 𝜂 a={\tilde{a}}/\eta which define the characteristic scales St − subscript St \operatorname{St}_{-} , St + subscript St \operatorname{St}_{+} and β eff subscript 𝛽 eff \beta_{\rm eff} as introduced in Eqs. ( 16 ) and ( 27 ) respectively; for details see Sec.  3.1 . The notion of being co-chiral or anti-chiral is made in terms of the flow helicity which is always taken in average positive in this paper. The parameter c 𝑐 c entering in the definition of β eff subscript 𝛽 eff \beta_{\rm eff} is obtained as c = ( ⟨ 𝛀 2 ⟩ / ⟨ 𝒖 2 ⟩ ) 1 / 2 = 0.1 𝑐 superscript delimited-⟨⟩ superscript 𝛀 2 delimited-⟨⟩ superscript 𝒖 2 1 2 0.1 c=(\langle\boldsymbol{\Omega}^{2}\rangle/\langle\boldsymbol{u}^{2}\rangle)^{1/2}=0.1 in DNS.

Helicoids type	$C_0$	$S$	$a$	${\mathrm{St}}_{-}$	${\mathrm{St}}_{+}$	${\beta }_{\mathrm{eff}}$
Anti-chiral	$-5$	$1$	$10$	$0.058$	$4.3$	$-0.1$
	$-1.6$	$0.1$	$30$	$0.013$	$1.3$	$-0.2$
Neutral	$0$	$1$	$-$	$1$	$3.3$	$0$
Co-chiral	$1.6$	$0.1$	$30$	$0.013$	$1.3$	$1.7$
	$5$	$1$	$10$	$0.058$	$4.3$	$0.5$

Equations

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}=-\nabla p+\nu \Delta \mathbf{u}+\mathbf{f},\nabla \cdot \mathbf{u}=0,$$

$$\mathbf{u}(\mathbf{x},t)=\sum _k[u_{\mathbf{k}}^{+}(t){\mathbf{h}}_{\mathbf{k}}^{+}+u_{\mathbf{k}}^{-}(t){\mathbf{h}}_{\mathbf{k}}^{-}]e^{-i\mathbf{k}\cdot \mathbf{x}},$$

$$E=\sum _{\mathbf{k}}|u_{\mathbf{k}}^{+}{|}^2+|u_{\mathbf{k}}^{-}{|}^2,H=\sum _{\mathbf{k}}k(|u_{\mathbf{k}}^{+}{|}^2-|u_{\mathbf{k}}^{-}{|}^2).$$

$$E(k)=E^{+}(k)+E^{-}(k),H(k)=k[E^{+}(k)-E^{-}(k)].$$

$$E^{+}(k)+E^{-}(k)\sim C_E{ϵ}^{2/3}{k}^{-5/3},E^{+}(k)-E^{-}(k)\sim C_{H}h{ϵ}^{-1/3}{k}^{-8/3},$$

$$E^{\pm }(k)\sim C_E{ϵ}^{2/3}{k}^{-5/3}\pm C_{H}h{ϵ}^{-1/3}{k}^{-8/3}.$$

$$\mathbf{f}(\mathbf{x},t)=\sum _{\mathbf{k}}{f}_{\mathbf{k}}^{+}(t){\mathbf{h}}_{\mathbf{k}}^{+}{e}^{-i\mathbf{k}\cdot \mathbf{x}},$$

$$⟨f_{\mathbf{k}}^{+}(t)f_{{\mathbf{k}}^{\prime }}^{+}(t^{\prime })⟩=D_0{k}^{1-d-y}\delta (t-t^{\prime }){\delta }_{\mathbf{k},{\mathbf{k}}^{\prime }},$$

$$ϵ(k)\sim \sum _{|\mathbf{q}|<k}|\mathbf{q}{|}^{1-d-y},h(k)\sim \sum _{|\mathbf{q}|<k}|\mathbf{q}{|}^{2-d-y}.$$

$$\{\begin{array}{l}{E}^{\pm }(k)\sim C_E^I{ϵ}^{2/3}{k}^{-5/3}\pm C_H^{I}h{ϵ}^{-1/3}{k}^{-8/3},y>2\\ E^{\pm }(k)\sim C_E^{II}{ϵ}^{2/3}{k}^{-5/3}\pm C_H^{II}{k}^{2-y}{ϵ}^{-1/3}{k}^{-8/3},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},y<1,\end{array}$$

$$\{\begin{array}{l}{R}^I(k)\sim k^{-1},y>2\\ R^{II}(k)\sim k^{1-y},1<y<2\\ R^{III}(k)\sim const.,y<1,\end{array}$$

$$\mathrm{Ku}=u_0{\tau }_0/{\eta }_0,$$

$$P_0(H)$$

$$H_0\equiv \frac{{\eta }_0}{u_0^2}⟨H{⟩}_{\mathrm{flow}}=\frac{8}{3}\sqrt{\frac{2}{\pi }}\frac{{\mu }^2-1}{{\mu }^2+1}.$$

$$\begin{array}{rl}\dot{\mathbf{v}}=& \frac{1}{{\tau }_{\mathrm{p}}}\left[\mathbf{u}(\mathbf{x}(t),t)-\mathbf{v}+\frac{2\tilde{a}}{9}{C}_0(\mathbf{\Omega }(\mathbf{x}(t),t)-\mathbf{\omega })\right]\\ \dot{\mathbf{\omega }}=& \frac{1}{{\tau }_{\mathrm{p}}}\left[\frac{10}{3}S(\mathbf{\Omega }(\mathbf{x}(t),t)-\mathbf{\omega })+\frac{5}{9\tilde{a}}{C}_0(\mathbf{u}(\mathbf{x}(t),t)-\mathbf{v})\right].\end{array}$$

$$\left(\begin{array}{c}{\dot{v}}_i\\ {\dot{\omega }}_i\end{array}\right)=\text{mathsfbi}D\left(\begin{array}{c}{u}_i-v_i\\ {\Omega }_i-{\omega }_i\end{array}\right),\text{mathsfbi}D=\frac{1}{\mathrm{St}}\left(\begin{array}{cc}1& \frac{2C_{0}a}{9}\\ \frac{5C_0}{9a}& \frac{10}{3}S\end{array}\right).$$

$$d_{\pm }$$

$${\mathbf{\xi }}_{\pm }$$

$$\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}\sim-\,\frac{\operatorname{St}}{27S-C_{0}^{2}}\Big{(}27S\mbox{Tr}[\mathsfbi{A}^{2}]-\frac{9aC_{0}}{5}\mbox{Tr}[\mathsfbi{A}\mathsfbi{V}]\Big{)}+\textit{o}(\operatorname{St})\,,$$

$$\Delta ={\left(\frac{\det \text{mathsfbi}A}{2}\right)}^2-{\left(\frac{\mathrm{tr}[\text{mathsfbi}{A}^2]}{6}\right)}^3,$$

$$P(r)\sim r^{D_2},r\ll 1.$$

$${\dot{\zeta }}_{-,i}$$

$${\dot{\zeta }}_{+,i}$$

$$\left(\begin{array}{c}{v}_i\\ {\omega }_i\end{array}\right)=\text{mathsfbi}X\left(\begin{array}{c}{\zeta }_{-,i}\\ {\zeta }_{+,i}\end{array}\right)\text{mbox}and\left(\begin{array}{c}{u}_i\\ {\Omega }_i\end{array}\right)=\text{mathsfbi}X\left(\begin{array}{c}{u}_{-,i}\\ u_{+,i}\end{array}\right),$$

$$\mathbf{\omega }-\mathbf{\Omega }=\frac{9}{2C_{0}a}(1-{\mathrm{St}}_{-})(\mathbf{u}-\mathbf{v}).$$

$$\dot{\mathbf{v}}=\frac{{\mathrm{St}}_{-}}{\mathrm{St}}(\mathbf{u}-\mathbf{v})+\frac{1}{3({\mathrm{St}}_{-}-{\mathrm{St}}_{+})}\left[3({\mathrm{St}}_{-}-1)\dot{\mathbf{u}}-\frac{2C_{0}a}{3}\dot{\mathbf{\Omega }}\right].$$

$$\dot{\mathbf{v}}=\frac{S{t}_{-}}{St}(\mathbf{u}-\mathbf{v})+\frac{9(S{t}_{-}-1)-2c{C}_{0}a}{9(S{t}_{-}-S{t}_{+})}{\mathrm{D}}_t\mathbf{u}.$$

$${\mathrm{St}}_{\mathrm{eff}}$$

$${\beta }_{\mathrm{eff}}$$

$$\dot{\mathbf{v}}=\frac{1}{St}(\mathbf{u}-\mathbf{v})+\beta {\mathrm{D}}_t\mathbf{u}.$$

$$D_{\mathrm{L}}=K+\sum _{i=1}^K{\lambda }_i/|{\lambda }_{K+1}|,$$

$$\mathbf{\omega }=\frac{9({\mathrm{St}}_{+}-1)}{2a{C}_0}\mathbf{v}.$$

Full text

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 (Frisch, 1995; Moffatt & Tsinober, 1992; Chen et al., 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; Kedia et al., 2016; Laing et al., 2015; Kerr, 2015). Helicity can be introduced in a flow by a stirring mechanism that breaks mirror symmetry and its effects on the turbulent energy cascade in three spatial dimensions have been widely studied since the pioneering work of Brissaud et al. (1973) (see Borue & Orszag (1997); Pelz et al. (1985); Kerr (1987); Kholmyansky et al. (1991); Kit et al. (1987) for other contributions). In geophysical flows, helicity plays an important role in the atmospheric Ekman layer, where there exist arguments supporting a turbulent helicity cascade in the logarithmic range of the boundary layer (Deusebio & Lindborg, 2014; Koprov et al., 2005; Kurgansky, 2017). Recent experimental advancements allowed 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 et al., 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 et al., 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 et al., 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, 1871; Gustavsson & Biferale, 2016), i.e. point-like particles whose dynamics is isotropic but breaks mirror symmetry. By using both direct numerical simulations (DNS) and a stochastic model for the Eulerian advecting velocity field (Gustavsson & Mehlig, 2016), we show that isotropic helicoids possess highly non-trivial preferential sampling of the underlying helical flow properties depending on the particle parameters. The paper is organized as follows. In Sec. 2 we describe the Eulerian part, discussing the different regimes for different helical injection power spectra and we present numerical simulations of the different regimes. In Sec. 3 we introduce the isotropic helicoids and their dynamical equations. We discuss the existence of two new scales of the Stokes number, $\operatorname{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 $\operatorname{St}$ is either much smaller than $\operatorname{St}_{+}$ or much larger than $\operatorname{St}_{-}$. We conclude the paper in Sec. 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:

[TABLE]

where $\nu$ is the kinematic viscosity and ${\boldsymbol{f}}$ is a parity-breaking external forcing with energy injection rate $\epsilon=\langle{\boldsymbol{u}}\boldsymbol{\cdot}{\boldsymbol{f}}\rangle$ and helicity injection rate $h=\langle{\boldsymbol{u}}\boldsymbol{\cdot}({\boldsymbol{\nabla}}\times{\boldsymbol{f}})+2{\boldsymbol{\Omega}}\boldsymbol{\cdot}{\boldsymbol{f}}\rangle$, where $2{\boldsymbol{\Omega}}={\boldsymbol{\nabla}}\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):

[TABLE]

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

[TABLE]

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 $\Delta k=2\pi/L$, such that the energy and helicity spectra become

[TABLE]

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 $\tau_{E}(r)\sim\epsilon^{-1/3}r^{2/3}$, we have for the semi-sum and semi-difference of the spectral components (Chen et al., 2003a):

[TABLE]

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

[TABLE]

It is known that for large-scale energy and helicity injection the Kolmogorov-like scaling (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 et al., 1977; Seoud & Vassilicos, 2007).

2.2 Multi-scale energy and helicity injections

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

[TABLE]

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

[TABLE]

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\pi/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_{\eta}$. Using this forcing, the energy and helicity injection rates up to the scale $k<k_{max}$ can be estimated as:

[TABLE]

By considering spherical symmetry, the sums in (8) 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, $\epsilon(k)\to const.$ and $h(k)\to const.$ when $k\to\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 (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, $\epsilon(k)\sim k^{1-y}$ and $h(k)\sim k^{2-y}$. The three regimes are summarized in Table 1.

As a result, the spectral properties (6) are valid only for regime (I), and we can summarize the scaling for all three different regimes as follows:

[TABLE]

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

[TABLE]

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: $\tau_{E}(k)\propto k^{-2/3}\epsilon(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 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 $\nu_{\alpha}\Delta^{\alpha}u$ as the viscous term, with $\alpha=2$. The external forcing ${\boldsymbol{f}}$ in (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 Fig. 1, we present four panels with the results for (a) the energy spectrum, (b) the helicity spectrum, (c) the energy flux and (d) the helicity flux as functions of $k$ and for the three representatives values of $y$. 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 (9) and energy fluxes (8) 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 (10). In particular, in Fig. 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 Fig. 2(b) 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 (9).

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}=\boldsymbol{\nabla}\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 $\eta_{0}$ and an exponential time-correlation function with decay rate, $\tau_{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

[TABLE]

the ratio between the Eulerian flow decorrelation time $\tau_{0}$ and the advecting time, $\eta_{0}/u_{0}$. The Kubo number can be seen as a dimensionless correlation time of the flow. If $\operatorname{Ku}$ tends to zero a white-noise flow is approached and if $\operatorname{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 $\operatorname{Ku}$ and qualitatively at large $\operatorname{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). Weighting the positive modes $\boldsymbol{h}^{+}_{\boldsymbol{\footnotesize k}}$ with a factor $\mu$ leads to flows where positive ($\mu>1$, $\langle H\rangle_{\rm flow}>0$) or negative ($\mu<1$, $\langle H\rangle_{\rm flow}<0$) helical structures are dominant. The resulting flow has the following exponential-like distribution of helicity (see Appendix A for details):

[TABLE]

where $K_{\nu}(x)$ is the modified Bessel function of the second kind and $H_{0}$ is the average dimensionless helicity

[TABLE]

Fig. 3 shows a comparison to $256^{3}$ DNS (fourth case in Table 2) using $H_{0}=0.85$ ($\mu\approx 1.5$) to make the shape of the distribution (12) 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{{\rm Re}_{\lambda}}$ for not too large ${\rm Re}_{\lambda}$ (Calzavarini et al., 2009). In our DNS we have ${\rm Re}_{\lambda}\sim 100$ and we therefore use $\eta_{0}\sim 10\eta_{\rm K}$ for the comparison. We observe that the distributions in Fig. 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.

3 Helical turbulent flows: suspensions of helicoidal particles

The helical flows described in Sec. 2 break parity invariance (chiral symmetry): in configurations dominated by positive helicity, as the flow in Fig. 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, 1871; 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 Kelvin (Kelvin, 1871) is illustrated in Fig. 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\degree$ (anti-chiral helicoid) or $-45\degree$ (co-chiral helicoid) with the great circle traversed clockwise, see Fig. 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 Fig. 4 where we show that isotropic helicoids move to different flow regions depending on their helicoidality also in our DNS of the NSE (1). In this section we use DNS, the stochastic model and theoretical approaches to analyse the motion of isotropic helicoids in helical turbulence.

3.1 Isotropic helicoids

The dynamics of an isotropic helicoid with position $\boldsymbol{x}$, velocity $\boldsymbol{v}$ and angular velocity $\boldsymbol{\omega}$ suspended in a fluid with velocity $\boldsymbol{u}$ and vorticity $2\boldsymbol{\Omega}=\boldsymbol{\nabla}\times\boldsymbol{u}$ is governed by the following equations (Kelvin, 1871; Happel & Brenner, 2012; Gustavsson & Biferale, 2016)

[TABLE]

Here dots denote time derivatives and ${\boldsymbol{u}}$ and ${\boldsymbol{\Omega}}$ 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, $\tau_{\rm p}$ is a relaxation time quantifying particle inertia. In the limit of $\tau_{\rm p}\to 0$ the particle approaches the dynamics of a tracer, $\boldsymbol{v}=\boldsymbol{u}$ and $\boldsymbol{\omega}=\boldsymbol{\Omega}$. 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 $\tau_{\rm 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 $\boldsymbol{\Omega}$, the particle accelerates either along the vorticity component, or opposite to it, see Fig. 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 $\eta$) for the point-particle approximation to be valid. The governing equations (14) 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 $\boldsymbol{\omega}$ 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 $\boldsymbol{\Omega}$, but is also affected by the strain rate of the flow (Jeffery, 1922).

Rescaling to dimensionless units $t^{\prime}=t/\tau_{\eta}$, $\boldsymbol{x}^{\prime}=\boldsymbol{x}/\eta$, ${{\boldsymbol{u}}}^{\prime}={\boldsymbol{u}}\tau_{\eta}/\eta$, ${{\boldsymbol{v}}}^{\prime}={\boldsymbol{v}}\tau_{\eta}/\eta$, ${{\boldsymbol{\omega}}}^{\prime}={\boldsymbol{\omega}}\tau_{\eta}$, and ${{\boldsymbol{\Omega}}}^{\prime}={\boldsymbol{\Omega}}\tau_{\eta}$ and dropping the primes in what follows, we can write the equations of motion for each pair of components $v_{i}$ and $\omega_{i}$ in dimensionless form:

[TABLE]

Here we have introduced the dimensionless size $a={\tilde{a}}/\eta$ and the Stokes number $\operatorname{St}=\tau_{\rm p}/\tau_{\eta}$. Interpreting the two-tensor $\mathsfbi{D}$ as a matrix, it has two eigenvalues $d_{\pm}$ and corresponding eigenvectors $\boldsymbol{\xi}_{\pm}$ given by

[TABLE]

These equations are well defined for all parameter values, but in the limit of isotropic particles, $C_{0}\to 0$, there is a complication. Taking the limit $C_{0}\to 0$ in Eqs. (16) and (17) 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 $\xi_{+}$ for $S<3/10$ and $\xi_{-}$ 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 $\operatorname{St}\sim 1$, see Table 3. This has been observed in simulations of inertial particles, where the most interesting dynamics occurs around values of $\operatorname{St}$ of order unity, see for example Fessler et al. (1994); Bec et al. (2007); Falkovich & Pumir (2007). The translational dynamics of isotropic helicoids on the other hand has two characteristic inertial scales $\operatorname{St}_{-}$ and $\operatorname{St}_{+}$ that depend on the helicoid parameters $C_{0}$ and $S$. These scales may be well separated in the meaning that $\operatorname{St}_{+}/\operatorname{St}_{-}$ can take arbitrarily large values. We therefore expect that isotropic helicoids may show significantly different behaviour depending on whether the Stokes number $\operatorname{St}$ is of the order of $\operatorname{St}_{-}$ or $\operatorname{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).

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 $\operatorname{St}$ ($\operatorname{St}\ll\operatorname{St}_{-}$)

[TABLE]

where $\mathsfbi{A}$ and $\mathsfbi{V}$ are matrices with elements $A_{ij}=\partial_{j}u_{i}$ and $V_{ij}=\partial_{j}\Omega_{i}$. Depending on the sign of $\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}$ trajectories of close-by particles may either converge ($\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}<0$) or diverge ($\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}>0$). It is expected that particles cluster in regions where $\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}<0$, i.e. where $27S\mbox{Tr}[\mathsfbi{A}^{2}]>9aC_{0}\mbox{Tr}[\mathsfbi{A}\mathsfbi{V}]/5$. For heavy spherical particles $C_{0}$ is zero and particles cluster in straining regions of the flow where $\mbox{Tr}[\mathsfbi{A}^{2}]>0$ (Maxey, 1987). For helicoids the structures in which particles with small values of $\operatorname{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\mbox{Tr}[\mathsfbi{A}\mathsfbi{V}]$, on the right-hand side of Eq. (18). As observed by Gustavsson & Biferale (2016), for a flow region with strong helical coherence, $\mathsfbi{V}\sim c\mathsfbi{A}$, particles cluster where $(27S-9acC_{0}/5)\mbox{Tr}[\mathsfbi{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 $\mbox{Tr}[\mathsfbi{A}^{2}]<0$, similar to light spherical particles. In order to quantify the preferential sampling of helical flow structures, we simulate the dynamics (15) for a number of parameters summarized in Table 4 using the flows described in Sec. 2.3. For each set of parameters, once the fluid reaches its statistically stationary state it is seeded with $2.4\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.

Fig. 4 illustrates that isotropic helicoids of opposing chirality preferentially sample different flow regions in DNS of a helical turbulent flow. Fig. 5a,b shows the fraction of particles in rotational regions of the flow. The data are plotted against $\operatorname{St}/\operatorname{St}_{-}$ in Fig. 5a and $\operatorname{St}/\operatorname{St}_{+}$ in Fig. 5b, i.e. against the inverse of the two eigenvalues (16) of the dynamics (when $C_{0}=0$ the data are only plotted against $\operatorname{St}/\operatorname{St}_{-}$ because for this case the preferential sampling cannot depend on $\operatorname{St}_{+}$). In rotational regions of a flow, the fluid gradient matrix $\mathsfbi{A}$ has complex eigenvalues, or equivalently, the sign of the discriminant

[TABLE]

is positive (Chong et al., 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 $\operatorname{St}$ they behave similar to light particles that oversample rotational regions where $\Delta>0$, while for larger values of $\operatorname{St}$ they instead behave as heavy inertial particles that oversample strain regions where $\Delta<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.

Fig. 5$\boldsymbol{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}{\boldsymbol{\Omega}}\rangle$, as functions of $\operatorname{St}/\operatorname{St}_{-}$ (Fig. 5$\boldsymbol{c}$) and $\operatorname{St}/\operatorname{St}_{+}$ (Fig. 5$\boldsymbol{d}$). Comparing Fig. 5a 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 $\operatorname{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 $\operatorname{St}$ helicoids with positive values of $C_{0}$ on average sample lower degree of helicity.

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 Fig. 4 which shows that the helicoids of opposite chirality distribute in different flow regions. It is also visible in Fig. 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$:

[TABLE]

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

3.4 Limiting cases

We now explain the main features of the data shown in Figs. 5 and 6 by analysing three limiting cases of the system parameters. The first two limiting cases, $\operatorname{St}\ll\operatorname{St}_{+}$ and $\operatorname{St}\gg\operatorname{St}_{-}$, apply to both DNS and the stochastic model, while the third limit is obtained for small values of $\operatorname{Ku}$ and therefore only applies to the stochastic model.

We use Eqs. (16) and (17) to write the dynamics (15) in its diagonal basis as

[TABLE]

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

[TABLE]

where the columns of the $2\times 2$ matrix $\mathsfbi{X}$ consist of the eigenvectors $\boldsymbol{\xi}_{-}$ and $\boldsymbol{\xi}_{+}$.

The dimensionless parameter groups in Eqs. (20) and (21) can be viewed as ratios of the time scales, $\tau_{\pm}/\tau_{\eta}\equiv\operatorname{St}/\operatorname{St}_{\pm}$, where $\tau_{\eta}$ is the Kolmogorov time of the flow and $\tau_{\pm}$ are two particle time scales ($\tau_{+}\leq\tau_{-}$) of the isotropic helicoid ($\tau_{+}=\tau_{-}=\tau_{\rm p}$ for a spherical particle). Eqs. (20) and (21) are only implicitly coupled through the trajectory dependence in $\boldsymbol{u}_{-}$ and $\boldsymbol{u}_{+}$. We therefore expect that the two limiting cases $\operatorname{St}\ll\operatorname{St}_{+}$ and $\operatorname{St}\gg\operatorname{St}_{-}$ can be taken, to a lowest-order approximation, in one of Eqs. (20) and (21) independent from the second equation. In summary, when $\operatorname{St}\ll\operatorname{St}_{+}$ ($\tau_{+}\ll\tau_{\eta}$), Eq. (21) becomes overdamped and the remaining equation (20) gives rise to strong preferential sampling when $\operatorname{St}/\operatorname{St}_{-}\sim 1$ ($\tau_{-}\sim\tau_{\eta}$) in analogue to the case $\operatorname{St}\sim 1$ ($\tau_{\rm p}\sim\tau_{\eta}$) for inertial spherical particles. In the second limit $\operatorname{St}\gg\operatorname{St}_{-}$ ($\tau_{-}\gg\tau_{\eta}$), Eq. (20) becomes underdamped and $\zeta_{-}$ can to lowest order be approximated by its mean value. The remaining equation (21) gives rise to strong preferential sampling when $\operatorname{St}/\operatorname{St}_{+}\sim 1$ ($\tau_{+}\sim\tau_{\eta}$). Below we discuss the two limiting cases in more detail.

3.4.1 Case $\operatorname{St}\ll\operatorname{St}_{+}$

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

[TABLE]

Using Eq.(23) and its time derivative to replace ${\boldsymbol{\omega}}$ and $\dot{\boldsymbol{\omega}}$ in Eq. (20), and reverting to the original coordinates, we obtain:

[TABLE]

Thus, a single equation determines the velocity of the isotropic helicoid in the limit $\operatorname{St}\ll\operatorname{St}_{+}$ and the angular velocity is given by Eq. (23). It can be noted that for the case of $C_{0}=0$, the constraint (23) becomes singular because velocity and angular velocity are uncoupled. However, Eq. (24) still shows the same results as those obtained by letting $C_{0}=0$ and $\operatorname{St}\ll\operatorname{St}_{+}$ in the original equation (15). When $C_{0}=0$ and $S>3/10$, Eq. (24) simplifies to $\dot{\boldsymbol{v}}=(\boldsymbol{u}-\boldsymbol{v})/\operatorname{St}$ while the rotational dynamics is overdamped, $\boldsymbol{\omega}-\boldsymbol{\Omega}=0$. When $C_{0}=0$ and $S<3/10$ on the other hand, the condition $\operatorname{St}\ll\operatorname{St}_{+}=1$ implies that the translational dynamics is overdamped and Eq. (24) simplifies to $\dot{\boldsymbol{v}}=(\boldsymbol{u}-\boldsymbol{v})/\operatorname{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 $\operatorname{St}\ll 1$.

The dynamics (24) can be further simplified using one or both of the following two assumptions. First, for small enough values of $\operatorname{St}/\operatorname{St}_{-}$, we approximate $\dot{\boldsymbol{u}}\sim{\rm D}_{t}{\boldsymbol{u}}$ and $\dot{\boldsymbol{\Omega}}\sim{\rm D}_{t}{\boldsymbol{\Omega}}$, where ${\rm D}_{t}\equiv\partial_{t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})$ are advective derivatives. Second, in a helical flow $\boldsymbol{\Omega}$ and $\boldsymbol{u}$ tend to be aligned and we may approximate ${\boldsymbol{\Omega}}\sim c{\boldsymbol{u}}$ with some proportionality constant $c$. Using these approximations, Eq.(24) simplifies to:

[TABLE]

If we define an effective Stokes number $\operatorname{St}_{\rm eff}$ and an effective density parameter $\beta_{\rm eff}$:

[TABLE]

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

[TABLE]

For spherical particles the density is characterized by $\beta$: $0\leq\beta<1$ corresponds to heavy particles more dense than the fluid, $\beta=1$ corresponds to neutrally buoyant particles and $1<\beta\leq 3$ corresponds to particles lighter than the fluid. We remark that $\beta_{\rm eff}$ in Eq. (27) is not constrained to the interval $0\leq\beta\leq 3$ as the case of spherical particles, $\beta_{\rm eff}$ may also take negative values as well as values larger than $3$.

*Comparison to DNS

*Although being a crude first-order approximation, Eq. (25) allows us to use what is known from spherical particles to explain the main features of the behaviour of helicoids in DNS for $\operatorname{St}\ll\operatorname{St}_{+}$ in Fig. 5a. The values of the effective density parameter $\beta_{\rm eff}$ in Eq. (27) are quoted in Table 4 for our parameters. Only the case $C_{0}=1.6$ has $\beta_{\rm 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 $\operatorname{St}_{\rm eff}=\operatorname{St}/\operatorname{St}_{-}$ is of order unity, which is consistent with the data in Fig. 5a. The cases of helicoids with $C_{0}=-1.6$ or $C_{0}=-5$ can be viewed as heavy particles because $\beta_{\rm eff}$ is close to zero. In these cases the helicoids preferentially sample strain regions of the fluid with strongest effect around $\operatorname{St}_{\rm eff}\sim O(1)$. Finally, the case $C_{0}=5$ has $\beta_{\rm 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 $\operatorname{St}_{\rm 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 Fig. 5a.

The approximation (25) allows us to also explain the observed clustering in Fig. 6a. The helicoids with $C_{0}=1.6$ show stronger clustering around $\operatorname{St}/\operatorname{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 $\operatorname{St}/\operatorname{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 $\operatorname{St}/\operatorname{St}_{-}$.

*Comparison to stochastic model

*Below we study in detail the validity of the approximations leading to (25). Since DNS is slow, it is hard to reach the steady state for helicoids with $\operatorname{St}\ll\operatorname{St}_{-}$ and to get good statistics in this limit. We therefore use the stochastic model to study the approximations.

Fig. 7a,b shows simulation results for isotropic helicoids with $\operatorname{St}\ll\operatorname{St}_{+}$ in the stochastic model described in Sec. 2.4 (solid lines). We choose parameters that are expected to correspond well with the DNS parameters in Fig. 5. We use a large Kubo number, $\operatorname{Ku}=10$, corresponding to the persistent flow limit where the dynamics most resembles the small scales in turbulence. As described in Sec. 2.4 we fix $H_{0}=0.85$ to match the distribution of flow helicity to that of the DNS and we take $\eta_{0}/\eta=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, $\tau_{\eta}\equiv\langle\mbox{Tr}(\mathsfbi{A}\mathsfbi{A}^{\rm T})\rangle^{-1/2}_{\rm flow}=\eta_{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 Fig. 7a and b to Fig. 5a and c we obtain a qualitative agreement with the DNS also for the isotropic helicoids. The general trends as functions of $\operatorname{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 $\Delta>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.

Fig. 7a,b also compares stochastic model simulations of the approximation (24) to the full simulation data of Eq. (15). We observe a quantitative agreement of the approximation in the expected limit $\operatorname{St}\ll\operatorname{St}_{+}$. Fig. 7c,d shows numerical simulations of (24) with $\dot{\boldsymbol{u}}$ and $\dot{\boldsymbol{\Omega}}$ replaced by ${\rm D}_{t}{\boldsymbol{u}}$ and ${\rm D}_{t}{\boldsymbol{\Omega}}$, which approach the results of Eq. (24) when $\operatorname{St}/\operatorname{St}_{-}\ll 1$ as expected. Finally, Fig. 7c,d also shows the approximation (25) using ${\boldsymbol{\Omega}}=c{\boldsymbol{u}}$ with $c=(\langle\boldsymbol{\Omega}^{2}\rangle/\langle\boldsymbol{u}^{2}\rangle)^{1/2}=\sqrt{5}/20$ for the stochastic model (rescaled using $\eta_{0}=10\eta$). We observe that the prediction using spherical particles reproduces the average sampling of helicity well for $\operatorname{St}\ll\operatorname{St}_{-}$ (Fig. 7d), while it does not work as well for the probability of being in vortex regions (Fig. 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}\pm 1.6$: for $C_{0}=1.6$ the helicoids do not oversample vortex regions to the degree that is predicted by Eq. (25) 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 Eq. (25).

Fig. 8 shows the positions of isotropic helicoids in a strongly helical flow in the stochastic model. Similar to the DNS in Fig. 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.

Fig. 9a shows stochastic model results for the correlation dimension $D_{2}$ with parameters corresponding to the DNS in Fig. 6. We observe qualitative agreement between Fig. 9a and the DNS.

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_{\rm L}$ (Frederickson et al., 1983). Denoting by $\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{9}$ the nine Lyapunov exponents of the equation system (14) together with $\dot{\boldsymbol{x}}=\boldsymbol{v}$, the Lyapunov dimension is

[TABLE]

where $K$ is the largest integer such that $\sum_{i=1}^{K}\lambda_{i}>0$. Fig. 9b shows numerical evaluation of the Lyapunov dimension for the stochastic model. For $D_{\rm L}<3$ it shows similar trends as the correlation dimension $D_{2}$. The values of $D_{\rm L}$ are somewhat higher than $D_{2}$, consistent with the fractal attractor in turbulence being a multifractal with $D_{\rm L}\geq D_{2}$ (Bec, 2005). Finally, we remark that the Lyapunov dimension defined in Eq. (29) 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 $\operatorname{St}$, the Lyapunov dimension is expected to approach the dimensionality $D$ of phase space. This is consistent with the data in Fig. 9b ($D=9$ for helicoids and $D=6$ for spherical particles).

Evaluation of the Lyapunov exponents also allows us to validate the approximation (18) 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\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}\rangle=\lambda_{1}+\lambda_{2}+\lambda_{3}$. Using this relation, we verify in Fig. 9c the small $\operatorname{St}$ approximation given by Eq. (18). We remark that the limit used to obtain Eq. (18) is the same as in this section ($\operatorname{St}\ll\operatorname{St}_{+}$), but requires first-order corrections in $\operatorname{St}/\operatorname{St}_{+}$ to the condition in Eq. (23). Moreover, we need to consider $\operatorname{St}/\operatorname{St}_{-}\ll 1$. Consistently with this limit, we find that the average of the approximation (18) tends to approach $\langle\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}\rangle$ as $\operatorname{St}/\operatorname{St}_{-}$ is reduced and that the two expressions agree approximately for $\operatorname{St}/\operatorname{St}_{-}\sim 0.1$.

3.4.2 Case $\operatorname{St}\gg\operatorname{St}_{-}$

We now consider the second limiting case, $\operatorname{St}\gg\operatorname{St}_{-}$ with general values of $\operatorname{St}_{+}$. In this limit $\zeta_{-,i}$ in Eq. (20) responds slowly to changes in the flow compared to $\zeta_{+,i}$. Due to the symmetries of the underlying flow, we expect the averages of $u_{-,i}$ and consequently of $\zeta_{-,i}$ to vanish. We therefore replace $\zeta_{-,i}$ by its vanishing average, $\zeta_{-,i}=0$, which gives the following constraint on $\boldsymbol{\omega}$

[TABLE]

Inserting this constraint and its time derivative into Eq. (21), gives the following equation for the velocity

[TABLE]

As for the first limiting case we can consider a helical flow with $\boldsymbol{\Omega}\sim c\boldsymbol{u}$ to obtain

[TABLE]

where $\beta_{\rm eff}$ is the parameter in Eq. (27) occurring in the limit of $\operatorname{St}\ll\operatorname{St}_{-}$. Thus, in the limit $\operatorname{St}\gg\operatorname{St}_{-}$ the equation of motion for helicoids in a helical flow is like a Stokes drag with effective Stokes number $\operatorname{St}/\operatorname{St}_{+}$ and with a rescaled amplitude of the fluid velocity.

*Comparison to DNS

*The limiting dynamics in Eq. (31) allows us to explain the DNS results in Fig. 5b. For helicoids with negative values of $C_{0}$, the coupling to the flow, $\beta_{\rm eff}\boldsymbol{u}$ in Eq. (31), is small, see Table 4. The particle motion is thus expected to be only weakly correlated to the underlying flow structures, which is consistent with the data: the particles with negative values of $C_{0}$ have approximately the same statistical properties as the flow (black dashed lines in Fig. 5b). The helicoids with positive values of $C_{0}$ on the other hand have $\beta_{\rm eff}\sim 1$ and are therefore expected to have preferential sampling similar to spherical particles with effective Stokes number $\operatorname{St}_{+}/\operatorname{St}$. This is what we observe for $C_{0}=1.6$, the shape and magnitude of the curve around $\operatorname{St}\sim\operatorname{St}_{+}$ is similar to that of spherical particles. However, the approximation (31) does not work as well for $C_{0}=5$: even though $\beta_{\rm 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 $\operatorname{St}\ll\operatorname{St}_{+}$, we conclude that the approximations (for example $\boldsymbol{\Omega}\sim c\boldsymbol{u}$) leading to Eq. (25) and Eq. (31) may not be so accurate for $C_{0}=5$ in DNS.

Using the approximation (31) to explain the fractal clustering observed in Fig. 6b, 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 $\operatorname{St}\sim\operatorname{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\to 0$. Indeed, as seen in the inset of Fig. 6b, 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 (31) 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 (25) mainly depend on the instantaneous flow and are therefore approximately the same for the two particles.

*Comparison to stochastic model

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

We compare stochastic model simulations of the approximation (30) [solid lines] to simulations of Eq. (14) in Fig. 10. When $\operatorname{St}\gg\operatorname{St}_{-}$ we observe a quantitative agreement for $C_{0}=\pm 1.6$ and qualitative agreement when $C_{0}=\pm 5$. Fig. 10 also shows that the approximation (31) based on $\boldsymbol{\Omega}=c\boldsymbol{u}$ (dash-dotted lines) does not work equally well. Eq. (31) reproduces that co-chiral helicoids have larger preferential sampling of straining regions around $\operatorname{St}\sim\operatorname{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 $\operatorname{St}\gg\operatorname{St}_{+}$ the preferential sampling for all parameter cases in Table 4 converges to that of the flow (not shown).

As seen in Fig. 9a the correlation dimension $D_{2}$ increases monotonously towards the spatial dimension $d=3$ after the peak at $\operatorname{St}\sim\operatorname{St}_{-}$. As a consequence, $D_{2}$ does not show any clustering around $\operatorname{St}\sim\operatorname{St}_{+}$. The deviations from the DNS data in Fig. 6b can be explained by the fact that the observed data are better resolved in the stochastic model: the correlation dimension around $\operatorname{St}\sim\operatorname{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 $\operatorname{Ku}$

In the stochastic model the properties of the flow can be modified by changing the value of the Kubo number (11). In general, this allows exploration of the robustness of results with respect to the nature of the flow. In the limit of small $\operatorname{Ku}$ we can solve the dynamics analytically in terms of the full set of model parameters $\operatorname{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 $\operatorname{Ku}$:

[TABLE]

where $\operatorname{{\overline{a}}}=\tilde{a}/\eta_{0}$ (velocity and position are made dimensionless in terms of $u_{0}$ and $\eta_{0}$). Fig. 11a shows the analytical solution for the mean helicity together with data for $\operatorname{Ku}=10$. In order to compensate for the different magnitudes of relevant time scales in flows with $\operatorname{Ku}=10$ and flows with small values of $\operatorname{Ku}$, the parameters $\operatorname{Ku}$ and $\operatorname{St}$ which depend on the correlation time of the flow have been rescaled in Eq. (32). 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}$). Fig. 11a shows that the small Kubo results have the same trends as those of the DNS in Fig. 5b and the stochastic model in Fig. 7b. This shows that the trends shown in Fig. 5 are robust, they do not depend on the particular nature of the underlying flow. Using this observation, we can use the theoretical solution of the stochastic model to get an estimate of the parameter dependence of preferential sampling of helicity for general values of the five model parameters. Two examples of this dependence are illustrated in Fig. 11b and c. Fig. 11b 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 Fig. 5b and d: co-chiral helicoids oversample helicity for small Stokes numbers, while anti-chiral helicoids oversample helicity for large Stokes numbers. Helicoids with small $|C_{0}|$ behave similar to neutral particles and undersample helicity. Similar trends are observed in a neutral flow (see Fig. 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 Fig. 11c: 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 Fig. 11c.

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 $\operatorname{Ku}$, similar to the expansions for the Lyapunov exponents of spherical particles (Gustavsson & Mehlig, 2011, 2016). The theory relies on the additional constraint $\operatorname{St}\ll\operatorname{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 (9) is qualitatively well reproduced by our DNS results as shown in Fig. 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 (14) 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 Figs. 5 and 6 have sensitive and complicated parameter dependence. Nevertheless, it is remarkable that the crude approximations (25) and (31) 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 $\operatorname{St}_{\rm eff}$ (26), $\beta_{\rm eff}$ (27) and $\operatorname{St}/\operatorname{St}_{+}$. The only term that is odd in $c$ or $C_{0}$ in these effective parameters, as well as in Eq. (18), 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\eta$), 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 $\beta_{\rm eff}$, as exemplified by comparing $|C_{0}|=5$ and $|C_{0}|=1.6$ in Fig. 5. By optimizing $\beta_{\rm eff}$ it is possible to find helicoids with smaller values of $a$ that behave as light particles, i.e. have $\beta_{\rm eff}>1$. We conclude by remarking that the construction in Fig. 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 research leading to these results has received funding from the European Union’s Seventh Framework Programme (FP7/2007-2013) under grant agreement No. 339032. K.G. acknowledges funding from the Knut and Alice Wallenberg Foundation, Dnar. KAW 2014.0048. The numerical computations used resources provided by C3SE and SNIC.

Appendix A Stochastic model for helical turbulence

To construct the incompressible, homogeneous and isotropic stochastic velocity field, $\boldsymbol{u}=\boldsymbol{\nabla}\times\boldsymbol{A}$, used in the article, we generate the components of the vector potential $\boldsymbol{A}$ as a Fourier sum

[TABLE]

Here $L=10$ is the system size (we use $L=10\eta_{0}$ in our simulations), the wave vector $\boldsymbol{k}$ is summed over the components $k_{j}=2\pi n_{j}/L$ with $n_{j}=-20,-19,\dots,+20$ and $j=1,2,3$. For each $\boldsymbol{k}$, the vector of Fourier coefficients, $\boldsymbol{a}_{\boldsymbol{\footnotesize k}}(t)$, has been expanded in terms of the eigenmodes $\boldsymbol{h}^{\pm}_{\boldsymbol{\footnotesize k}}$ of the curl operator, weighted by a factor $\mu$ to give a bias to positive helical modes if $\mu>1$, and to negative helical modes if $0\leq\mu<1$. The coefficients $a_{i,\boldsymbol{\footnotesize k}}(t)$ are complex random Gaussian numbers fulfilling the condition $a^{*}_{i,\boldsymbol{\footnotesize k}}=a_{i,-\boldsymbol{\footnotesize k}}$ and having the statistics

[TABLE]

The exponential time correlation in (34) is generated from an underlying Ornstein-Uhlenbeck processes

[TABLE]

where $\delta t$ is the time step of the simulation and $b_{i,\boldsymbol{\footnotesize k}}(t)$ are independent random Gaussian numbers that are white noise in time with statistics

[TABLE]

The Gaussian cutoff for large $k$ in Eq. (33) ensures a Gaussian spatial correlation function with correlation length $\eta_{0}$. When $L\gg\eta_{0}$, Eq. (33) implies the correlation function

[TABLE]

From this correlation function the statistics of $\boldsymbol{u}$ and its spatial derivatives follows. To obtain the distribution $P_{0}(H)$ of helicity $H=2\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\Omega}$ for the stochastic flow in Eq. (33), we start from the joint distribution of $\boldsymbol{u}$ and $\boldsymbol{\Omega}$

[TABLE]

where $\boldsymbol{X}=(u_{1},u_{2},u_{3},\Omega_{1},\Omega_{2},\Omega_{3}){{}^{\rm T}}$ and $\mathsfbi{C}$ is the corresponding covariance matrix (velocity is made dimensionless in terms of $u_{0}$ and position in terms of $\eta_{0}$)

[TABLE]

obtained from Eq. (37) with $H_{0}\equiv\sqrt{2/\pi}8(\mu^{2}-1)/(3(\mu^{2}+1))$. After a change of coordinates $\Omega_{z}=(H/2-\Omega_{x}u_{x}-\Omega_{y}u_{y})/u_{z}$ and integration over $\Omega_{x}$ and $\Omega_{y}$ the remaining joint distribution of $H$, $u_{x}$, $u_{y}$ and $u_{z}$ depends only on $H$ and the combination $\sqrt{u_{x}^{2}+u_{y}^{2}+u_{z}^{2}}$. Changing to spherical coordinates in $\boldsymbol{u}$-space and integrating them away gives the final distribution of helicity, Eq. (12):

[TABLE]

where $K_{\nu}(x)$ is the modified Bessel function of the second kind. The average helicity of the flow is determined from the helicity bias $\mu$ as follows:

[TABLE]

Appendix B Expansion around deterministic trajectories

Here we outline the series expansion used to calculate the mean helicity (32) of isotropic helicoids for small values of $\operatorname{Ku}$ in the statistical model. The expansion follows the method introduced in Gustavsson & Mehlig (2011) and reviewed in Gustavsson & Mehlig (2016). We want to expand the dynamics in the dimensionless equations of motion (Eq. (15) together with $\dot{\boldsymbol{r}}=\operatorname{Ku}\boldsymbol{v}$) around the deterministic solution $\boldsymbol{r}^{(\rm d)}$ obtained without flow, i.e. when $\boldsymbol{u}=\boldsymbol{\Omega}=\boldsymbol{0}$. Since we in this work only consider homogeneous steady-state statistics, we can put all initial conditions to zero for simplicity. We therefore expand around the simple deterministic solution $\boldsymbol{r}^{(\rm d)}=\boldsymbol{0}$. An implicit solution to the dynamics in Eq. (15) together with $\dot{\boldsymbol{r}}=\operatorname{Ku}\boldsymbol{v}$ can be found by first solving the diagonal equations for $\boldsymbol{\zeta}$ in Eqs. (20) and (21), then transforming back to $\boldsymbol{v}$ and finally integrating to obtain $\boldsymbol{r}$. We find that the following is an exact implicit solution to the dynamics:

[TABLE]

where

[TABLE]

A series expansion of the flow velocity (and spatial derivatives thereof) around the deterministic trajectory $\boldsymbol{r}^{(\rm d)}=0$ gives

[TABLE]

Eq. (50) is an expansion of the flow velocity in terms of the displacement from $\boldsymbol{r}=0$. Recursively substituting $\boldsymbol{u}(\boldsymbol{r}_{t},t)$ (and derivatives thereof) from Eq. (50), and $\boldsymbol{r}_{t}$ from Eq. (48) into Eq. (50), we obtain an increasingly refined approximation of the flow evaluated along the true trajectory $\boldsymbol{r}_{t}$. Since $\boldsymbol{r}_{t}$ is of order $\operatorname{Ku}$ in Eq. (48), we can use $\operatorname{Ku}$ to keep track of the order of $\boldsymbol{r}_{t}$ in the expansion. Truncating the recursive expansion of Eq. (50) at some order in $\operatorname{Ku}$, one obtains the approximate expression for $\boldsymbol{u}$ along a trajectory $\boldsymbol{r}_{t}$ to this order in $\operatorname{Ku}$. To evaluate an approximation for the helicity along a particle trajectory, we do a similar expansion for $\boldsymbol{\Omega}(\boldsymbol{r}_{t},t)$ to form $H=2\boldsymbol{u}(\boldsymbol{r}_{t},t)\boldsymbol{\cdot}\boldsymbol{\Omega}(\boldsymbol{r}_{t},t)$. Finally, we evaluate the steady-state average $\langle H\rangle=2\langle\boldsymbol{u}(\boldsymbol{r}_{t},t)\boldsymbol{\cdot}\boldsymbol{\Omega}(\boldsymbol{r}_{t},t)\rangle$, where the average is taken over an ensemble of trajectories and can be explicitly evaluated for the stochastic model in terms of the known Eulerian correlation function in Eq. (37), see Gustavsson & Mehlig (2011, 2016) for more details. As a result, we obtain the expression in Eq. (32).