Gyrotactic phytoplankton in laminar and turbulent flows: a dynamical systems approach

TL;DR

This paper reviews how gyrotactic phytoplankton, influenced by fluid flows, exhibit complex distribution patterns like trapping and clustering, explained through dynamical systems theory.

Contribution

It provides a comprehensive overview of recent theoretical, numerical, and experimental findings on gyrotactic phytoplankton in various flow regimes, highlighting new insights into their distribution behaviors.

Findings

Gyrotactic trapping occurs in nonlinear shear flows.

Fractal clustering of phytoplankton is observed in turbulent flows.

Dynamical systems theory effectively explains these phenomena.

Abstract

Gyrotactic algae are bottom heavy, motile cells whose swimming direction is determined by a balance between a buoyancy torque directing them upwards and fluid velocity gradients. Gyrotaxis has, in recent years, become a paradigmatic model for phytoplankton motility in flows. The essential attractiveness of this peculiar form of motility is the availability of a mechanistic description which, despite its simplicity, revealed predictive, rich in phenomenology, easily complemented to include the effects of shape, feed-back on the fluid and stochasticity (e.g. in cell orientation). In this review we consider recent theoretical, numerical and experimental results to discuss how, depending on flow properties, gyrotaxis can produce inhomogeneous phytoplankton distributions on a wide range of scales, from millimeters to kilometers, in both laminar and turbulent flows. In particular, we focus on the phenomenon of gyrotactic trapping in nonlinear shear flows and in fractal clustering in turbulent flows. We shall demonstrate the usefulness of ideas and tools borrowed from dynamical systems theory in explaining and interpreting these phenomena.

Equations

$$\dot{\mathbf{p}}=\frac{1}{2B}\left[\hat{\mathbf{z}}-(\hat{\mathbf{z}}\cdot \mathbf{p})\mathbf{p}\right]+\frac{1}{2}\mathbf{\omega }\times \mathbf{p}$$

$$\dot{\mathbf{x}}=\mathbf{u}+v_s\mathbf{p}.$$

$$\begin{array}[]{l}\mathrm{p}^{eq}_{x}=4UBx/L^{2}\\ \mathrm{p}^{eq}_{z}=\left(1-16B^{2}U^{2}x^{2}/L^{4}\right)^{1/2}\end{array}$$

$$\dot{\mathbf{p}}=-\frac{1}{v_0}\left[\mathbf{A}-(\mathbf{A}\cdot \mathbf{p})\mathbf{p}\right]+\frac{1}{2}\mathbf{\omega }\times \mathbf{p}+{\mathbf{\Gamma }}_r$$

$${\mathbf{p}}^{eq}=\left(\frac{-\gamma \mathbf{r}}{\sqrt{1+(\gamma r{)}^2}},\frac{1}{\sqrt{1+(\gamma r{)}^2}}\right)$$

$$n_s(r)=\mathcal{N}\exp \left(-\frac{\gamma r^2}{2v_{s}B{F}_3^2(\lambda )}\right)$$

$${\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u}=-\nabla p+\nu \Delta \mathbf{u}+F\cos (z/L)\hat{\mathbf{x}}$$

$$\dot{X}$$

$$\dot{Y}$$

$$\dot{Z}$$

$${\dot{\mathrm{p}}}_x$$

$${\dot{\mathrm{p}}}_y$$

$${\dot{\mathrm{p}}}_z$$

$$\mathcal{C}(\mathbf{p},Z)$$

$$\mathcal{H}(\mathbf{p},Z)$$

$$\dot{\theta }$$

$$\dot{Z}$$

$$\cos \theta ={\mathrm{p}}_x=\frac{\Psi (2\Phi \Psi \cos z-\sin z)}{1+4{\Phi }^2{\Psi }^2}.$$

$$\mathcal{P}(T)=\frac{L_B}{(4\pi D_z{T}^3{)}^{1/2}}{e}^{-\frac{(V_{d}T-L_B/2{)}^2}{4D_{z}T}}.$$

$$\dot{\Gamma }=\sum _{i=1}^3⟨\frac{\partial \dot{x_i}}{\partial x_i}+\frac{\partial \dot{p_i}}{\partial p_i}⟩=-\frac{2}{\Psi }⟨p_z⟩,$$

$${\mathbf{p}}_{eq}\sim \left(\Psi {\omega }_y,-\Psi {\omega }_x,1\right).$$

$$\nabla \cdot \mathbf{v}=\Phi \nabla \cdot {\mathbf{p}}_{eq}=-\Psi \Phi {\nabla }^2{u}_z,$$

$$⟨u_{z}n⟩=\Psi \Phi {\int }_{-\infty }^{t}ds⟨u_z(\mathbf{x}(t),t){\nabla }^2{u}_z(\mathbf{x}(s),s)⟩.$$

$$\nabla \cdot \mathbf{a}=\sum _{ij}\left({\hat{S}}_{ij}^2-{\hat{\Omega }}_{ij}^2\right),$$

$$R^{\ast }\simeq \eta \Phi \Delta \Psi .$$

Full text

Microorganisms in a flow and Fractal clustering and Motility

I Introduction

Biological and geophysical fluids host a sea of microorganisms many of which are motile. An often overlooked aspect of the life of such microorganisms is that the fluids where they are suspended are not still but flowing. For instance, in fresh water and marine environments microorganisms, such as unicellular algae or bacteria, are exposed to turbulent motion Guasto et al. (2012) and turbulence is believed to have been one of the the main factors in shaping the huge variety of form and strategies of such aquatic microbes Margalef (1978). In (photo-)bio-reactors microorganisms are grown in continuously stirred tanks Rodolfi et al. (2009); Croze et al. (2013). Bacteria composing the microbioma of mammals undergo flowing fluids in several organs that they inhabit such as, e.g., the gut or renal tissues Duan et al. (2008); Kim et al. (2012). Microbes are also exposed to a variety of water flows in the soil Ebrahimi and Or (2014).

In the presence of flows, microorganisms are at the mercy of the velocity field which transports them possibly, for motile species, in addition to their intrinsic swimming, and modifies their orientation by action of velocity gradients. Flow can thus impact the motility of microorganism, their spatial distribution, interaction with surfaces, response to nutrients or other chemicals for motile species Rusconi and Stocker (2015); Chilukuri et al. (2014); Locsei and Pedley (2009); Stocker and Seymour (2012), and nutrient uptake also for non motile ones Karp-Boss et al. (1996). Moreover, flow influences the encounter rate of microorganisms Rothschild and Osborn (1988), affecting their reproduction and competition Kiørboe and MacKenzie (1995); Pigolotti et al. (2013). In summary, flows are key in shaping microbial responses and ecology Kiørboe (2008).

In this brief review we aim at emphasizing some interesting phenomena that can emerge due to the modification of the microorganisms’ swimming direction by velocity gradients, which affects both the individual motion of microorganisms and their spatial distribution in dilute suspensions. More specifically, we focus on the case study of gyrotactic phytoplankton.

Phytoplankton is responsible for about half of the photosynthetic activity on Earth. It is composed by thousands of species many of which are able of swim. Motility confers phytoplankton the ability to reach well-lit waters near the sea surface during daylight and migrate into deeper water, richer of nutrients, during the night Lampert (1989). For several species, the upward migration is guided by a stabilizing torque, induced, e.g., by bottom heaviness, which biases cell’s swimming direction upwards, and is opposed by hydrodynamic shear, which exerts a viscous torque tending to overturn the cell. When the swimming direction results from the competition between the cell’s stabilizing torque and the shear-induced viscous torque, we speak of gyrotaxis Kessler (1985); Pedley and Kessler (1987, 1992). In both laminar and turbulent flows gyrotaxis can promote heterogeneous spatial distribution of phytoplankton cells.

Unicellular algae of the genus Chlamydomonas aggregate in the center (walls) of downwelling (upwelling) vertical pipe flows Kessler (1985), a similar phenomenon can be induced by phototaxis in horizontal pipe flows Garcia et al. (2013); Martin et al. (2016). Gyrotaxis not only alters cells’ spatial distribution but also the cells’ dispersion properties Thorn and Bearon (2010); Bees and Croze (2010); Bearon et al. (2012), which can be important for photobioreactors Croze et al. (2013). When the fluid acceleration is not negligible with respect to the gravitational one, the stabilizing torque biases the motion in a position dependent direction given by superposition of fluid and gravitational acceleration. This causes, for instance, cell focusing toward the rotation axis, when this is directed along the vertical De Lillo et al. (2014); Cencini et al. (2016). Remarkably, also when the shear-induced viscous torque is not balanced by the stabilizing one, interesting phenomena can happen. In this condition, cells overturn due to the shear induced torque and start tumbling without directed motion. In inhomogeneous shear flows, this tumbling motion can trap cells in regions of high shear, a phenomenon, discovered in microfluidic experiments Durham et al. (2009), which can explain the formation of high phytoplankton concentrations in very thin layers as observed in coastal oceans Durham and Stocker (2012). This gyrotactic trapping has an interesting interpretation from a dynamical systems point of view Santamaria et al. (2014). In more complex laminar flows, such as Tayolor-Green vortices, the combination of gyrotactic motility and flow can give rise to small scale aggregation and complex trajectories Durham et al. (2011). Such effects become even more striking in turbulent flows of moderate intensity, where gyrotaxis can generate small-scale fractal clusters which are dynamically formed and dissolved, both when the fluid acceleration can be neglected and when it cannot Durham et al. (2013); De Lillo et al. (2014), see also Fouxon and Leshansky (2015); Gustavsson et al. (2016). The properties of such fractal clusters is also influenced by cell morphology Gustavsson et al. (2016); Zhan et al. (2014). The physics underlying the formation of such fractal clusters can be easily understood using dynamical systems concepts: cells swimming in turbulence can be described in terms of a chaotic, dissipative system, whose trajectories naturally evolve onto (multi - )fractal sets Ott (1993).

Here, we review some of the above phenomena emphasizing their interpretation within the framework of dynamical systems theory. This point of view is chosen on the basis of the background of the Authors. This might have biased the choice of some topics and it is inevitable that some relevant works on the subject have been not properly discussed. This choice also left out many interesting phenomena arising in dense suspensions of gyrotactic organisms such as bioconvection Pedley et al. (1988); Ghorai and Hill (2007) or complex rheological effects Rafaï et al. (2010). Before presenting the organization of the matter, we would like to mention some interesting phenomena induced by the interplay between motility and flow that arise in other kinds of microorganisms, such as bacteria, and that are connected with those observed in gyrotactic phytoplankton.

For instance, owing to their elongated shape, bacterial cells swimming in a pipe (Poiseuille) flow are affected by Jeffery orbits Jeffery (1922) which can induce upstream swimming in low shear regions and tumbling in high shear ones. This was predicted on the basis of mathematical models in Zöttl and Stark (2012, 2013) and observed in microfluidic experiments Rusconi et al. (2014), which revealed also the accumulation of cells in high shear regions due to trapping induced by Jeffery orbits. This trapping causes spatial inhomogeneity and can reduce the efficiency of chemotaxis Rusconi et al. (2014), as also confirmed by mathematical analysis Bearon and Hazel (2015). Motility also modifies the transport properties with respect to passively advected particles Khurana et al. (2011). More in general, flow effects on motility are expected to impact chemotaxis and other kinds of taxis Stocker and Seymour (2012). Although such effects are largely unexplored, we mention here the numerical study of chemotaxis in the presence of a turbulent flow Taylor and Stocker (2012). The details of cell morphology, such as the chirality of the flagella, can also induce further directional biases in the presence of velocity gradients Marcos et al. (2012). Flows can also impact the interaction of microorganisms with surfaces Chilukuri et al. (2014), e.g., inducing accumulation Berke et al. (2008) and upstream swimming Kaya and Koser (2012) close to the walls of pipe flows. Fractal clustering, similar to that found in gyrotactic cells, was also found in models of bacteria swimming in cellular and chaotic flows with or without taxis Torney and Neufeld (2007, 2008).

The material is organized as follows. In Sect. II we describe the equations which have become the standard model of gyrotactic motility and their modification when fluid acceleration is important. Moreover, we briefly discuss the experimental validation of that modification and the generic phenomenology one can derive from these models. Section III is devoted to the phenomenon of gyrotactic trapping in inhomogeneous shear flows. In particular, we will consider the case of the Kolmogorov flow and show how, borrowing ideas from conservative dynamical systems, gyrotactic trapping can be ascribed to the presence of effective barriers to transport. We will also briefly discuss how the presence of small-scale turbulence destabilizes the trapping. In Sect. IV, we shall focus on the effects of turbulence on gyrotactic motility, showing how fractal clustering emerges due to the chaotic dissipative character of the dynamics. Moreover, we will link this phenomenon to the well known clustering of inertial particles in turbulence, showing that there are many conceptual analogies. We will also briefly consider the case of a population in which cells have different characteristics, closer to what one may observe in the ocean, to discuss how clustering can be revealed in these cases. Finally, Sect. V is devoted to conclusions.

II Mathematical models

Gyrotaxis is observed in several species of unicellular algae such as, e.g., biflagellate spheroidal algae belonging to the genus Chlamydomonas and Dunaliella , and also in some monoflagellate species such as e.g. Heterosigma akashiwo Kessler (1985); Guasto et al. (2010); Durham et al. (2009, 2013). Having as a reference Chlamydomonas, we can consider the cell close to spherical with a diameter of about $10\mu{\rm m}$ and swimming speeds around $100\mu{\rm m}\,s^{-1}$ Guasto et al. (2010). Most phytoplankton cells have density very close to that of water and thus can be considered as neutrally buoyant. A characteristic of most cells displaying gyrotaxis is to have an inhomogeneous distribution of mass leading to a displacement of the cell center of gravity with respect to its center of symmetry. In particular, their center of mass is displaced opposite to their direction of swimming, so they are generally defined as bottom-heavy Kessler (1985).

The mathematical model for gyrotactic algae was introduced by Kessler Kessler (1985) (see also Pedley and Kessler (1990, 1992)) on the basis of the observation that bottom-heavy swimming micro-organisms focus in the center of a pipe when the fluid flows downwards. The swimming direction ${\bf p}$ results from the competition between gravity-buoyant torque, due to bottom-heaviness, and the shear-induced viscous torque and evolves according to

[TABLE]

where ${\bm{\omega}}={\bm{\nabla}}\times{\bf u}$ is the vorticity at the position of the cell, $B=\nu\alpha_{\perp}/(2hg)$ is a characteristic orientation time which depends on viscosity $\nu$, gravity ${\bf g}=-g\hat{\bf z}$, on the displacement $h$ of the cell center of mass relative to the geometrical center and on $\alpha_{\perp}$ the dimensionless resistance coefficient for rotation about an axis perpendicular to $\bf p$. For a sphere $\alpha_{\perp}=6$ and thus $B=3\nu/(hg)$. The third term on the rhs is the rotation due to local rate of strain $\hat{S}_{ij}=\frac{1}{2}(\partial_{j}u_{i}+\partial_{i}u_{j})$ and it is controlled by the shape factor $\alpha=(l^{2}-d^{2})/(l^{2}+d^{2})$ measuring the elongation $l$ with respect to the width $d$ of the cell Jeffery (1922). Prolate (respectively oblate) particles have $\alpha>0$ (respectively $\alpha<0$) and the particular cases of rods, disks and spheres are described by $\alpha=1,-1,0$. As a consequence, this term vanishes for spheres. The stochastic term ${\bf\Gamma}_{r}$ represents rotational diffusion of the swimming direction as a results of irregularities in the cell propulsion Hill and Häder (1997), indeed cells are too large to be subjected to Brownian rotation. Equation (1) has been written assuming the general case of ellipsoidal cells, although in this review we will mainly consider spherical cells, and we will thus assume $\alpha=0$ unless otherwise specified.

Owing to their small size and small density mismatch with the fluid, gyrotactic cells can be represented as point-like and neutrally buoyant particles transported by the fluid velocity ${\bf u}({\bf x},t)$ with a superimposed swimming velocity of intensity $v_{s}$ along the direction ${\bf p}$

[TABLE]

We remark that in writing the model defined by Eqs. (1)-(2) many details have been neglected including the unsteadiness of swimming due to flagella beating, cell-–cell interactions, the feedback of cell motion on the surrounding fluid, and (see below) effects due to fluid acceleration. The model can thus be appropriate only for dilute suspensions where interactions between cells and flow modifications are expected to be negligible. In spite of its simplicity, however, the gyrotactic model (1-2) is able to predict remarkable features observed in experiments, such as the focusing observed in pipe flows Kessler (1985).

To illustrate the phenomenon of focusing, we consider a two-dimensional, laminar flow in a vertical pipe, in which the velocity field has only the vertical component $w(x)=U(1-(x/L)^{2})$ and the coordinate $x\in[-L:L]$ varies in the wall-normal direction. Neglecting the stochastic term, with this choice it is easy to find the equilibrium, quasi stationary solution of Eq. (1):

[TABLE]

The equilibrium solution exists for $BU/L\leq 1/4$, i.e. for small stability number $\Psi\equiv BU/L$, which is the ratio of the second to the first term in the RHS of (1). When used in (2), the stationary solution (3) shows that for downwelling flow ($U<0$) the horizontal swimming direction is toward the symmetry axis (since $\mathrm{p}_{x}\propto-x$). This leads to the gyrotactic focusing, i.e. an accumulation of motile cells at the center of the pipe Kessler (1985). For upwelling flows ($U>0$) the sign reverses and accumulation happens at the pipe’s walls. When the stochastic term is included one realizes that cells focus in a region with a finite width controlled by the noise strength Pedley and Kessler (1990, 1992).

In analogy to the stability numbers, from (2) we also introduce a swimming number defined as $\Phi=v_{s}/U$, measuring the swimming speed with respect to the fluid velocity. More complex dynamics of gyrotactic swimmers in laminar flows, with or without stochastic behavior, have been investigated Thorn and Bearon (2010).

II.1 Effects of fluid acceleration

Recently, the gyrotactic model has been generalized by taking into account the effect of fluid acceleration in the orientation of the cell De Lillo et al. (2014). By introducing the total acceleration on the cell ${\bf A}={\bf g}-{\bf a}$, where ${\bf a}=d{\bf u}/dt$, Eq. (1) becomes

[TABLE]

where $v_{0}=Bg$ is the reorientation speed.

The general model (4) has been validated by laboratory experiments of a suspension of gyrotactic cells in a cylindrical vessel of radius $R$ which rotates with constant angular velocity $\Omega$ Cencini et al. (2016). When rotation is sufficiently fast, centripetal acceleration overcomes gravity and the cells are expected to swim towards the axis of the cylinder. Neglecting the stochastic term, the deterministic motion given by (1-2) in the velocity field produced by solid body rotation ${\bf u}=(-\Omega y,\Omega x,0)$ can be easily obtained in cylindrical coordinates ${\bf x}=({\bf r},z)$ under the hypothesis of local equilibrium in the swimming direction, i.e. by assuming $\dot{\bf p}=0$ locally. In this limit one obtains Cencini et al. (2016)

[TABLE]

where $\gamma=\Omega^{2}/g$. For $\gamma r\ll 1$, from (2) and (5) we have $\dot{r}=-\gamma v_{s}r$ which implies that the cell position relaxes exponentially towards the rotation axis $r=0$ as $r(t)=r(0)e^{-\gamma v_{s}t}$.

Figure 1 shows a sketch of the experimental setup together with two pictures of the asymptotic cell distributions for different angular velocities. The experiments show that, at late times, the asymptotic cell distribution is over a finite radius around the axis, and this has been interpreted as due to stochastic effects, modeled as rotational diffusion, which deviate the swimming direction with respect to the equilibrium one (5).

The stationary distribution in presence of rotational diffusivity can be obtained within the so-called Generalized Taylor Dispersion theory Frankel and Brenner (1989); Bearon et al. (2011, 2012). The basic idea is to reduce the Fokker-Planck equation for the probability ${\cal P}({\bf x},{\bf p},t)$ to an advection-diffusion equation for the population density $n({\bf x},t)=\int d{\bf p}{\cal P}$ which includes an effective drift and diffusivity tensor which can be analytically derived in the approximation of fast orientation of the swimming direction Bearon et al. (2011, 2012). The final result (see Cencini et al. (2016) for details) is, for $\gamma r\ll 1$, a Gaussian distribution of the radial population density

[TABLE]

where $F_{3}(\lambda)$ a dimensionless function of the parameter $\lambda=1/(2Bd_{r})$ ($d_{r}$ is the rotational diffusivity), which can be analytically expressed as series, and the coefficient $\mathcal{N}$ can be written in terms of the total number of cells $N_{s}=\int drn_{s}(r)$ as $\mathcal{N}=N_{s}\gamma/(2\pi Hv_{s}BF_{3}^{2}(\lambda))$ ($H$ is the vertical size of the vessel). Figure 2 shows the evolution of the experimental radial population density measured at different times, together with the asymptotic theoretical density Eq.(6) for the case $\gamma=0.20\,mm^{-1}$ with parameters $v_{s}=0.1\,mm$, $d_{r}=0.067\,rads^{-1}$ taken from the literature and $B=7.5\,s$ as fitting parameter. The inset of the figure shows the average radius at stationarity as a function of the rotation angular velocity. As one can see the full line, obtained on the basis of Eq.(6) (once the effect of background nonmotile impurities is taken into account, see Cencini et al. (2016) for details), is in very good agreement with the experimental data demonstrating the validity of the model (4).

III Gyrotactic trapping in laminar and turbulent Kolmogorov flow

The mechanisms of cell accumulation described above rely on the assumption that $\Psi<1$, since they are based on the existence of a quasi-static solution of Eq.(1). However, when $\Psi>1$, the second term dominates and the cell tumbles with no persistent direction.

Laboratory experiments Durham et al. (2009) demonstrated that in an inhomogeneous shear flow with vorticity varying along the gravity direction, gyrotactic cells can get trapped where $\Psi>1$ locally. Such gyrotactic trapping has been proposed as a possible mechanism for the formation of thin phytoplankton layers (TPLs) often observed in ocean coastal areas Dekshenieks et al. (2001); Cheriton et al. (2007); Churnside and Donaghay (2009); Steinbuck et al. (2009). TPLs are regions of high vertical concentration of phytoplankton, centimeters to one meter thick, which extend horizontally up to kilometers and last from hours to a few days (see Durham and Stocker (2012) for a review). Though TPLs can be formed by non-motile or motile species characterized by different swimming styles, Heterosigma Akashiwo – a toxic, gyrotactic algae – is known to form harmful thin layers Ryan et al. (2011) and gyrotactic trapping could be an effective explanation. For the sake of completeness, we mention that other mechanisms, not discussed here, of phytoplankton cells accumulation in shear flows can emerge when the cell shape is elongated, see e.g. Barry et al. (2015).

III.1 The Kolmogorov Flow (KF)

We focus on the Kolmogorov Flow (KF), a periodic shear flow model, much studied for the transition to turbulence Sivashinsky (1985); She (1987); Borue and Orszag (1996), characterized by an inhomogeneous distribution of vorticity, necessary for gyrotactic trapping. The KF solves the incompressible Navier-Stokes equations subjected to a monochromatic body force:

[TABLE]

where $p$ is the pressure, density is taken to unity $\rho=1$, and $\hat{\bm{x}}$ denotes the unit vector in the horizontal direction. The physical domain is a cube of side $L_{B}=2\pi L$ with periodic boundary conditions in all directions. Equation (7) admits a stationary solution, the laminar KF $\bm{u}=U\cos(z/L)\hat{\bm{x}}$ with $U=L^{2}F/\nu$, which becomes unstable to large scale perturbations when the Reynolds number, $Re=UL/\nu$, exceeds the critical value $Re_{c}=\sqrt{2}$. The first instability is two-dimensional, at increasing $Re$ three-dimensional motion develops and the flow eventually becomes turbulent Borue and Orszag (1996); Musacchio and Boffetta (2014). Also in the turbulent KF the time averaged velocity $\bar{\bm{u}}$ remains monochromatic, $\overline{\bm{u}}=U\cos(z/L)\hat{\bm{x}}$, but with a different amplitude $U<L^{2}F/\nu$ Musacchio and Boffetta (2014). This latter property is useful to study the effect of turbulent fluctuations, which will be briefly discussed at the end of the section.

III.2 Gyrotactic Cells in Laminar Kolmogorov Flow

After non-dimensionalization by measuring lengths, velocities and times in terms of $L$, $U$ and $L/U$, for the KF Eq. (2) reads

[TABLE]

where $\Phi=v_{s}/U$ is the ratio of the swimming speed over the flow velocity. While, Eq. (1) becomes

[TABLE]

where $\Psi=BU/L$ measures the cells’ stability. The box size in dimensionless units is $L_{B}/L=2\pi$. Capital letters denote the coordinates in non-dimensional units.

Considering that $|\mathrm{\bf p}|=1$, and noticing that $X$ and $Y$ do not enter the dynamics of the other variables, the dynamics is effectively three dimensional. Moreover, using Eqs. (10)-(12) and Eqs. (10)-(11) one can easily check that the system admits the two following conserved quantities, respectively:

[TABLE]

Thus the dynamics has only one degree of freedom and it is integrable (see Zöttl and Stark (2012, 2013) for a similar approach for bacteria swimming in a shear flow).

For vertically migrating cells, the conservation of $\mathcal{C}$ in Eq. (14) implies that $p_{y}$ should vanish exponentially as $Z$ increases. This observation tells us that asymptotically the motion is described by the two-dimensional dynamics

[TABLE]

where $(\mathrm{p}_{x},\mathrm{p}_{z})=(\cos\theta,\sin\theta)$. Strictly speaking the reduction to the above two dimensional system applies only for vertically migrating cells, however the following derivations extend also to the general three-dimensional dynamics also for non-vertically migrating cells, see Santamaria et al. (2014) for details. Here below we will denote with $Z$ the unrestricted vertical coordinate and with $z=Z\mod 2\pi$ its restriction to the periodic cell, while $\theta\in[-\pi,\pi]$.111We remark that $\mathcal{H}$ (15) is not periodic in $Z$: when $Z\to Z\pm 2\pi\,n$, we have $\mathcal{H}(\theta,z)\to\mathcal{H}(\theta,z)e^{\pm\pi n/(\Phi\Psi)}$. We also notice that $\mathcal{H}(\theta,z)$ plays a similar role to the Hamiltonian and though the system is not Hamiltonian it can be made so by a non-canonical change of variables Santamaria et al. (2014).

For $\Psi<1$, Eqs. (16)-(17) do not admit fixed points and $Z$ grows in time. The exponential dependence on $Z$ in (15) and the conservation of $\mathcal{H}$ implies that the term in square brackets in (15) must decrease exponentially with $Z$. Therefore, for large $Z$, the swimming angle depends on $z$ (i.e. restricted on the torus) as given by

[TABLE]

Thus the vertical velocity will change with height and cells will accumulate where it is minimal, i.e. for $Z=n\,\pi-\arctan[1/(2\Phi\Psi)]$, for any integer $n$. Around these positions one observes ephemeral layers of high density of cells, even though cells are trapped asymptotically . This is confirmed in Fig. 3a showing the time evolution of the vertical probability density distribution (PDF), $\rho(Z,t)$, resulting from an initially uniform distribution in $Z\in[0:2\pi]$ for $\Psi<1$. The transient accumulations last longer for smaller values of $\Phi$.

Since $|\cos\theta|\leq 1$, after some algebra, one can realize that Eq. (18) corresponds to a well defined relation between $\theta$ and $z$ on the torus $[-\pi:\pi]\times[0:2\pi]$ either for fast swimmers $\Phi\geq\Phi_{c}=1/2$ or if $\Phi<\Phi_{c}$ whenever $\Psi<\Psi_{c}=(1-4\Phi^{2})^{-1/2}$. Considering first the case $\Phi<\Phi_{c}$, since $\Psi_{c}>1$ we have two possible instances. For $1<\Psi<\Psi_{c}$ (Fig. 3b), we observe the coexistence of upward migrating cells that, as for the $\Psi<1$ case, asymptotically satisfy Eq. (18), and of cells remaining vertically confined with high concentration in thin regions unrelated to the minima of $\mathrm{p}_{z}$. Conversely, for $\Psi>\Psi_{c}$ (Fig. 3c), we observe that all cells are asymptotically vertically trapped and organized in thin layers.

The origin of these two behaviors can be traced back to the properties of the phase-space dynamics on the $(\theta,z)$ torus. We start noticing that for $\Psi>1$, Eqs. (16)-(17) possess two hyperbolic (at $(\theta^{*},z^{*})=(0,2\pi-\arcsin\Psi^{-1})$ and $(\pi,\pi-\arcsin\Psi^{-1})$) and two elliptic fixed points (at $(\theta^{*},z^{*})=(0,\pi+\arcsin\Psi^{-1})$ and $(\pi,\arcsin\Psi^{-1})$). However, depending on $\Psi$ being smaller or larger than $\Psi_{c}$ the separatrices, corresponding to the orbits emerging from the hyperbolic fixed points, qualitatively change (see Fig. 4). $1<\Psi<\Psi_{c}$ (Fig. 4a), the hyperbolic point is the vertex of a slip-knot containing the elliptic fixed point. Orbits whose initial condition falls into this loop cannot escape it, those starting outside it migrate vertically asymptotically following (18). We thus have only partial trapping. Conversely, for $\Psi>\Psi_{c}$ (Fig. 4b), the separatrices roll up around the torus in the $\theta$ direction becoming impenetrable barriers to vertical transport and thus trapping the cells. As clear from Fig. 4, the compression between the elliptic and hyperbolic points is responsible for the higher concentration layers observed in Fig. 3b,c. When $\Phi>\Phi_{c}$ only purely vertically migrating cells (for $\Psi<1$) or partial trapping (for $\Psi>1$) is possible. Figure 3d summarizes the various regime in the parameter space $(\Phi,\Psi)$.

III.3 Gyrotactic Cells in Turbulent Kolmogorov Flow

Fluctuations of velocity and/or vorticity, due to turbulence, allow a swimmer to escape from the homoclinic loops in the partial trapping regime or from the separatrices in the trapping one (Fig. 4). As a result, layers are less intense and, more importantly, they become transient. Below, we briefly discuss the transient dynamics of such layers, for further details see Santamaria et al. (2014). Such a phenomenology is indeed observed in field experiments with simultaneous measurement of biological and physical properties, which have shown that while thin layers are weakly affected by turbulence of moderate intensity, stronger turbulence dissolves them Wang and Goodman (2010); Sullivan et al. (2010). As discussed in Santamaria et al. (2014) a similar effect can be induced, in the laminar KF, by the presence of rotational diffusion, but this effect is typically much smaller than that due to turbulence.

Increasing the Reynolds number the Kolmogorov flow becomes unstable and eventually turbulence for large values of $Re$. Nonetheless, since the (time-averaged) mean flow remain monochromatic, velocity and vorticity fields can be still decomposed in the time-averaged fields with superimposed fluctuations, i.e. $\bm{u}^{\prime}$ and $\bm{\omega}^{\prime}$ as $\bm{u}=U\cos(z/L)\hat{\bm{x}}+\bm{u}^{\prime}({\bm{x}},t)$ and $\omega=-(U/L)\sin(z/L)\hat{\bm{y}}+\bm{\omega}^{\prime}({\bm{x}},t)$. However, even at relatively low $Re$, fluctuations are non-negligible with respect to the mean flow, indeed it has been found that $u^{\prime}_{\rm rms}/U\simeq 0.5$Musacchio and Boffetta (2014). Conversely, in real oceans fluctuations are typically much smaller than the mean flow being depleted by, e.g., stratification Thorpe (2007). For this reason, here, while solving Eq. (7) by means of direct numerical simulations (see Santamaria et al. (2014) for details), gyrotactic swimmers are evolved by modulating the fluctuations $\bm{u}^{\prime}$ and $\bm{\omega}^{\prime}$ with a multiplicative factor $\gamma<1$, so to control their intensity. Besides the possibility to control fluctuation intensity, another advantage of such an approach is that the statistical properties of the turbulent fluctuations do not change with $\gamma$ as they would by introducing stratification.

Numerical simulations show that as soon as turbulent fluctuations are considered (i.e. $\gamma>0$), even if $\Phi$ and $\Psi$ are chosen in the trapping region (black circle in Fig. 3d), vertical migration becomes possible as clear by comparing Figs. 3e-g with Fig. 3c. At moderate values of the turbulent intensity (Fig. 3e), velocity and vorticity fluctuations allow cells to escape from the trapping regions by moving them to regions of lower shear, where upward directed swimming is possible. Then cells get trapped around another high shear region. We recall that the periodic layer structure is inherited from the flow periodicity. As a result, the average vertical cell velocity, $\langle v_{z}\rangle$, which was zero in the absence of turbulent fluctuations, becomes positive (Fig. 5a). However, very intense turbulence rotates the cell swimming direction randomly and, moreover, fluctuations of the vertical velocity also mix cells. As a consequence, the average vertical motion $\langle v_{z}\rangle$ decreases for large values of $\gamma$. An intermediate turbulence intensity maximizes the vertical migration velocity.

From the point of view of the single cell dynamics, the above phenomenology means that, for low intensity fluctuations, a cell is trapped for a finite time until fluctuations of the vertical velocity or of the vorticity makes the cell able to escape its trapped trajectory and to swim upwards until it gets trapped again. At high turbulent intensity trapping is less and less effective, even transiently, and the motion becomes basically diffusive in the vertical direction due to turbulent diffusion.

The above observation suggests to model the vertical dynamics of gyrotactic swimmers in terms of a diffusive process with drift $V_{d}$ (due to the average vertical migration speed, i.e. $V_{d}=\langle v_{z}\rangle$) and diffusion constant $D_{z}$, both depending on turbulent fluctuations. With this modelization in mind one can study the exit times statistics – a standard problem in stochastic processes Redner (2001) – from a layer, i.e. of the times $T$ needed for a swimmer to travel a vertical distance $L_{B}/2$, separating two consecutive layers. For diffusion with drift, the probability density function of the exit time $T$ is given by the inverse Gaussian function, which reads

[TABLE]

The above expression provides a prediction for the exit time PDF that can be directly tested against the measured one. For the drift velocity we have $V_{d}=\langle v_{z}\rangle$, which is measured in DNS. The diffusion constant $D_{z}$ can be estimated by measuring $\langle T^{2}\rangle$ in the DNS and noticing that in Eq. (19) $\langle T^{2}\rangle=L_{B}(D_{z}+LV_{d}/2)/V_{d}^{3}$.

Figure 5b shows the comparison between measured exit-time PDF $p(T)$ and the inverse Gaussian prediction (19), with $D_{z}$ and $V_{d}$ obtained as discussed above. The prediction results rather accurate for the right tail (long exit times) for all turbulent intensities $\gamma$, while the left tail is fairly well captured only for large values of $\gamma$. The latter deviations can be interpreted as the result of the suppression of fast escapes due to gyrotactic trapping, which is more effective in the limit $\gamma\to 0$. Conversely, long escape times allow trajectories to sum-up many uncorrelated contributions, thus recovering a diffusive behavior, which explains the good agreement on the right tail. The average exit time of single trajectories, $T_{e}$, is given by $T_{e}=L_{B}/(2\langle v_{z}\rangle)$ (as implied by the argument of the exponential in (19). While for $D_{z}$ the estimation is more difficult because the turbulent diffusivity is influenced by swimming, especially for low $\gamma$ values (see Santamaria et al. (2014) for a discussion).

The persistence time of the layer $T_{p}$, can be heuristically estimated as the time needed for $90\%$ of the cells to escape the layer, $\int_{0}^{T_{p}}p(T)dT\approx 0.9$. Using (19) one finds that $T_{p}$ is in the order of a few (typically $\sim 2-3$) $T_{e}$ depending on the value of the vertical diffusivity $D_{z}$. Ignoring some difficulties in estimating $D_{z}$ for realistic oceanic flows, if we consider average swimming speed $\langle v_{z}\rangle$ in the range $0.2-0.6\,v_{s}$, as suggested by Fig. 5a (with $v_{s}\approx 100-200\mu m/s$) and typical lengths $L_{B}$ of the order of a few centimeters, we obtain an estimation of $T_{p}$ from a few hours to a few days, which is akin to values found in TPLs observed on the field Sullivan et al. (2010).

IV Gyrotactic clustering in homogeneous turbulence

As discussed in previous sections, inhomogeneous distributions can arise in dilute suspensions of gyrotactic phytoplankton in laminar flow, either due to directed motility for stable cells ($\Psi<1$) Kessler (1985) or due to trapping of unstable cells ($\Psi>1$) in high shear regions Durham et al. (2009), phenomena well reproduced by the mechanistic model defined by Eqs. (1)-(2). Clustering of cells has been numerically predicted also in a Taylor-Green vortical flow both for gyrotactic swimmers Durham et al. (2011) and for elongated, non-gyrotactic motile cells Torney and Neufeld (2007). It is rather difficult to infer, from these results, the effects of unsteadiness of the flow or of the presence, typical of turbulence, of many active scales and intense fluctuations.

While data on large scale phytoplankton patchiness abound Martin (2003); Mackas et al. (1985), small-scale inhomogeneities have not been widely studied until recently, when in-situ observations by high resolution submarine cameras became feasible Malkiel et al. (1999); Gallager et al. (2004). Although data are still scarce, they provide evidence of clustered distribution (patchiness) of phytoplankton cells on scales comparable to the dissipative scales of turbulence Malkiel et al. (1999). In particular, motile cells exhibit a more intense patchiness on these scales. One possible interpretation is that swimming can interact with turbulence to produce inhomogeneous distributions at small scales.

The above hypothesis seems, apparently, at odds with the mixing properties of turbulent flow, which typically smooth out inhomogeneities, e.g., of tracer particles. However, it is a now common observation in fluid dynamic research that inertial particles, i.e. finite-size impurities with a density different from that of the advecting fluid, distribute inhomogeneously in an incompressible flow, typically in the form of small-scale fractal aggregates. It can be shown, for example by making use of a standard description of the force acting on small spherical particles Maxey and Riley (1983), that particles heavier than the fluid are ejected from vortices, while lighter particles (such as bubbles) are attracted to the center of vortices thus forming quasi-one-dimensional clusters in the core of vortex filaments Squires and Eaton (1991); Calzavarini et al. (2008); Bec et al. (2007), a phenomenon dubbed preferential concentration. It is based on those observations that one can understand how gyrotaxis can indeed give rise to non-trivial cell distributions in generic turbulent flows, even without a definite large-scale mean flow.

As we will see in the following, in order for small-scale clustering to occur, the swimming speed, $v_{s}$, and the gyrotactic orientation time, $B$, of a cell must be of the order of the typical, small-scale speeds and turn-over times of turbulence, i.e. of the Kolmogorov velocity $u_{\eta}\approx(\nu\epsilon)^{1/4}$ and time $\tau_{\eta}\approx(\nu/\epsilon)^{1/2}$ ($\nu$ being the fluid viscosity) Frisch (1995). In the ocean, typical turbulence intensities, measured in terms of the energy dissipation rate $\epsilon$, fall in the range $\epsilon\approx 10^{-8}\div 10^{-4}{\rm m^{2}/s^{3}}$ Thorpe (2007). These values correspond to $u_{\eta}\sim 3\times 10^{-4}\div 3\times 10^{-3}{\rm m/s}$ and $\tau_{\eta}\sim 10\div 10^{-1}{\rm s}$. Given the typical values of $B\sim 1\div 10{\rm s}$ and $v_{s}\sim 100\div 300\mu{\rm m/s}$ Harvey et al. (2015); Guasto et al. (2010), one can see that both $\Phi=v_{s}/u_{\eta}$ and $\Psi=B/\tau_{\eta}$ can be of order 1.

IV.1 Fractal clustering of gyrotactic phytoplankton in turbulence

Oceanic turbulence in the bulk of the mixed layer can be considered as homogeneous and isotropic Thorpe (2007). Direct numerical simulations of homogeneous and isotropic turbulent flows seeded with cells following the dynamics (1)-(2), relevant for most oceanic applications where fluid accelerations are much smaller than gravity, showed indeed that intense clustering (as clear from Fig. 6a) can occur even in the absence of coherent large-scale structures or of a mean flow Durham et al. (2013); Zhan et al. (2014). In analogy to the case of inertial particles Bec (2003), this effect can be explained in the framework of dynamical systems theory. Equations (1)-(2) describe a dissipative, potentially chaotic (depending on the advecting flow) system. It has five, in general non-null, Lyapunov exponents222Five and not six dimensions because ${\bf p}$ is a unit vector. $\lambda_{i}$ whose sum gives the rate of change $\dot{\Gamma}$ of an infinitesimal phase-space volume Ott (1993). The latter is in turn given by the average divergence in phase space of Eqs. (1)-(2), namely

[TABLE]

where $B$ has been made non-dimensional in terms of the stability parameter $\Psi=B/\tau_{\eta}$. Though the above average cannot be done analytically in general, it should be noticed that, in the limit of small $\Psi$ (and therefore fast orientation), $p_{z}\to 1$ and the sign of (20) is negative. This is the signature of the dissipative nature of the system and it has the consequence that an initially uniform distribution of cells will converge onto a fractal set of dimension ${\mathcal{D}}<5$ in phase space and, consequently, produce a distribution of particles in physical space with dimension $D=\min\{{\mathcal{D}},3\}$ Falconer (1986). If $\mathcal{D}<3$ the dimension in real space will be $D=\mathcal{D}<3$ signalling the presence of observable clustering, indeed a homogeneous distribution would give $D=3$.

Figure 6b shows the fractal dimension $D$ as a function of $\Psi$ for different values of the swimming parameter $\Phi$. The particular definition of $D$ used here is the correlation dimension $D_{2}$, which represents the scaling exponent of the $r$-dependence of the probability $p_{2}(r)\sim r^{D_{2}}$ to find a pair of particles with separation less than $r$ Grassberger and Procaccia (1983). This is, of course, only one possible definition of $D$ as cell distributions will in general be multifractal Paladin and Vulpiani (1987). The qualitative behavior of $D_{2}$ as a function of the stability parameter $\Psi$ can be understood considering the limits of fast and slow orientation. In the limit $\Psi\to 0$, cells rapidly orient upwards, so that one can write $\dot{\bf x}={\bf u}({\bf x},t)+\Phi\hat{\bf z}$, in non-dimensional units. The latter expression has zero divergence and leads to no accumulation, the dynamics collapses onto the real space, where it preserves the volume Ott (1993). In the opposite limit $\Psi\to\infty$ cells are unstable and randomly tumble due to the local vorticity fluctuations. Moreover, in this limit, as clear from Eq. (20) $\dot{\Gamma}\to 0$, meaning that the dynamics preserve phase-space volumes so that one expects $\mathcal{D}\to 5$ and thus $D\to 3$. Again this means no clustering. Any clustering must therefore happen for intermediate values of $\Psi$, as indeed shown in Fig. 6b.

The above phenomenological argument suggests the possibility to obtain a quantitative theory for clustering in the limit of fast orientation ($\Psi\ll 1$). In order to make the derivation clear it is useful to better define the parameter range in which we are working: to this aim we refer to the theory developed in Ref. Fouxon and Leshansky (2015). First of all we notice that the orientation dynamics (1) is fully determined by the history of the velocity gradients along the trajectory of the particle. In turbulence the typical correlation time of velocity gradients along a tracer path is of the order of $\tau_{\eta}$. However, for swimming particles such correlation time depends on the swimming speed. Indeed when $\Phi>1$, i.e. $v_{s}>u_{\eta}$, the correlation time will be of the order of $\eta/v_{s}=\tau_{\eta}/\Phi$, being the Kolmogorov length $\eta$ the typical scale of the gradients. Hence the correlation time of the vorticity along the swimming particle path will be $t_{cor}=\min\{\tau_{\eta},\tau_{\eta}/\Phi\}$. As detailed in Ref. Fouxon and Leshansky (2015), from Eq. (1), the above observation and the fact that the vorticity magnitude is order $1/\tau_{\eta}$ one can derive that ${\bf p}$ is weakly affected by turbulence under two circumstances: either when $\Psi\ll t_{cor}/\tau_{\eta}=\min\{1,1/\Phi\}$, or when $\Psi\gg t_{cor}/\tau_{\eta}$ and $\Psi/\Phi\ll 1$. The latter regime corresponds to the case in which the gradients seen by the particle are almost uncorrelated and the vertical polarization results from the central limit theorem. In this case precise predictions on the PDF of the orientation can be obtained Fouxon and Leshansky (2015). In the former regime, instead, the vorticity changes slowly with respect to the orientation dynamics, and local equilibrium is a good approximation. Therefore, similarly to what was done in Sect. II to explain cell focusing, one can impose $\dot{\bf p}=0$ and expand (1) to first order in $\Psi$ obtaining the following expression for the instantaneous equilibrium orientation,

[TABLE]

In Ref. Fouxon and Leshansky (2015) the reader can find a detailed discussion of the limit of validity of the flow description (21) in terms of the appearence of singularities in the gradients of the particle orientation, which are analogous to caustics in inertial particles Wilkinson and Mehlig (2005). We now discuss the consequence of Eq. (21).

First, in this limit the equations of motion reduce to $\dot{\bf x}={\bf v}={\bf u}({\bf x},t)+\Phi{\bf p}_{eq}$. Thus, the swimmers behave as tracers transported by an effective velocity field $\bf v$ with divergence

[TABLE]

Thus the effective velocity field is the sum of an incompressible part (the fluid velocity) and a compressible one proportional to $\Phi\Psi$. This situation is in close analogy with that found when considering the small-Stokes-number limit for inertial particles Falkovich et al. (2002); Balkovsky et al. (2001); Fouxon (2012). Following closely those derivations, one can show that $3-D\sim(\Phi\Psi)^{2}$ as indeed observed in the inset of Fig. 6b, where the left branches of the curves approximately rescale when plotted as a function of $\Phi\Psi$ and are well approximated by a parabola. It is worth observing that the $(\Phi\Psi)^{2}$ dependence for the codimension of the fractal clusters is based on (22) and is valid when both $\Psi\ll 1$ and $\Phi\ll 1$. However, it was then shown Fouxon and Leshansky (2015) (see also Gustavsson et al. (2016)) that the same dependence should be expected also when $\Phi\gg 1$ provided $\Psi\ll 1/\Phi$, though with a different prefactor. These derivations are based on an estimation of the Kaplan-Yorke dimension Fouxon and Leshansky (2015) and perturbative results based on a small-Kubo expansion Gustavsson et al. (2016). The analogy with weakly-compressible flows allows one to extend this argument to the generalized fractal co-dimension of any order $q$, which is predicted to depend linearly on $q$ Fouxon (2012).

Secondly, trajectories will tend to concentrate in regions where the divergence of the effective flow is negative. Therefore, from the expression for $\dot{\Gamma}$ we can predict that cells will preferentially sample regions where $\nabla^{2}u_{z}>0$. By considerations of isotropy of the velocity field and a positive energy dissipation at small scales, one can conclude that such regions are also, on average, regions where $u_{z}<0$ Durham et al. (2013). In other words, fast re-orienting gyrotactic swimmers preferentially sample downwelling regions of the flow, even when the background flow is turbulent and statistically isotropic.333Clearly, the symmetry breaking happens at the level of the dynamics (1). This conclusion is confirmed in Fig. 6c, showing that the average vertical fluid speeds at cells’ positions is negative and, moreover, as shown in the inset rescale with the swimming parameter $\Phi$.

The above argument to explain the prefential sampling of downwelling flow is based on an Eulerian average and not on the average on the particle position. In Fouxon and Leshansky (2015) a different, more precise, derivation was obtained based on the observation that denoting with $n({\bf x},t)$ the density of swimming particles, the average flow vertical velocity at the particle positions is $\langle u_{z}n\rangle$. Then one observes that in the limit in which they can be described as tracers in a compressible flow, $n$ evolves as $\dot{n}=-\bm{\nabla}\cdot{\bf v}=\Psi\Phi\nabla^{2}u_{z}$, where (22) is used. By the formal solution of the above equation one can write

[TABLE]

Again, to show that $\langle u_{z}n\rangle<0$, the Authors of Fouxon and Leshansky (2015) assume that, at equal time, $\langle u_{z}\nabla^{2}u_{z}\rangle=-\langle(\nabla u_{z})^{2}\rangle<0$. Moreover, in the regime in which $\Phi\gg 1$ and $\Psi\ll 1/\Phi$ the negative sign can be proved computing explicitly the correlation in (23) as an Eulerian correlation, because, thanks to the fast swimming speed, the turbulent velocity field can be assumed as frozen on the time scales involved. We also mention that the preferential sampling of downwelling regions was also derived, following a completely different approach, in Ref. Gustavsson et al. (2016), where sharp analytical results were obtained for random flows.

In Refs. Zhan et al. (2014) the model (1) was studied, by means of numerical simulations, for generic values of $\alpha$ thus taking into account the effect of the cell’s shape. Theoretical predictions as a function of $\alpha$ were obtained in Refs. Gustavsson et al. (2016), some of which were recently numerically confirmed in turbulent flows Refs. Borgnino et al. (2018). For elongated, ellipsoidal cells, besides vorticity also the strain rate can indeed influence the cell orientation. Remarkably, it was found, both by theoretical considerations Gustavsson et al. (2016) and by numerical simulations Gustavsson et al. (2016); Borgnino et al. (2018) that when the cells’ aspect ratio increases, i.e. they become more rod-like, and for large swimming velocities the preferential sampling reverses and cells spend more time in upwelling regions. Here we will focus for simplicity on the case of rods. The effect can be appreciated in Fig. 7(a), showing the average of $u_{z}$ computed on the positions of rod-like swimmers in a DNS of homogeneous turbulence Borgnino et al. (2018). Clearly there is a critical swimming speed $\Phi_{c}(\Psi)$ above which they experience positive, instead of negative, average vertical fluid speeds. If this analysis is repeated at varying Re, the observed value of $\Phi_{C}$ increases with Re (Fig. 7(b), inset) and appears to scale with $u_{rms}$ (which is a large scale quantity), as can be seen in Fig. 7(b) (main panel) where a large-scale swimming parameter $\Phi_{L}=v_{s}/u_{rms}$ has been introduced. The latter results lead one to conclude that the inversion in the preferential sampling is controlled by velocity correlations, while the small-$\Phi$ dinamics is controlled by velocity-gradient correlations. In developed turbulence, the two quantities have different characteristic times, namely $\tau_{\eta}$ for the gradients and the large-eddy turnover time $\tau_{L}$ for velocities, with $\tau_{L}/\tau_{\eta}\sim{\rm Re}^{1/2}$, so that the separation of the two dynamics increases at increasing Re. Differences between the elongated and the spherical case were found also for fractal clustering, both in stochastic models Gustavsson et al. (2016) and in DNS of homogeneous turbulence Zhan et al. (2014); Borgnino et al. (2018) It was found that the intensity of clustering, i.e. the extent of the deviation of $D_{2}$ from the homogeneous value $3$, decreases for such more elongated cells. However, elongated cells display a tendency to accumulate even in the absence of gyrotaxis ($\Psi\rightarrow\infty$). A theoretical framework for this phenomenon was provided in the context of stochastic models Gustavsson et al. (2016).

IV.2 Effects of turbulent accelerations on gyrotactic clustering

In the case of intense turbulence, swimmers may experience extreme accelerations La Porta et al. (2001), which in principle require the use of model (4) instead of (2). As seen from the experimental results in Fig. 1, when centrifugal acceleration exceeds gravity swimmers are expected to concentrate towards the center of a stationary vortex. Again one can wonder what would happen when unsteady, turbulent flows are considered. Direct numerical simulations model (1)-(4) in homogeneous and isotropic turbulence De Lillo et al. (2014) show that indeed this dynamics is relevant also in more general flows. Figure 8 shows the correlation dimension of cell clusters in the case where fluid acceleration is explicitly considered. Comparing with Fig. 6b, the main qualitative change due to fluid acceleration is observed for fast orientations ($\Psi\ll 1$): at increasing turbulence, and thus the weight of fluid acceleration, $D_{2}$ deviates more and more from the uniform value $3$ also when $\Psi\to 0$. In the limit of negligible gravity, $D_{2}$ is monotonous in $\Psi$ and one has maximum clustering at $\Psi\rightarrow 0$.

When fluid acceleration is very intense and gravity negligible, the dynamics is better described in terms of a stability number based on the typical acceleration encountered $\Psi_{a}=\omega_{rms}v_{0}/a_{rms}$. With the same ideas used in the previous case, it is easy to realize that when $\Psi_{a}$ is small, cells will orient along the local acceleration so that an effective velocity field can be written as ${\bf v}=u+\Phi\hat{\bf a}$ (with $\hat{\bf a}={\bf a}/|{\bf a}|$). Again one has that swimmers behave as tracers in a compressible velocity field as $\bm{\nabla}\cdot\hat{\bf a}\neq 0$ and they will concentrate in regions characterized by $\bm{\nabla}\cdot\hat{\bf a}<0$. Numerical simulations lead to conclude that the signs of $\bm{\nabla}\cdot\hat{\bf a}$ and $\bm{\nabla}\cdot{\bf a}$ are strongly correlated. The request that the latter quantity be negative is tantamount to requiring that the flow is locally dominated by vorticity, since

[TABLE]

with $\hat{S}_{ij}=\frac{1}{2}(\sigma_{ij}+\sigma_{ji})$ and $\hat{\Omega}_{ij}=\frac{1}{2}(\sigma_{ij}-\sigma_{ji})$ being the symmetric (rate of strain) and antisymmetric (vorticity) part of the velocity gradient $\sigma_{ij}=\partial_{j}u_{i}$, respectively. At increasing $Re_{\lambda}$, therefore, gyrotactic swimmers are more and more attracted towards regions of high vorticity. This is confirmed in the inset of Fig.8, showing the average mean square vorticity measured at the position of cells is plotted. As one can see, swimmers experience larger vorticities than tracers (which would remain homogeneous and therefore be subject to the Eulerian average $\langle\omega^{2}\rangle_{E}$) or purely gravitactic swimmers, which also show very small deviations from the Eulerian value. Remarkably, the above argument is essentially the same one can use to explain the clustering of light inertial particles in turbulent vortex filaments Balkovsky et al. (2001); Calzavarini et al. (2008).

IV.3 Clustering of polydisperse populations

In the above discussion, fractal clustering has been numerically demonstrated for monodisperse suspensions, i.e. when all cells have identical stability and swimming parameters. Natural colonies are of course characterized by a distribution of $\Phi$ and $\Psi$, it is thus natural to wonder whether fractal clustering can be observed in realistic conditions. We discuss here the robustness of clustering in two cases: when the swimming parameters are Gaussian distributed within the swimmers population and when two sub-populations are present with different parameters, e.g. representing two hypothetical strains. These two cases will be referred as Gaussian and bimodal cases.

When considering two sub-population of swimmers, we necessarily have to extend the previous definition of $D_{2}$ to quantify a cross correlation dimension Bec et al. (2005); Borgnino et al. (2017); the latter can be defined via the probability of finding two swimmers, characterized by two values of parameters $(\Phi_{1},\Psi_{1})$ and $(\Phi_{2},\Psi_{2})$, at a distance smaller than $r$, namely $P_{12}(r)\propto r^{D^{(12)}_{2}}$. Clearly, for a monodisperse population (i.e. $\Phi_{1}=\Phi_{2}$ and $\Psi_{1}=\Psi_{2}$), one recovers the correlation dimension introduced in the previous section. In order to understand what we should expect for such quantity, it is useful, for simplicity, to consider the case in which the two sub-populations are characterized by the same swimming number $\Phi=\Phi_{1}=\Phi_{2}$ and a small mismatch in the value of the stability parameter $\Delta\Psi=\Psi_{1}-\Psi_{2}$. Using non-dimensional units, from Eq.(2) one can write that the separation ${\bf R}$ between two swimmers evolves according to $\dot{\bf R}=\Delta{\bf u}({\bf R})+\Phi\Delta{\bf p}$. At very small separation, i.e. below the Kolmogorov length $\eta$, the velocity field is smooth and we can write $\Delta{\bf u(R)}\sim u_{\eta}(R/\eta)$. Thus from the balance of the two terms in the above equation a characteristic scale emerges Borgnino et al. (2017):

[TABLE]

A similar argument can be done for swimmers with the same stability parameter and different swimming velocities where one finds $R^{*}\simeq\eta\Psi\Delta\Phi$ Borgnino et al. (2017). For $R<R^{*}$ the separation dynamics is dominated by the parameter mismatch and one expect poor correlation between the swimmers, i.e. homogeneous distribution.

For $R>R^{*}$, conversely, the dynamics is dominated by the velocity field and it is almost indistinguishable from that of two cells with the same stability parameter. In this latter case, one expects to observe fractal clustering. In practice, this means that the cross correlation dimension will depend on $r$ and $D_{2}^{(12)}(r)\approx 3$ for $r<R^{*}$ and $D_{12}(r)<3$, signaling fractal clustering, above $R^{*}$. This is very similar to what happens to inertial particles with slightly different Stokes times Bec et al. (2005). This scenario is confirmed in Fig. 9, which shows the cross probability $P_{12}(r)$ that quantifies the case of two sub-population with different values of $\Psi$. Curves obtained with different $\Delta\Psi$ but the same $\overline{\Psi}=\frac{1}{2}(\Psi_{1}+\Psi_{2})$ display the same behavior at very small scale, where the scaling exponent is close to $3$, thus the sub-populations considered see each other as uniformly distributed or, in other words, they have uncorrelated distributions. Conversely, at larger scales, a nontrivial power-law behaviour is observed with a exponent close to that one of a homogeneous population with stability number $\overline{\Psi}$. Finally, the inset of Fig. 9 confirms the linear scaling of the characteristic scale $R^{*}$ discussed above by showing a remarkable collapse of the different crossover scales when plotted as a function of the combination $\Phi\Delta\Psi$. Similar results, both concerning the cross probability and the characteristic scale, can be obtained when a bimodal distribution with two different swimming numbers $\Phi_{1}$ and $\Phi_{2}=\Phi_{1}-\Delta\Phi$ and same $\Psi$ is considered.

We now consider a more realistic polydisperse suspension, whose stability number $\Psi$ (or equivalently $\Phi$) is Gaussian distributed with mean value $\overline{\Psi}$ and standard deviation $\sigma_{\Psi}$. In particular, assuming a $\sim 30\%$ of relative variation in gyrotactic parameters, which is rather realistic Sengupta et al. (2017); Harvey et al. (2015), is fractal clustering still observable? Here, it is useful to consider the cumulative probability, $P(r)$, of having two cells separated by a distance less than $r$ integrated on all the possible pairs within the population. In the same way of the bimodal distribution, as shown in Fig. 10, $P(r)$ is characterized by two different scaling regions: at very small scales all curves recover a uniform scaling $r^{3}$, while for $r\gtrsim\eta$ fractal clustering can be observed, as signaled by a power-law behavior of $P(r)$ with a scaling exponent very close to that expected for a monodisperse population with $D_{2}=D_{2}(\overline{\Psi})\simeq 2$, but for a weak dependence on $\sigma_{\Psi}$. The exponent for the smallest $\sigma_{\Psi}$ is $\simeq 2$ while it grows to above $2.3$ for the population with the largest value of $\sigma_{\Psi}$. Also, similarly with the bimodal case, the transition between the two different scaling ranges shifts to larger scales as $\sigma_{\Psi}$ is increased. Remarkably, also the case with largest standard deviation ($\sigma_{\Psi}/\overline{\Psi}\simeq 0.29$) shows a strongly inhomogeneous distribution at large scale indicating the presence of fractal clustering even in a population with a high variability in the gyrotactic parameters. This result provides a strong indication that fractal clustering should be detectable also in experimental suspensions.

V Conclusions and perspectives

The range of flow effects on motile microorganisms here reviewed is by no means exhaustive. We mostly focused on those aspects which can be well characterized and understood by using tools and ideas from dynamical systems. To this aim we have specialized our discussion on gyrotactic motility showing that ideas from integrable (conservative) dynamical systems can explain the phenomenon of gyrotactic trapping in terms of barriers to transport induced by the separatrices emerging from hyperbolic fixed point of the dynamics. Similar ideas can, and have been used, e.g., to explain the behavior of bacterial trajectories in shear flows Zöttl and Stark (2012, 2013). On the other hand the physics of dissipative, chaotic dynamical systems has been used to successfully explain micro-patchiness of gyrotactic phytoplankton in turbulent flows. Again similar ideas can be extended to other kind of motile microorganisms Torney and Neufeld (2007) and also to artificial micro-swimmers, such as phoretic colloids Schmidt et al. (2016); Shukla et al. (2017).

More in general we think that the kind of approach here reviewed can be helpful, possibly by also accounting for stochastic effects (which can also alter the swimming direction via, e.g., rotational Brownian motion or intrinsic stochasticity of the propulsion mechanisms) for understanding most of the effects arising from the interplay between flow and motility in dilute suspensions of living or artificial micro-swimmers. In particular, we believe there are many directions that still need to be explored or fully understood.

For instance, the dynamics of nonspherical gyrotactic organisms Zhan et al. (2014); Gustavsson et al. (2016); Borgnino et al. (2018) certainly deserves further study in view of its richer phenomenology. Still in the context of gyrotactic motility, some recent experiments have shown that gyrotactic cells can, to some extent, display some level of adaptation when exposed to frequent reorientation, such as those induced by turbulence, for a long time Sengupta et al. (2017). It would thus be very interesting to model and characterize such behaviors. Also the motility of other kinds of phytoplankton cells have been shown to be influenced by shear flows, for instance, at high shear rate, Dunaliella primolecta swims in the direction of local flow vorticity Chengala et al. (2013). Interestingly, this behavior seems to be the result of active shear-induced response and this opens several interesting questions both in the direction of understanding the origin of such adaptation and the consequences for their spatial distribution.

Another direction of interest is to investigate the combined effect of flow and directed motility as due to chemotaxis. As briefly discussed in the introduction, the presence of shear can trap bacteria and deplete their chemotactic efficiency Rusconi et al. (2014). This effect is mostly due to Jeffery orbits which align elongated cells along the shear. How and to what extent similar phenomena alter the motility in unsteady, turbulent flows, such as those encountered in the oceans, is largely unknown. In the presence of turbulence, also the chemical field is advected and mixed by the flow. What is the effect on chemotaxis? Some attempt in this direction Taylor and Stocker (2012) has shown that the chemotactic response can be altered by turbulence also without considering the direct impact of the flow on motility. It would thus be very interesting to include and quantify such effects. This is particularly interesting in view of the fact that aquatic bacteria often display different motility strategies with respect to enteric ones Stocker and Seymour (2012), which may be an indication of adaptation to the exposure to turbulent motion for such microbes.

Acknowledgements.

We thank M. Barry, E. Climent, W. M. Durham, M. Franchino F. Santamaria, R. Stocker, for their contribution to the papers that form the basis of this review. We acknowledge support from the COST Action MP1305 “Flowing Matter”. G.B., M.B. and F.D. acknowledge support from the “Departments of Excellence” (L. 232/2016) grant, funded by the Italian Ministry of Education, University and Research (MIUR), Italy.

Authors contribution statement

All the authors contributed equally to the writing of this review.