# Loopy Lévy flights enhance tracer diffusion in active suspensions

Brownian motion is widely used as a paradigmatic model of diffusion in equilibrium media throughout the physical, chemical, and biological sciences. However, many real world systems, particularly biological ones, are intrinsically out-of-equilibrium due to the energy-dissipating active processes underlying their mechanical and dynamical features Needleman2017 (). The diffusion process followed by a passive tracer in prototypical active media such as suspensions of active colloids or swimming microorganisms Koch:2011aa () indeed differs significantly from Brownian motion, manifest in a greatly enhanced diffusion coefficient Wu2000 (); Leptos2009 (); Mino:2011aa (); Kurtuldu:2011aa (); Mino:2013aa (); Jepson:2013aa (); Jeanneret2016 (); Kurihara2017 (), non-Gaussian tails of the displacement statistics Kurtuldu:2011aa (); Jeanneret2016 (); Kurihara2017 (), and crossover phenomena Jeanneret2016 (); Kurihara2017 () from non-Gaussian to Gaussian scaling. While such characteristic features have been extensively observed in experiments, there is so far no comprehensive theory explaining how they emerge from the microscopic active dynamics. Here we present a theoretical framework of the enhanced tracer diffusion in an active medium from its microscopic dynamics by coarse-graining the hydrodynamic interactions between the tracer and the active particles as a stochastic process. The tracer is shown to follow a non-Markovian coloured Poisson process that accounts quantitatively for all empirical observations. The theory predicts in particular a long-lived Lévy flight regime Hughes1981 () of the tracer motion with a non-monotonic crossover between two different power-law exponents. The duration of this regime can be tuned by the swimmer density, thus suggesting that the optimal foraging strategy of swimming microorganisms might crucially depend on the density in order to exploit the Lévy flights of nutrients Bartumeus2002 (). Our framework not only provides the first validation of the celebrated Lévy flight model Hughes1981 () from a physical microscopic dynamics, but can also be applied to address important conceptual questions, such as the thermodynamics of active systems KanazawaPRL2015 (), and practical ones regarding, e.g., the interaction of swimming microorganisms with nutrients and other small particles like degraded plastic Goldstein:2015aa () and the design of artificial nanoscale machines Bechinger:2016aa ().

A passive tracer immersed in a fluid at equilibrium moves randomly due to its collisions with the surrounding fluid molecules. Understanding how the observed stochastic process followed by the tracer relates to the statistical mechanics of the surrounding fluid, as accomplished in the seminal works by Einstein, Smoluchowski, and Langevin Gardiner1985 (), has provided deep insight into the connection between molecular transport and equilibrium thermodynamics , which has been widely exploited to describe soft matter and other complex physical systems Coffey2004 (). However, when either artificial self-propelled colloids or biological swimming micro-organisms, such as bacteria like Escherichia coli or algae like Volvox and Chlamydomonas reinhardtii Koch:2011aa (), are also suspended, the diffusion of the tracer changes dramatically due to the active stirring of the fluid exerted by the self-propelled particles. Indeed, the active diffusion of the tracer experimentally exhibits the following unique features that can no longer be explained as a Brownian motion: (i) the tracer exhibits dancing loopy trajectories Leptos2009 (); Jepson:2013aa (); Jeanneret2016 (); (ii) its mean square displacement (MSD) exhibits a crossover between superdiffusion with () for short times and normal diffusion () for long times, where the effective diffusion coefficient is greatly enhanced compared with that of the equilibrium case, , revealing for both three and quasi-two dimensional systems a linear dependence on the density of swimmers : , with the coefficient being system dependent Wu2000 (); Leptos2009 (); Mino:2011aa (); Mino:2013aa (); Jepson:2013aa (); Jeanneret2016 (); (iii) the probability density function (PDF) of position displacements in a given time interval exhibits strong non-Gaussian features manifest as power-law tails Kurtuldu:2011aa (); Kurihara2017 (); (iv) eventually reverts to a Gaussian shape for large Jeanneret2016 (); Kurihara2017 (); (v) the associated non-Gaussian parameter exhibits a scaling regime for large times Kurihara2017 ().

Developing a single theory that captures all features (i)–(v) has been a major challenge, due to the multi-particle hydrodynamic, and thus long-range, interactions of the tracer with the swimmers underlying its transport. While the loop-like motion (i) results from an individual scattering event of the tracer in the dipolar flow field of a single swimmer Dunkel:2010aa (); Lin:2011aa (); Mino:2013aa (), and the linear form of (ii) has been explained phenomenologically based on the active flux of the swimmers, which is defined as the product of their number density and characteristic swimming speed Mino:2011aa (); Mino:2013aa (); Jepson:2013aa (); Lin:2011aa (); Morozov:2014aa (), the statistical observations (iii)–(v) could so far not be explained consistently. The power-law tails in and their convergence to Gaussian scaling for long observation times, which is expected based on the central limit theorem (CLT) arguments Pushkin:2013aa (); Thiffeault2015 (), have been reproduced in Zaid2011 (); Zaid2016 (); Kurihara2017 () assuming a static force distribution akin to the Holtsmark theory of gravitating particles (see SI for a review) Holtsmark1919 (). However, this approach neglects any dynamics of the swimmers and is thus not sufficient to capture the enhanced diffusion observed in experiments Wu2000 (); Leptos2009 (); Mino:2011aa (); Kurtuldu:2011aa (); Mino:2013aa (); Jepson:2013aa (); Jeanneret2016 (); Kurihara2017 (). Here, we present a derivation of the stochastic process underlying the enhanced diffusion of the tracer from microscopic dynamics that is valid at all timescales. The resulting process captures all characteristic features (i)–(v) and is in excellent quantitative agreement with simulation data.

We consider a three-dimensional system composed of active particles and a passive tracer suspended in a viscous fluid in a cubic box of linear length (Fig. 2a). The active particles (swimmers) are assumed self-propelled and their motion is driven by unidirectional forces with constant amplitude Lauga2009 (). In general, the dynamics of such multi-particle systems where interactions are mediated by the fluid environment in the form of hydrodynamic forces are complex and analytically intractable. However, suspensions of micro-organisms often considered in experiments (see above) are generally characterized by (a) low Reynolds number swimming and (b) a low density of swimmers (dilute conditions, see below). In particular the dilute condition (b) allows us to neglect the mutual hydrodynamic interactions of the swimmers Drescher2011 (), thus leading to the overdamped equations of motion:

(1) |

where is the force on the tracer generated by a single swimmer and is a viscous coefficient for the passive particle. In Eq. (1), and denote the positions of the -th active particle and the passive particle, respectively, and is the constant amplitude of the active self-propulsion force with unit vector specifying the swimming direction.

The low Reynolds number condition (a) further yields a closed form expression for as the solution of the Stokes equation. For force- and torque-free swimmers, the leading-order term of this solution in a far-field expansion is a dipole force Lauga2009 (). Therefore, is specified without loss of generality as the stresslet

(2) |

for , with the difference vector , the dipole strength and the system specific cut-off . The dipole strength parameter specifies the universal features of the far-flow hydrodynamic field Lauga2009 (): denotes pusher swimmers whose flow lines are oriented outward along the direction of its velocity vector and inward laterally (e.g., E. coli Lauga2009 ()); denotes instead puller swimmers whose flow lines are oriented in the opposite directions (e.g., Chlamydomonas Lauga2009 ()). The length parameter separates the far-flow field from any near-flow field hydrodynamic contributions and hard-core interactions, and thus determines when the approximation (2) is valid. The model has one further length parameter , which can be related to the typical length scale of the swimmer Lauga2009 (). All these parameters can be determined experimentally: for E. coli and Mino:2011aa (); Drescher2011 (); for Chlamydomonas and Drescher2010 (); Zaid2011 (); Zaid2016 (); Kurihara2017 (), satisfying .

Conversely, for the interaction force is not universal but system-specific. Nevertheless, all swimmer-tracer interactions in this regime can be accurately captured by using arguments based on the CLT (see below), which do not require a detailed form of . Therefore, in our present implementation we set for simplicity in this regime, but we remark that our qualitative results are independent of this specific choice and that this approximation has also been validated experimentally Kurihara2017 ().

Numerical simulations of Eqs. (1,2) in dilute conditions reproduce all features (i)–(v) (see Fig. 2c and Fig. 3). Crucially, an inspection of time-series trajectories under such dilute conditions shows that the dynamics of the tracer can be resolved as a sequence of individual scattering events, where only two-body tracer-swimmer interactions are relevant (Fig. 2a. See also SI). In the kinetic theory of gases the dilute condition is ensured by requiring , where is the number density and the range of interparticle interactions. In our system, however, the hydrodynamic force is long-range such that an interaction range cannot be well defined. Nevertheless, we can identify with the maximum geometrical length parameter as and define dilute conditions accordingly. The motivation behind this definition of is two-fold. Firstly, with this definition the condition is indeed realized in experiments that exhibit the features (i)–(v) (e.g., in Kurihara2017 () and in Mino:2011aa ()). Secondly, this parameter regime allows for a self-consistent description of the tracer dynamics, in which the dipole interaction governs the displacement statistics on short and intermediate time scales, while the statistics for longer times reverts to a Gaussian form due to the CLT. To this end, we note that in a dilute system, at every instant in time, swimmers on average have , where is introduced as the impact parameter of a binary swimmer-tracer interaction (see Fig. 2b). The tracer statistics will then be governed by three distinct dynamical regimes: (1.) For short times , the tracer experiences the long-range forces of effectively static swimmers as in the Holtzmark theory (“Holtzmark regime”). Accordingly, . (2.) For times , the tracer is scattered by the moving swimmers in a sequence of binary interactions that we call “scatterings” (“scattering regime”). Since the system is dilute, the swimmer-tracer interaction is governed by the far-flow field as in Eq. (2) in this regime. The time scale (see Fig. 2b) thus estimates the time necessary for a swimmer to arrive close enough to the tracer to interact via hard-core and near-field hydrodynamic interactions. (3.) For , these interaction events thus represent “collisions” and the tracer is displaced by an accumulation of such collisions with the swimmers such that the CLT applies (“CLT regime”).

Remarkably, these three regimes are captured by a coarse-grained description of the tracer dynamics in terms of the Langevin equation

(3) |

where counts the number of scattering events up to a fixed time and is the force shape function (FSF) describing the force exerted on the tracer during each scattering. The transition from the fully deterministic dynamics of Eqs. (1, 2) to a stochastic description by Eq. (3) is realised by assuming to be a Poisson process with intensity . The FSF is characterized by the set of impact parameter and injection angles (denoted by in Fig. 2b) and is obtained in analytical form by solving the binary swimmer–tracer scattering problem. Remarkably, for (valid on average in dilute conditions and in the Holtzmark and scattering regimes, see above) an analytical approximation of the FSF can be obtained using a Picard iteration up to 2nd order and subsequent Taylor expansion (see SI), which is in excellent agreement with numerics (Fig. 2c). The FSF is centred at the scattering time point , which is set by the condition . The intensity can be shown to satisfy directly from its microscopic dynamics (see SI), and the total intensity diverges in the large system size limit as . Consequently, is a coloured Poisson noise with infinite intensity, also known as generalized Campbell’s process Campbell1909 (). This infinite intensity is a typical singular character of a Lévy process, and is a physical consequence of the long-range hydrodynamic interactions, which cause an infinite number of small scatterings at possibly infinite impact parameter.

We now calculate the statistics resulting from Eq. (3) for the tracer displacement in the -direction in agreement with the experimental protocols typically used Leptos2009 (); Kurtuldu:2011aa (); Jeanneret2016 (); Kurihara2017 (). Using functional techniques (see SI), we derive the displacement PDF and obtain its scaling behaviours in the regimes (1.)–(3.) as

(4) |

with power-law exponents , and positive constant that depends on (see SI). For finite a truncation appears in the power-law tails, which is realistic because any physical system must accommodate finite cutoffs Mantegna1994 ().

The Holtsmark regime (1.) yields the same scaling behaviour as found in other purely static approaches Zaid2011 (); Zaid2016 (); Kurihara2017 (). For , these theories are not applicable because the rearrangements of the active swimmers are no longer negligible. In fact, the force induced during binary scatterings has to be fully taken into account, and the stochastic process is then non-Markovian (see SI). In the scattering regime (2.), the description nevertheless becomes effectively Markovian as the FSF can be approximated effectively as a Dirac -function (see SI). In this regime the coloured Poisson model (3) is equivalent to a compound Poisson process with jump-length distribution prescribed as a truncated power-law with scaling behaviour . Our results thus validate the well known Lévy flight model Hughes1981 () as an approximate description of the tracer dynamics at this timescale. A related jump-diffusion process underlying the enhanced diffusion has previously been proposed in Jeanneret2016 () on purely phenomenological grounds. In the collision regime (3.), collisions become dominant over scatterings. Since collisional impact has finite cutoff, accumulation of a sufficient number of collisions leads to the Gaussian tail as a consequence of the central limit theorem. We note that the detailed form of for is renormalized into the variance , though it is irrelevant to the power-law tail (e.g., and ). Overall, the striking non-monotonic behaviour of the scaling exponents of the tracer displacement statistics as predicted by Eq. (4) is in excellent agreement with simulation results (Fig. 3a).

Using the exact expression for (SI), the tracer MSD for different concentrations of active swimmers can be shown to collapse onto the same universal curve upon rescaling by , which reveals a crossover from ballistic motion for short timescales to normal diffusion for longer ones (see Fig. 3b). This yields in the ballistic regime: . Likewise, in the diffusive regime we find: . Eq. (3) thus predicts the enhancement of the tracer diffusive motion at all time scales. While in the diffusive regime it captures the linear dependence of the diffusion coefficient on the active flux that has been validated experimentally Wu2000 (); Leptos2009 (); Mino:2011aa (); Mino:2013aa (); Jepson:2013aa (); Jeanneret2016 (), in the ballistic regime it predicts a dependence of the on .

Likewise, we can also calculate the non-Gaussian parameter for the tracer displacement statistics . Upon rescaling by , simulation data for different concentrations of active swimmers indeed collapse onto the same curve, that scales as for long observation times (see Fig. 3c). This scaling prediction is confirmed by the model that yields (see SI).

The tracer is also subjected to thermal fluctuations that can be included in the model (3) as Brownian motion (with variance , the Boltzmann constant and the temperature of the bath without swimmers). While the scaling predictions (4) for the tracer displacement distribution and that for the are preserved, thermal noise breaks the data collapse for both the and the and shifts the by a term with (see SI for details).

Our approach can be extended straightforwardly to any dimensional systems. Further, more general hydrodynamic force fields can be readily incorporated than Eq. (2), such as higher-order terms in the far-field expansion (e.g., the quadrupole). Let us assume the scaling relations for the far-flow hydrodynamic interaction force (or equivalently, ) and for the net tracer displacement during a scattering event (see SI). As a general recipe, our theory predicts power-law crossover with and , considering that for -dimensional active flows. This suggests that, while is universally determined by the power-law exponent of the hydrodynamic interaction force, depends crucially on its detailed form (i.e., the existence of loops).

The enhanced diffusion process as determined by (i)–(v) is a linear-in-time diffusion process with non-Gaussian displacement statistics. Such processes represent a ubiquitous feature of non-equilibrium diffusion phenomena Wang:2012aa (), but yet they pose considerable challenges already for purely phenomenological modelling approaches Chechkin:2017aa (). The fact that the stochastic process Eq. (3) captures all its characteristic features is thus remarkable, even more so since it is derived through an exact coarse-graining of microscopic dynamics. In contrast to many models that solve this “Brownian yet non-Gaussian” diffusion conundrum, Eq. (3) is easily interpretable physically due to the well-defined force pulse nature. Models of this type might thus be applicable much more generally. A striking prediction of our theory is the intermediate fatter tail of the tracer displacement distribution in the scattering regime compared with the Holtsmark regime, which is counter-intuitive since the tails are typically expected to become thinner monotonically for longer times due to the CLT. Such signature should be experimentally detectable. Moreover, our results make predictions on the possible foraging behaviour of real swimming micro-organisms like Chlamydomonas. On the one hand, diverges for , which implies that the power-law tail persists for longer timescales the more dilute the system is. On the other hand, Lévy flights have been shown to increase encounter probabilities in the stochastic search for sparsely and randomly distributed revisitable targets Viswanathan1999 (); Humphries:2012aa (). These results thus suggest that for swimming micro-organisms an optimal foraging strategy should significantly depend on the density. For large densities the displacements of nutrients are primarily Gaussian (i.e. spatially localized), which make an active search like intermittent Lévy-Brownian strategies Benichou2011 () more efficient for the forager. On the contrary, at low densities such displacements are power-law distributed, such that it might be advantageous for the forager to simply wait for a nutrient to come close dragged by the other swimmers Bartumeus2002 (). Finally, we remark that superimposing additional force fields might lead to novel mechanisms to control and exploit enhanced diffusion in artificial devices.

## References

- (1) Needleman, D. & Dogic, Z. Active matter at the interface between materials science and cell biology. Nat. Rev. Mat. 2, 17048 (2017).
- (2) Koch, D. L. & Subramanian, G. Collective hydrodynamics of swimming microorganisms: Living fluids. Annu. Rev. Fluid Mech. 43, 637–659 (2011).
- (3) Wu, X.-L. & Libchaber, A. Particle diffusion in a quasi-two-dimensional bacterial bath. Phys. Rev. Lett. 84, 3017 (2000).
- (4) Leptos, K. C., Guasto, J. S., Gollub, J. P., Pesci, A. I. & Goldstein, R. E. Dynamics of enhanced tracer diffusion in suspensions of swimming eukaryotic microorganisms. Phys. Rev. Lett. 103, 198103 (2009).
- (5) Miño, G. et al. Enhanced diffusion due to active swimmers at a solid surface. Phys. Rev. Lett. 106, 048102 (2011).
- (6) Kurtuldu, H., Guasto, J. S., Johnson, K. A. & Gollub, J. P. Enhancement of biomixing by swimming algal cells in two-dimensional films. Proc. Natl. Acad. Sci. 108, 10391–10395 (2011).
- (7) Mino, G. L., Dunstan, J., Rousselet, A., Clément, E. & Soto, R. Induced diffusion of tracers in a bacterial suspension: theory and experiments. J. Fluid Mech. 729, 423–444 (2013).
- (8) Jepson, A., Martinez, V. A., Schwarz-Linek, J., Morozov, A. & Poon, W. C. K. Enhanced diffusion of nonswimmers in a three-dimensional bath of motile bacteria. Phys. Rev. E 88, 041002 (2013).
- (9) Jeanneret, R., Pushkin, D. O., Kantsler, V. & Polin, M. Entrainment dominates the interaction of microalgae with micron-sized objects. Nat. Commun. 7, 12518 (2016).
- (10) Kurihara, T., Aridome, M., Ayade, H., Zaid, I. & Mizuno, D. Non-Gaussian limit fluctuations in active swimmer suspensions. Phys. Rev. E 95, 030601 (2017).
- (11) Hughes, B. D., Shlesinger, M. F. & Montroll, E. W. Random walks with self-similar clusters. Proc. Natl. Acad. Sci. 78, 3287–3291 (1981).
- (12) Bartumeus, F., Catalan, J., Fulco, U., Lyra, M. & Viswanathan, G. Optimizing the encounter rate in biological interactions: Lévy versus Brownian strategies. Phys. Rev. Lett. 88, 097901 (2002).
- (13) Kanazawa, K., Sano, T. G., Sagawa, T. & Hayakawa, H. Minimal model of stochastic athermal systems: Origin of non-Gaussian noise. Phys. Rev. Lett. 114, 090601 (2015).
- (14) Goldstein, R. E. Green algae as model organisms for biological fluid dynamics. Annu. Rev. Fluid Mech. 47, 343–375 (2015).
- (15) Bechinger, C. et al. Active particles in complex and crowded environments. Rev. Mod. Phys. 88, 045006 (2016).
- (16) Gardiner, C. Stochastic methods. Springer Series in Synergetics (Springer-Verlag, Berlin, 2009) (1985).
- (17) Coffey, W. T. & Kalmykov, Y. P. The Langevin equation: with applications to stochastic problems in physics, chemistry and electrical engineering (World Scientific, 2004).
- (18) Dunkel, J., Putz, V. B., Zaid, I. M. & Yeomans, J. M. Swimmer-tracer scattering at low Reynolds number. Soft Matter 6, 4268–4276 (2010).
- (19) Lin, Z., Thiffeault, J.-L. & Childress, S. Stirring by squirmers. J. Fluid Mech. 669, 167–177 (2011).
- (20) Morozov, A. & Marenduzzo, D. Enhanced diffusion of tracer particles in dilute bacterial suspensions. Soft Matter 10, 2748–2758 (2014).
- (21) Pushkin, D. O. & Yeomans, J. M. Fluid mixing by curved trajectories of microswimmers. Phys. Rev. Lett. 111, 188101 (2013).
- (22) Thiffeault, J.-L. Distribution of particle displacements due to swimming microorganisms. Phys. Rev. E 92, 023023 (2015).
- (23) Zaid, I. M., Dunkel, J. & Yeomans, J. M. Lévy fluctuations and mixing in dilute suspensions of algae and bacteria. J. Royal Soc. Interface 8, 1314–1331 (2011).
- (24) Zaid, I. & Mizuno, D. Analytical limit distributions from random power-law interactions. Phys. Rev. Lett. 117, 030602 (2016).
- (25) Holtsmark, J. Über die verbreiterung von Spektrallinien. Ann. Phys. (Berl.) 363, 577–630 (1919).
- (26) Lauga, E. & Powers, T. R. The hydrodynamics of swimming microorganisms. Rep. Prog. Phys. 72, 096601 (2009).
- (27) Drescher, K., Dunkel, J., Cisneros, L. H., Ganguly, S. & Goldstein, R. E. Fluid dynamics and noise in bacterial cell–cell and cell–surface scattering. Proc. Natl. Acad. Sci. 108, 10940–10945 (2011).
- (28) Drescher, K., Goldstein, R. E., Michel, N., Polin, M. & Tuval, I. Direct measurement of the flow field around swimming microorganisms. Phys. Rev. Lett. 105, 168101 (2010).
- (29) Campbell, N. The study of discontinuous phenomena. In Proc. Camb. Philos. Soc., vol. 15, 117–136 (1909).
- (30) Mantegna, R. N. & Stanley, H. E. Stochastic process with ultraslow convergence to a Gaussian: the truncated Lévy flight. Phys. Rev. Lett. 73, 2946 (1994).
- (31) Wang, B., Kuo, J., Bae, S. C. & Granick, S. When Brownian diffusion is not Gaussian. Nature Mater. 11, 481–485 (2012).
- (32) Chechkin, A. V., Seno, F., Metzler, R. & Sokolov, I. M. Brownian yet non-Gaussian diffusion: From superstatistics to subordination of diffusing diffusivities. Phys. Rev. X 7, 021002 (2017).
- (33) Viswanathan, G. M. et al. Optimizing the success of random searches. Nature 401, 911 (1999).
- (34) Humphries, N. E., Weimerskirch, H., Queiroz, N., Southall, E. J. & Sims, D. W. Foraging success of biological Lévy flights recorded in situ. Proc. Natl. Acad. Sci. 109, 7169–7174 (2012).
- (35) Bénichou, O., Loverdo, C., Moreau, M. & Voituriez, R. Intermittent search strategies. Rev. Mod. Phys. 83, 81 (2011).

###### Acknowledgements.

We appreciate D. Mizuno, H. Takayasu, M. Takayasu, H. Hayakawa, and F. van Wijland for fruitful discussions. This work was supported by Grant-in-Aid for JSPS Fellows (Grant No. 16J05315), JSPS KAKENHI (Grant Nos. 16K16016 and 18K13519), the Research Fellowship granted by the Royal Commission for the Exhibition of 1851, and Atoms program granted by the Yukawa Institute for Theoretical Physics. The numerical calculations were carried out on XC40 at Yukawa Institute for Theoretical Physics in Kyoto University.## Author contributions

KK conceived the original idea. KK, TGS, AC and AB designed research. KK and AC performed analytical calculations. TGS performed numerical simulations. KK, AC, and AB wrote the paper.

## Competing interests

The authors declare no competing financial interests.

## Data availability

We will provide all the numerical codes and plot files upon requests.