Preferential sampling and small-scale clustering of gyrotactic microswimmers in turbulence

0.45 MB·English

by  K. Gustavsson

Preferential sampling and small-scale clustering of gyrotactic microswimmers in
turbulence
∗
K. Gustavsson1,2), F. Berglund1), P.R. Jonsson3), and B. Mehlig1)
1)Department of Physics, Gothenburg University, SE-41296 Gothenburg, Sweden
2)Department of Physics and INFN, University of Rome ’Tor Vergata’, 00133 Rome, Italy and
3)Department of Biological and Environmental Sciences - Tja¨rn¨o, SE-45296 Str¨omstad,Sweden
Recentstudiesshowthatsphericalmotilemicro-organismsinturbulencesubjecttogravitational
torques gather in down-welling regions of the turbulent flow. By analysing a statistical model we
analyticallycomputehowshapeaffectsthedynamics,preferentialsampling,andsmall-scalespatial
clustering. Wefindthatoblongorganismsmayspendmoretimeinup-wellingregionsoftheflow,and
6
1 that all organisms are biased to regions of positive fluid-velocity gradients in the upward direction.
0 We analyse small-scale spatial clustering and find that oblong particles may either cluster more or
2 less than spherical ones, depending on the strength of the gravitational torques.

y
PACSnumbers: 05.40.-a,47.63.Gd,47.27.-i,92.20.jf
a
M
Patchiness in suspensions of micro-organisms is fre- ganisms respond not only to turbulent vorticity but also
9
quentlyobservedonarangeofspatialscales. Theunder- to turbulent strain [12–15]. This causes passive rods
2
lying mechanisms differ, depending on the properties of to exhibit intricate orientational patterns on the sur-

] the micro-organisms, and upon the spatial scale. Patch- face of turbulent and other complex flows [16–18]. Also,
n
iness can be caused by density stratification and ver- shape strongly affects the trajectories of active particles
y
tical shears [1], by predator-prey cycles, or by interac- inmodelflows[19–21],andrecentDNSindicatethatpro-
d
- tionsbetweentheorganismsandwater-columngradients late gyrotactic organisms cluster less than spherical ones
u
–inlight,chemistry,turbulence,andinhydrostaticpres- when gyrotaxis is strong [22]. Second, where do the or-
l
f sure[2]. Patchinessisimportantbecausemanybiological ganismsgointurbulence? Aretherecircumstanceswhere
.
s processes (mating, feeding, predation) rely on individual the organisms may not gather in down-welling regions,
c
encounters [3], and the encounter rate is strongly influ- or where other mechanisms of preferential sampling may
i
s enced by small-scale number-density fluctuations. apply? Third, the fact that orgamisms tend to gather in
y
certain regions of the flow (preferential sampling) does
h Gravitaxismaycausesuchinhomogeneitiesinthespa-
p tial distribution of motile micro-organisms. Density- or not explain which mechanisms actually cause them to
[ get in contact. To determine these one must follow the
drag-asymmetries of the body give rise to torques affect-

dynamics of two organisms that are initially very close
3 ing the swimming direction [4–6]. When the effects of
together, and determine whether they tend to approach
v gyrotactic torques and fluid-velocity gradients balance,
6 inhomogeneities may form in the spatial distribution, as further or move apart. We refer to the resulting small-
8 scalespatialfluctuationsinthenumberdensityas‘small-
shown by the micro-alga Chlamydomonas nivalis swim-
3 scale clustering’.
ming up against a down-welling pipe flow. The micro-
2
Statistical model. To answer these questions we use a
0 algae gather in the centre of the pipe where the down-
simplified model [7, 11, 22] for the translation and ro-
. welling velocity is largest [7]. Gyrotaxis may trap motile
1
tation of small axisymmetric active particles subject to
0 organismsinmacroscopicsheargradients[8,9],andfluc-
turbulence and gyrotaxis:
5 tuating vorticity may cause patchiness [10]. This is con-
1 firmed by recent direct numerical simulations (DNS) of r˙ v =u(r,t)+v n and n˙ =ω(r,t) n. (1)
v: motile, spherical micro-organisms in turbulence [11] re- ≡ s ∧
i vealingthattheorganismsaremorelikelytobefoundin Dots denote derivatives w.r.t. time t, r is the particle
X
down-welling regions of the turbulent flow, they ‘prefer- position, and u is the flow velocity. Each particle swims
ar entially sample’ such regions. withconstantspeedvs inthedirectionnofitssymmetry
axis (n =1). The angular velocity of the particle is
These results raise three fundamental questions that | |
we address and answer in this Letter. First, how does ω(r,t)=(gˆ n)/(2B)+Ω(r,t)+Λn S(r,t)n . (2)
∧ ∧
shape affect the dynamics in turbulence of motile micro-
The first term on the r.h.s describes g(cid:2)yrotaxis.(cid:3) The
organisms subject to gyrotaxis? In Ref. [11] the organ-
unit vector gˆ points in the direction e of gravity, and
isms were assumed to be spherical. Non-spherical or- − z
B is the reorientation time [7, 11]. It depends on the
mass distribution within the particle, and on its shape
through hydrodynamic resistance. The other terms on
the r.h.s. of Eq. (2) represent the effect of the turbu-
∗Version accepted for publication (postprint) in Phys. Rev. Lett.
116,108104(2016) lent velocity gradients upon the particle orientation [12]:
2
Ω=( u)/2, and S is the symmetric part of the ma- terms of np for p < q. The lowest-order solution in Ku
∇∧
trixAoffluid-velocitygradients. TheparameterΛchar- is just n0 = gˆ. This yields a lowest-order deterministic
−
acterises particle shape: Λ = 0 for spheres, and Λ = 1 approximation for the particle position at time t:
0
for infinitely thin rods. Eq. (2) disregards turbulent ac-
r (t)=r Φgˆt . (4)
celerations. In most marine conditions this is an excel- d0et 0 00− 0
lent approximation [23]. We model the dissipative range
Third, we expand Eqs. (1,2) in terms of deviations from
of turbulence by incompressible, homogeneous, isotropic
thedeterministictrajectoryδr (t) r (t) r (t)pre-
Gaussian random functions with typical length- , time- 0 0 ≡ 0 0 − d0et 0
cisely as described in Ref. [24]. In the fourth and final
, and speed-scales η, τ, u [24]. This neglects inertial-
0 stepweaverageoverthefluid-velocityfluctuationsinthe
range properties which may become important for parti-
statisticalmodel. IntheremainderofthisLetterwesum-
clesthatarelargerthantheKolmogorovlength[25]. We
marise the results obtained in this way.
notethatthedissipative-rangeturbulentfluctuationsare
Preferential sampling. Consider the steady-state aver-
universal [26], but they are not Gaussian. We comment
ages of the z-component u of the fluid velocity and of
z
on this difference between turbulence and the statistical
itsgradient,A ,bothevaluatedattheparticleposition.
zz
model below.
Analyticalresultsfortheseaveragesarederivedtolowest
There are four dimensionless parameters: the shape
orderinKuintheSupplementalMaterial[31],Eqs.(S15)
parameter Λ, the reorientation time Ψ = B/τ, the
and (S16). These expressions are plotted in Fig. 1. Here
swimming speed Φ = v τ/η, and the Kubo number
s we quote only limiting results. For small Φ we have
Ku = u τ/η. We vary the parameters independently,
0
keeping B constant as Λ is changed. Ku is a dimension- A η KuΦ2 d(1−Λ)+2(Λ+2) Ψ(4Ψ+1), (5a)
less measure of the correlation time of the flow. h zzi∞u0∼ d (2Ψ+1)2
Ourchoiceofthedimensionlessparametersisdictated u d(1 Λ)+2 Ψ
z
by the method (explained below). DNS employ differ- h i∞ KuΦ − , (5b)
u ∼− d 2Ψ+1
0
ent de-dimensionalisations [11]: B by the Kolmogorov
time τK 1/ Tr AAT , and vs by the corresponding d is the spatial dimension. For large Φ we find
≡ h i∞
Kolmogorov speed u . Our dimensionless parameters
p K η Ku d+1 π
tΦranslate Φto/tKhous.e uWseedeixnptehcetDthNaSt atsheΨDstNaStis∼tiKcaul-Ψmoadnedl hAzzi∞u0 ∼ Φ 2d (1−Λ)r2 , (6a)
DNS
∼ u Ku (d(Λ 1)+2Λ)
results become independent of Ku at large Ku and qual- h zi∞ − . (6b)
u ∼ Φ 2d
itatively agree with DNS results [24]. 0
Typical values of B, vs are given in Ref. [5]: B 1– What can we learn from these analytical results?
∼
5s and vs 0.1–1mm/s. Typical ocean dissipation Eqs. (5a) and (6a) show that the particles collect in
∼
rates are ε 1–102mm2/s3 for surface water [27], giv- the sinks of the transversal flow-velocity field, Tr A
ing Kolmogo∼rov times, lengths, and speeds in the range A < 0. This is because gyrotaxis breaks up⊥-dow≡n
zz
τK 0.1–1s, η 0.3–1mm, and uK 1–3mm/s. These −symmetry: whenΨissmalltheparticlesswimessentially
∼ ∼ ∼
estimates yield ΨDNS ∼ 1–50 and ΦDNS ∼ 0.03–1. In upwards (in the ez-direction), and gather in transversal
Ref. [28] smaller dissipation rates, ε 10−4mm2/s3, are sinks irrespective of their shape and swimming speed.
∼
quoted for the very deep sea. This extends the ranges to Simulations (Fig. 1a) confirm the theory.
ΨDNS 0.01–50 and ΦDNS 0.03–10. Motivated by Kessler’s study in pipe flows [7] the au-
∼ ∼
Method. Eqs.(1,2)canbesolvedbyiterativelyrefining thors of Ref. [11] concluded that spherical particles pref-
approximations for the path a particle takes through the erentiallysampledown-wellingregionsalsointurbulence.
flow [24, 29, 30]. This results in expansions of steady- This is not in contradiction with the result discussed
state averages in powers of Ku and allows to determine above because particles may preferentially sample dif-
how the remaining parameters (Φ, Ψ, and Λ) affect pref- ferent observables. In fact Eqs. (5b) and (6b) explain
erential sampling and small-scale clustering. The details that spherical particles are biased towards down-welling
of this calculation are given in the Supplemental Ma- regions (as observed in DNS [11]), in addition to sinks in
terial [31]. Here we outline the essential steps. First, the transversal flow. But (6b) also shows that elongated
to consistently track the orders in the expansion we de- particles[Λ>d/(d+2)]preferentiallysampleup-welling
dimensionalise t0=t/τ,r0=r/η,u0=u/u0. Second, we regions for large enough Φ. This is seen in Fig. 1c which
expand the dynamics of the vector n in powers of Ku: shows u (Eq. (S15) in the Supplemental Material).
z
h i∞
for Ku = 0.1 and d = 2 as a function of Φ. Also shown
n(t)= ∞ n (t)Kuq . (3) are results of statistical-model simulations, in excellent
0 q 0
agreement with theory. Fig. 1d shows that the same
q=0
X conclusions hold in three spatial dimensions for Ku=1.
Inserting this ansatz into (1,2) and identifying terms of RodssampleupwellingregionswhenΦislargerthan(ap-
order Kuq yields equations for n that can be solved in proximately) unity.
q
3
0.03 0.1
u0 a b consider the dependence of uz on the reorientation
/ h i∞
η time B. Fig. 1f shows that the statistical-model result
∞ 0.02
i becomes independent of Ku for large Ku, and that it
A 0.05
⊥ 0.01 reproduces the DNS results fairly well, it explains the
Tr B/τ -dependence of u of the DNS results up to a
−h 0 0 prefaKctoroforderunityh. zTih∞isfactorisduetothefactthat
10−2 100 102 10−2 100 102
fully-developedturbulentvelocityfluctuationsinthedis-
c d
u0 0.02 0.1 sipative range differ from those in the statistical model:
/
∞ they are not Gaussian, more persistent, and the proba-
i0.01 0.05
z bility of straining regions to occur is higher [24].
u
−h 0 0 Small-scale clustering. Which mechanisms cause two
particles caught in the same flow region to actually col-
−0.0110−2 100 102−0.0150−2 100 102 lide? This is a two-particle problem, only indirectly
Φ Φ
related to preferential sampling. Fluctuations in the
e f separations between nearby particles are determined by
vs 0.3
/ the dynamics of the particle-velocity gradients ∂v /∂r .
∞0.05 i j
i 0.2 Small-scale clustering occurs where v < 0. We have
uz ∇·
h computed v to lowest order in Ku. The result is
− 0.1 h∇· i∞
quite lengthy [Eq. (S32) in the Supplemental Material].
0
0 For small Φ the full expression simplifies to:
10−2 100 102 10−1 100 101 102
Φ B/τK v η/u Ku(ΦΨ)2B (Λ) for Φ 1, (7)
0 d
h∇· i∞ ∼− (cid:28)
FpaIGrt.ic1l:es.PRreefseurletnstoiaflssimamulpaltiniogn.saofhTstra⊥tiAstii∞caalsmsoadmepl:lesdymbybothlse. withBd(Λ)≡[(d+2)(d+4)−2d(d+4)Λ+(4+2d+d2)Λ2 /d.
Since B (Λ) > 0, Eq. (7) implies small-scale clustering.
Eq. (5a): lines. Parameters: Ku = 0.1, Ψ = 0.1, Λ = 0 ( ), d
(cid:3)
Ψ = 1, Λ = 0 ( ), Ψ = 0.1, Λ = 1 ( ) Ψ = 1,Λ = 1 ( ). Forsphericalparticles(Λ=0)thequadraticdependence
b The same but for Ku = 1 (theory only for large Φ). c of v on ΦΨ was derived in Ref. [11] (and also in
uz as sampled by the particles for Ku=0.1. d Same but Rehf.∇[3·2])i:∞expanding Eqs. (1,2) for B τ gives
hforiK∞u = 1. e u /v as a function of Φ for Ku = 0.1. f (cid:28)
z s
(huz),i∞K/uvs=ve5rs(ush),BKi/∞uτK=fo1r0Λ( =).0HaonlldowKuma=rk1er(s s)h,oKwuD=NS2 ∇·v ∼vsB −(1+Λ)∂z2uz +(1−Λ)(∂z2uz−∆uz) . (8)
datafromFig.3dinRef.[11](atReλ =64). Alldataarefor Substituting Λ(cid:2) = 0 yields Eq. (6) of Ref. [11], and(cid:3)aver-
values of Φ from the small-Φ plateau observed in DNS [11],
agingEq.(8)alongparticlepathsresultsinEq.(7). The
andalsointhestatisticalmodel(forKu=0.1thisplateauis
factorv B in(8)correspondstoonefactorofΦΨin(7).
shown in panel e, for Ku=1 in Fig. S1 in the Supplemental s
The second factor of ΦΨ comes from averaging the ve-
Material). Panels a,c,e are for d=2, panels b,d,f for d=3.
locity derivatives in Eq. (8). We note that Tr A does
⊥
not figure in Eq. (8): preferential sampling of sinks in
the flow-velocity field perpendicular to gravity does not
It is remarkable that the shapes of the curves at Ku=
contributetosmall-scaleclustering,showingthatthetwo
0.1 are very similar to those at large Ku. This means
effects are distinct [24].
that the small-Ku theory qualitatively explains what is
Expanding the full result (S32) for large Φ gives:
observed in the statistical-model simulations at large Ku
and in DNS for spherical particles [11].
v η/u Ku ΦΨ2E (Λ) for Φ 1, (9)
0 d
In the limit of large Φ particles swim rapidly upwards h∇· i∞ ∼− (cid:29)
and experience the flow as a white-noise signal (as in withE (Λ) π/2(d+1)(d+3)(Λ 1)2/d. Forspherical
d
rapid gravitational settling [30]). This limit is universal, particles the≡ΦΨ2-dependence was−derived in Ref. [32].
particles in any homogeneous, isotropic and incompress- p
Let us now analyse the shape dependence of Eq. (7).
ibleflowshowpreferentialsamplingaccordingtoEq.(6).
The Λ-dependence of B (Λ) explains that rods (Λ = 1)
d
Thismeansthatthesmall-Kutheoryshoulddescribere-
cluster less than spheres (Λ = 0), consistent with the
sults of statistical-model simulations at Ku = 1 quanti-
DNS results reported in Ref. [22]. But when gyrotaxis
tatively for large Φ. This is confirmed by Figs. 1b,d.
isweak,spheresareessentiallyrandomlyoriented,unlike
For small Φ DNS [11] show that the average of uz is neighbouring rods that are aligned by turbulent shears.
proportional to Φ for small Φ, so that uz /Φ is con- In this limit motile rods must cluster more than spheres.
h i∞
stant. Eq. (5b) shows this behaviour, in good agreement This is demonstrated below, but it is not captured by
with simulations (Fig. 1e). Eqs. (7) and (9) which must fail for large Ψ because the
We conclude with a quantitative comparison of limit Ψ is singular, and due to the occurrence of
→ ∞
statistical-model and DNS results [11]. As an example singularitiesinthedynamicsofthegradientsofnatlarge
4
2
butfinitevaluesofΨ(SupplementalMaterial). Thefirst
Rodscluster
caveat also applies to Eqs. (5) and (6). L 100 less
Fractal dimension. DNSshow[11]thatthesmall-scale −d 1.8
spatial patterns of motile gyrotactic organisms are frac- d dL
=10−2 Rods
t[3a3l.].TWhiesamnaaylyssuebtshtaenftriaaclltyalenphatatnecrensthfeoirrfienncitoeubntuetrsrmataelsl ∆L 1.6 cmluosrteer
a b
Ku,intwodimensions. Weexpectqualitativelythesame 10−4
10−2 100 0.1 0.2 0.3 0.4
resultinthreedimensions. Thefractalpatternsarechar- Φ Ψ
acterised by ‘Lyapunov exponents’ λ and λ
1 2
FIG. 2: a Fractal dimension deficit ∆ d d for d=2,
L L
λ1 limt−1ln R(t) and λ1+λ2 limt−1ln A(t) . (10) Ku=0.1. Numericalsimulationsofstatis≡tical−modelforΨ=
≡t R(0) ≡t (0) 0.1,Λ=0( ),Ψ=1,Λ=0( ),Ψ=0.1,Λ=1( ),andΨ=
→∞ →∞ A
1,Λ = 1 ( ). Theory [Eqs. (11), (S31), and (S32)] including
These exponents quantify the expansion (contraction) the Ku4-contribution to Eq. (S31): solid lines. Asymptote
rates of the distance R(t) between two initially nearby (ΦΨ)2: dashed lines. b Numerical simulations, d for d=2∝,
L
particles, and of the area element (t) spanned by the Φ=1,Ku=1,Λ=0( ),0.2( ),0.4( ),0.6( ),0.8( ),1
A
separation vectors between three nearby particles. The ( ).
fractal Lyapunov dimension is defined by [24, 34]
dL 1 λ1/λ2, (11) often found in regions where uz is negative, explaining
≡ −
the behaviour found in DNS [11]. But our calculations
assuming λ > 0 and λ +λ < 0. When d < 2 fractal also predict that rod-like particles preferentially sample
1 1 2 L
clustering occurs. To evaluate d we use λ +λ = up-wellingregionsofhomogeneousisotropicflowssuchas
L 1 2
v ,Eq.(S32),andcomputeλ toorderKu4. Therehs∇ult· turbulence,providedthattheyswimfastenough. Third,
1
i∞
is lengthy, in the Supplemental Material [31] we quote we have analytically computed how the degree of small-
the result to order Ku2, Eq. (S31). To this order it is scale spatial clustering depends on particle shape. This
independentofΨandΛ. ForsmallvaluesofΦEq.(S31) is important because small-scale fractal clustering may
simplifies to λ τ Ku2(1 3Φ2) for d = 2. Together enhance particle-encounter rates. We find a transition
1
∼ −
with (7) this implies ∆ d d Φ2Ψ2 consistent that we predict to be observable in DNS as well: when
L L
≡ − ∼
with the results of Refs. [11, 32] for spherical particles. gyrotaxisisstrong(smallΨ)oblongparticlesclusterless
Fig.2ashowstheanalyticalresultford asafunctionof thansphericalones,whileatlargeΨtheoppositeistrue.
L
Φ. Itisingoodagreementwithnumericalsimulationsof Ourcalculationsalsoshowthatsingularitiesinthemo-
the statistical model (d = 2) for Ku = 0.1. We see that tion of nearby micro-organisms occur, much like ‘caus-
sphericalorganismsclustermorethanrods. Asexplained tics’ for heavy particles in turbulence [35–38]. We pre-
above this is expected for strong gyrotaxis. dict that such singularities must also be observed in the
Butwhentheeffectofthegravitationaltorqueissmall DNS of gyrotactic microswimmers in turbulence. It is
thenprolateorganismsclustermore: intheabsenceofgy- of interest to estimate how often the singularities occur
rotaxis, rotational symmetry ensures that active spheri- because their effect may modify the predictions of phe-
cal particles remain uniformly distributed, but rod-like nomenological models for encounter rates [39].
particles show fractal clustering. Panel b in Fig. 2 The analytical results obtained in this Letter were de-
demonstrates this cross-over. It shows d for Φ = 1, rived for small Ku (or large Φ). But we have shown
L
Ku = 1 as a function of Ψ. We arrive at qualitatively thatouranalyticalresultsandthecorrespondingmecha-
similar conclusions by numerically computing the frac- nismsqualitativelyexplainwhatisobservedinDNS,and
tal correlation dimension d . But the numerical values explain also the results of statistical-model simulations
2
found for d differ from d This shows that the spatial at large values of Ku. We find fairly good quantitative
2 L
distribution is multifractal [24]. agreementbetweenourstatistical-modelcalculationsand
Conclusions. First, our statistical-model calculations DNS results for fully developed turbulence. To achieve
explain how the dynamics of gyrotactic motile micro- even better quantitative agreement with the DNS would
organisms depends on the dimensionless parameters of require to account for the universal non-Gaussian small-
the problem: Λ (shape), Φ (swimming speed), and Ψ scale fluctuations of fully developed turbulence [26].
(reorientation time). Second, we find that the particles But the fluctuations of the unsteady ocean are nei-
tend to preferentially sample positive values of A , cor- ther fully-developed turbulent, nor are they Gaussian.
zz
responding to sinks in the transversal flow, regardless of Therefore the fact that the much simpler Gaussian sta-
shape. We predict that this must also be observed in tistical model explains the dynamics observed in DNS
DNS, it is simply a consequence of the fact that gravity of fully developed turbulence [11] shows that the ana-
breaks the symmetry of the problem. At the same time lytical theory (and the underlying mechanisms) describe
our calculations show that spherical particles are more robust behaviour, that must be taken into account in
5
the analysis of patchiness and encounter rates of motile B. Mehlig, and E. Variano, Phys. Fluids 27, 035101
micro-organisms in the ocean. (2015).
Acknowledgments. We thank the authors of Ref. [11] [16] M.Wilkinson,V.Bezuglyy,andB.Mehlig,Phys.Fluids
21, 043304 (2009).
for permission to reproduce their data in Fig. 3d. This
[17] V. Bezuglyy, B. Mehlig, and M. Wilkinson, Europhys.
work was supported by Vetenskapsr˚adet (VR), by a
Lett. 89, 34003 (2010).
Linnaeus-grant from VR and Formas, by the G¨oran
[18] M. Wilkinson, V. Bezuglyy, and B. Mehlig, J. Fluid
Gustafsson Foundation for Research in Natural Sciences Mech. 667, 158 (2011).
and Medicine, and by the grant Bottlenecks for parti- [19] C. Torney, and Z. Neufeld, Phys. Rev. Lett. 99, 78101
cle growth in turbulent aerosols from the Knut and Alice (2007).
Wallenberg Foundation, Dnr. KAW 2014.0048. K.G. [20] W. M. Durham, E. Climent, and R. Stocker, Phys. Rev.
Lett. 106, 238102 (2011).
acknowledges partial funding from the European Re-
[21] N.KhuranaandN.T.Ouellette,Phys.Fluids24,091902
search Council under the European Communitys Sev-
(2012).
enth Framework Programme, ERC Grant Agreement N.
[22] C. Zhan, G. Sardina, E. Lushi, and L. Brandt, J. Fluid
339032. Thenumericalcomputationsusedresourcespro-
Mech. 793, 22 (2014).
vided by the Centre for Scientific and Technical Com- [23] F. De Lillo, M. Cencini, W. M. Durham, M. Barry,
puting at Chalmers University of Technology in Gothen- R.Stocker,E.Climent,andG.Boffetta,Phys.Rev.Lett.
burg (Sweden) and the Swedish National Infrastructure 112, 044502 (2014).
for Computing. [24] K. Gustavsson and B. Mehlig, arxiv:1412.4374 (2014).
[25] H. P´ecseli, J. Trulsen, and Ø. Fiksen, Progress in
Oceanography 101, 14 (2012).
[26] J.Schumacher,J.Scheel,D.Krasnov,D.Donzis,V.Yal-
hot, and K. Sreenivasan, PNAS 11, 10961 (2014).
[27] H. Yamazaki and K. D. Squires, Mar. Ecol. Prog. Series
[1] P. J. S. Franks and J. S. Jaffe, J. Marine Sys. 69, 254
144, 299 (1996).
(2008).
[28] A.E. Gargett, Ann. Rev. Fluid Mech. 21, 419 (1989).
[2] C. L. Folt and C. W. Burns, Trends Ecol. Evol. 300, 14
[29] K.GustavssonandB.Mehlig,Europhys.Lett.96,60012
(1999).
(2011).
[3] T. Kiørboe, A mechanistic approach to plankton ecology
[30] K. Gustavsson, S. Vajedi, and B. Mehlig, Phys. Rev.
(Princeton University Press, New Jersey, USA, 2008).
Lett. 112, 214501 (2014).
[4] A. M. Roberts, J. Exp. Biology 53, 687 (1970).
[31] See Supplemental Material at [http://link.aps.org/ 405
[5] P. R. Jonsson, Mar. Ecol. Prog. Ser. 52, 39 (1989).
supplemental/10.1103/PhysRevLett.000.000000] for de-
[6] A. M. Roberts and F. M. Deacon, J. Fluid. Mech. 452,
tailsofthecalculations,andforaSupplementalFigure.).
405 (2002).
[32] I. Fouxon, A. M. Leshansky, Phys. Rev. E 92, 013017
[7] J. O. Kessler, Nature 313, 218 (1985).
(2015).
[8] W. M. Durham, J. O. Kessler, and R. Stocker, Science
[33] B.Andersson,K.Gustavsson,B.Mehlig,andM.Wilkin-
323, 1067 (2009).
son, Europhys. Lett. 80, 69001 (2007).
[9] F. Santamaria, De Lillo. F., M. Cencini, and Boffetta,
[34] J. Kaplan and J. A. Yorke, Springer Lecture Notes in
Phys. Fluids 26, 111901 (2014).
Mathematics 730, 204 (1979).
[10] J. G. Mitchell, A. Okubo, and J. A. Fuhrmann, Limnol-
[35] M. Wilkinson and B. Mehlig, Europhys. Lett. 71, 186
ogy and Oceanography 35, 123 (1990).
(2005).
[11] W. M. Durham, E. Climent, M. Barry, F. De Lillo,
[36] G. Falkovich, A. Fouxon, and G. Stepanov, Nature 419,
G. Boffetta, M. Cencini, and R. Stocker, Nature Comm.
151 (2002).
4, 2148 (2013).
[37] K. Gustavsson, E. Meneguz, M. Reeks, and B. Mehlig,
[12] G. B. Jeffery, Proc. R. Soc. A 102, 161 (1922).
New J. Phys. 14, 115017 (2012).
[13] S.Parsa,E.Calzavarini,F.Toschi,andG.A.Voth,Phys.
[38] K. Gustavsson and B. Mehlig, Phys. Rev. E 87, 023016
Rev. Lett. 109, 134501 (2012).
(2013).
[14] K. Gustavsson, J. Einarsson, and B. Mehlig, Phys. Rev.
[39] B.J.RothschildandT.R.Osborn,J.PlanktonRes.10,
Lett. 112, 014501 (2014).
465 (1988).
[15] M. Byron, J. Einarsson, K. Gustavsson, G. A. Voth,
Supplemental material for ‘A model for preferential sampling and small-scale
clustering of gyrotactic micro-swimmers in turbulence’
K. Gustavsson1,2), F. Berglund1), P.R. Jonsson3), and B. Mehlig1)
1)Department of Physics, Gothenburg University, SE-41296 Gothenburg, Sweden
2)Department of Physics and INFN, University of Rome ’Tor Vergata’, 00133 Rome, Italy and
3)Department of Biological and Environmental Sciences - Tja¨rn¨o, SE-45296 Str¨omstad,Sweden
In this Supplemental Material we outline the steps that are needed to derive Eqs. (5,6) in the Letter [S1]. This
Supplemental Material also describes how we computed an analytical expression for the Lyapunov dimension shown
6
in Fig. 2a of the Letter, and it contains a Supplemental Figure (Fig. S1).
1
0
2
PREFERENTIAL SAMPLING
y
a
M The results described in the Letter [S1] rest on the perturbation method described in Refs. [S2–S4]. In this Section
we summarise how we have applied this method to the equations of motion (1,2) in the Letter, in order to compute

9 theaverageofthefluid-velocityfieldanditsgradientsalongparticlepaths. Theseaveragesdeterminetowhichextent
2
the particle preferentially sample certain fluid-velocity configurations.

In the dimensionless variables t =t/τ,r =r/η,u =u/u ,ω =ωτ the equations of motion read:
] 0 0 0 0 0
n
y dr0
=Kuu +Φn (S1)
d dt 0
0
- dn
u =ω n (S2)
0
l dt ∧
f 0
s. ω0 =−(n∧gˆ)/(2Ψ)+KuΩ0(r0(t0),t0)+KuΛn∧ S0(r0(t0),t0)n . (S3)
c
i The dimensionless parameters Ku = u0τ/η, Ψ = B/τ, and Φ = vsτ/η are(cid:2)chosen so tha(cid:3)t all terms containing the
s
fluid velocity u or its gradients are proportional to Ku. This implies as we shall see that a perturbation expansion
y
h about the deterministic (u=0) solution is equivalent to a perturbation expansion in the Kubo number Ku.
p To expand Eqs. (S1) and (S2) in the Kubo number we first require a Ku-expansion of the orientation vector n:
[

∞
3 n(t)= n (t)Kuq, (S4)
0 q 0
v
q=0
6 X
8 where n (t) are time-dependent expansion coefficients. Inserting this ansatz into Eq. (S2) and identifying terms of
q 0
3 order Kuq we find that the expansion coefficients must satisfy:
2
0 q
d 1
01. dt0nq(t0)= 2Ψ(cid:16)Xj=0(cid:0)nq−j(t0)·gˆ(cid:1)nj(t0)−δq,0gˆ(cid:17)+B0(r0(t0),t0)nq−1(t0)
5 q 1q j 1
− − −
1 nj(t0) B0(r0(t0),t0)nk(t0) nq j k 1(t0). (S5)
: − · − − −
v j=0 k=0
X X (cid:0) (cid:1)
i
X Here B0 = O0 +ΛS0, and S0 and O0 are the symmetric and antisymmetric parts of the matrix A0 of fluid-velocity
r gradients. Eq. (S5) can be solved for nq in terms of integrals over solutions np with p<q. The lowest-order solution
a
in Ku, n0(t0), is deterministic. It does not depend on B, and is thus independent of the stochastic fluid-velocity
gradients. This deterministic solution n (t) depends on the initial orientation in a complicated way. But in the limit
0 0
of large times the dependence on the initial orientation must drop out (provided that Ψ = 0). Using this fact we
6
obtain n (t )= gˆ. We have verified explicitly that the dependence on the initial orientation drops out for all
0 0
→∞ −
averages evaluated in the Letter, by asymptotic expansion of n (t) for large values of t and by numerical evaluation
0 0
of the exact solution for n (t).
0 0
Using n (t)= gˆ we find the following solution for the expansion coefficients n (t):
0 0 q 0
−
t0
nq(t0)=−δq,0gˆ+ dt01 e(t01−t0)/(2Ψ)f(t01)+(e(t01−t0)/Ψ−e(t01−t0)/(2Ψ))(f(t01)·gˆ)gˆ , (S6)
Z0 h i
2
with
q 1 q 1q j 1
1 − − − −
f(t0)= nq j(t0) gˆ nj(t0)+B0(r0(t0),t0)nq 1(t0) nj(t0) B0(r0(t0),t0)nk(t0))nq j k 1(t0).
2Ψ − · − − · − − −
j=1 j=0 k=0
X(cid:0) (cid:1) X X (cid:0)
(S7)
Now we insert the resulting solution for n(t) into the equation of motion (S1). This yields:
0
dr0 =Kuu r (t),t +Φ ∞ n (t)Kuq . (S8)
0 0 0 0 q 0
dt
0 q=0
(cid:0) (cid:1) X
This equation has the implicit solution
t0 ∞
r0(t0)=rd0et(t0)+Ku dt01[u0 r0(t01),t01 +Φ nq(t01)Kuq−1]. (S9)
Z0 q=1
(cid:0) (cid:1) X
Here r (t) denotes the deterministic solution, obtained by letting Ku=0 in Eq. (S8) and using n (t)= gˆ:
d0et 0 0 0 −
r (t)=r (0) KuΦgˆt . (S10)
d0et 0 0 − 0
Now consider an actual particle path r (t). The difference δr (t) r (t) r (t) between this path and its
0 0 0 0 ≡ 0 0 − d0et 0
deterministic approximation is obtained from Eq. (S9). We see that δr (t) Ku. This ensures that δr (t) is small if
0 0 0 0
∼
Kuissmallenough,anditallowsustoexpandthefluidvelocityu r (t),t inδr aroundthedeterministictrajectory
0 0 0 0 0
r (t):
d0et 0 (cid:0) (cid:1)
u0 r0(t0),t0 =u0 rd0et(t0),t0 + δr0(t0)·∇0 u0 rd0et(t0),t0 +... . (S11)
Inserting this expansion and(cid:0)the corr(cid:1)espond(cid:0)ing expan(cid:1)sion(cid:2) for B0 r0(t(cid:3)0),t(cid:0)0 in (S7)(cid:1)into Eq. (S9) gives an expression
forδr (t)intermsofu r (t),t anditsderivatives,analogoustotheexpressionsdiscussedfortheinertial-particle
0 0 0 d0et 0 0 (cid:0) (cid:1)
problem in Ref. [S4].
(cid:0) (cid:1)
Finally,toevaluatetheaverage u weiterativelysubstituteδr (t)intoEq.(S11),andcutofftheresultingseries
0 0
h i∞
at the desired order in Ku. To first order in Ku we find for example:
t0
u0 r0(t0),t0 =u0 rd0et(t0),t0 +Ku dt01 u0 r0(t01),t01 +Φn1(t01) ·∇0 u0 rd0et(t0),t0 , (S12)
(cid:0) (cid:1) (cid:0) (cid:1) Z0 n(cid:2) (cid:0) (cid:1) (cid:3) o (cid:0) (cid:1)
with
t01
n1(t01)= dt02e(t02−t01)/(2Ψ) −B0 r0(t02),t02 gˆ+ gˆ·B0 r0(t02),t02 gˆ gˆ . (S13)
Z0 n (cid:0) (cid:1) (cid:2) (cid:0) (cid:1) (cid:3) o
Up to this point the expressions are valid for fluid-velocity fields with arbitrary distribution. Now we make use of the
fact that the fluid-velocity field is assumed to be Gaussian in the statistical model with, zero mean, Gaussian spatial
correlations, and with exponential time correlations [S4]:
1
hu0i(r01,t01)u0j(r02,t02)i= d(d 1) δij(d−1−(r01−r02)2)+(r01−r02)i(r01−r02)j e−(r01−r02)2/2−|t01−t02|. (S14)
−
(cid:2) (cid:3)
The form of the prefactor follows from isotropy, homogeneity, incompressibility, and normalisation u(0,0)2 = 1.
0
h i
Taking the average of Eqs. (S12) and (S13) we find:
KuΦ C ΨC 1 4Ψ2C +(4Ψ+1)C 2Ψ+1
u = gˆ 0 + Φ Φ 0 , (S15)
h 0i∞ − d 2Φ2 √2Φ3 F √2Φ − 4√2Φ3Ψ F √8ΦΨ
n h i h io
with C Λ + 1 + (d + 1)(Λ 1)x2 and [x] √πex2erfc[x]. In the same way we obtain the average of the
x
≡ − F ≡
fluid-velocity gradient matrix A0:
KuΦδ dgˆgˆ 4Ψ+1 ΨC 1 (2Ψ+1) 4Ψ2C +(4Ψ+1)C 2Ψ+1
ij i j Φ Φ 0
hA0i∞ = d 1−−d nC04Φ3Ψ + √2Φ4 Fh√2Φi− (cid:2) 8√2Φ4Ψ2 (cid:3)Fh√8ΦΨio. (S16)
3
The z-components of Eqs. (S15) and (S16) are plotted in Fig. 1 in the Letter [S1]. The asymptotic behaviour for
small values of Φ is obtained by series expansion in Φ:
d(1 Λ)+2(Λ+2) Ψ(4Ψ+1)
hA0zzi∞ ∼KuΦ2 − d (2Ψ+1)2 , (S17)
d(1 Λ)+2 Ψ
hu0zi∞ ∼− KuΦ −d 2Ψ+1. (S18)
These equations correspond to Eqs. (5a) and (5b) in the Letter.
A corresponding series expansion for large Φ of Eqs. (S15) and (S16) gives
Kud+1 π
A (1 Λ) , (S19)
h 0zzi∞ ∼ Φ 2d − 2
r
Ku(d(Λ 1)+2Λ)
hu0zi∞ ∼ Φ −2d . (S20)
These are Eqs. (6a) and (6b) in the Letter. We note that this result can also be found by an expansion for large Φ
of the original equations for an arbitrary homogeneous, isotropic and incompressible flow, for arbitrary values of Ku.
Because this limit is universal it is not necessary to assume that Ku is small.
4
SMALL-SCALE CLUSTERING
The small-scale spatial patterns of gyrotactic micro-swimmers are fractal [S5]. Fractal small-scale clustering is
quantified by fractal dimensions (see Ref. [S4] for a review). There are different fractal dimensions, characterising
different geometrical aspects of the spatial distribution. In the Letter we show results for the fractal Lyapunov
dimension in d = 2 spatial dimensions. Eq. (11) in the Letter shows that this dimension is defined in terms of the
Lyapunov exponents:
1 R(t) 1 (t)
0 0 0 0
λ τ lim ln and (λ +λ )τ lim lnA . (S21)
1 1 2
≡t0→∞t0 R0(0) ≡t0→∞t0 A0(0)
The exponent λ gives the rate which the separation R(t) between two initially nearby micro-swimmers decreases
1 0 0
or increases exponentially at long times. The sum λ +λ determines the rate at which the infinitesimal area A(t)
1 2 0 0
of the parallelogram spanned by the separation vectors between three initially nearby micro-swimmers decreases or
increases at long times. The problem of computing small-scale fractal clustering in d dimensions is a d+1-particle
problem.
We solve this problem using the method described in Ref. [S4]. Linearising the equation of motion Eq. (S1) we see
that the separation between two initially nearby micro-swimmers evolves according to
dR
0 =Ku(Rˆ Z0Rˆ)R0, (S22)
dt ·
0
where Rˆ denotes the unit vector along the separation direction between the two micro-swimmers, and Z0 =A0+ΦY0
is the particle-velocity gradient matrix. Here Y0 is the matrix of orientation gradients with elements Yi0j = ∂ni/∂rj0.
It follows from Eq. (S22) that the exponent λ is given by
1
λ1τ =Ku Rˆ Z0Rˆ . (S23)
h · i∞
In two spatial dimensions the sum of λ and λ is simply
1 2
(λ1+λ2)τ =Ku TrZ0 . (S24)
h i∞
To evaluate these averages we require the time evolution of Rˆ and Y0. This follows from linearisation of Eqs. (S1)
and (S2):
dRˆ
i =Ku[Z Rˆ (Rˆ Z Rˆ )Rˆ ], (S25)
dt i0j j − j j0k k i
0
dYi0j 1
= [gˆ Y n +n gˆ Y ]+Ku[(∂ B )n +B Y 2Λ(n S Y )n Λ(n (∂ S )n )n
dt 2Ψ k k0j i k k i0j j i0k k i0k k0j − k k0l l0j i− k j k0l l i
0
Λ(n S n )Y Y (A +ΦY )], (S26)
− k k0l l i0j − i0k 0kj k0j
where repeated indices are summed over. We solve Eq. (S26) in terms of a series expansion
∞
Y0(t0)= Y0q(t0)Kuq . (S27)
q=0
X
Inserting this expression into Eq. (S26) we find the following expression for the expansion coefficients:
t0
Yq;ij(t0)= dt01 e(t01−t0)/(2Ψ)Fij(t01)+(e(t01−t0)/Ψ−e(t01−t0)/(2Ψ))gˆigˆkFkj(t01) , (S28)
Z0 h i
where we the initial condition was put to zero because the corresponding terms do not contribute to steady-state
averages. The coefficients F (t) are given by
ij 0
q 1
1 −
F (t)= [gˆ Y n +n gˆ Y +δ [(∂ B (r (t),t))n Λn (∂ S (r (t),t)n )n
ij 0 2Ψ k r0;kj q−r;i q−r;k k r0;ij q,1 j i0k 0 0 0 k− k j k0l 0 0 0 l i
r=0
X
+B (r (t),t)Y 2Λn S (r (t),t)Y n Λn S (r (t),t)n Y Y A (r (t),t) ΦY Y ]. (S29)
i0k 0 0 0 k0j − k k0l 0 0 0 l0j i− k k0l 0 0 0 l i0j − i0k 0kj 0 0 0 − i0k k0j
5
Eq. (S25) can be solved implicitly as
t0
Rˆ (t)=Ku dt Z (t )Rˆ (t ) Rˆ (t )Z (t )Rˆ (t ) Rˆ (t ) . (S30)
i 0 01 i0j 01 j 01 − j 01 j0k 01 k 01 i 01
Z0 n (cid:2) (cid:3) o
By substituting Rˆ and Y0 in the right-hand side of Eqs. (S29) and (S30), Eqs. (S29) and (S30) can be recursively
iteratedtoanydesiredorderinKu. Expandingthefluidvelocityandgradientsofthefluidvelocityasintheprevious
Section gives equations for Rˆ and Y0 in terms of the fluid velocity evaluated along the deterministic trajectories.
Averaging gives the Lyapunov exponents to order Ku2
Ku2 1
λ τ = (d 3)Φ2+(3 d+(15 4d+d2)Φ2) (Rˆ gˆ)2 2(1+5Φ2) (Rˆ gˆ)4 Φ+ (3 d)Φ2 (S31)
1 d(d 1)Φ5 − − − h · i− h · i √2 −
− nh i h
+(3 2d+d2)Φ4 (3 d+(18 5d d2)Φ2+(9 2d+d2)Φ4) (Rˆ gˆ)2 +(2+12Φ2+6Φ4) (Rˆ gˆ)4 [1/(√2Φ)] ,
− − − − − − h · i h · i F
Ku2Ψ2 i o
(λ +λ )τ = (Λ+1) C +Φ2(d(Λ 1) 2) Φ C2 +2(d+1)(Λ 1)2Φ4 1/(√2Φ) /√2 . (S32)
1 2 dΦ3 Φ − − − Φ − F
(cid:8) (cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3) (cid:9)
Thesteady-statemoments (Rˆ gˆ)2p inEq.(S31)characterisetheanisotropyofthespatialpatterns. Themoments
h · i∞
are given by Eq. (8) derived in Ref. [S3] for particle pairs settling in a turbulent aerosol, upon replacing G in that
equation by Φ. The fact that the moments agree is a coincidence, only valid to lowest order in Ku. While
−
the problems are superficially related (to lowest order particles move on a deterministic trajectory r (t) through
d0et 0
turbulence)thedetailsarequitedifferent. UsingtheaboveexpressionsfortheLyapunovexponents(withλ extended
1
to order Ku4) we obtain the curves plotted in Fig. 2a in the Letter.
Series expansions of Eq. (S32) for small and large values of Φ give Eqs. (7) and (9) in the Letter
v η/u Ku(ΦΨ)2B (Λ) forΦ 1, (S33)
0 d
h∇· i∞ ∼− (cid:28)
v η/u Ku ΦΨ2E (Λ) forΦ 1. (S34)
0 d
h∇· i∞ ∼− (cid:29)
The shape-dependent prefactors B (Λ) and E (Λ) are found to be:
d d
B (Λ) [(d+2)(d+4) 2d(d+4)Λ+(4+2d+d2)Λ2 /d (S35)
d
≡ −
Ed(Λ) π/2(d+1)(d+3)(Λ 1)2/d. (cid:3) (S36)
≡ −
p
The corresponding limits for λ are more subtle due to the complicated dependence of the moments (Rˆ gˆ)2p on
1
Φ. For small values of Φ we have (Rˆ gˆ)2p (2p 1)!!/(2pp!) [S3] and consequently h · i∞
h · i∞ ∼ −
d2+d 3
λ τ 2 − Ku2Φ2. (S37)
1
∼− d(d 1)
−
This equation is quoted in the Letter for d=2.
The results derived here are obtained for general values of Φ, Ψ, and Λ. It is important to note that the limit of
large Ψ is singular. In this limit the dependence on the initial orientation n(t=0) does not decay fast enough for it
to be disregarded. This explains, as pointed out in the Letter, that the theory must fail for very large Ψ.
WealsonotethatEqs.(S31)and(S32)arederivedassumingthattheelementsofY0 aresmallenough,allowingusto
expandEq.(S26)inKu. ThismayfailwhensingularitiesoccurinthedynamicsofY0. ElementsofYmayrepeatedly
diverge to . These singularities are analogous to caustics in inertial-particle dynamics [S4] and correspond to
−∞
instances where the n-field becomes multi-valued so that nearby particles can have very different orientations. The
perturbation expansion leading to Eqs. (S31) and (S32) is expected to diverge when caustics are frequent, precisely
as in the inertial-particle problem (this point is discussed in detail in Ref. [S4]).