Distribution of particle displacementsdue to swimming microorganisms

# Distribution of particle displacements due to swimming microorganisms

Jean-Luc Thiffeault Department of Mathematics, University of Wisconsin – Madison, 480 Lincoln Dr., Madison, WI 53706, USA
###### Abstract

The experiments of Leptos et al.. [Phys. Rev. Lett. 103, 198103 (2009)] show that the displacements of small particles affected by swimming microorganisms achieve a non-Gaussian distribution, which nevertheless scales diffusively — the ‘diffusive scaling.’ We use a simple model where the particles undergo repeated ‘kicks’ due to the swimmers to explain the shape of the distribution as a function of the volume fraction of swimmers. The net displacement is determined by the inverse Fourier transform of a single-swimmer characteristic function. The only adjustable parameter is the strength of the stresslet term in our spherical squirmer model. We give a criterion for convergence to a Gaussian distribution in terms of moments of the drift function, and show that the experimentally-observed diffusive scaling is a transient related to the slow crossover of the fourth moment from a ballistic to a linear regime with path length. We also present a simple model, with logarithmic drift function, that can be solved analytically.

microswimmers; swimming microorganisms; effective diffusivity; particle transport

## I Introduction

The study of microswimming has exploded in recent years with the advent of precise, well-controlled experiments. (See for instance the reviews of Pedley and Kessler (1992) and Lauga and Powers (2009).) This has uncovered a plethora of fascinating behavior, for example the complex interaction of microswimmers with boundaries Rothschild (1963); Winet et al. (1984); Cosson et al. (2003); Lauga et al. (2006); Berke et al. (2008); Drescher et al. (2009), or the collective suspension instability (swirls and jets) at high concentrations of ‘pushers,’ organisms whose propulsion mechanism is at the rear Dombrowski et al. (2004); Hernandez-Ortiz et al. (2005); Underhill et al. (2008); Underhill and Graham (2011); Saintillan and Shelley (2007); Sokolov et al. (2009); Saintillan and Shelley (2012).

Another fruitful research direction is biogenic mixing, or biomixing for short. Does the motion of swimmers influence the effective diffusivity of passive scalars advected by the fluid, such as the nutrients the organisms depend on? This has been proposed as a mixing mechanism in the ocean Huntley and Zhou (2004); Dewar et al. (2006); Kunze et al. (2006); Katija and Dabiri (2009); Dabiri (2010); Leshansky and Pismen (2010); Thiffeault and Childress (2010); Lorke and Probst (2010); Katija (2012), though its effectiveness is still very much open to debate Visser (2007); Gregg and Horne (2009); Kunze (2011); Noss and Lorke (2014). Biomixing has also been studied in suspensions of small organisms Ishikawa et al. (2010); Kurtuldu et al. (2011); Miño et al. (2011); Zaid et al. (2011).

The main ingredient in formulating a theory for the enhanced diffusion due to swimming organisms is the drift caused by the swimmer Maxwell (1869); Darwin (1953); Lighthill (1956). Katija and Dabiri (2009) and Thiffeault and Childress (2010) proposed that the enhanced diffusivity is due to the repeated displacements induced by a swimmer on a particle of fluid. Thiffeault and Childress (2010) and Lin et al. (2011) formulated a probabilistic model where, given the drift caused by one swimmer, an effective diffusivity could be computed. This model has been tested in physical and numerical experiments Jepson et al. (2013); Morozov and Marenduzzo (2014); Kasyap et al. (2014) and modified to include curved trajectories Pushkin and Yeomans (2013) and confined environments Pushkin and Yeomans (2014). Miño et al. Miño et al. (2011, 2013) observe that effective diffusivity is inversely related to swimming efficiency, and find increased diffusivity near solid surfaces, both theoretically and experimentally. The drift caused by individual microswimmers has also been studied in its own right Dunkel et al. (2010); Pushkin et al. (2013). Pushkin and Yeomans (2013) also found an analytical expression for stresslet displacements, valid in the far field.

The studies mentioned above have typically been concerned with the effective diffusivity induced by the swimmers, but one can also ask more detailed questions about the distribution of displacements of fluid particles. Wu and Libchaber (2000) studied the displacement of spheres larger than the swimming organisms. More recently, Leptos et al. (2009) studied the microscopic algae Chlamydomonas reinhardtii. They used spheres that are much smaller than the organisms, so their distributions can be taken to be close to the displacements of idealized fluid particles. The probability density function (pdf) of tracer displacements was found to be strongly non-Gaussian, though the distributions scaled ‘diffusively’: they collapsed onto each other if rescaled by their standard deviation.

Several papers have dealt with these non-Gaussian distributions. Zaid et al. (2011) examine the velocity fluctuations due to swimmers modeled as regularized point stresslets, and obtain strongly non-Gaussian tails. The non-Gaussianity in their case is due to the divergence of the stresslet near the singularity, which indicates large displacements. While the broad outline of this mechanism is surely correct, examining this singular limit is questionable: it is never valid to evaluate terms such as the stresslet in the singular limit, since the swimmer’s body necessarily regularizes the velocity. In addition, no direct comparison to experiments is offered beyond a comment that the data ‘resemble the measurements of Leptos et al. (2009).’ Pushkin and Yeomans (2014) extended this work to confined environments, and we will contrast their results to ours. As we will show here, the non-Gaussianity arises from the rarity of interaction events — the system is very far from the Gaussian limit. Note also that Eckhardt and Zammert (2012) have fitted the distributions of Leptos et al. (2009) very well to a continuous-time random walk model, but this does not suggest a mechanism and requires fitting different parameters at each concentration.

What causes the non-Gaussian form of the displacement distribution? As was pointed out by Pushkin and Yeomans (2014), the experiments are run for a very short time. Let us quantify what is meant by ‘short.’ Leptos et al. (2009) define a ‘sphere of influence’ of radius  around a particle: swimmers outside that sphere do not significantly displace the particle. If swimmers with number density  moves a distance  in random directions, the expected number of ‘interactions’ with a target particle is roughly

 nλπR2eff∼0.4.

Here we took  and , which are the largest values used in the experiments, and  as estimated in Leptos et al. (2009). Hence, a typical fluid particle feels very few near-encounters with any swimmer. In order for the central limit theorem to apply, the net displacement must be the sum of many independent displacements, and this is clearly not the case here for the larger values of the displacement. We thus expect a Gaussian core (due to the many small displacements a particle feels) but non-Gaussian tails (due to the rarity of large displacements), which is exactly what was observed in the experiments.

Here, we present a calculation that quantitatively predicts essentially all the details of the distributions obtained by Leptos et al. (2009). The underlying model is not new, being based on the particle-displacement picture of Thiffeault and Childress (2010) and Lin et al. (2011). However, the analysis is new: we show how to combine multiple displacements to obtain the probability density function due to multiple swimmers, and take the appropriate infinite-volume limit. As we go, we discuss the mathematical assumptions that are required. Upon comparing with experiments, we find the agreement to be excellent, in spite of the differences between our model swimmer and the experiments. Only a single parameter needs to be fitted: the dimensionless stresslet strength, .

The paper is organized as follows. In Section II we derive the probability density of displacements based on the drift function of a single swimmer, in the infinite-volume limit. We use numerical simulations of a model swimmer (of the squirmer type Lighthill (1952); Blake (1971); Ishikawa et al. (2006); Ishikawa and Pedley (2007); Drescher et al. (2009)) in Section III to obtain a distribution of displacements which we match to the experiments of Leptos et al.  In Section IV we give a different interpretation of the main formula of Section II in terms of ‘interactions’ between swimmers and a fluid particle. This alternative form can be used to obtain some analytic results, in particular when the drift function is logarithmic. We examine in Section V the long-time (or long swimming path) asymptotics of the model, and find what features of the drift function affect the convergence to Gaussian. In Section VI we address the ‘diffusive scaling’ observed in the experiments, and show that it is a transient phenomenon. Finally, we discuss our results as well as future directions in Section VII.

## Ii Distribution of displacements

The setting of our problem is a large volume  that contains a number of swimmers , also typically large. The swimmers move independently of each other in randomly directions. In the dilute limit that we consider, the velocity field of one swimmer is not significantly affected by the others. A random fluid particle (not too near the edges of the domain), will be displaced by the cumulative action of the swimmers. If we follow the displacements of a large number of well-separated fluid particles, which we treat as independent, we can obtain the full pdf of displacements. Our goal is to derive the exact pdf of displacements from a simple probabilistic model. Our starting point is the model described by Thiffeault and Childress (2010) and improved by Lin et al. (2011), which captures the important features observed in experiments.

For simplicity, we assume the swimmers move along straight paths at a fixed speed . The velocity field induced at point  by a swimmer is , with the time dependence reflecting the motion of the swimmer. The main ingredient in the model is the finite-path drift function  for a fluid particle, initially at , affected by a single swimmer:

 Δλ(η)=∫λ/U0u(x(s)−Us)ds,˙x=u(x−Ut),x(0)=η. (1)

Here  is the swimming distance. To obtain  we must solve the differential equation  for each initial condition . Assuming homogeneity and isotropy, we obtain the probability density of displacements Pushkin and Yeomans (2014),

 pR1λ(r)=1Ωrd−1∫Vδ(r−Δλ(η))dVηV (2)

where  is the area of the unit sphere in  dimensions: , . Here  is a random variable that gives the displacement of the particle from its initial position after being affected by a single swimmer with path length . We denote by the pdf of . Because of the isotropy assumption, only the magnitude  enters (2).

Before we continue with finding the pdf for multiple swimmers, let us investigate how the variance of displacements evolves. The second moment of  is

 ⟨(R1λ)2⟩=∫Vr2pR1λ(r)dVr=∫VΔ2λ(η)dVηV. (3)

This typically goes to zero as , since a single swimmer in an infinite volume shouldn’t give any fluctuations on average. We write  for the random particle displacement due to  swimmers; the second moment of  is

 ⟨(RNλ)2⟩=N⟨(R1λ)2⟩=n∫VΔ2λ(η)dVη (4)

with  the number density of swimmers. This is nonzero (and might diverge) in the limit , reflecting the cumulative effect of multiple swimmers. Note that this expression is exact, within the problem assumptions: it doesn’t even require  to be large.

The expression (4) will lead to diffusive behavior if the integral grows linearly in  (or if the swimmers change direction Lin et al. (2011), which we shall not treat here). Surprisingly, it has been found to do so in two distinct ways. In the first, exemplified by bodies in inviscid flow Thiffeault and Childress (2010); Lin et al. (2011), the support of  grows linearly with , but the displacements themselves become independent of  when  is large. The intuition is that the swimmer pushes particles a finite distance as it encounters them. As we wait longer, the volume of such displaced particles grows linearly in , but once particles are displaced they are left behind and suffer no further displacement. This diffusive behavior is thus appropriate for very localized interactions, where the only displaced particles are very near the axis of swimming. This tends to occur in inviscid flow, or for spherical ‘treadmillers’ in viscous flow. See Fig. 1 for an illustration.

The second situation in which (4) shows diffusive behavior even for straight swimming paths is when the far-field velocity has the form of a stresslet, as is the case for a force-free swimmer in a viscous fluid. This diffusive behavior was observed in Lin et al. (2011) but it was Pushkin and Yeomans (2013) who provided a full explanation. For a stresslet swimmer, the main contributions to (4) come from  of order , so it is appropriate to use a point singularity model swimmer for large . In that case the drift function has the scaling , where  is a dimensionless variable and the function  is independent of  for large  Pushkin and Yeomans (2013). Inserting this form in (4), we find

 ∫Δ2λ(η)dVη=∫(λ−2D2(ζ))(λ3dVζ)∼λ. (5)

The integral of  converges despite having singularities Pushkin and Yeomans (2013). We thus see that the integral in (4) grows linearly in  for very different reasons than our first case: here the volume of particles affected by the swimmer grows as  (particles are affected further and further away), but they are displaced less (since they are further away, see Fig. 1). Any truncation of the integral in (5) (because of finite volume effect) will lead to a decrease in the diffusivity, a possible origin for the decrease in diffusivity with path length observed in Jepson et al. (2013). Note also that the reorientation mechanism discussed by Lin et al. (2011) is not necessary in this case to achieve the diffusive behavior, as pointed out by Pushkin and Yeomans (2014).

Having addressed the growth of the variance, we continue with finding the pdf of displacements for multiple swimmers. We write  for a single coordinate of  (which coordinate is immaterial, because of isotropy). From (2) with  we can compute , the marginal distribution for one coordinate:

 pX1λ(x)=∫∞−∞pR1λ(r)dy=∫V∫∞−∞12πrδ(r−Δλ(η))dydVηV. (6)

Since , the -function will capture two values of , and with the Jacobian included we obtain

 pX1λ(x)=1π∫V1√Δ2λ(η)−x2[Δλ(η)>|x|]dVηV, (7)

where  is an indicator function: it is  if  is satisfied, otherwise.

The marginal distribution in the three-dimensional case proceeds the same way from (2) with :

 pX1λ(x)=∫∞−∞pR1λ(r)dydz=∫V∫∞−∞∫∞−∞14πr2δ(r−Δλ(η))dydzdVηV. (8)

Again with  the -function captures two values of , and with the Jacobian included we obtain

 pX1λ(x)=12π∫V∫∞−∞1Δλ(η)1√Δ2λ(η)−x2−y2[Δ2λ(η)>x2+y2]dydVηV. (9)

Now we integrate over  to get

 pX1λ(x)=12∫V1Δλ(η)[Δλ(η)>|x|]dVηV (10)

which is the three-dimensional analogue of (7). The integrand of (10) has an intuitive interpretation. The indicator function says that a displacement in a random direction must at least be larger than  to project to a value . The factor of  in the denominator then tells us that large displacements in a random direction are less likely to project to a value . (The two-dimensional form (7) has essentially the same interpretation, with a different weight.)

In order to sum the displacements due to multiple swimmers, we need the characteristic function of , defined by

 ⟨eikX1λ⟩=∫∞−∞pX1λ(x)eikxdx. (11)

For the two-dimensional pdf (7), we have

 ⟨eikX1λ⟩=∫VJ0(kΔλ(η))dVηV (12)

where  is a Bessel function of the first kind. For the three-dimensional pdf (10), the characteristic function is

 ⟨eikX1λ⟩=∫Vsinc(kΔλ(η))dVηV (13)

where  for , and .111Beware that this function is sometimes defined as , most notably by Matlab. The expression (13) appears in Pushkin and Yeomans (2014), except here we compute it directly from a spatial integral rather than from the pdf of . The main difference will come in the way we take the limit below, which will allow us to study the number density dependence directly.

We define

 γ(x)\raisebox0.2967pt:={1−J0(x),d=2;1−sincx,d=3, (14)

We have , , so  as . For large argument, . We can then write the two cases (12)–(13) for the characteristic function together as

 ⟨eikX1λ⟩=1−(vλ/V)Γλ(k) (15)

where

 Γλ(k)\raisebox0.2967pt:=1vλ∫Vγ(kΔλ(η))dVη. (16)

Here  is the volume ‘carved out’ by a swimmer moving a distance :

 vλ=λσ (17)

with  the cross-sectional area of the swimmer in the direction of motion.

Since we are summing independent particle displacements, the probability distribution of the sum is the convolution of  one-swimmer distributions. Using the Fourier transform convolution property, the characteristic function for  swimmers is thus . From (15),

 ⟨eikX1λ⟩N=(1−vλΓλ(k)/V)nV, (18)

where we used , with  the number density of swimmers. We will need the following simple result:

###### Proposition 1.

Let  as  for an integer ; then

 (1−εy(ε))1/ε=exp(−M∑m=1εm−1ym(ε)m)(1+o(ε0)),ε→0. (19)

See Appendix A for a short proof.

Let’s examine the assumption of Proposition 1 for  applied to (18), with  and . For , the assumption of Proposition 1 requires

 Γλ(k)∼o(V1/2),V→∞. (20)

A stronger divergence with  means using a larger  in Proposition 1, but we shall not need to consider this here. Note that it is not possible for  to diverge faster than , since  is bounded. In order for  to diverge as , the displacement must be nonzero as , an unlikely situation that can be ruled out.

Assuming that (20) is satisfied, we use Proposition 1 with  to make the large-volume approximation

 ⟨eikX1λ⟩N=(1−vλΓλ(k)/V)nV∼exp(−nvλΓλ(k)),V→∞. (21)

If the integral  is convergent as  we have achieved a volume-independent form for the characteristic function, and hence for the distribution of  for a fixed swimmer density. We define the quantity

 νλ\raisebox0.2967pt:=nvλ=λ/ℓmfp (22)

where  is the swimmer mean free path. Since is the volume carved out by a single swimmer moving a distance  (Eq. (17)), is the expected number of swimmers that will ‘hit’ a given fluid particle.

A comment is in order about evaluating (16) numerically: if we take  to , then , and thus , which then leads to  in (21). This is negligible as long as the number of swimmers  is moderately large. In practice, this means that  only needs to be large enough that the argument of the decaying exponential in (21) is of order one, that is

 νλΓλ(kmax)∼O(1). (23)

Wavenumbers  do not contribute to (21). (We are assuming monotonicity of  for , which will hold for our case.) Note that (23) implies that we need larger wavenumbers for smaller densities : a typical fluid particle then encounters very few swimmers, and the distribution should be far from Gaussian.

Now that we’ve computed the characteristic function for  swimmers (21), we finally recover the pdf of  for  swimmers as the inverse Fourier transform

 pXλ(x)=12π∫∞−∞exp(−νλΓλ(k))e−ikxdk, (24)

where we dropped the superscript  from  since the number of swimmers no longer enters the expression directly.

## Iii Comparing to experiments

We now compare the theory discussed in the previous sections to the experiments of Leptos et al.  in particular the observed dependence of the distribution on the number density . (Another aspect of their experiments, the ‘diffusive scaling’ of the distributions, will be discussed in Section VI.) In their experiments they use the microorganism C. reinhardtii, an alga of the ‘puller’ type, since its two flagella are frontal. This organism has a roughly spherical body with radius . They observe a distribution of swimming speeds with a strong peak around . They place fluorescent microspheres of about a micron in radius in the fluid, and optically measure their displacement as the organisms move. The volume fraction of organisms varies from  (pure fluid) to .

They measure the displacement of the microspheres along a reference direction, arbitrarily called  (the system is assumed isotropic). Observing many microspheres allows them to compute the pdf of tracer displacements , which we’ve denoted . Thus, is the probability of observing a particle displacement  after a path length . (They write their density , where  are the same as our .)

At zero volume fraction (), the pdf  is Gaussian, due solely to thermal noise. For higher number densities, Leptos et al. see exponential tails appear and the Gaussian core broaden. The distribution is well-fitted by the sum of a Gaussian and an exponential:

 pXλ(x)=1−f√2πδ2ge−x2/2δ2g+f2δee−|x|/δe. (25)

They observe the scalings  and , where  and  depend on . The dependence on  is referred to as the ‘diffusive scaling’ and will be discussed in Section VI. Exploiting this scaling, Eckhardt and Zammert (2012) have fitted these distributions very well to a continuous-time random walk model, but this does not suggest a mechanism.

We shall use a model swimmer of the squirmer type Lighthill (1952); Blake (1971); Ishikawa et al. (2006); Ishikawa and Pedley (2007); Drescher et al. (2009), with axisymmetric streamfunction Lin et al. (2011)

 Ψsf(ρ,z)=12ρ2U{−1+ℓ3(ρ2+z2)3/2+32βℓ2z(ρ2+z2)3/2(ℓ2ρ2+z2−1)} (26)

in a frame moving at speed . Here  is the swimming direction and  is the distance from the  axis. To mimic C. reinhardtii, we use  and . (Leptos et al. (2009) observe a distribution of velocities but the peak is near .) We take  for the relative stresslet strength, which gives a swimmer of the puller type, just like C. reinhardtii. The contour lines of the axisymmetric streamfunction (26) are depicted in Fig. 2. The parameter  is the only one that was fitted (visually) to give good agreement later.

First we compute the drift function  for a single swimmer moving along the  axis. The model swimmer is axially symmetric, so  can be written in terms of  and , the perpendicular distance to the swimming axis. We take , since the time is  in Fig. 2(a) of Leptos et al., and our swimmer moves at speed . We compute  for a large grid of  and  values, using the analytic far-field stresslet form for the displacement Dunkel et al. (2010); Miño et al. (2013); Pushkin and Yeomans (2013) when far away from the swimmer’s path.

From the drift function  we now want to compute  defined by (16). To estimate how large a  value we will need, we start from the smallest volume fraction in the experiments, . For spherical swimmers of radius  (with cross-sectional area ), this gives a number density of . We thus get . The criterion (23) then tells us that we need  large enough that .

Figure 3 shows the numerically-computed  for several values of , with  the broadest curve. We can see from the figure that choosing  will ensure that  is large enough. As  gets larger, decreases, reflecting the trend towards the central limit theorem (which corresponds to the small- expansion of , see Section V). Note also that  tends to become independent of  as  gets larger.

To obtain  and compare to Leptos et al., we must now take the inverse Fourier transform of , as dictated by (24). This is straightforward using Matlab’s ifft routine. The ‘period’ (domain in ) is controlled by the spacing of the grid, so we make sure the grid is fine enough to give us the largest values of  required. We also convolve with a Gaussian distribution of half-width  to mimic thermal noise. This follows from the value  measured by Leptos et al. for the diffusivity of the microspheres. The value of  is consistent with the Stokes–Einstein equation for the diffusivity of thermally-agitated small spheres in a fluid.

The results are plotted in Fig. 4 and compared to the data of Fig. 2(a) of Leptos et al. (2009). The agreement is excellent: we remind the reader that we adjusted only one parameter, . This parameter was visually adjusted to the  data in Fig. 4, since the larger concentration is most sensitive to ; a more careful fit is unnecessary given the uncertainties in both model and data. (The inset shows the sensitivity of the fit to .) All the other physical quantities were gleaned from Leptos et al.  What is most remarkable about the agreement in Fig. 4 is that it was obtained using a model swimmer, the spherical squirmer, which is not expected to be such a good model for C. reinhardtii. The real organisms are strongly time-dependent, for instance, and do not move in a perfect straight line. Nevertheless the model captures very well the pdf of displacements, in particular the volume fraction dependence. The model swimmer slightly underpredicts the tails, but since the tails are associated to large displacements they depend on the near-field details of the swimmer, so it is not surprising that our model swimmer should deviate from the data.

In Figure 4 we compare our results to the phenomenological fit of Eckhardt and Zammert (2012) based on continuous-time random walks: their fit is better in the tails, but our models disagree immediately after the data runs out. Our model has the realistic feature that the distribution is cut off at the path length , since it is extremely unlikely that a particle had two close encounters with a swimmer at these low volume fractions.

A possible explanation as to why the squirmer model does so well was provided by Pushkin and Yeomans (2014). They used numerical simulations of squirmers (with a larger value  that leads to a trapped volume) to show that the tails of distribution scale as , which is the asymptotic form of the stresslet displacement distribution. Figure 5 shows that our computations have a similar tail, though we emphasize here that our agreement with the experiments of Leptos et al. (2009) is quantitative and correctly reproduces the volume fraction dependence. We also point out that though the trend in Fig. 5 follows , the slope changes gradually and does not have a clear power law (the log scale means the deviations are quite large). The inset in Fig. 5 is a numerical simulation that includes only the singularity in the stresslet displacement, , as assumed in the analysis of Pushkin and Yeomans (2014). Though the  tails are eventually achieved, they have far lower probability than needed to explain the numerics. Pushkin and Yeomans’s use of the far-field stresslet form to predict the tails is thus questionable, at least for short path lengths.

For the effective diffusivity, Leptos et al. (2009) give the formula , with  and . Elsewhere in their paper they also give  for the diffusivity of the microspheres in the absence of swimmers, but their fitting procedure changes the intercept slightly. (Here we used , but the difference is minute.) Figure 6 shows the numerically-computed effective diffusivity for our squirmer model, as a function of . This curve is as in Lin et al. (2011), Fig. 6(a), except that we corrected the integrals in the far field using the analytic expression of Pushkin and Yeomans (2013), which gives a more accurate result. The Figure also shows the fit

 Deff−D0Unℓ4≃0.266+34πβ2, (27)

which is fairly good over the whole range (keeping in mind that this is a logarithmic plot, so the discrepancy at moderate  are of the order of ). Here the value  is the diffusivity due to spheres in inviscid flow (, see Thiffeault and Childress (2010)), and is the large- analytic expression Pushkin and Yeomans (2013) for stresslets. From the data in Fig. 6 we find , significantly larger than Leptos et al. (2009), as can be seen in Fig. 6. The solid line is their fit, the dashed is our model prediction for . The overestimate is likely due to the method of fitting to the squared displacement: their Fig. 3(a) clearly shows a change in slope with time, and the early times tend to be steeper, which would increase the effective diffusivity. Note also that their Fig. 3(a) has a much longer temporal range than their PDFs, going all the way to  (compared to ), raising the possibility that particles were lost by moving out of the focal plane.

## Iv The ‘interaction’ viewpoint

Equation (24) gives the exact solution for the distribution of uncorrelated displacements due to swimmers of number density . In this section we derive an alternative form, in terms of an infinite series, which is often useful and provides an elegant interpretation for (24).

The displacement  typically decays rapidly away from the swimmer, so that it may often be taken to vanish outside a specified ‘interaction volume’ . Then from (16), since , we have

 Γλ(k)=1vλ∫˜Vλγ(kΔλ(η))dVη=˜Vλvλ(1−˜Γλ(k)) (28)

where

 ˜Γλ(k)=1˜Vλ∫˜Vλ(1−γ(kΔλ(η)))dVη. (29)

Define ; we insert (28) into (24) and Taylor expand the exponential to obtain

 pXλ(x)=∞∑m=0~νmλm!e−~νλ12π∫∞−∞˜Γmλ(k)e−ikxdk. (30)

The factor  is a Poisson distribution for the number of ‘interactions’ between swimmers and a particle: it measures the probability of finding  swimmers inside the volume . The inverse transform in (30) gives the -fold convolution of the single-swimmer displacement pdf. This was the basis for the model used in Thiffeault and Childress (2010); Lin et al. (2011) and in an earlier version of this paper Thiffeault (2014). We have thus shown that formula (24) is the natural infinite-volume limit of the interaction picture.

Formula (30) is very useful in many instances, such as when  is small, in which case only a few terms are needed in (30) for a very accurate representation. Note that the first term of the sum in (30) is a -function, which corresponds to particles that are outside the interaction volume . This singular behavior disappears after is convolved with a Gaussian distribution associated with molecular noise.

Let us apply (30) to a specific example. A model for cylinders and spheres of radius  traveling along the  axis in an inviscid fluid Thiffeault and Childress (2010); Lin et al. (2011) is the log model,

 Δλ(η)={Cln+(ℓ/ρ),if 0≤z≤λ,0,otherwise, (31)

where  is the perpendicular distance to the swimming direction and . The logarithmic form comes from the stagnation points on the surface of the swimmer, which dominate transport in this inviscid limit. This model is also appropriate for a spherical ‘treadmiller’ swimmer in viscous flow. The drift function (31) resembles Fig. 1.

For the form (31) the interaction volume  is the same as , the volume carved out during the swimmer’s motion (Eq. (17)). By changing integration variable from  to  in (16) we can carry out the integrals explicitly to obtain (see Appendix B)

 ˜Γλ(k)={(1+(Ck)2)−1/2,(% cylinders);(Ck/2)−1arctan(Ck/2),(spheres). (32)

This is independent of , even for short paths (but note that (31) is not a good model for ).

Furthermore, for  we can also explicitly obtain the convolutions that arise in (30) to find the full distribution,

 pXλ(x)=e−νλ(δ(x)+∞∑m=1νmλm!1C√πΓ(m/2)(|x|/2C)(m−1)/2K(m−1)/2(|x|/C)), (33)

where  are modified Bessel functions of the second kind, and  is the Gamma function (not to be confused with  above). Equation (33) is a very good approximation to the distribution of displacements due to inviscid cylinders. Unfortunately no exact form is known for spheres: we must numerically evaluate (24) with (32) or use asymptotic methods (see Section V).

## V Long paths: Large-deviation theory

In Section IV we derived an alternative form of our master equation (24) as an expansion in an ‘interaction’ volume. Here we look at another way to evaluate the inverse Fourier transform in (24), using large-deviation theory Gärtner (1977); Ellis (1984, 1985); Touchette (2009). In essence, large-deviation theory is valid in the limit when a particle encounters many swimmers, so that  is large (in practice ‘large’ often means order one for a reasonable approximation). This includes the central limit theorem (Gaussian form) as a special case. In this section we provide a criterion for how much time is needed before Gaussian behavior is observed, which can help guide future experiments.

Earlier we used the characteristic function (21). Here it is more convenient to work with the moment-generating function, which in our case can be obtained simply by letting . The moment-generating function of the distribution is then

 ⟨esXλ⟩=exp(−νλΓλ(−is))=exp(νλΛ(s))

where  was defined by Eq. (22), and

 Λ(s)\raisebox0.2967pt:=1νλln⟨esXλ⟩=−Γλ(−is) (34)

is the scaled cumulant-generating function. As its name implies, this function has the property that its derivatives at  give the cumulants of  scaled by , for example

 Λ′′(0)=ν−1λ⟨X2λ⟩,Λ′′′′(0)=ν−1λ(⟨X4λ⟩−3⟨X2λ⟩2), (35)

where we left out the vanishing odd moments. We left out the  subscript on  since we assume that it becomes independent of  for large .

If  is differentiable over some interval of interest,  satisfies a large-deviation principle Gärtner (1977); Ellis (1984, 1985); Touchette (2009),

 pXλ(x)∼e−νλI(x/νλ)+o(νλ),νλ≫1, (36)

where  is the rate function, which is the Legendre–Fenchel transformation of :

 I(X)=sups∈R{sX−Λ(s)}. (37)

The large-deviation principle is in essence an application of the method of steepest descent for large .

The scaled cumulant-generating function  is always convex, which guarantees a unique solution to (37). The rate function  is also convex, with a global minimum at . This means that for small  we can use the Taylor expansion

 I(X)=12I′′(0)X2+14!I′′′′(0)X4+O(X6) (38)

to write

 pXλ(x)∼e−12I′′(0)x2/νλ,x≪cνλ,νλ≫1, (39)

with . This is a Gaussian approximation with variance , which can be shown to agree with (4) after multiplying by . To recover a Gaussian distribution over an appreciable range of  (say, a standard deviation) we insert  in the condition  to find the Gaussian criterion

 νλ≫112|I′′′′(0)|(I′′(0))2=112|Λ′′′′(0)|(Λ′′(0))2. (40)

After using  to find the cumulants, we can rewrite this as

 (41)

where  is the volume of one swimmer. When (40) or (41) is satisfied, we can expect that the distribution will be Gaussian (except in the far tails). (The constant prefactor in (41) is only valid for  or .) The criterion (41) can be interpreted as the minimum volume fraction  required to observe Gaussian behavior, roughly within a standard deviation of the mean. We note that, at small swimmer volume fraction, a long time (i.e., path length ) is required to achieve the Gaussian form. Figure 7 highlights this: the solid curve is  from Eq. (41) for the squirmer model in Section III, with parameter values appropriate for the experiments of Leptos et al. (2009). Their experiments had , so they are in the slowly-decreasing region of Fig. 7, before more rapid  convergence sets in after . It is thus not surprising that Gaussian tails were not observed in the experiments.

As an illustration of the large-deviation approach, we consider again the inviscid cylinder and sphere results (32). We have then respectively

 Λ(s)={(1−(Cs)2)−1/2−1,(cylinders);(Cs/2)−1arctanh(Cs/2)−1,(spheres). (42)

We can see from (37) that the singularities in (42) ( for cylinders, for spheres) immediately lead to  and  as , respectively, corresponding to exponential tails in (36) independent of . These are the displacements of particles that come near the stagnation points at the surface of the cylinder or sphere Lin et al. (2011). We can also use (42) to compute the constant on the right-hand side of (40): (cylinders) and  (spheres), which are both of order unity. This reflects the fact that the drift function  is very localized, so convergence to Gaussian is tied directly to the volume carved out by the swimmers. For swimmers with a longer-range velocity field, such as squirmers, the constant is much larger, as reflected by the large difference between the solid (squirmers) and the dashed (inviscid spheres) curves in Fig. 7.

For inviscid cylinders the Legendre–Fenchel transform (37) can be done explicitly to find (with )

 I(X)=1−√3πα(12−α2X−2)−1/2+12√πα((πα−4)α−2X2+13)1/2 (43)

where  is defined by

 α3(X)=6(√(9πX4)2+48X6−9πX4). (44)

For spheres (37) must be solved numerically for each , which is straightforward since this is a one-dimensional problem with a unique solution. The function  for both cylinders and spheres is plotted in Fig. 8.

## Vi The diffusive scaling

One of the most remarkable property of the pdfs found by Leptos et al. is the diffusive scaling. This is illustrated in Fig. 9: the unrescaled displacement pdfs are shown in Fig. 9; the same pdfs are shown again in Fig. 9, but rescaled by their standard deviation. The pdfs collapse onto a single curve (the shortest path length collapses more poorly). Figure 9 was obtained in the same manner as Fig. 4, using our probabilistic approach. Hence, the diffusive scaling is also present in our model, as it was in the direct simulations of Lin et al. (2011) for a similar range of path lengths. In Fig. 9 we left out thermal diffusion completely, which shows that it is not needed for the diffusive scaling to emerge.

Here we have the luxury of going much further in time and to examine the probability of larger displacements, since we are simply carrying out integrals and not running a statistically-limited experiment or simulation. (The numerical integrals are of course limited by resolution.) Figure 10 shows much longer runs (maximum  compared to  in the experiments). We see that, though the diffusive scaling holds in the core (as it must, since the core is Gaussian), the tails are narrowing, consistent with convergence to a Gaussian distribution but breaking the diffusive scaling. We now explain why the diffusive scaling appears to hold for some time, but eventually breaks down.

To understand the origin of the diffusive scaling, let us first examine how the integrated moments of  change with . Figure 11 shows the evolution of the spatial integrals of  and  for our squirmer model. For short , the moment of  grows as . This is a typical ‘ballistic’ regime: it occurs because for short times the integrals are dominated by fluid particles that are displaced proportionately to the swimmer’s path length. These particles are typically very close to the swimmer, and get dragged along for a while. This regime is visible for  in Fig. 11.

As  becomes larger, the particles initially near the swimmer are left behind, and thus undergo only a finite displacement even as  increases. Eventually, for  the scenario illustrated in Fig. 1 takes over and leads to linear growth of the moment with . This can be seen in Fig. 11 (triangles) for , though the scaling already looks fairly linear at . For  the moment also eventually grows linearly with , but the mechanism is different: the larger power downplays the far-field stresslet effect, and the near-field dominates. The linear growth is thus due to a corresponding linear growth of the support of  as in Fig. 1. This can be seen in Fig. 11 (dots) for , as indicated by a dashed line (see Appendix B for the computation of this asymptotic form). The crucial fact is that for  the crossover from  to  takes much longer than for . This is because the larger power weighs the largest displacements (with ) more heavily, so they dominate for longer before becoming too rare. This crossover is at the heart of the diffusive scaling, as we now show.

Let us assume that the distribution  does satisfy a diffusive scaling, such that is independent of . From (24), after changing integration variable to ,

 ~pXλ(~x)=√λpXλ(√λ~x)=12π∫∞−∞exp(−νλΓλ(~k/√λ))e−i~k~xd~k. (45)

Hence, a diffusive scaling law requires that be independent of . Using this scaling in (16), we have

 νλΓλ(~k/√λ)=n∫Vγ(Δλ(η)~k/√λ)dVη. (46)

We Taylor expand  (for ):

 νλΓλ(~k/√λ)/n=16~k2λ−1∫VΔ2λ(η)dVη+1120~k4λ−2∫VΔ4λ(η)dVη+O(~k6). (47)

The first term recovers the Gaussian approximation; the second is the first correction to Gaussian. Again this must be independent of  to obtain a diffusive scaling, so we need

 ∫VΔ2λ(η)dVη∼λ,∫VΔ4λ(η)dVη∼λ2, (48)

and clearly in general we would need each even moment  to scale as . However, we’ve already seen that all the moments typically eventually scale linearly with , so there can be no diffusive scaling. Because there is a transition from a power larger than  () to one less that  (), observe that in Fig. 11 there is a range of  (roughly ) where  is tangent to the  curve, as indicated by the line segment. In that range the distribution will appear to have a reasonably good diffusive scaling, consistent with Fig. 9. But, as we saw in Fig. 10, the diffusive scaling does not persist for larger . It is a coincidence that the range of  used in the experiments of Leptos et al. (2009) were exactly in that intermediate regime.

## Vii Discussion

In this paper, we showed how to use the single-swimmer drift function to fully derive the probability distribution of particle displacements. We took the limit of infinite volume and discussed the underlying assumptions, such as the need for the function  in (20) to not diverge too quickly with volume. In typical cases, the function becomes independent of volume as we make  large, but it is possible for the integral to diverge with , as may occur for example in sedimentation problems. If the divergence is rapid enough a larger value of  would need to be used when applying Proposition 1, potentially leading to interesting new distributions. Whether this can happen in practice is a topic for future investigation.

An intriguing question is: why does the squirmer model do so well? As was observed previously Lin et al. (2011); Pushkin and Yeomans (2014), it reproduces the pdf very well in the core and part of the tails (Fig. 4). However, the high precision of our calculation reveals that the experiments have slightly ‘fatter’ tails. This means that the specific details of the organisms only begin to matter when considering rather large displacements. In future work, we shall attempt to determine what is the dominant cause of large displacements in the near-field for a more realistic model of C. reinhardtii. The large displacements could arise, for instance, from the strong time-dependence of the swimming organism, or from particles ‘sticking’ to the no-slip body of the organism or to stagnation points.

We have not discussed at all the role of reorientation, that is, running-and-tumbling or orientation diffusion. Pushkin and Yeomans (2013) showed that some curvature in the paths does not influence the diffusivity very much, so it is likely not a very important factor here. In experiments involving different organisms it could matter, especially if the swimmer carries a volume of trapped fluid.

One glaring absence from the present theory is any asymmetry between pushers and pullers. This suggests that correlations between swimmers must be taken into account to see this asymmetry emerge. These correlations begin to matter as swimmer densities are increased. However, how to incorporate these correlations into a model similar to the one presented here is a challenge.

###### Acknowledgements.
The author thanks Bruno Eckhardt and Stefan Zammert for helpful discussions and for providing the digitized data from Leptos et al.  The paper also benefited from discussions with Raymond Goldstein, Eric Lauga, Kyriacos Leptos, Peter Mueller, Tim Pedley, Saverio Spagnolie, and Benedek Valko. Much of this work was completed while the author was a visiting fellow of Trinity College, Cambridge. This research was supported by NSF grant DMS-1109315.

## References

• Pedley and Kessler (1992) T. J. Pedley and J. O. Kessler, Annu. Rev. Fluid Mech. 24, 313 (1992).
• Lauga and Powers (2009) E. Lauga and T. R. Powers, Rep. Prog. Phys. 72, 096601 (2009).
• Rothschild (1963) A. J. Rothschild, Nature 198, 1221 (1963).
• Winet et al. (1984) H. Winet, G. S. Bernstein,  and J. Head, J. Reprod. Fert. 70, 511 (1984).
• Cosson et al. (2003) J. Cosson, P. Huitorel,  and C. Gagnon, Cell Motil. Cytoskel. 54, 56 (2003).
• Lauga et al. (2006) E. Lauga, W. R. DiLuzio, G. M. Whitesides,  and H. A. Stone, Biophys. J. 90, 400 (2006).
• Berke et al. (2008) A. P. Berke, L. Turner, H. C. Berg,  and E. Lauga, Phys. Rev. Lett. 101, 038102 (2008).
• Drescher et al. (2009) K. Drescher, K. C. Leptos, I. Tuval, T. Ishikawa, T. J. Pedley,  and R. E. Goldstein, Phys. Rev. Lett. 102, 168101 (2009).
• Dombrowski et al. (2004) C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein,  and J. O. Kessler, Phys. Rev. Lett. 93, 098103 (2004).
• Hernandez-Ortiz et al. (2005) J. P. Hernandez-Ortiz, C. G. Stoltz,  and M. D. Graham, Phys. Rev. Lett. 95, 204501 (2005).
• Underhill et al. (2008) P. T. Underhill, J. P. Hernandez-Ortiz,  and M. D. Graham, Phys. Rev. Lett. 100, 248101 (2008).
• Underhill and Graham (2011) P. T. Underhill and M. D. Graham, Phys. Fluids 23, 121902 (2011).
• Saintillan and Shelley (2007) D. Saintillan and M. J. Shelley, Phys. Rev. Lett. 99, 058102 (2007).
• Sokolov et al. (2009) A. Sokolov, R. E. Goldstein, F. I. Feldchtein,  and I. S. Aranson, Phys. Rev. E 80, 031903 (2009).
• Saintillan and Shelley (2012) D. Saintillan and M. J. Shelley, J. Roy. Soc. Interface 9, 571 (2012).
• Huntley and Zhou (2004) M. E. Huntley and M. Zhou, Mar. Ecol. Prog. Ser. 273, 65 (2004).
• Dewar et al. (2006) W. K. Dewar, R. J. Bingham, R. L. Iverson, D. P. Nowacek, L. C. St. Laurent,  and P. H. Wiebe, J. Mar. Res. 64, 541 (2006).
• Kunze et al. (2006) E. Kunze, J. F. Dower, I. Beveridge, R. Dewey,  and K. P. Bartlett, Science 313, 1768 (2006).
• Katija and Dabiri (2009) K. Katija and J. O. Dabiri, Nature 460, 624 (2009).
• Dabiri (2010) J. O. Dabiri, Geophys. Res. Lett. 37, L11602 (2010).
• Leshansky and Pismen (2010) A. M. Leshansky and L. M. Pismen, Phys. Rev. E 82, 025301 (2010).
• Thiffeault and Childress (2010) J.-L. Thiffeault and S. Childress, Phys. Lett. A 374, 3487 (2010)arXiv:0911.5511 .
• Lorke and Probst (2010) A. Lorke and W. N. Probst, Limnol. Oceanogr. 55, 354 (2010).
• Katija (2012) K. Katija, J. Exp. Biol. 215, 1040 (2012).
• Visser (2007) A. W. Visser, Science 316, 838 (2007).
• Gregg and Horne (2009) M. C. Gregg and J. K. Horne, J. Phys. Ocean. 39, 1097 (2009).
• Kunze (2011) E. Kunze, J. Mar. Res. 69, 591 (2011).
• Noss and Lorke (2014) C. Noss and A. Lorke, Limnol. Oceanogr. 59, 724 (2014).
• Ishikawa et al. (2010) T. Ishikawa, J. T. Locsei,  and T. J. Pedley, Phys. Rev. E 82, 021408 (2010).
• Kurtuldu et al. (2011) H. Kurtuldu, J. S. Guasto, K. A. Johnson,  and J. P. Gollub, Proc. Natl. Acad. Sci. USA 108, 10391 (2011).
• Miño et al. (2011) G. L. Miño, T. E. Mallouk, T. Darnige, M. Hoyos, J. Dauchet, J. Dunstan, R. Soto, Y. Wang, A. Rousselet,  and E. Clément, Phys. Rev. Lett. 106, 048102 (2011).
• Zaid et al. (2011) I. M. Zaid, J. Dunkel,  and J. M. Yeomans, J. Roy. Soc. Interface 8, 1314 (2011).
• Maxwell (1869) J. C. Maxwell, Proc. London Math. Soc. s1-3, 82 (1869).
• Darwin (1953) C. G. Darwin, Proc. Camb. Phil. Soc. 49, 342 (1953).
• Lighthill (1956) M. J. Lighthill, J. Fluid Mech. 1, 31 (1956), corrigendum 2, 311–312.
• Lin et al. (2011) Z. Lin, J.-L. Thiffeault,  and S. Childress, J. Fluid Mech. 669, 167 (2011)http://arxiv.org/abs/1007.1740 .
• Jepson et al. (2013) A. Jepson, V. A. Martinez, J. Schwarz-Linek, A. Morozov,  and W. C. K. Poon, Phys. Rev. E 88, 041002 (2013).
• Morozov and Marenduzzo (2014) A. Morozov and D. Marenduzzo, Soft Matter 10, 2748 (2014).
• Kasyap et al. (2014) T. V. Kasyap, D. L. Koch,  and M. Wu, Phys. Fluids 26, 081901 (2014).
• Pushkin and Yeomans (2013) D. O. Pushkin and J. M. Yeomans, Phys. Rev. Lett. 111, 188101 (2013).
• Pushkin and Yeomans (2014) D. O. Pushkin and J. M. Yeomans, J. Stat. Mech.: Theory Exp. 2014, P04030 (2014).
• Miño et al. (2013) G. L. Miño, J. Dunstan, A. Rousselet, E. Clément,  and R. Soto, J. Fluid Mech. 729, 423 (2013).
• Dunkel et al. (2010) J. Dunkel, V. B. Putz, I. M. Zaid,  and J. M. Yeomans, Soft Matter 6, 4268 (2010).
• Pushkin et al. (2013) D. O. Pushkin, H. Shum,  and J. M. Yeomans, J. Fluid Mech. 726, 5 (2013).
• Wu and Libchaber (2000) X.-L. Wu and A. Libchaber, Phys. Rev. Lett. 84, 3017 (2000).
• Leptos et al. (2009) K. C. Leptos, J. S. Guasto, J. P. Gollub, A. I. Pesci,  and R. E. Goldstein, Phys. Rev. Lett. 103, 198103 (2009).
• Eckhardt and Zammert (2012) B. Eckhardt and S. Zammert, Eur. Phys. J. E 35, 96 (2012).
• Lighthill (1952) M. J. Lighthill, Comm. Pure Appl. Math. 5, 109 (1952).
• Blake (1971) J. R. Blake, J. Fluid Mech. 46, 199 (1971).
• Ishikawa et al. (2006) T. Ishikawa, M. P. Simmonds,  and T. J. Pedley, J. Fluid Mech. 568, 119 (2006).
• Ishikawa and Pedley (2007) T. Ishikawa and T. J. Pedley, J. Fluid Mech. 588, 437 (2007).
• Thiffeault (2014) J.-L. Thiffeault, “Short-time distribution of particle displacements due to swimming microorganisms,”  (2014), arXiv:1408.4781v1.
• Gärtner (1977) J. Gärtner, Theory Probab. Appl. 22, 24 (1977).
• Ellis (1984) R. S. Ellis, Ann. Prob. 12, 1 (1984).
• Ellis (1985) R. S. Ellis, Entropy, Large Deviations, and Statistical Mechanics (Springer-Verlag, New York, 1985).
• Touchette (2009) H. Touchette, Physics Reports 478, 1 (2009).
• Thiffeault (2010) J.-L. Thiffeault, in Proceedings of the 2010 Summer Program in Geophysical Fluid Dynamics, edited by N. J. Balmforth (Woods Hole Oceanographic Institution, Woods Hole, MA, 2010) http://www.whoi.edu/main//gfd/proceedings-volumes/2010.

## Appendix A Proof of Proposition 1

###### Proof.

Observe that  as . Writing , we expand the exponent as a convergent Taylor series:

 (1−εy)1/ε =exp(−ε−1∞∑m=1(εy)mm)(converges since εy∼o(ε1/(M+1))) =exp(−ε−1(M∑m=1(εy)mm+O((εy)M+1))) =exp(−ε−1M∑m=1(εy)mm)exp(O(εMyM+1)) =exp(−ε−1M∑m=1(εy)mm)(1+o(ε0)).\qed

## Appendix B The log model

A reasonable model for objects in an inviscid fluid Thiffeault and Childress (2010); Lin et al. (2011) is the drift function

 Δλ(ρ,z)={˜Δ(ρ),if 0≤x≤λ,0,otherwise, (49)

where  is the perpendicular distance to the swimming direction. The integral (16) then simplifies to

 Γλ(k)=1vλ∫Vγ(kΔλ(η))dVη=(d−1)ℓd−1∫∞0γ(k˜Δ(ρ))ρd−2dρ. (50)

Assume a monotonic relationship between  and ; we can then change the integration variable to :

 Γλ(k)=(d−1)ℓd−1∫∞0γ(k˜Δ)ρd−2(˜Δ)d˜Δ|˜Δ′(ρ(˜Δ))|.