# Accumulation and alignment of elongated gyrotactic swimmers in turbulence

Zehua Liu  Department of Engineering Mechanics, School of Aerospace
Engineering, Tsinghua University, Beijing 100084, China  Linfeng Jiang
<EMAIL_ADDRESS>Center for Combustion Energy, Key Laboratory
for Thermal Science and Power Engineering of Ministry of Education, Department
of Energy and Power Engineering, Tsinghua University, Beijing 100084, China
Chao Sun  Center for Combustion Energy, Key Laboratory for Thermal Science and
Power Engineering of Ministry of Education, Department of Energy and Power
Engineering, Tsinghua University, Beijing 100084, China  Department of
Engineering Mechanics, School of Aerospace Engineering, Tsinghua University,
Beijing 100084, China

###### Abstract

We study the dynamics of gyrotactic swimmers in turbulence, whose orientation
is governed by gravitational torque and local fluid velocity gradient. The
gyrotaxis strength is measured by the ratio of the Kolmogorov time scale to
the reorientation time scale due to gravity, and a large value of this ratio
means the gyrotaxis is strong. By means of direct numerical simulations, we
investigate the effects of swimming velocity and gyrotactic stability on
spatial accumulation and alignment. Three-dimensional Voronoï analysis is used
to study the spatial distribution and time evolution of the particle
concentration. We study spatial distribution by examing the overall
preferential sampling and where clusters and voids (subsets of particles that
have small and large Voronoï volumes respectively) form. Compared with the
ensemble particles, the preferential sampling of clusters and voids is found
to be more pronounced. The clustering of fast swimmers lasts much longer than
slower swimmers when the gyrotaxis is strong and intermediate, but an opposite
trend emerges when the gyrotaxis is weak. In addition, we study the
preferential alignment with the Lagrangian stretching direction, with which
passive slender rods have been known to align. We show that the Lagrangian
alignment is reduced by the swimming velocity when the gyrotaxis is weak,
while the Lagrangian alignment is enhanced for the regime in which gyrotaxis
is strong.

##  I Introduction

Understanding the dynamics and transport of motile particles in a flow is
essential for the ecology of aquatic environments. Marine bacteria and
phytoplankton live in a world of chemical and physical gradients, and many
species actively propel themselves to exploit such heterogeneity [1, 2].

The swimming direction of a micro-organism is influenced by the viscous torque
exerted on cells by fluid shear, and by the response to external gradients and
biases like gradients of dissolved organic matter and gravity. According to
the type of biases, the resulting directed motility is termed chemotaxis [3],
gyrotaxis [4], and phototaxis [5, 6].

In this study, we consider small, elongated, gyrotactic, swimming particles in
homogenous isotropic turbulence (HIT). Many motile phytoplankton species are
gyrotactic, i.e., their swimming direction results from the competition
between shear-induced viscous torque and the stabilizing torque due to bottom-
heaviness [7]. The stabilizing torque biases their swimming in the vertical
direction, inducing an accumulation of cells in well-lit regions near the
surface during daylight. Moreover, both single phytoplankton cells and
multicellular phytoplankton chains can have elongated shapes, which makes the
study of gyrotactic active particles with prolate shapes necessary [8].

The velocity gradients of flows exert viscous torques on microorganisms, which
can modify their swimming direction and transportation. Flow can thus affect
the spatial distribution of microorganisms, which then determines their
encounter rates with prey, predators, and conspecifics [9]. The viscous
torques experienced by microorganisms have a shape dependence. While vorticity
generates a torque for both spherical and elongated cells, the strain of rate
only affects the swimming direction of elongated cells. The additional source
in torques for elongated cells causes their different spatial distribution
from spherical ones. In laminar flows, spherical cells accumulate in the
center of downwelling vertical pipe flows [4, 10], become trapped in a thin
layer of high shear in horizontal shear flows [11, 12]. Gyrotactic swimmers
are found to accumulate due to gravitactic focusing and wall accumulation in
channel Poiseuille flows [13, 14]. In turbulent flows, gyrotactic motility can
generate small-scale clusters that are dynamically evolving, both when the
effect of the fluid acceleration is included or not [15, 16, 17]. Although
both spherical and elongated cells can form clusters in turbulent flows, the
extent of the clustering and where cells accumulate in the water column are
very different. The clustering of elongated gyrotactic swimmers is shown to be
generally weaker than spherical ones, except when the gyrotaxis is very weak
[18, 19]. While spherical gyrotactic cells collect in downwelling regions
[15], elongated counterparts preferentially visit either downwelling or
upwelling regions, depending on the swimming velocity and the gyrotaxis
strength [20, 21, 22]. To date, the spatial structures and temporal evolution
of clusters of gyrotactic swimmers are still open questions.

Owing to the incompressibility of the fluid, clusters are formed due to the
intrinsic swimming of gyrotactic swimmers, with their swimming direction
converging within specific regions of the flow. Therefore, the analysis for
the alignment is important, which can also influence ocean light climate by
altering light scattering [23]. To describe the alignment of elongated
gyrotactic swimmers, one needs a reference frame. In an Eulerian approach, the
alignment with the eigenvectors of the strain rate and the fluid velocity has
been studied [18, 19, 24]. It has been found that a simpler picture emerges if
one describes the local alignment in the reference formed in the Lagrangian
framework [25, 26, 27]. The physical image of the Lagrangian reference frame
is suggested as below. An infinitesimal, spherical fluid element deforms, over
time, into a tri-axial ellipsoid as it is stretched by the turbulent flow,
whose dynamic evolution can be naturally described by following its
trajectory. After a sphere is distorted into an ellipsoid, the orientations of
the three principal axes form an orthogonal coordinate system, which is in
accordance with the eigenvectors of the left Cauthy-Green strain tensor [25].
The Lagrangian direction is defined to be the orientation of the longest
principal axis of this ellipsoid, implying the strongest stretching along this
direction. However, the orientation behavior of active particles with respect
to the Lagrangian stretching direction has not been addressed.

In this study, we focus on the relationship between the clustering of
gyrotactic particles and the flow, using the three-dimensional Voronoï
analysis. The Voronoï method allows us to track the concentration of particles
in the Lagrangian framework. We examine the p.d.f.s of the Voronoï volume, the
locations where clusters form, and the lifetime of the clusters. The remaining
paper is organized as follows. In Sec. II, we briefly describe the background,
explaining the model equations for gyrotactic self-propelled particles and the
numerical method we use. The definitions of three-dimensional Voronoï
tessellation and Lagrangian stretching are also discussed in this section. In
Sec. III, we present results on the accumulation and alignment. The
conclusions, Sec. IV, summarize the main finding of this study.

##  II Method

###  II.1 Model equations

The turbulent flow is governed by the incompressible Navier-Stokes equations
(NSE) driven by an external random large-scale statistically homogeneous and
isotropic force with constant in-time global energy input:

$\displaystyle\frac{\partial{\bm{u}}}{\partial t}+\bm{u}\cdot{\nabla}\bm{u}$
$\displaystyle=$ $\displaystyle-\rho^{-1}{\nabla}p+\nu\
\nabla^{2}\bm{u}+\bm{f},$ (1) $\displaystyle{\nabla}\cdot\bm{u}$
$\displaystyle=$ $\displaystyle 0.$ (2)

Here, $\bm{u}(\bm{x},t)$ is the velocity field, $p(\bm{x},t)$ is the pressure,
$\nu$ is the kinematic viscosity, and $\rho$ is the liquid density. The vector
$\bm{f}$ refers to the external large-scale body force needed to sustain
statistically stationary turbulence. The force $\bm{f}$ to stir turbulence is
modulated by means of a sum of sine waves with small wavenumbers and applied
at each position and at each time step. The phases of sine waves are evolved
in time by means of a stochastic process to ensure a homogeneous and isotropic
stirring [28]. The turbulent intensity of the flow is characterised by a
dimensionless parameter, the Taylor-Reyolds number,
$Re_{\lambda}=u_{rms}\lambda/\nu$, where $u_{rms}=\sqrt{\langle
u_{i}u_{i}\rangle^{v,t}/3}$ is the root-mean-square (r.m.s.) velocity of the
flow ($\langle\ldots\rangle^{v,t}$ denotes here volume and time average),
$\lambda=\sqrt{15\nu u_{rms}^{2}/\epsilon}$ is the Taylor length scale and
$\epsilon=(\nu/2)\Sigma_{i,j}\langle(\nabla_{i}u_{j}+\nabla_{j}u_{i})^{2}\rangle^{v,t}$
is the mean global energy dissipation rate.

A gyrotactic swimmer is modeled as an axisymmetric ellipsoid that swims along
with its long axis direction $\bm{p}$ with a constant swimming speed $v_{s}$.
The total particle velocity is the sum of the swimming velocity $v_{s}\bm{p}$
and the flow velocity $\bm{u}$ at the particle location. Assuming the swimmers
are small compared to the dissipative scale of the turbulent flow and
neutrally-buoyant, particle trajectories, $\bm{x}(t)$, and swimming
directions, $\bm{p}$ can be integrated using the following equations [29, 4,
7]:

$\displaystyle\frac{\mathrm{d}\bm{x}}{\mathrm{d}t}$ $\displaystyle=$
$\displaystyle\bm{u}(\bm{x}(t),t)+v_{s}\bm{p},$ (3)
$\displaystyle\frac{\mathrm{d}\bm{p}}{\mathrm{d}t}$ $\displaystyle=$
$\displaystyle\frac{1}{2B}[\bm{z}-(\bm{z}\cdot\bm{p})\bm{p}]+\bm{\Omega}\bm{p}+\Lambda\left(\bm{S}\bm{p}-(\bm{p}\cdot\bm{S}\bm{p})\bm{p}\right),\quad$
(4)

In the above equations, $\bm{z}$ is the vertically upward direction,
$\bm{\Omega}=\frac{1}{2}(\nabla\bm{u}-(\nabla\bm{u})^{\mathrm{T}})$ and
$\bm{S}=\frac{1}{2}(\nabla\bm{u}+(\nabla\bm{u})^{\mathrm{T}})$ are rotation
and strain rate tensor, respectively, $\Lambda=(\alpha^{2}-1)/(\alpha^{2}+1)$
is the particle eccentricity, $\alpha$ is the particle aspect ratio, which is
defined as the ratio the length $l$ to the diameter $d$ of the particle. $B$
is the gyrotactic reorientation timescale, which measures how long it takes
for a perturbed cell to return to a vertical direction in a quiescent flow.
According to the bottom heaviness, the reorientation time scale can be
computed using $B=\nu\alpha_{\perp}/2gh$, where $\nu$ is the kinematic
viscosity of the fluid, $\alpha_{\perp}$ is a shape-dependent dimensionless
coefficient that measures the viscous resistance of the tumbling rotation of a
particle, $h$ is the distance between the center of gravity and buoyancy, and
$g$ is gravity[22, 7]. Therefore, a strong gravity effect or equivalently a
strong bottom-heaviness-induced gravitational torque is corresponding to a
small $B$. The three terms on the right-hand side of (4) reflect that the
tumbling rate of an elongated gyrotactic swimmer is governed by the
gravitational torque, and the torque of the ambient flow through the fluid’s
vorticity and strain of rate, respectively. In what follows, only particles
with $\alpha=20$ will be considered, corresponding to $\Lambda=0.995$. This
choice of particle aspect ratio is based on two considerations. First, the
individual cell of many phytoplankton species is elongated, often with large
aspect ratios [30]. Secondly, the chain formation of phytoplankton when
daughter cells remain attached to one another after cell division leads to
large aspect ratios. Considering the chain formed by individual spheric
particles, $\Lambda$ rapidly saturates with the chain length. Chains of 2, 4,
8, 16, and 32 cells are corresponding to
$\Lambda=0.600,0.882,0.969,0.992,0.998$ [22]. This large aspect ratio
($\Lambda$ approaches 1) has also been adopted by many previous studies to
investigate the clustering behavior of prolate particles [21, 22]. The
suspension is assumed to be diluted, so the interactions between particles are
neglected. For studies about particle interaction, we refer to Ref. 31, 32,
33.

For elongated gyrotactic swimmers herein, two dimensionless parameters control
the particle’s fate:

$\displaystyle\Phi=v_{s}/u_{\eta},$ (5) $\displaystyle\Psi=\tau_{\eta}/B.$ (6)

First, the swimming number $\Phi$ in (5) measures the swimming velocity
$v_{s}$ relative to the Kolmogorov velocity $u_{\eta}=(\nu\epsilon)^{1/4}$.
Second, the stability number $\Psi$ in (6) is defined to be the ratio of
Kolmogorov time $\tau_{\eta}=(\nu/\epsilon)^{1/2}$ to the gyrotactic
reorientation time $B$, as a predictor of a particle’s ability to maintain
vertical migration. A small gyrotactic reorientation times $B$ implies that
the particle can quickly orient upward again after a perturbation. Therefore,
$\Psi$ parametrizes the importance of directional swimming with respect to the
ambient flow overturning. $\Psi\gg 1$ means strong stability, producing a
swimming direction that is very stable and uniformly upwards. Conversely,
$\Psi\ll 1$ reflects that the gyrotactic effect is negligible, and the
swimming direction is completely determined by the ambient flow.

Figure 1: Visualization of a three-dimensional Voronoï tessellation. Particles
are represented by red points. Each Voronoï cell is a polyhedron assigned to a
particle at its center.

###  II.2 Three-dimensional Voronoï analysis

Previous studies have quantified the clustering of gyrotactic active particles
by employing different approaches, such as pair correlation function [18],
clustering index [21], and correlation dimension [15, 16]. In the above
studies, these methods do not provide information on the exact locations where
the strongest clustering occurs or the Lagrangian evolution of the clusters.
An alternative mathematical tool, the Voronoï tessellation, can be used to
study the clustering, which provides more information about the clustering
[34]. As shown in Figure 1, the Voronoï diagram is a unique decomposition of
the three-dimensional space into independent cells associated with each
particle. One Voronoï cell is defined as the ensemble of points that is the
closest to the particle position compared to any other. From the definition of
the Voronoï diagrams, the volume of a Voronoï cell is the inverse of the local
particle concentration. Voronoï volumes can be used to measure the
concentration of each individual particle present in the field, and this
feature makes Lagrangian tracking of concentration along particle trajectories
possible [35]. Ferenc and Néda [36] have found that a $\Gamma$-distribution
can adequately describe the probability density functions (p.d.f.s) of the
Voronoï volumes normalized by the mean volume for randomly distributed
particles. In the three-dimensional case, the $\Gamma$-distribution has the
following prefactor and exponent:

$f(x)=\frac{3125}{24}x^{4}\text{exp}(-5x).$ (7)

Here $x$ is the Voronoï volume normalized by the mean volume. Particles that
have a p.d.f. that deviates the $\Gamma$-distribution, indicate that the
combined effect of flow and characteristics of particles can affect the
spatial distribution. For those particles at the border of the domain, there
may not be enough other particles surrounding them. These Voronoï cells are
ill-defined, and thus not considered for the analyses.

###  II.3 Lagrangian stretching direction

In this part, we briefly discuss how the Lagrangian stretching direction is
obtained. For more details about this Lagrangian reference frame, we refer to
Ref. 37, 25, 27. In order to characterize the deformation experienced by the
fluid element, we first define a deformation tensor $\bm{F}(t)$ as the
solution of the differential equation:

$\frac{d}{dt}\bm{F}(t)=\bm{A}(\bm{x},t)\bm{F}(t).$ (8)

Here $\bm{A}(\bm{x},t)$ is the instantaneous velocity gradient. We obtain
$\bm{F}(t)$ by integrating (8) using the velocity gradient along the particle
trajectories, with initial condition $\bm{F}(0)=\mathbf{I}$ (identity matrix).
For a non-active particle, the particle velocity is always the same as the
background fluid velocity at the particle position. In this case, the physical
image of $\bm{F}(t)$ is suggested as below. Consider a fluid element that
moves with a non-active particle, and this can be done because they have the
same velocity and thus the same trajectory. This infinitesimal initial
spherical fluid element is stretched, over time, into an ellipsoid, due to the
deformation effect exerted by the fluid velocity gradient. And the deformation
tensor characterizes the fluid element deformation in a Lagrangian way [38,
39].

The velocity of a motile particle, $\bm{u}(\bm{x},t)+v_{s}\bm{p}$, is
different from the fluid element velocity at the particle position
$\bm{u}(\bm{x},t)$. We remark that this time there is no real fluid element
following a swimming particle. However, by integrating (8) the velocity
gradient $\bm{A}(\bm{x},t)$ along a swimmer trajectory, an imaginary initial
spherical deformable element is used to follow the swimming particle and
record the stretching of the ambient flow. Although there is a slight
difference between the meaning of the deformation tensor between non-active
and motile particles, in both these two cases, the deformation tensor is a
useful intermediate quantity to form a Lagrangian reference frame at the
particle position and reflect the accumulative stretching on particles exerted
by the ambient flow.

The left Cauthy-Green tensor $\bm{C}^{(L)}$ is then obtained from the inner
product of $\bm{F}$ with itself:

$\bm{C}^{(L)}=\bm{FF}^{T}.$ (9)

The eigenvalues of $\bm{C}^{(L)}$ are denoted with $\Lambda_{i}(i=1,2,3)$, and
eigenvectors $\bm{e}_{Li}$ correspondingly. In the case of an incompressible
flow, $Tr\bm{A}=0$ and the eigenvalues of the left Cauchy–Green tensor satisfy
$\Lambda_{1}(t)\Lambda_{2}(t)\Lambda_{3}(t)=1$. The largest eigenvalue,
$\Lambda_{1}>1$, indicates extension, and the corresponding eigenvector
$\bm{e}_{L1}$ is called the Lagrangian stretching direction. The smallest
eigenvalue, $\Lambda_{3}<1$, indicates contraction in the $\bm{e}_{L3}$
direction. As for the intermediate eigenvalue $\Lambda_{2}$, it can either be
greater or less than $1$, reflecting either stretching or compressing in the
$\bm{e}_{L2}$ direction.

###  II.4 Numerical simulations

The NSE here is numerically simulated by using a Lattice Boltzmann Method
code, the ch4-project [40], which has been extensively employed in the studies
of Lagrangian tracers and point-like particle dynamics in turbulence [41, 42,
43]. More specifically, recent studies of rotational dynamics of inertialess
ansiotropic particles [41] and gyrotactic particles [43] in turbulence have
validated the particle dynamics model and the code through comparison with
experimental results. This code uses a tri-linear scheme for the Lagrangian-
Eulerian frame interpolations. We refer to Ref. 40 for the details of the
numerical method. For the present simulations, a resolution of $256^{3}$ grid
is used for the HIT flow at $Re_{\lambda}\approx 60$. The boundary conditions
are periodic in all directions of the three-dimensional cubic simulation
domain. As indicated in Figure 2, a wide range of particle parameters has been
explored, with $\Phi=1,~{}6.3,~{}\text{and}~{}16$ and $\Psi\in[0.01,16]$. For
each parameter set, $5\times 10^{4}$ Lagrangian trajectories are seeded in the
flow and the duration is about 20 large-eddy turnover times. The deformation
tensor $\bm{F}(t)$ is integrated from the DNS data using a fourth-order Runge-
Kutta scheme according to (8). The random initial conditions are adopted for
the particle positions and orientations. All analyses are performed after
particles have reached a statistically stationary state.

Figure 2: Explored parameter space of $\Psi$ and $\Phi$ for our numerical
study of gyrotactic swimmers.

##  III Results and discussion

We aim to address the following problems: (i) How the clustering is affected
by the swimming number $\Phi$ and the stability number $\Psi$? In particular,
how to characterize the extent of clustering based on the three-dimensional
Voronoï analysis? (ii) How the clustering is related to the flow structures?
This question has two perspectives. From the space perspective, what are the
regions that particles accumulate in? From the time perspective, how long do
particles exist in an aggregated state?

###  III.1 Particle accumulation

Figure 3: The p.d.f.s for the normalized Voronoï cell volumes of fast
($\Phi=16$), moderate ($\Phi=6.3$), and slow swimmers ($\Phi=1.0$) at
$\Psi=0.01$, as compared with the distributions obtained from non-active
particles (dashed line) and randomly distributed particles (7). All three
types of motile particles show clustering, whereas non-active particles remain
randomly distributed. Figure 4: The normalized Voronoï volume p.d.f.s for (a)
slow swimmers ($\Phi=1.0$), (b) moderate swimmers ($\Phi=6.3$), (c) fast
swimmers ($\Phi=16$). The results of different $\Psi$ are plotted for each
type of particle. No preferential concentration is observed in the gravity-
dominated case ($\Psi=16$), whereas rod-like slow and moderate swimmers
maintain observable clustering even when the gravity effect is nearly
negligible ($\Psi=0.01$). At intermediate values of $\Psi$, the clustering is
most pronounced.

As shown in Figure 3, the p.d.f.s of the Voronoï volume ($\mathscr{V}$)
normalized by the global mean value ($\mathscr{\bar{V}}$),
$\mathscr{V}/\mathscr{\bar{V}}$, is plotted for particles of different $\Phi$
at $\Psi=0.01$ to investigate the effect of swimming velocity on the
clustering. The gyrotactic effect can be considered to be negligible at
$\Psi=0.01$. For non-active particles ($\Phi=0$), no preferential
concentration is expected, with the p.d.f. agreeing well with that of randomly
distributed particles. In contrast, for active particles, we observe different
behavior compared to the randomly distributed particles. For all three types
of active particles at $\Psi=0.01$, the possibilities of attaining either
small or large normalized Voronoï volumes are larger than that of randomly
distributed particles, which are corresponding to clusters and voids in the
domain. For slow swimmers, the effect of the swimming velocity on the particle
dynamics is not significant, which leads to a slight increase in the
probability of small Voronoï volumes. However, when the particle dynamics is
dominated by the swimming velocity ($\Phi=16$), particle accumulation becomes
less evident as the particles may tend not to be trapped due to high swimming
velocity. The strongest clustering is observed for moderate swimmers
($\Phi=6.3$). As studied in Ref. 19, anisotropic swimmers without gyrotaxis
accumulate in regions of lower fluid vorticity. In the swimming number range
they covered, $\Phi\in[0,5]$, the extent of accumulation monotonically
increases with the particle swimming velocity. This is consistent with the
variation of the p.d.f.s at $\Psi=0.01$ from $\Phi=1$ to $\Phi=6.3$.

Figure 5:  The normalized standard deviation of Voronoï cell volumes,
$\sigma/\sigma_{\Gamma}$, as a function of $\Psi$ at $\Phi$ of 1.0, 6.3, and
16.

Figure 4 shows the p.d.f.s of the normalized Voronoï volume for different
$\Psi$ at $\Phi$ of 1.0, 6.3, and 16. The gravity-induced torque further
affects the particle dynamic so as to change the particle accumulation. It
turns that, for all swimming numbers, gyrotaxis enhances the probability of
finding either small or large Voronoï volumes and the influence becomes
highest when the flow time scale and the gyrotactic reorientation time are
comparable ($\Psi=1.0$). More importantly, this phenomenon even exists when
the swimming velocity is dominated ($\Phi=16$), whereas no particle
accumulation is observed for non-gyrotactic fast swimmers as shown in Figure
3. Furthermore, the overall extent of particle accumulation increases with the
swimming number. When the gyrotactic effect is dominated ($\Psi=16$), the
particles are compelled to point upward, which leads to random distributions
similar to $\Gamma$-distribution (7).

Figure 6: Several largest clusters and voids in a snapshot of particles at
$\Phi=16$ and $\Psi=0.40$, with the Voronoï-cell faces colored uniformly per
cluster. The particles outside these clusters and voids have been omitted from
the figure. Panel (a) and (b) shows the clusters and voids, respectively.
Figure 7: The mean vertical fluid velocity at the positions of (a) overall
particles, (b) clusters, (c) voids.

The extent of particle accumulation can be further quantified by the standard
deviation of normalized Voronoï volumes, $\sigma/\sigma_{\Gamma}$, where
$\sigma_{\Gamma}$ is the standard deviation of Voronoï volumes for randomly
distributed particles. This quantity has been adopted by other studies on
inertial particle clustering in turbulence [44, 34]. A value of
$\sigma/\sigma_{\Gamma}=1$ represents that no clustering can be observed, and
high values of this term reveal the existence of high and low concentration
events. We remark that $\sigma/\sigma_{\Gamma}$ is dependent on the number of
particles in the domain. To test this particle number dependence, we
calculated $\sigma/\sigma_{\Gamma}$ for a reduced number of particles by
randomly removing some particles. The value of the indicator
$\sigma/\sigma_{\Gamma}$ can be biased by this subsampling procedure, and it
is reduced by nearly $20\%$ when $80\%$ of the $50000$ particles in total
belonging to the group of $\Psi=1$ and $\Phi=16$ are randomly removed from the
domain. However, $\sigma/\sigma_{\Gamma}$ remains a good indicator for a fixed
number of particles, with larger values indicating more clustering. As seen in
Figure 5, for the range of explored stability numbers (spanning from 0.01 to
16), the standard deviation of the normalized Voronoï volumes of all three
types of swimmers shows a non-monotonic trend, with the indicator value
peaking at $\Psi=O(1)$. Under the intermediate gyrotaxis, fast swimmers
exhibit the strongest clustering, reaching a value of 3.9. The clustering
result has a consistent trend with that of Ref. 18, which evaluated the
clustering extent by the slope of the radial distribution function. Similar
clustering extent trend has also been presented through the correlation
dimension $D_{2}$ [21], which rules the small-scale scaling behaviour of the
probability to find a pair of cells with distance less than $r$, i.e.,
$p_{2}(r)\sim r^{D_{2}}$ when $r\to 0$.

###  III.2 Clusters and preferential sampling

In this section, we aim to study how the spatial distribution of particles is
related to the turbulent structures. The particle parameters affect not only
the extent of the clustering but, importantly, also the locations of the
clustering. It has been found that rod-like gyrotactic swimmers preferentially
visit either downwelling or upwelling regions depending on the $\Psi$ and
$\Phi$ [21, 22]. Therefore, we investigate this issue by examing the vertical
fluid velocity sampled by particles and their connections with the Voronoï
cell volumes. In order to identify where clusters and voids form, we denote
$\mathscr{V}_{f}$ as the Voronoï volume which a fraction $f$ of Voronoï cells
are smaller than, i.e. $P(\mathscr{V}<\mathscr{V}_{f})=f$. Here, normalized
Voronoï volumes smaller than $\mathscr{V}_{0.05}$ are defined as belonging to
clusters. Equivalently, Voronoï cells larger than $\mathscr{V}_{0.95}$ are
defined as belonging to voids. In other words, Voronoï cells whose volume is
below the threshold defined by $\mathscr{V}_{0.05}$ belong to clusters, while
Voronoï cells above $\mathscr{V}_{0.95}$ belong to voids. Figure 6 shows an
example of several identified largest clusters and voids in one snapshot from
the simulation. It appears that these Voronoï cells that possess extreme sizes
tend to be connected in groups of various sizes and shapes. The sizes of voids
are naturally larger than that of the clusters, due to the larger sizes of
connected Voronoï cells. And the elongation of the shapes of voids in the
vertical direction is observed. It is now tempting to ask these questions:
where do the overall particles ($\mathscr{V}_{0}<\mathscr{V}<\mathscr{V}_{1}$)
preferentially sample? where do clusters ($\mathscr{V}<\mathscr{V}_{0.05}$)
and voids ($\mathscr{V}>\mathscr{V}_{0.95}$) preferentially form?

First, we focus on the overall preferential sampling. We use the non-
dimensional vertical fluid velocity (normalized by the Kolmogorov velocity)
averaged over the positions of all particles, $\left\langle
u_{z}\right\rangle/u_{\eta}$, to show this preferential sampling of fluid
velocities, where $\langle\ldots\rangle$ denotes the ensemble average over the
locations of overall particles ($\mathscr{V}_{0}<\mathscr{V}<\mathscr{V}_{1}$)
and over time. Similarly, in what follows, $\langle\ldots\rangle_{c,t}$ and
$\langle\ldots\rangle_{v,t}$ denote the time average over the locations of
clusters ($\mathscr{V}<\mathscr{V}_{0.05}$) and voids
($\mathscr{V}>\mathscr{V}_{0.95}$) respectively. In Figure 7(a), the different
curves are plotted for $\left\langle u_{z}\right\rangle/u_{\eta}$ of three
types of swimmers at different values of $\Phi$ as a function of the stability
number $\Psi$. There is a clear difference among slow, moderate, and fast
swimmers. Slow swimmers, at all $\Psi$ in the range we have covered, tend to
stay in the downwelling flow, whereas fast swimmers consistently
preferentially visit upwelling regions. Moderate swimmers, preferentially
sample regions with a positive vertical velocity at small $\Psi$, but regions
with a negative vertical velocity at large $\Psi$. Field observation shows
that the dinoflagellate Cochlodinium polykrikoides tends to be in long chain
groups of $\geq 6$ cells before sunrise, whereas the reverse is observed after
that with the densities of short chain groups of $\leq 5$ cells increases
[45]. The study notes that this is a ecological strategy through which the
population can use the sufficient propulsion produced by more cells to swim
fast to move upward before sunrise. After arrival in the surface layers, the
shape of short chain provides a high area-to-volume ratio for them to capture
energy efficiently. The tendency of fast swimmers to preferentially sample in
upwelling regions in Figure 7(a) implies that a long chain can help the
vertical upward migration not only through fast swimming speed but also
through a turbulent upwelling water column.

Here, we answer the question of where clusters and voids form. Figure 7(b, c)
show the average non-dimensional vertical fluid velocity sampled by clusters
$\left\langle u_{z}\right\rangle_{c,t}/u_{\eta}$ and voids $\left\langle
u_{z}\right\rangle_{v,t}/u_{\eta}$ respectively. The curves in Figure 7(b) for
clusters roughly show a similar trend with the overall case, reflecting that
the cluster locations qualitatively agree with locations that overall
particles preferentially visit. The opposite trend is observed for voids in
Figure 7(c), reflecting that they form at regions where the overall particles
tend to avoid. From a quantitative perspective, the mean velocity sampled by
the clusters is larger than that of overall particles for most particle
parameters (namely $\Phi$ and $\Psi$). This demonstrates that the clusters
preferentially form at higher vertical velocity flow regions compared with the
overall preferential sampling. Particles that preferentially sample in
upwelling regions or downwelling regions, respectively, form clusters in
upflows or downflows with more extreme vertical velocity. The voids also tend
to form at higher velocity vertical flows compared with the overall particles,
though the vertical flow is in the opposite direction to that of the overall
particles.

###  III.3 Voronoï Lagrangian autocorrelation

Figure 8: Lagrangian Voronoï analysis for fast ($\Phi=16$), moderate
($\Phi=6.3$), slow swimmers ($\Phi=1.0$) and non-active particles (dashed
line). (a) Temporal autocorrelation functions of Voronoï volumes at $\Psi=1$.
(b) Decorrelation time.

In this section, we study the temporal evolution of the Voronoï volumes.
Figure 8(a) shows the temporal autocorrelation function $C_{V}(\tau)$ of the
associated Voronoï volumes for three types of swimmers at a fixed $\Psi=1$,
which is defined in (10) as the correlation coefficient between the Voronoï
volumes at two different time as a function of the time lag $\tau$:

$C_{V}(\tau)=\frac{\left\langle
V^{\prime}(t)V^{\prime}(t+\tau)\right\rangle_{o}}{\left\\{\left\langle
V^{\prime}(t)^{2}\right\rangle_{o}\left\langle
V^{\prime}(t+\tau)^{2}\right\rangle_{o}\right\\}^{1/2}},$ (10)

where $\langle\ldots\rangle_{o}$ represents the average taken over for the
overall particles at a specified time step. Note that the mean value of the
Voronoï volumes at the corresponding time have been subtracted, so the
fluctuating component of the Voronoï volumes $V^{\prime}=V-\left\langle
V\right\rangle_{o}$ is used for the calculation of the temporal
autocorrelation function. The results for non-active particles are also shown
for comparison. Compared with the non-active particles, the $C_{V}(\tau)$ of
the slow and moderate swimmers decreases relatively slowly, reflecting the
trapping of particles by turbulence prolonged the time needed for the
decorrelation of the Voronoï volumes. We observe that, for fast swimmers,
$C_{V}$ decorrelates much slower as compared to the slow swimmers, which is
related to the result that fast swimmers have a high possibility to accumulate
as shown in Figure 4.

Figure 8(b) shows the decorrelation time scales of the Voronoï volume
($\tau_{V}$) as a function of the stability number ($\Psi$) for fast,
moderate, slow swimmers and non-active particles. The decorrelation time for
Voronoï volume $\tau_{V}$ is defined to be the time when the autocorrelation
function decreased to 0.5, i.e. $C_{V}(\tau_{V})=0.5$. We remark that the
spatial distributions of non-active particles are independent of $\Psi$ since
they simply follow the fluid velocity at the particle position. We observe
that $\tau_{V}$ for all three types of swimmers shows non-monotonic dependence
on $\Psi$, with a peak around $\Psi=1\sim 2.5$. We remark the decorrelation
process is slow at $\Psi=16$, with $\tau_{V}$ approximating the peak value.
The stable upward swimming direction that helps maintain the spatial
distribution may partially explain this large $\tau_{V}$ for particles at
large $\Psi$. Using the results of non-active particles as a baseline, we can
study the swimming velocity effects on the decorrelation time. The
decorrelation time ($\tau_{V}$) for fast swimmers is the largest among the
three types of swimmers at large $\Psi$, while is the smallest at small
$\Psi$. In the gravity torque dominated region ($\Psi\gg 1$), the uniformly
stable upward swimming direction combined with a larger swimming velocity
$\Phi$ (relatively smaller fluid velocity effects) results in a more stable
spatial distribution and increasing decorrelation time. In contrast, in the
turbulence torque dominated region ($\Psi\ll 1$), the random swimming
direction combined with a larger swimming velocity $\Phi$ facilitates the
change of distance between particles and accelerates the decorrelation
process.

Figure 9: Alignments of particles with the Lagrangian stretching direction
$\bm{e}_{L1}$ as a function of the stability number $\Psi$ for slow, moderate,
and fast swimmers at $\alpha=20$. Particles with greater swimming velocity
exhibit stronger alignment when $\Psi\gg 1$, whereas the opposite trend is
observed when $\Psi\ll 1$.

###  III.4 Alignment with Lagrangian fluid stretching

The quantitative measurement of this Lagrangian alignment is obtained via the
ensemble and temporal average of the cosine of the angle between the particle
swimming direction and the Lagrangian stretching direction,
$\left\langle|\bm{e}_{L1}\cdot\bm{p}|\right\rangle$.
$\left\langle|\bm{e}_{L1}\cdot\bm{p}|\right\rangle$ takes the value 1 in the
case of a perfect parallel alignment along the $\bm{e}_{L1}$ direction, 0 in
the perfect perpendicular case, and 0.5 if the particles point randomly with
respect to $\bm{e}_{L1}$.

We look at the mean Lagrangian alignment of particles of different stability
numbers $\Psi$ for three different swimming numbers $\Phi$. Figure 9 shows the
measurements for slow ($\Phi=1$), moderate ($\Phi=6.3$) and fast swimmers
($\Phi=16$). Slow swimmers with weak gyrotaxis ($\Psi\ll 1$) approach perfect
alignment with the Lagrangian stretching direction, with
$\left\langle|\bm{e}_{L1}\cdot\bm{p}|\right\rangle$ being close to 1. This
strong alignment is consistent with the results of passive rods [25],
reflecting the alignment behavior of slow swimmers shows only small deviations
from that of non-active particles. In contrast, when $\Psi\gg 1$,
$\left\langle|\bm{e}_{L1}\cdot\bm{p}|\right\rangle$ of slow swimmers reaches
0.5, implying that the particle swimming direction and the Lagrangian
stretching direction point randomly relative to each other. This is consistent
with our expectations since the influence of fluid stretching is almost
negligible compared with the dominated gravity when $\Psi\gg 1$. We now
consider the effect of the swimming number $\Phi$ on the Lagrangian alignment,
using the alignment of slow swimmers as a baseline. The dependence of
alignment on $\Phi$ shows opposite trends in the turbulence ($\Psi\ll 1$) and
the gravity torque ($\Psi\gg 1$) dominated regions. The alignment decreases
with the swimming velocity in the small $\Psi$ regime, whereas it becomes
stronger in the large $\Psi$ regime as $\Phi$ increases.

In order to make sense of the dependence on $\Phi$ when $\Psi\ll 1$, we
measure the principal axes of length $l_{i}=\sqrt{\Lambda_{i}}(i=1,2,3)$ of an
ellipsoid deformed from a sphere due to Lagrangian stretching. We
quantitatively measure the Lagrangian stretching strength via the p.d.f.s of
the geometrical ratio $l_{1}/l_{2}$ for three types of swimmers, see Figure
10(a). An important feature of the p.d.f. of the slow swimmer is that the
highest probability occurs at the largest $l_{1}/l_{2}$ among three swimmers.
With increasing swimming velocity, the geometrical ratio has a higher
probability to be in the small range. It has been found that the alignment
between $\bm{p}$ and $\bm{e}_{L1}$ increases monotonically with the
geometrical aspect ratio $l_{1}/l_{2}$ from a statistical perspective[25].
Therefore, the geometrical aspect ratio $l_{1}/l_{2}$, to some extent, can
reflect the strength of this Lagrangian stretching. To understand the weakened
Lagrangian stretching due to swimming, we conjecture the following physical
image. The Lagrangian velocity gradient experienced by particles with a large
swimming velocity decorrelates faster than that of non-active particles. As a
result, the typical time scale for the process of Lagrangian stretching
becomes shorter. Therefore, the alignment with the Lagrangian stretching
direction decreases as the swimming velocity increases at $\Psi\ll 1$, as
shown in Figure 9.

Figure 10: (a) The p.d.f.s of the geometrical ratios $l_{1}/l_{2}$ for slow,
moderate, and fast swimmers at $\Psi=0.01$. The geometrical ratio p.d.f.s are
qualitatively similar for $\Psi$ in the range of [0.01,16], thus only the
results at $\Psi=0.01$ are plotted here. For slow swimmers, the geometrical
ratio tends to be more extreme, reflecting a stronger Lagrangian stretching
effect. (b) The mean of the absolute value of the vertical component of the
Lagrangian stretching direction.

The analysis above is for $\Psi\ll 1$. As we have examined, three swimmers
give qualitatively similar geometrical ratio ($l_{1}/l_{2}$) p.d.f.s for
$\Psi\in[0.01,16]$. These p.d.f.s for other $\Psi$ are not shown in Figure
10(a), since these curves nearly collapse with the shown one. Therefore, the
conclusion that fast swimmers experience weaker Lagrangian stretching is
robust in the large $\Psi$ range. However, we see in Figure 9 that fast
swimmers under strong gyrotaxis ($\Psi\gg 1$) show better alignment than
moderate and slow swimmers, although they have encountered weaker fluid
stretching. This deviates from the usual expectations. Figure 10(b) shows the
mean of the vertical component of $\bm{e}_{L1}$,
$\left\langle|\bm{e}_{L1}\cdot\bm{z}|\right\rangle$, as a function of $\Psi$
for different $\Phi$. It turns out that the Lagrangian stretching direction
can be biased by gravity. The Lagrangian stretching direction of fast swimmers
shows the largest extent to align with the vertical direction. Under the
strong gyrotaxis effect, the swimming direction of particles also approaches
the vertically upward direction $\bm{z}$. Thus, the counterintuitive strong
alignment for fast swimmers at large $\Psi$ can be understood. However, the
bias of $\bm{e}_{L1}$ in the vertical direction is not revealed. The
Lagrangian stretching direction is obtained by integrating the velocity
gradient $\bm{A}(\bm{x},t)$ (8), which does not include the gravity effect
explicitly. One possible reason could be the biased sampling of particle
positions $\bm{x}$ in the flow, which was found to be upwelling or downwelling
regions [20, 22]. Even if there is no preferential accumulation when $\Psi\gg
1$, the uniformly upward swimming direction of particles may prolong their
time staying in structures with large vertical velocity. In these regions the
flow may tend to stretch particles in the vertical direction, thus causing the
bias of $\bm{e}_{L1}$ in the vertical direction. More studies are needed to
fully understand this issue.

##  IV Conclusions

In this paper, we consider simulations of the dynamics of elongated gyrotactic
swimmers in turbulent flows. We carry out the study of clustering of
gyrotactic swimmers, using the three-dimensional Voronoï analysis. Three
different types of swimmers (fast, moderate, and slow) of different $\Psi$ are
covered. At small and intermediate $\Psi$, a higher probability of finding
depleted regions (large Voronoï volumes) and concentrated regions (small
Voronoï volumes) is identified for active particles, implying the clustering
behavior. We use the standard deviations of the normalized Voronoï volumes to
measure the clustering extent quantitatively. At intermediate $\Psi$, fast
swimmers show stronger clustering than the other two types of swimmers, which
is consistent with the trend based on the radial distribution function
analysis. While the gyrotaxis-dominated ($\Psi\gg 1$) swimmers exhibit no
clustering, non-gyrotactic particles ($\Psi\ll 1$) are amenable to concentrate
preferentially, with a non-monotonic dependence on $\Phi$.

Next, we study the connections between the Voronoï volume and the local flow
properties from the spatial and time perspective. By comparing the p.d.f. of
the overall particles and that of particles belonging to clusters and voids,
the locations that overall, aggregated, and sparse particles preferentially
sample are further distinguished. Clusters and voids preferentially form at
locations with more extreme vertical velocity compared with the overall
preferential sampling. The preferential vertical flow regions of clusters and
voids are roughly the opposite. From the Lagrangian autocorrelation of the
Voronoï volumes, the lifetime of the clusters is examined. The clustering of
fast swimmers lasts longer than that of moderate and slow swimmers at large
$\Psi$, but the opposite dependence on swimming velocity emerges at small
$\Psi$. In the future, more insights may be gained by following individual
clusters of gyrotactic swimmers, focusing on how they are formed and destroyed
by turbulent fluctuations.

Finally, we study their Lagrangian alignment, focusing on the Lagrangian
stretching directions. We have demonstrated how the motility can weaken this
alignment, by visualizing stretching as the process of deforming a sphere into
an ellipsoid. With faster velocity, the deformation is reduced, reflecting a
weakened Lagrangian stretching. However, when gravity effects become more
dominated, the Lagrangian alignment shows a counterintuitive increasing trend
with a particle moving velocity. We found, under this case, the Lagrangian
stretching direction can also be biased by gravity. Gravity biases both the
particle swimming direction and the Lagrangian stretching direction in the
vertical direction. Therefore, the remarkable alignment trend with the
particle moving velocity at large $\Psi$ can be understood.

Acknowledgments We thank Enrico Calzavarini and Z. Wang for their comments on
the manuscript. This work is financially supported by the Natural Science
Foundation of China under Grant No. 11988102, 91852202, the National Training
Program of Innovation and Entrepreneurship of China for Undergraduates
(202110003017), and Tencent Foundation through the XPLORER PRIZE.

## References

  * Stocker [2012] R. Stocker, “Marine microbes see a sea of gradients,” Science 338, 628–633 (2012).
  * Qiu  _et al._ [2022] J. Qiu, N. Mousavi, K. Gustavsson, C. Xu, B. Mehlig, and L. Zhao, “Navigation of micro-swimmers in steady flow: the importance of symmetries,” J. Fluid Mech. 932, 1–21 (2022).
  * Chakraborty  _et al._ [2018] S. Chakraborty, F. Ivancic, M. Solovchuk, and T. Sheu, “Stability and dynamics of a chemotaxis system with deformed free-surface in a shallow chamber,” Phys. Fluids 30 (2018).
  * Kessler [1985] J. O. Kessler, “Hydrodynamic focusing of motile algal cells,” Nature 313, 218–220 (1985).
  * Panda [2020] M. K. Panda, “Effects of anisotropic scattering on the onset of phototactic bioconvection with diffuse and collimated irradiation,” Phys. Fluids 32 (2020).
  * Dervaux, Capellazzi Resta, and Brunet [2017] J. Dervaux, M. Capellazzi Resta, and P. Brunet, “Light-controlled flows in active fluids,” Nat. Phys. 13, 306–312 (2017).
  * Pedley and Kessler [1992] T. J. Pedley and J. O. Kessler, “Hydrodynamic phenomena in suspensions of swimming microorganisms,” Annu. Rev. Fluid Mech. 24, 313–358 (1992).
  * Smayda [2010] T. Smayda, “Adaptations and selection of harmful and other dinoflagellate species in upwelling systems. 2\. motility and migratory behaviour,” Prog. Oceanogr. 85, 71–91 (2010).
  * Wheeler  _et al._ [2019] J. D. Wheeler, E. Secchi, R. Rusconi, and R. Stocker, “Not just going with the flow: the effects of fluid flow on bacteria and plankton,” Annu. Rev. Cell Dev. Biol. 35, 213–237 (2019).
  * Jiang and Chen [2020] W. Jiang and G. Chen, “Dispersion of gyrotactic micro-organisms in pipe flows,” J. Fluid Mech. 889 (2020).
  * Durham, Kessler, and Stocker [2009] W. M. Durham, J. O. Kessler, and R. Stocker, “Disruption of vertical motility by shear triggers formation of thin phytoplankton layers,” Science 323, 1067–1070 (2009).
  * Cencini, Boffetta, and Borgnino [2019] M. Cencini, G. Boffetta, and F. Borgnino, M.and De L., “Gyrotactic phytoplankton in laminar and turbulent flows: a dynamical systems approach,” Eur. Phys. J. E 42, 1–15 (2019).
  * Wang  _et al._ [2021] B. Wang, W. Jiang, G. Chen, L. Tao, and Z. Li, “Vertical distribution and longitudinal dispersion of gyrotactic microorganisms in a horizontal plane Poiseuille flow,” Phys. Rev. Fluids 6, 1–23 (2021).
  * Wang, Jiang, and Chen [2022] B. Wang, W. Jiang, and G. Chen, “Cross-channel distribution and streamwise dispersion of micro-swimmers in a vertical channel flow: A study on the effects of shear, particle shape, and convective inertial torque,” Phys. Fluids 34, 011904 (2022).
  * Durham  _et al._ [2013] W. M. Durham, E. Climent, M. Barry, F. De Lillo, G. Boffetta, M. Cencini, and R. Stocker, “Turbulence drives microscale patches of motile phytoplankton,” Nat. Commun. 4, 1–7 (2013).
  * De Lillo  _et al._ [2014] F. De Lillo, M. Cencini, W. M. Durham, M. Barry, R. Stocker, E. Climent, and G. Boffetta, “Turbulent fluid acceleration generates clusters of gyrotactic microorganisms,” Phys. Rev. Lett. 112, 044502 (2014).
  * Fouxon and Leshansky [2015] I. Fouxon and A. Leshansky, “Phytoplankton’s motion in turbulent ocean,” Phys. Rev. E 92, 013017 (2015).
  * Zhan  _et al._ [2014] C. Zhan, G. Sardina, E. Lushi, and L. Brandt, “Accumulation of motile elongated micro-organisms in turbulence,” J. Fluid Mech. 739, 22–36 (2014).
  * Pujara, Koehl, and Variano [2018] N. Pujara, M. A. R. Koehl, and E. A. Variano, “Rotations and accumulation of ellipsoidal microswimmers in isotropic turbulence,” J. Fluid Mech. 838, 356 (2018).
  * Gustavsson  _et al._ [2016] K. Gustavsson, F. Berglund, P. R. Jonsson, and B. Mehlig, “Preferential sampling and small-scale clustering of gyrotactic microswimmers in turbulence,” Phys. Rev. Lett. 116, 108104 (2016).
  * Borgnino  _et al._ [2018] M. Borgnino, G. Boffetta, F. De Lillo, and M. Cencini, “Gyrotactic swimmers in turbulence: shape effects and role of the large-scale flow,” J. Fluid Mech. 856, R1 (2018).
  * Lovecchio  _et al._ [2019] S. Lovecchio, E. Climent, R. Stocker, and W. M. Durham, “Chain formation can enhance the vertical migration of phytoplankton through turbulence,” Sci. Adv. 5, eaaw7879 (2019).
  * Seymour  _et al._ [2011] J. R. Seymour, M. Luhar, W. M. Durham, J. G. Mitchell, A. Macke, and R. Stocker, “Microbial alignment in flow changes ocean light climate,” Proc. Natl. Acad. Sci. USA 108, 3860–3864 (2011).
  * Borgnino  _et al._ [2019] M. Borgnino, K. Gustavsson, F. De Lillo, G. Boffetta, M. Cencini, and B. Mehlig, “Alignment of nonspherical active particles in chaotic flows,” Phys. Rev. Lett. 123, 138003 (2019).
  * Ni, Ouellette, and Voth [2014] R. Ni, N. T. Ouellette, and G. A. Voth, “Alignment of vorticity and rods with lagrangian fluid stretching in turbulence,” J. Fluid Mech. 743 (2014).
  * Hejazi, Mehlig, and Voth [2017] B. Hejazi, B. Mehlig, and G. A. Voth, “Emergent scar lines in chaotic advection of passive directors,” Phys. Rev. Fluids. 2, 124501 (2017).
  * Cui  _et al._ [2020] Z. Cui, A. Dubey, L. Zhao, and B. Mehlig, “Alignment statistics of rods with the lagrangian stretching direction in a channel flow,” J. Fluid Mech. 901, A16 (2020).
  * Perlekar  _et al._ [2012] P. Perlekar, L. Biferale, M. Sbragaglia, S. Srivastava, and F. Toschi, “Droplet size distribution in homogeneous isotropic turbulence,” Phys. Fluids 24 (2012).
  * Jeffery [1922] G. B. Jeffery, “The motion of ellipsoidal particles immersed in a viscous fluid,” Proc. R. Soc. Lond. A 102, 161–179 (1922).
  * Karp-Boss, Azevedo, and Boss [2007] L. Karp-Boss, L. Azevedo, and E. Boss, “Lisst-100 measurements of phytoplankton size distribution: evaluation of the effects of cell shape,” Limnol. Oceanogr. Methods 5, 396–406 (2007).
  * Bechinger  _et al._ [2016] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, “Active particles in complex and crowded environments,” Rev. Mod. Phys. 88 (2016).
  * Benfenati and Coscia [2013] A. Benfenati and V. Coscia, “Nonlinear microscale interactions in the kinetic theory of active particles,” Appl. Math. Lett. 26, 979–983 (2013).
  * Sun [2006] R. Sun, “The coagulation and re-separation process of particles in a slow viscous flow,” J. Hydrodyn. 18, 107–110 (2006).
  * Tagawa  _et al._ [2012] Y. Tagawa, J. M. Mercado, V. N. Prakash, E. Calzavarini, C. Sun, and D. Lohse, “Three-dimensional Lagrangian Voronoï analysis for clustering of particles and bubbles in turbulence,” J. Fluid Mech. 693, 201–215 (2012).
  * Monchaux, Bourgoin, and Cartellier [2012] R. Monchaux, M. Bourgoin, and A. Cartellier, “Analyzing preferential concentration and clustering of inertial particles in turbulence,” Int. J. Multiph. Flow . 40, 1–18 (2012).
  * Ferenc and Néda [2007] J. S. Ferenc and Z. Néda, “On the size distribution of poisson voronoi cells,” Physica A 385, 518–526 (2007).
  * Chadwick [1999] P. Chadwick, “Continuum mechanics: Concise theory and problems. Dover publications,” (1999).
  * Parsa  _et al._ [2011] S. Parsa, J. S. Guasto, M. Kishore, N. T. Ouellette, J. P. Gollub, and G. A. Voth, “Rotation and alignment of rods in two-dimensional chaotic flow,” Phys. Fluids 23, 043302 (2011).
  * Wilkinson, Bezuglyy, and Mehlig [2011] M. Wilkinson, V. Bezuglyy, and B. Mehlig, “Emergent order in rheoscopic swirls,” J. Fluid Mech. 667, 158–187 (2011).
  * Calzavarini [2019] E. Calzavarini, “Eulerian–lagrangian fluid dynamics platform: The ch4-project,” Softw. Imp. 1, 100002 (2019).
  * Jiang, Calzavarini, and Sun [2020] L. Jiang, E. Calzavarini, and C. Sun, “Rotation of anisotropic particles in Rayleigh–Bénard turbulence,” J. Fluid Mech. 901 (2020).
  * Calzavarini, Jiang, and Sun [2020] E. Calzavarini, L. Jiang, and C. Sun, “Anisotropic particles in two-dimensional convective turbulence,” Phys. Fluids 32 (2020).
  * Jiang  _et al._ [2021] L. Jiang, C. Wang, S. Liu, C. Sun, and E. Calzavarini, “Rotational dynamics of bottom-heavy rods in turbulence from experiments and numerical simulations,” Theor. App. Mech. Lett. , 100227 (2021).
  * Monchaux, Bourgoin, and Cartellier [2010] R. Monchaux, M. Bourgoin, and A. Cartellier, “Preferential concentration of heavy particles: A voronoï analysis,” Phys. Fluids 22, 103304 (2010).
  * Jong Gyu Park  _et al._ [2001] Jong Gyu Park, Min Kyung Jeong, Jin Ae Lee, K. J. Cho, and O. S. Kwon, “Diurnal vertical migration of a harmful dinoflagellate, Cochlodinium polykrikoides (Dinophyceae), during a red tide in coastal waters of Namhae Island, Korea,” Phycologia 40, 292–297 (2001).