Interplay between Brownian and hydrodynamic tracer diffusion in suspensions of swimming microorganisms

TL;DR

This study investigates how thermal Brownian motion and hydrodynamic advection influence tracer diffusion in suspensions of microswimmers, revealing that Brownian effects are negligible at typical biological swimming speeds, contrary to prior beliefs.

Contribution

The paper demonstrates through simulations and theory that Brownian motion's impact on active diffusion is minimal at realistic swimming speeds, challenging previous assumptions.

Findings

Brownian motion has negligible effect on active diffusion at typical biological speeds.

Active and Brownian diffusion effects are additive for moderate and high swimmer speeds.

Reinterpretation of experimental data on tracer diffusion in bacterial suspensions is suggested.

Abstract

The general problem of tracer diffusion in non-equilibrium baths is important in a wide range of systems, from the cellular level to geographical lengthscales. In this paper, we revisit the archetypical example of such a system: a collection of small passive particles immersed in a dilute suspension of non-interacting dipolar microswimmers, representing bacteria or algae. In particular, we consider the interplay between thermal (Brownian) diffusion and hydrodynamic (active) diffusion due to the persistent advection of tracers by microswimmer flow fields. Previously, it has been argued that even a moderate amount of Brownian diffusion is sufficient to significantly reduce the persistence time of tracer advection, leading to a significantly reduced value of the effective active diffusion coefficient $D_A$ compared to the non-Brownian case. Here, we show by large-scale simulations and kinetic theory that this effect is in fact only practically relevant for microswimmers that effectively remain stationary while still stirring up the surrounding fluid, so-called \emph{shakers}. In contrast, for moderate and high values of the swimming speed $v_s$, relevant for biological microswimmer suspensions, the effect of Brownian motion on $D_A$ is negligible, leading to the effects of advection by microswimmers and Brownian motion being additive. This conclusion contrasts with previous results from the literature, and encourages a reinterpretation of recent experimental measurements of $D_A$ for tracer particles of varying size in bacterial suspensions.

Equations

$$D_A\approx \frac{7{\kappa }^{2}n}{2048\lambda \epsilon +336\pi v_s}.$$

$$\mathbf{U}(\mathbf{r},t)=\sum _{i=1}^N{\mathbf{u}}_s(\mathbf{r};{\mathbf{r}}_i,{\mathbf{p}}_i),$$

$$D_A=\frac{1}{3}{\int }_0^{\infty }⟨{\dot{\mathbf{r}}}_T(t)\cdot {\dot{\mathbf{r}}}_T(0)⟩dt=\frac{1}{3}{\int }_0^{\infty }⟨\mathbf{U}({\mathbf{r}}_T,t)\cdot \mathbf{U}({\mathbf{r}}_T,0)⟩dt\equiv \frac{1}{3}{\int }_0^{\infty }{C}_T(t)dt,$$

$${\dot{\mathbf{r}}}_T=\mathbf{U}({\mathbf{r}}_T)+\sqrt{2D_0}\mathbf{\eta },$$

$$D_{\mathrm{Tot}}=D_0+\frac{1}{3}{\int }_0^{\infty }{C}_T(t)dt=D_0+D_A,$$

$$\mathrm{Pe}\equiv \frac{D_A(D_0=0)}{D_0},$$

$$\xi \equiv \frac{D_A(D_0)}{D_A(D_0=0)}.$$

$$L\equiv \frac{v_s}{\epsilon \lambda }.$$

$${\mathbf{u}}_s(\mathbf{r};{\mathbf{r}}_i,{\mathbf{p}}_i)=\frac{\kappa }{8\pi }\left[3\frac{({\mathbf{p}}_i\cdot {\mathbf{r}}^{\prime }{)}^2{\mathbf{r}}^{\prime }+{\epsilon }^2({\mathbf{p}}_i\cdot {\mathbf{r}}^{\prime }){\mathbf{p}}_i}{(r^{\prime 2}+{\epsilon }^2{)}^{5/2}}-\frac{{\mathbf{r}}^{\prime }}{(r^{\prime 2}+{\epsilon }^2{)}^{3/2}}\right],$$

$$C_U(t)=\frac{{\kappa }^{2}n}{15{\pi }^2\epsilon }{\int }_0^{\infty }{A}^2(\zeta )e^{-\tau }\mathcal{F}(L\zeta \tau )d\zeta ,$$

$$A(x)=\frac{1}{2}{x}^2{K}_2(x),$$

$$\mathcal{F}(x)=15\frac{(5x^2-12)\sin x-x(x^2-12)\cos x}{x^5},$$

$$D_A=\frac{{\kappa }^{2}n}{45{\pi }^2\lambda \epsilon }{\int }_0^{\infty }{A}^2(\zeta )\mathcal{G}(L\zeta )d\zeta ,$$

$$\mathcal{G}(x)=\frac{5}{2}\frac{3x+2x^3-3(1+x^2)\arctan x}{x^5},$$

$${\dot{\mathbf{r}}}_i=v_s{\mathbf{p}}_i+\sqrt{2D_0}\mathbf{\eta },$$

$$C_U(t;D_0)=\frac{{\kappa }^{2}n}{15{\pi }^2\epsilon }{\int }_0^{\infty }{A}^2(\zeta )e^{-(1+\tilde{D}{\zeta }^2)\tau }\mathcal{F}(L\zeta \tau )d\zeta ,$$

$$D_A(D_0)=\frac{{\kappa }^{2}n}{45{\pi }^2\lambda \epsilon }{\int }_0^{\infty }\frac{A^2(\zeta )}{1+\tilde{D}{\zeta }^2}\mathcal{G}\left(\frac{L\zeta }{1+\tilde{D}{\zeta }^2}\right)d\zeta .$$

$$C_U(t)=\underset{t^{\prime }\to \infty }{\lim }\frac{1}{V}\int d\mathbf{r}\overline{U^{\alpha }\left(\mathbf{r},t^{\prime }\right)U^{\alpha }\left(\mathbf{r},t+t^{\prime }\right)},$$

$$U^{\alpha }\left(\mathbf{r},t\right)=\sum _{i=1}^N{u}_s^{\alpha }(\mathbf{r};{\mathbf{r}}_i(t),{\mathbf{p}}_i(t)).$$

$$\displaystyle\partial_{t}F_{N}+\sum_{i=1}^{N}\Big{(}v_{s}p_{i}^{\alpha}\nabla_{i}^{\alpha}-D_{0}\nabla_{i}^{2}\Big{)}F_{N}=-N\lambda F_{N}+\frac{\lambda}{4\pi}\sum_{i=1}^{N}\int\mathrm{d}{\bm{p}}_{i}F_{N},$$

$${\partial }_{t}h+v_s{p}^{\alpha }{\nabla }^{\alpha }h-D_0{\nabla }^{2}h+\lambda h-\frac{\lambda }{4\pi }\int d\mathbf{p}h=\chi (\mathbf{r},\mathbf{p},t),$$

$$⟨\chi (\mathbf{r},\mathbf{p},t)⟩=0,$$

$$\displaystyle\langle\chi({\bm{r}},{\bm{p}},t)\chi({\bm{r}}^{\prime},{\bm{p}}^{\prime},t^{\prime})\rangle=\frac{n}{4\pi}\delta(t-t^{\prime})\Bigg{[}2\lambda\delta({\bm{r}}-{\bm{r}}^{\prime})\left(\delta({\bm{p}}-{\bm{p}}^{\prime})-\frac{1}{4\pi}\right)$$

$$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-D_{0}\delta({\bm{p}}-{\bm{p}}^{\prime})\Big{(}\nabla^{2}+\nabla^{\prime 2}\Big{)}\delta({\bm{r}}-{\bm{r}}^{\prime})\Bigg{]}.$$

$$C_U(t)=\underset{t^{\prime }\to \infty }{\lim }\frac{1}{V}\int d\mathbf{r}\int d{\mathbf{r}}^{\prime }d{\mathbf{r}}^{\prime \prime }d{\mathbf{p}}^{\prime }d{\mathbf{p}}^{\prime \prime }{u}_s^{\alpha }(\mathbf{r};{\mathbf{r}}^{\prime },{\mathbf{p}}^{\prime })u_s^{\alpha }(\mathbf{r};{\mathbf{r}}^{\prime \prime },{\mathbf{p}}^{\prime \prime })$$

$$\times ⟨h({\mathbf{r}}^{\prime },{\mathbf{p}}^{\prime },t)h({\mathbf{r}}^{\prime \prime },{\mathbf{p}}^{\prime \prime },t+t^{\prime }){⟩}_{\chi },$$

$$\hat{h}(\mathbf{k},\mathbf{p},s)={\int }_0^{\infty }dt{e}^{-st}\int d\mathbf{r}{e}^{i\mathbf{k}\cdot \mathbf{r}}h(\mathbf{r},\mathbf{p},t),$$

$$\hat{h}(\mathbf{k},\mathbf{p},s)=\frac{\hat{\chi }(\mathbf{k},\mathbf{p},s)}{\sigma (\mathbf{k},\mathbf{p},s)}+\frac{\lambda }{4\pi \sigma (\mathbf{k},\mathbf{p},s)}\frac{\int d{\mathbf{p}}^{\prime }\frac{\hat{\chi }(\mathbf{k},{\mathbf{p}}^{\prime },s)}{\sigma (\mathbf{k},{\mathbf{p}}^{\prime },s)}}{1-\frac{\lambda }{4\pi }\int \frac{d{\mathbf{p}}^{\prime }}{\sigma (\mathbf{k},{\mathbf{p}}^{\prime },s)}}.$$

$$C_U(t)=\frac{n{\kappa }^2}{16{\pi }^4}\underset{t^{\prime }\to \infty }{\lim }{\mathcal{L}}_{s_1,t^{\prime }}^{-1}{\mathcal{L}}_{s_2,t^{\prime }+t}^{-1}\int d\mathbf{k}\frac{A^2(kϵ)}{k^4}\left(\lambda +D_0{k}^2\right)$$

$$\times \int d\mathbf{p}(\mathbf{k}\cdot \mathbf{p}{)}^2\left[1-\frac{(\mathbf{k}\cdot \mathbf{p}{)}^2}{k^2}\right]\frac{1}{s_1+s_2}$$

$$\times \frac{1}{\lambda +s_1+D_0{k}^2+i{v}_s(\mathbf{k}\cdot \mathbf{p})}\frac{1}{\lambda +s_2+D_0{k}^2-i{v}_s(\mathbf{k}\cdot \mathbf{p})}.$$

Full text

1 Introduction

Understanding the mass transport of colloidal and molecular species in non-equilibrium environments is crucial for various processes, ranging from active intracellular transport (Mogre et al., 2020) to the dispersion of nutrients in world oceans (Katija, 2012). Apart from its practical importance, the transport properties of tracer particles in generic “active baths” has attracted much interest from a statistical physics perspective, where they can be viewed as a minimal example of particles driven by external, non-equilibrium noise (Argun et al., 2016; Park et al., 2020). Beyond the level of tracer particles driven by generic non-equilibrium noise, the archetypical example of a tracer particle in an active bath is a collection of point-like tracers being advected by a set of microswimmers such as bacteria or algae (Lauga & Powers, 2009). When swimming through a viscous fluid, these swimmers create long-ranged flow fields that advect the tracers, leading to tracer dynamics that is ballistic at short times and diffusive over timescales longer than the autocorrelation time of the local flow field (Lin et al., 2011). Realisations of this system have been extensively studied both experimentally, typically in suspensions of E. coli bacteria (Wu & Libchaber, 2000; Drescher et al., 2011; Jepson et al., 2013; Kim & Breuer, 2004; Koumakis et al., 2013; Miño et al., 2013, 2011; Patteson et al., 2016; Peng et al., 2016; Semeraro et al., 2018) or Chlamydomonas algae (Leptos et al., 2009; Ortlieb et al., 2019; Yang et al., 2016; von Rüling et al., 2021), and theoretically, with microswimmers typically being modelled either as force dipoles acting on the surrounding fluid (Morozov & Marenduzzo, 2014; Pushkin et al., 2013; Pushkin & Yeomans, 2013; Nordanger et al., 2022), as spherical “squirmers” with an imposed slip velocity along their body (Lin et al., 2011; Thiffeault & Childress, 2010; Thiffeault, 2015), or as needle-shaped “slender swimmers” with imposed stresses along their body lengths (Krishnamurthy & Subramanian, 2015; Saintillan & Shelley, 2012). While the details of these three microswimmer models differ, the results regarding enhanced tracer diffusion are largely generic and consistent with experimental results, which have shown the swimmer-induced, hydrodynamic diffusivity $D_{A}$ to scale linearly with microswimmer density $n$ in the dilute limit where swimmer-swimmer correlations can be neglected (Lin et al., 2011; Thiffeault & Childress, 2010; Miño et al., 2013). In this limit, a fruitful way of calculating $D_{A}$ is to consider the net displacement due to binary swimmer-tracer scattering events (Pushkin & Yeomans, 2013; Morozov & Marenduzzo, 2014); two examples of resulting (deterministic) tracer trajectories for scattering events are shown in Fig. 1. For a non-tumbling swimmer, starting and ending at $x=\pm\infty$, the resulting loop is closed, leading to a vanishing tracer net displacement $\Delta$ (Fig. 1a). For tumbling swimmers with a finite persistence length, the trajectory is however punctuated mid-way through the tracer loop, leading to significantly larger values of $\Delta$ (Fig. 1b). The resulting value of $D_{A}$ due to a large set of such scattering events can then be obtained by explicitly summing over all possible sets of scattering parameters.

In spite of the dependence of $\Delta$ on the microswimmer tumbling rate $\lambda$ for single scattering events such as that in Fig. 1, Pushkin & Yeomans (2013) showed that, in the limit of large swimming speeds $v_{s}$, $D_{A}$ is in fact independent of $\lambda$ when summed over all possible swimmer trajectories. This result was later generalised by Škultéty et al. (2020) to arbitrary swimming speeds, leading to the following approximate expression for $D_{A}$:

[TABLE]

Here, $\kappa$ is the magnitude of the microswimmer dipole (in units of volume over time) and $\varepsilon$ is a characteristic linear size of the microswimmer, which we take to be equal to the short-range regularisation length of the dipolar flow field; unless stated otherwise, we will in the following use $\kappa$ and $\varepsilon$ to non-dimensionalise the numerical data.

Rather than using the scattering approach outlined above, Eq. (1) was derived by formulating a kinetic theory for the spatiotemporal correlations of the disturbance velocity field ${\bm{U}}$ created by the swimmers. Due to the linearity of Stokes flow, we can write ${\bm{U}}$ as the superposition of the individual swimmer flow fields ${\bm{u}}_{s}$:

[TABLE]

where ${\bm{r}}_{i}$ and ${\bm{p}}_{i}$ is, respectively, the position and orientation of swimmer $i$. Knowing the statistical properties of ${\bm{U}}$, $D_{A}$ can readily be calculated via the Green-Kubo relation

[TABLE]

where, in the second equality, we have assumed point-like tracers advected by the disturbance flow, so that $\dot{{\bm{r}}}_{T}={\bm{U}}({\bm{r}}_{T})$, and the third equality defines the velocity autocorrelation function $C_{T}(t)$ in the co-moving tracer frame. While yielding identical results for dilute suspensions as the scattering approach discussed above, kinetic theories are however more readily extended to accommodate the effect of swimmer-swimmer correlations due to the mutual advection and reorientation of swimmers (Škultéty et al., 2020). Importantly, these interactions break the symmetry between rear-actuated (pusher) swimmers, such as most bacteria, and front-actuated (puller) ones, such as Chlamydomonas, leading to a super-linear scaling of $D_{A}$ with $n$ for pushers and a corresponding sub-linear scaling for pullers (Stenhammar et al., 2017).

Equation (1) shows two qualitatively different regimes for high and low $v_{s}$: For $v_{s}\rightarrow 0$ – the so-called shaker limit – the dominant mechanism controlling the decorrelation of $C_{T}(t)$ is tumbling of the swimmer. For fast swimmers, with $v_{s}\gg\lambda\varepsilon$, the decorrelation of the fluid velocity is instead dominated by the swimmer self-propulsion and thus independent of $\lambda$, and Eq. (1) reduces to the expression derived by Pushkin & Yeomans (2013). A third, somewhat less explored, mechanism affecting $C_{T}(t)$ is Brownian translational diffusion of the tracer: even though Brownian diffusion does not affect the statistics of the flow field ${\bm{U}}({\bm{r}},t)$ as measured in the lab frame, the positional noise of the tracer particle will cause it to cross the streamlines of the disturbance flow, thus perturbing its trajectory compared to the athermal case shown in Fig. 1 and leading to a lower $D_{A}$. Since the Brownian diffusion constant $D_{0}$ depends inversely on the tracer radius per the Stokes-Einstein relation, the magnitude of this effect is expected to be significant primarily for small tracer particles, and it has been hypothesised to explain the non-monotonic size dependence of enhanced tracer diffusion seen in experiments with colloids in E. coli suspensions, where Patteson et al. (2016) observed a maximum in $D_{A}$ for a tracer radius of approximately 5 $\mu$m. The effect of Brownian motion on enhanced tracer diffusion was furthermore theoretically analysed by Kasyap et al. (2014) for a model of slender-body swimmers, showing that $D_{A}$ is a non-monotonic function of $D_{0}$, with $D_{A}$ first showing a small increase for intermediate $D_{0}$, before falling below the athermal value as $D_{0}$ grows larger. In a more recent study of finite-size spherical tracers in microswimmer suspensions, Dyer & Ball (2021) numerically analysed the combined effect of thermal fluctuations and near-field flows on the size-dependent tracer dynamics, finding a similar non-monotonic behaviour as observed experimentally. In this Paper, we will revisit the problem of the interplay between Brownian and hydrodynamic diffusion for the simple case of point-like tracers immersed in a dilute suspension of microswimmers described via a regularised dipolar flow field. Using kinetic theory and large-scale lattice Boltzmann simulations of E. coli-like suspensions, we show that the effect of Brownian diffusion on active diffusion is only practically relevant whenever $v_{s}<\lambda\varepsilon$, which corresponds to extremely slow (or frequently tumbling) swimmers. For biologically relevant values of $v_{s}$ and $\lambda$, swimming is instead the dominant decorrelation mechanism, so that $D_{A}$ becomes independent of both $D_{0}$ and $\lambda$. In contrast to previous studies, our results thus indicate that the effect of Brownian motion on the enhanced diffusion is in fact negligible for most microswimmer realisations, and thus that the experimentally observed non-monotonic size dependence of $D_{A}$ on tracer size reported by Patteson et al. (2016) has other explanations.

2 Model and Method

We consider a collection of $N$ non-interacting microswimmers at number density $n=N/V$ moving through a three-dimensional viscous fluid of viscosity $\mu$. Each microswimmer is composed of two equal and opposite point forces of magnitude $F$ separated by a length $\ell$ and swims with a constant speed $v_{s}$. The resulting reduced hydrodynamic dipole strength is $\kappa=F\ell/\mu$. The swimming direction ${\bm{p}}_{i}$ furthermore relaxes through Poisson-distributed random tumbles with uncorrelated directions occurring with average frequency $\lambda$.

The position ${\bm{r}}_{T}$ of a point-like tracer obeys the equation of motion

[TABLE]

where $\bm{\eta}$ is a unit-variance white noise, $\delta$-correlated in time, and $D_{0}$ is the Brownian diffusion constant. Thus, the effect of Brownian motion is fully contained in the tracer dynamics, while we assume the effects of thermal fluctuations on the fluid and on the pairwise swimmer-tracer dynamics to be subdominant. The fluid disturbance velocity ${\bm{U}}({\bm{r}}_{T})$ due to the presence of all swimmers is given by Eq. (2) and can, in principle, be explicitly summed up on each timestep. However, to avoid the (prohibitively costly) pairwise summation over all swimmers and tracers, we instead numerically solve for the flow field using an efficient point-force implementation of the lattice Boltzmann (LB) method described previously (Bárdfalvy et al., 2020; Nash et al., 2008). Our simulations comprise a system with periodic boundaries and a size of $100^{3}$ lattice units. In LB units, set by the LB lattice spacing $\Delta l$ and time step $\Delta t$, the microswimmer density was kept constant at $n=0.01$ unless otherwise stated, corresponding to $N=10^{4}$ microswimmers. Furthermore, $N_{T}=5\times 10^{4}$ point tracers were included for statistical averaging. In addition to the LB simulations, in Section 3 we will furthermore extend the kinetic theory developed previously by Škultéty et al. (2020) to the case of microswimmers undergoing Brownian diffusion and show that, in the limit of non-interacting swimmers, the derived expression for $D_{A}$ is equivalent to that for Brownian tracers in a suspension of non-Brownian microswimmers.

Inserting Eq. (4) into the Green-Kubo relation (3) yields

[TABLE]

where $C_{T}$ was defined in Eq. (3). Thus, to obtain $D_{A}$, we numerically evaluate the time correlation of the disturbance velocity measured in the co-moving tracer frame. Since the tracer position ${\bm{r}}_{T}$ depends on $D_{0}$, $C_{T}(t)$, and thus $D_{A}$, will depend implicitly on $D_{0}$.

To characterise the system, we will use three dimensionless quantities. Firstly, we define the Péclet number, which measures the relative importance of active and thermal forces, as

[TABLE]

where $D_{A}(D_{0}=0)$ is the active diffusivity of a tracer immersed in an equivalent microswimmer suspension but in the absence of Brownian tracer motion. It should be noted that our definition of Pe is qualitatively different from that of Kasyap et al. (2014), who, instead of $D_{A}(D_{0}=0)$, use the swimming speed $v_{s}$ to characterise the active forces. We however argue that the activity experienced by the tracers depend on the magnitude of the velocity fields generated by the swimmers, and is thus dependent on $\kappa$ and encoded in $D_{A}(D_{0}=0)$, while $v_{s}$ is instead a measure of the swimmers’ self-propulsion. For experimental realisations of microswimmers, $\kappa$ and $v_{s}$ are directly proportional to each other; however, the specific relation between between them will nevertheless be specific to each type (or species) of swimmer, and decoupling them conveniently enables us to study separately the effects of self-propulsion and fluid forcing, as we demonstrate further below.

Secondly, we measure the change in active diffusion due to Brownian motion through the quantity

[TABLE]

In the limit $\mathrm{Pe}\rightarrow\infty$, where Brownian motion becomes negligible, we thus expect that $\xi\rightarrow 1$ as $D_{A}$ approaches its non-Brownian value. Finally, in accordance with Škultéty et al. (2020), we account for the effect of microswimmer self-propulsion using the reduced swimmer persistence length $L$, defined by

[TABLE]

3 Kinetic theory

In this Section, we will outline the main steps in the derivation of $D_{A}$ for a suspension of Brownian tracer particles immersed in a dilute microswimmer suspension, whose dynamics are governed by Eq. (4). Just as in our previous works (Škultéty et al., 2020; Stenhammar et al., 2017) , we describe the flow field measured at ${\bm{r}}$ due to a swimmer with position ${\bm{r}}_{i}$ and orientation ${\bm{p}}_{i}$ by a regularised dipolar flow field ${\bm{u}}_{s}({\bm{r}})$

[TABLE]

where ${\bm{r}}^{\prime}={\bm{r}}-{\bm{r}}_{i}$, $r^{\prime}=|{\bm{r}}^{\prime}|$, and $\varepsilon$ is the regularisation length. Our starting point is the derivation of Škultéty et al. (2020), where we formulated and solved a kinetic theory describing the fluctuations of the velocity field ${\bm{U}}({\bm{r}},t)$ due to a superposition of single-swimmer flow fields. In the limit of non-interacting swimmers, which is the case that we consider here, the temporal correlations of the steady-state velocity field ${\bm{U}}$ measured in the lab frame, $C_{U}(t)\equiv\langle{\bm{U}}({\bm{r}},t)\cdot{\bm{U}}({\bm{r}},0)\rangle$, is given by

[TABLE]

where

[TABLE]

with $K_{2}$ being the modified Bessel function of the second kind, and

[TABLE]

defined such that $\mathcal{F}(0)=1$. We furthermore used the dimensionless variables $L=v_{s}/(\lambda\varepsilon)$, $\tau=t\lambda$, and $\zeta=k\varepsilon$, where $k=|{\bm{k}}|$ is the wavevector magnitude. Equation (10) can equivalently be expressed in closed form in terms of elliptic integrals, as given by Eq. (72) of Škultéty et al. (2020). Equation (10) contains decorrelation of the flow field due to two separate mechanisms: exponential decay of $C_{U}(t)$ due to tumbling, and a more complex, oscillatory behaviour due to swimming, encoded in the function $\mathcal{F}$.

To obtain the hydrodynamic diffusivity of a passive tracer, Škultéty et al. (2020) used a stationary tracer approximation, implying that the tracer advection by the swimmer flow field is negligible compared to the self-propulsion of the swimmer. This implies that ${\bm{r}}_{T}$ remains effectively constant over the time it takes for ${\bm{U}}$ to relax, such that $C_{T}(t)=\langle{\bm{U}}({\bm{r}}_{T}[t],t)\cdot{\bm{U}}({\bm{r}}_{T}[t=0],0)\rangle\approx\langle{\bm{U}}({\bm{r}}_{T}[t=0],t)\cdot{\bm{U}}({\bm{r}}_{T}[t=0],0)\rangle=C_{U}(t)$, where, in the last equality, we have made the additional assumption that the tracers are homogeneously distributed in space so that they sample an unweighted spatial average of the flow field. Thus, if we can replace the correlation function $C_{T}$ in the co-moving tracer frame with that in the stationary lab frame, $C_{U}$, we can insert Eq. (10) into the Green-Kubo relation (3) and integrate over time to yield the following expression for $D_{A}$, identical to Eq. (85) of Škultéty et al. (2020):

[TABLE]

where

[TABLE]

defined such that $\mathcal{G}(0)=1$. By matching the asymptotic behaviours for $L\rightarrow 0$ and $L\rightarrow\infty$, Equation (13) can furthermore be approximated by the simple expression given in Eq. (1)

The generalisation of Eq. (13) to the case of Brownian tracers might seem straightforward, but unfortunately is not: Since the approximation $C_{T}(t)=C_{U}(t)$ amounts to the tracer remaining effectively stationary during a swimmer-tracer scattering event, this approximation will, by construction, not capture any effects on $D_{A}$ coming from Brownian diffusion across streamlines. This can easily be realised by noticing that $C_{U}(t)$ is solely a property of the swimmer suspension, and will be strictly unaffected by the tracer dynamics; thus, $D_{A}$ in Eq. (13) remains unaffected by the inclusion of tracer diffusion. Instead of going beyond the stationary tracer approximation, we circumvent this problem by noticing that, in a suspension of non-interacting microswimmers, the dynamics of a Brownian tracer will be statistically identical to that of a non-Brownian tracer in a suspension of Brownian swimmers with the same translational diffusivity $D_{0}$. This is because the single-swimmer flow field in Eq. (9) solely depends on the separation vector ${\bm{r}}_{T}-{\bm{r}}_{i}$, implying that the flow field experienced by a diffusing tracer (i.e., noise acting on ${\bm{r}}_{T}$) is identical to that experienced by a non-Brownian tracer sampling the flow field from a swimmer with the same noise instead applied to ${\bm{r}}_{i}$. In a non-interacting microswimmer suspension this equivalence is exact as long as the noise has zero mean and identical spectral properties, and we verify it numerically in Fig. 3. It however breaks down as soon as swimmer-swimmer correlations become significant since swimmer diffusion will then affect the magnitude of such correlations, which tracer diffusion will not.

Thus, as outlined in Appendix A, we instead calculate $C_{U}(t)$ for the case of a suspension of diffusing swimmers, described by the dynamics

[TABLE]

including the same tumbling dynamics as before. This yields the following generalisation of Eq. (10):

[TABLE]

where we have additionally defined the non-dimensional diffusivity $\tilde{D}=D_{0}/(\lambda\varepsilon^{2})$. Since the effect of Brownian diffusion is now fully incorporated into the properties of $C_{U}$, we again use the stationary-tracer approximation and insert this expression into the Green-Kubo relation (3), leading to the following expression for $D_{A}$ in the presence of Brownian diffusion:

[TABLE]

In Section 4, we numerically evaluate Eqs. (16) and (17) and compare the results with direct numerical simulations of microswimmer suspensions.

4 Results and Discussion

In Fig. 2 we begin by verifying Eq. (13) for the hydrodynamic diffusion coefficient $D_{A}$ in the limit $D_{0}=0$. In Fig. 2a, we demonstrate the $\lambda^{-1}$ dependence of $D_{A}$ in the shaker limit $v_{s}=0$, while Fig. 2b shows its more complex dependence on $v_{s}$ for constant $\lambda$. Apart from numerically verifying the kinetic theory expression (13), these results illustrate how $D_{A}$ decreases abruptly due to the temporal decorrelation of the flow field ${\bm{U}}({\bm{r}}_{T},t)$ induced respectively by tumbling and swimmer self-propulsion. The slight deviation between the theoretical curve and simulation results at small $v_{s}$ in Fig. 2b is likely due to the specific form for the short-range regularisation, which becomes important as $v_{s}\rightarrow 0$. In our derivation of Eq. (13) we use the regularised flow field (9), based on the regularisation first introduced by Cortez et al. (2005). In the LB simulations we instead use a numerical interpolation scheme based on a regularisation of the $\delta$ function (Peskin, 2002) acting separately on the two point forces that make up each microswimmer. Unlike the expression in (9), this numerical regularisation does not allow a direct mapping (or adjustment) of the regularisation length $\varepsilon$. We thus do not expect perfect agreement between kinetic theory and simulation in the low-$v_{s}$ regime where the short-range regularisation becomes important, and therefore treat $\varepsilon$ as a fitting parameter when comparing data from LB simulations with kinetic theory predictions. However, we find that the fitted value of $\varepsilon$ only varies slightly ($\varepsilon\in[1.9\Delta l,2.5\Delta l]$) for the values of $v_{s}$ used throughout this work, in good accordance with the fact that the regularised $\delta$ function is interpolated over a support of $2\Delta l$ in each Cartesian direction; for a more in-depth discussion of the effect of the interpolation scheme on the tracer dynamics, see de Graaf & Stenhammar (2017).

In Fig. 3, we study the additional effect of varying the Brownian diffusion coefficient $D_{0}$, as encoded in Eq. (17). From the data in Fig. 3a, it is clear that, for large enough $D_{0}$, the active diffusivity $D_{A}$ decreases compared to its non-Brownian value. To enable an easier analysis of the effect of varying swimming speed, in panel (b) we present the same data instead plotted as a function of the reduced variables $\xi$ and Pe. For shakers with $L=0$ (blue curve in Fig. 3b), $D_{A}$ is reduced compared to its non-Brownian value ($\xi<1$) as soon as $\mathrm{Pe}<1$, reaching a value as low as $\xi=0.2$ for $\mathrm{Pe}\approx 10^{-3}$. For finite values of $v_{s}$, this effect on $D_{A}$ however occurs for gradually lower values of Pe; for the fastest swimmers considered here, with $L=4.0$, no significant reduction of $D_{A}$ is observed even for Pe as low as $10^{-3}$. Instead, we observe a small but significant increase in the active diffusion compared to its non-Brownian value, in accordance with what was previously observed for slender swimmers by Kasyap et al. (2014); we discuss this effect further below. Crucially, a reduced persistence length $L=4$ nevertheless corresponds to relatively slow swimming from a biological perspective: According to the approximate calculation in Škultéty et al. (2020), $L$ for E. coli bacteria lies somewhere in the range between 5 and 20, indicating that the effect of Brownian motion on active diffusion is likely negligible in suspensions of swimming bacteria due to their fast self-propulsion. In our LB simulations, studying values higher than $L\approx 4$ is challenging, as these large swimming speeds both require very large systems to avoid significant finite-size effects and yields artifacts due to the effect of finite Reynolds number in the swimmer-tracer scattering dynamics (de Graaf & Stenhammar, 2017). We nevertheless numerically studied $\xi(\mathrm{Pe})$ using Eq. (17) for larger values of $L$, verifying that both the peak and the subsequent decrease in $\xi$ continues to move to even lower values of Pe as $L$ is increased.

To further understand the mechanism behind the reduction in $D_{A}$ with $D_{0}$, we consider the two autocorrelation functions $C_{T}(t)$ and $C_{U}(t)$, which respectively measure the fluid autocorrelation in the co-moving tracer frame and in the lab frame. Figure 4 shows these correlation functions for $L=0$ and $L=2.4$, with the top row corresponding to LB results for $C_{T}$ and the bottom row to kinetic theory results from Eq. (16) for $C_{U}$. First, we notice that the two sets of curves are very similar, implying that the stationary tracer approximation $C_{T}\approx C_{U}$ is indeed accurate. Secondly, we notice that the decay of the correlation function is significantly faster for swimmers than for shakers, again illustrating that self-propulsion acts an efficient decorrelation mechanism for ${\bm{U}}$. The effect of finite $D_{0}$ for shakers (left column) is simply to decrease the relaxation time of the exponential decay, in accordance with the $L=0$ limit of Eq. (16). For swimmers, the situation is more complex: For short times, the flow field decays faster with decreasing Pe, while the long-time tail of $C_{T}$ and $C_{U}$ instead becomes somewhat more extended with decreasing Pe. For fast enough swimmers, the latter effect leads to the local maximum at $\xi>1$ for intermediate Pe observed in Fig. 3b for the two highest $L$. Finally, we note that the equal-time fluid velocity variance $\langle U^{2}({\bm{r}}_{T})\rangle$, corresponding to the $t=0$ values of $C_{T}$ and $C_{U}$, is independent of Pe. This means that, regardless of the ratio between diffusive and active motion, the tracer particles sample the overall flow field homogeneously. This fact is non-trivial, since dry active particles that move autonomously on a solid substrate are known to preferentially sample regions where they move slowly (Stenhammar et al., 2016). Our results thus highlight that this generic mechanism is absent for inertialess point tracers advected by an incompressible fluid. However, for microswimmer systems where entrainment due to tracers being captured by the near-field flows of passing swimmers (Jeanneret et al., 2016), we would expect the fluid flow sampled by tracers to be significantly different from the average flow field in the system.

In Fig. 5, we study the dependence of the suppression of active diffusion on the microswimmer density $n$ in the shaker limit $L=0$. At first sight, this dependence might appear trivial, since $D_{A}$ is well-known to be linearly dependent on $n$ (Jepson et al., 2013; Miño et al., 2013) in the limit of non-interacting swimmers, a fact which is unaffected by Brownian motion as shown by Eq. (17). Since $\xi$ measures the ratio between the Brownian and non-Brownian values of $D_{A}$, one would naïvely expect $\xi$ to be independent of $n$. However, since the Péclet number itself, as defined in Eq. (6), increases with $n$ for constant $D_{0}$, the relative effect of Brownian motion on $D_{A}$ is in fact a complex function of $n$ even for noninteracting swimmers, as shown in Fig. 5. More specifically, according to Fig. 5b the suppression of active diffusion becomes more significant with increasing microswimmer density. In physical units, the highest concentration considered ($n\varepsilon^{3}\approx 0.08$) approximately corresponds to a bacterial concentration of $10^{9}$ mL*-1*, which is somewhat higher than the highest E. coli concentration considered by Jepson et al. (2013) but still within the range of concentrations where swimmer-swimmer correlations are reasonably small (Stenhammar et al., 2017). In summary, our results thus show that, for Brownian diffusion to have any measurable influence on the hydrodynamic diffusion, it is necessary to create a system with a relatively high density of very slow microswimmers; as we discuss in Section 5, this set of parameters is likely not achievable for suspensions of biological microswimmers.

5 Summary and conclusions

In this study we have demonstrated a number of theoretical and computational results regarding the effect of Brownian diffusion on the swimmer-induced hydrodynamic diffusion of tracer particles in a suspension of dipolar microswimmers. Our key finding is that the effect of Brownian diffusivity $D_{0}$ on the activity-induced, hydrodynamic diffusivity $D_{A}$ is only significant when the Péclet number as defined by Eq. (6) is below unity, meaning that Brownian diffusivity needs to dominate over the hydrodynamic one. However, the necessary requirement $\mathrm{Pe}<1$ is only sufficient in the shaker limit $v_{s}\rightarrow 0$: for swimmers with persistence lenghts larger than the organism size, significantly lower values of Pe are required to perturb the tracer trajectories sufficiently to affect $D_{A}$. This conclusion is analogous to the independence of $D_{A}$ on the tumbling rate $\lambda$ for large $v_{s}$ illustrated in Eq. (1): Whenever $v_{s}$ is large, the decorrelation of ${\bm{U}}$ by swimming will dominate over the decorrelation due to tumbling and translational diffusion, and the dependence on $\lambda$ and $D_{0}$ will thus vanish in the limit $v_{s}\rightarrow\infty$. While this effect is expected, what is perhaps surprising is the rather moderate values of $L=v_{s}/(\varepsilon\lambda)$ necessary to render the coupling between $D_{0}$ and $D_{A}$ negligible, as illustrated in Fig. 3b. To put these values into perspective, we use the conservative estimate $L=5$ for E. coli. By virtue of Fig. 3b, we require that $\mathrm{Pe}\leq 10^{-3}$ for Brownian motion to have an effect of $\sim 5$ percent on $D_{A}$. Using as an example the minimum value $D_{A}\approx${10}^{-2}\text{,}{\mathrm{\SIUnitSymbolMicro m}}^{2}\text{,}{\mathrm{s}}^{-1} measured by Jepson *et al.* ([2013](#bib.bib9)) in a 3-dimensional *E. coli* suspension, this requirement thus implies that $D_{0}\geq$10\text{\,}{\mathrm{\SIUnitSymbolMicro m}}^{2}\text{\,}{\mathrm{s}}^{-1}, which by virtue of the Stokes-Einstein relation corresponds to a tracer radius of $R_{0}\sim 20$ nm. While this is significantly smaller than used in typical measurements on colloidal tracers (Miño et al., 2013; Patteson et al., 2016; Leptos et al., 2009), this diffusion coefficient is close to the value of $D_{0}$ measured for dextran in E. coli suspensions by Kim & Breuer (2004). It is also fully feasible to realise such low Péclet numbers for micron-sized spheres by instead decreasing the bacterial density to very low values; however, measuring the correction to $D_{A}$ for $\mathrm{Pe}\sim 10^{-3}$ represents a major difficulty, since it amounts to measuring a $\sim$ 5 percent deviation of an effective diffusivity which is itself a thousand times smaller than the Brownian diffusion. It would thus require an extremely accurate determination of $D_{0}$, which then needs to be subtracted from the total measured diffusion constant to determine $D_{A}$. Obtaining this accuracy in a colloidal suspension would be very challenging due to particle polydispersity, interactions with boundaries, temperature gradients, and other system-specific complications. We thus conclude that, for typical $L$ values relevant for biological microswimmers, our results imply that the effect of Brownian motion on $D_{A}$ is likely negligible for all practical purposes. To experimentally observe the reduction in $D_{A}$, one would instead need to study a system of dipolar shakers, which stir up the surrounding fluid without self-propelling. While this is a somewhat exotic type of system, it could potentially be realised by anchoring molecular motors or biological microswimmers to a surface. In a biological setting, the shaker limit furthermore resembles previously developed models of enzymes anchored to lipid bilayers that induce dipolar flows through cyclical conformation changes (Hosaka et al., 2020).

Our results are qualitatively consistent with the previous theoretical results obtained by Kasyap et al. (2014), including the non-monotonic behaviour of $D_{A}$ with respect to $D_{0}$ at high swimming speed. Their results however differ in three important ways. First of all, their results consider a “slender swimmer” model, where the fluid is forced via a stress applied along a line representing the swimmer body, rather than by two point forces. While this model should lead to a dipolar flow in the far field, the near-field differences are significant, as illustrated in Fig. 8 of Kasyap et al. (2014). Furthermore, they consider only the fast-swimming limit where, according to Eq. (1), $D_{A}$ for dipolar microswimmers is independent of $\lambda$. Finally, they parametrise their model in a qualitatively different way than us: as discussed above, they adopt a definition of $\mathrm{Pe}$ based on $v_{s}$ rather than on the tracer diffusivity as in Eq. (6). In this description, $v_{s}$ is furthermore directly coupled to $\kappa$, so that changing Pe simultaneously changes the activity of the bath (via $\kappa$) and the swimming speed $v_{s}$. These differences makes it difficult to compare directly with our results, as we consider the effects of fluid advection and self-propulsion separately via Pe and $L$.

Thus, even though Brownian motion is unlikely to provide a significant dependence of $D_{A}$ on $R_{0}$ for tracers immersed in suspensions of biological microswimmers, there are several other mechanisms that need to be studied to explain the non-monotonic dependence observed experimentally (Patteson et al., 2016) and computationally (Dyer & Ball, 2021). First of all, the effect of tracer entrainment by the near field flows of the swimmer is strongly dependent on the size ratio between the swimmer and the tracer (Jeanneret et al., 2016), although we expect this term to be small for micron-sized tracer particles in E. coli suspensions. Secondly, the finite size of the tracer will change the equation of motion (4) into the Faxén equation that takes into account the non-linearity of the flow field (Kim & Karrila, 1991), an effect which was implicitly included in the wavelet Monte Carlo simulations by Dyer & Ball (2021) and, together with tracer entrainment, is a significant explanation of their observed $R_{0}$ dependence of $D_{A}$. Finally, non-hydrodynamic interactions such as direct collisions, electrostatic interactions, and artifacts due to container walls are likely to depend in a non-trivial manner on the tracer size for each system in question. We thus conclude by noting that significant further experimental and theoretical work is necessary to disentangle the system-specific properties from the generic properties of tracer dynamics in microswimmer suspensions.

Acknowledgement

Discussions with Cesare Nardini are kindly acknowledged. The computations were enabled by resources provided by LUNARC. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

Funding

This work was financed through the Knut and Alice Wallenberg Foundation (project grant KAW 2014.0052). JS acknowledges financial support from the Swedish Research Council (Project No. 2019-03718).

Competing Interests

The authors report no competing interests.

Appendix A Kinetic theory for $D_{A}(D_{0})$.

In this Appendix, we demonstrate the main steps in the derivation of Eq. (16). The analysis follows closely a similar derivation presented in Škultéty et al. (2020), albeit with two major differences. First, in addition to that work, we include the effect of microswimmer Brownian diffusivity, as discussed in Section 3. Second, we consider the case of non-interacting microswimmers, which significantly simplifies the analysis. Due to the similarity with the derivation in Škultéty et al. (2020), we here present the key steps of the derivation and refer the interested reader to that paper for technical details.

The quantity of interest is the fluid velocity autocorrelation function $C_{U}(t)$, formally defined as

[TABLE]

where $V$ is the volume of the system, and the overbar denotes an average over the stochastic history, i.e., the history of tumble events, and the long-time limit ensures independence of the initial state of the system. Here and in the following, the superscript indices denote Cartesian components of vectors. For a given state of the system, the instantaneous fluid velocity $U^{\alpha}$ at a position ${\bm{r}}$ and time $t$ is readily obtained as a superposition of individual velocity fields generated by the microswimmers

[TABLE]

Here, ${\bm{r}}_{i}(t)$ gives the instantaneous position of particle $i$, while the unit vector ${\bm{p}}_{i}(t)$ gives its instantaneous orientation; the index $i=1\dots N$ enumerates the particles, where $N$ is the total number of microswimmers. In the following, we assume ${\bm{u}}_{s}$ to be given by the regularised hydrodynamic dipole (9).

Time evolution of the suspension comprises spatial motion of the microswimmers according to Eq. (15) and their random re-orientation with rate $\lambda$, as discussed in Section 1. These dynamics can be alternatively described by the master equation

[TABLE]

where $F_{N}=F_{N}({{\bm{r}}}_{1},\dots,{{\bm{r}}}_{N},{{\bm{p}}}_{1},\dots,{{\bm{p}}}_{N},t)$ is the $N$-particle probability distribution function. Here, ${\bm{\nabla}}_{i}$ denotes spatial derivatives with respect to the coordinates of particle $i$. As shown in Škultéty et al. (2020), the same dynamics can be conveniently encoded by the following equation:

[TABLE]

where $h=h({\bm{r}},{\bm{p}},t)$ is an auxiliary field related to the Klimontovich correlation function (Klimontovich & Silin, 1962), and $\chi$ is a noise term with the following properties:

[TABLE]

The advantage of this representation lies in its direct relation to the phase-space density correlation function (Klimontovich & Silin, 1962), which allows us to express the fluid velocity autocorrelation function as

[TABLE]

where the angular brackets denote the average of the (fictitious) noise $\chi$ (see Škultéty et al. (2020) for details).   The linear equation (21) is readily solved by introducing the Fourier–Laplace transform $\hat{h}$ of the auxiliary field $h$:

[TABLE]

which yields

[TABLE]

Here, $\sigma({\bm{k}},{\bm{p}},s)=s+\lambda+D_{0}k^{2}+iv_{s}({\bm{k}}\cdot{\bm{p}})$, $\hat{\chi}({\bm{k}},{\bm{p}},s)$ is the Fourier–Laplace transform of the noise, and we have dropped the initial condition $\hat{h}({\bm{k}},{\bm{p}},t=0)$ which does not contribute in the large-$t$ limit.

As demonstrated by Škultéty et al. (2020), only the first term in Eq. (26) contributes to $C_{U}(t)$. Performing the Fourier–Laplace transform in Eq. (24) and combining it with Eq. (26) yields

[TABLE]

Performing the inverse Laplace transforms, denoted symbolically by ${\mathcal{L}}^{-1}$, and integrating over $\bm{p}$ and the orientation of $\bm{k}$, we finally arrive at Eq. (16) in the main text.