The role of correlations in the collective behaviour of microswimmer suspensions

TL;DR

This paper investigates the collective behavior of microswimmers in fluid, revealing correlations below turbulence transition, developing a new kinetic theory, and validating it with large-scale simulations and tracer diffusivity measurements.

Contribution

It introduces a novel kinetic theory capturing correlations in microswimmer suspensions, validated by large-scale simulations and experimental measures.

Findings

Swimmers exhibit correlated motion below turbulence transition.

The kinetic theory accurately predicts the asymmetry between pusher and puller swimmers.

Tracer diffusivity shows non-trivial density dependence.

Abstract

In this Letter, we study the collective behaviour of a large number of self-propelled microswimmers immersed in a fluid. Using unprecedently large-scale lattice Boltzmann simulations, we reproduce the transition to bacterial turbulence. We show that, even well below the transition, swimmers move in a correlated fashion that cannot be described by a mean-field approach. We develop a novel kinetic theory that captures these correlations and is non-perturbative in the swimmer density. To provide an experimentally accessible measure of correlations, we calculate the diffusivity of passive tracers and reveal its non-trivial density dependence. The theory is in quantitative agreement with the lattice Boltzmann simulations and captures the asymmetry between pusher and puller swimmers below the transition to turbulence.

Equations

$${\dot{r}}_i^{\alpha }=v_s{p}_i^{\alpha }+U^{\alpha }({\mathbf{r}}_i),{\dot{p}}_i^{\alpha }={\mathbb{P}}_i^{\alpha \beta }{\nabla }_i^{\gamma }{U}^{\beta }({\mathbf{r}}_i)p_i^{\gamma },$$

$$k^2{u}_{\mathbf{k}}^{\alpha }(\mathbf{p})=-i\kappa A(kϵ)(\mathbf{k}\cdot \mathbf{p})[p_{\alpha }-k^{-2}{k}_{\alpha }(\mathbf{k}\cdot \mathbf{p})].$$

$${\partial }_{t}g+L_f^{(1)}[g]+L_f^{(2)}[g]=C,$$

$$C\equiv 3{\left(\frac{n}{4\pi }\right)}^2[p_1^{\alpha }{p}_1^{\beta }{\nabla }^{\alpha }{u}^{\beta }({\mathbf{r}}_1-{\mathbf{r}}_2,{\mathbf{p}}_2)$$

$$+p_2^{\alpha }{p}_2^{\beta }{\nabla }^{\alpha }{u}^{\beta }({\mathbf{r}}_2-{\mathbf{r}}_1,{\mathbf{p}}_1)],$$

$$L_f[h]$$

$${\partial }_t\delta f+L_f[\delta f]=\eta .$$

$${\widetilde{\delta U}}_{\mathbf{k}}^{\alpha }(\omega )=\frac{1}{C_0(k,\omega )}\int d^2\mathbf{p}\frac{{\widetilde{\delta f}}_{\mathbf{k}}(\mathbf{p},t=0)u_{\mathbf{k}}^{\alpha }(\mathbf{p})}{-i\omega +\lambda }$$

$$\frac{⟨U^2⟩}{v_0^2}\approx \frac{21\rho {\kappa }_n^2}{2048{ϵ}_n}\left[1\pm \frac{\rho (2{\rho }_c\mp \rho )}{2{\rho }_c({\rho }_c\mp \rho )}\right],$$

$$D_h=\frac{1}{3}{\int }_0^{\infty }ds\mathcal{C}({\mathbf{x}}_s,s,{\mathbf{x}}_0,0).$$

$$\frac{D_h}{D_h^{\mathrm{free}}}\approx \frac{1}{2}\left[1+\frac{{\rho }_c}{{\rho }_c\mp \rho }\pm \left(\frac{\rho }{{\rho }_c}+\frac{\rho {\rho }_c}{({\rho }_c\mp \rho {)}^2}\right)\right].$$

Full text

A suspension of particles that can extract energy from their surroundings and transform it into self-propulsion is an archetypal example of active matter Ramaswamy (2010); Marchetti et al. (2013). Such systems do not obey the principle of detailed balance Cates (2012) at a single-particle level and their behaviour often differs significantly from that of passive suspensions at the same conditions Poo (2013). Experiments on self-propelled particles like bacteria Dombrowski et al. (2004); Sokolov et al. (2007); Gachelin et al. (2014), sperm cells Creppy et al. (2015), mixtures of microtubules and molecular motors Sanchez et al. (2012), vibrated granular rods Blair et al. (2003), “Quincke rollers” Bricard et al. (2013), and “colloidal surfers” Palacci et al. (2013) reveal the existence of non-equilibrium steady states with non-zero macroscopic fluxes in these systems. One of the most striking examples is the phenomenon of “bacterial turbulence” Sokolov and Aranson (2012); Dunkel et al. (2013); Gachelin et al. (2014); Wensink et al. (2012), whereby a suspension of swimming bacteria at sufficient density forms a state with large-scale, coherent fluid motion.

Previous analytical Saintillan and Shelley (2008); Subramanian and Koch (2009); Hohenegger and Shelley (2010); Saintillan and Shelley (2013) and numerical Saintillan and Shelley (2007); Wolgemuth (2008); Saintillan and Shelley (2011); Lushi and Peskin (2013); Krishnamurthy and Subramanian (2015) studies recognise long-range hydrodynamic interactions between swimmers as a key ingredient of their collective motion; in the absence of external forces and torques, these interactions can be described by a dipolar field Lauga and Powers (2009). The main observation of the previous studies is that bacterial turbulence only emerges in suspensions of pushers (dipolar swimmers that expel fluid along their long direction), and is absent for pullers (dipolar swimmers that do the opposite). This conclusion has been corroborated by mean-field kinetic theories that consider the dynamics of a single swimmer in an average hydrodynamic field produced by other particles Saintillan and Shelley (2008); Subramanian and Koch (2009); Hohenegger and Shelley (2010); Saintillan and Shelley (2013); Krishnamurthy and Subramanian (2015).

Below the transition to collective motion, previous theoretical studies view motile suspensions as random and featureless, only acquiring non-trivial properties above the transition. In this Letter we demonstrate that due to the long-range nature of (unscreened) hydrodynamic interactions between swimmers, pre-transitional suspensions develop very strong correlations that dominate their dynamical properties (see also Qian et al. (2017)). Strong spatial, temporal and orientational correlations between swimmers act as precursors to bacterial turbulence and are essential to understanding the transition. We therefore develop a novel kinetic theory that goes beyond mean-field and compare its predictions to particle-resolved lattice Boltzmann (LB) simulations of up to 4 million hydrodynamically interacting microswimmers.

To illustrate the significance of swimmer-swimmer correlations, we consider the advection of passive tracer particles immersed in a microswimmer suspension. The effective tracer diffusion constant is an observable sensitive to the dynamical state of the system, and it has been extensively studied experimentally in suspensions of bacteria Wu and Libchaber (2000); Miño et al. (2011); Jepson et al. (2013) and algae Leptos et al. (2009); Jeanneret et al. (2016); Kurtuldu et al. (2011). For low densities of swimmers, it has been predicted to scale linearly with the swimmer density and to be identical for pushers and pullers of equal dipolar strengths Lin et al. (2011); Miño et al. (2011); Jepson et al. (2013); Pushkin and Yeomans (2013); Kasyap et al. (2014); Morozov and Marenduzzo (2014). Here, we use LB simulations and kinetic theory to show how correlations break the pusher-puller symmetry and result in non-linear scaling of the enhanced diffusivity with swimmer density. Our analysis suggests that these correlations become significant even for densities as low as $10\%$ of the critical density, with the latter being estimated at a volume fraction of $\sim 2\%$ for E. coli-like parameters (see below).

Model description. We consider a 3-dimensional suspension of $N$ microswimmers immersed in a fluid of volume $V$ at number density $n=N/V$. Each swimmer exerts a force on the fluid, $-F\mathbf{p}$, representing the flagellum, and an equal and opposite force, $F\mathbf{p}$, applied a distance $l$ from the propulsive force, representing the cell body; here $\mathbf{p}$ is the swimmer orientation (see Fig. 1a). The body is modelled as a sphere with hydrodynamic radius $a$, while the velocity scale of this model is defined as $v_{0}\equiv F/(\mu l)$, where $\mu$ is the viscosity of the fluid. The dipolar strength of each swimmer is given as $\kappa=Fl/\mu$, with $\kappa>0$ representing pushers and $\kappa<0$ pullers, and the non-dimensional swimmer density is defined as $\rho=nl^{3}$. To facilitate comparison with experimental results, we also define the reduced density $\rho_{\mathrm{body}}=(4\pi/3)a^{3}n$, where $a$ is estimated as described in SI , which gives an estimate of the volume fraction based on the bacterial body volume. Each swimmer $i$ moves according to the following equations of motion:

[TABLE]

where $\mathbb{P}_{i}^{\alpha\beta}=\delta_{\alpha\beta}-p_{i}^{\alpha}p_{i}^{\beta}$, $\mathbf{U}(\mathbf{r}_{i})$ is the fluid velocity at the position of swimmer $i$, and $v_{s}$ is the swimming speed (see SI ); Greek indices denote Cartesian components. In addition to being rotated by the fluid, the orientation of each swimmer is randomized with average tumbling frequency $\lambda$ 111The rotational dynamics of Eq. (1) corresponds to flow-aligning swimmers, and is valid in the limit $l\gg a$; the generalisation to finite aspect ratios is straightforward and does not change the main conclusions drawn from the model..

Lattice Boltzmann simulations. Large-scale numerical simulations of up to $\mathcal{O}(10^{6})$ hydrodynamically interacting microswimmers in cubic boxes with periodic boundary conditions were performed using a D3Q15 lattice Boltzmann (LB) algorithm. The swimmers are described using the point-force implementation developed by Nash et al. Nash et al. (2008, 2010), which accurately captures the full far-field interactions between the particles, while neglecting short-range hydrodynamics, lubrication effects and non-hydrodynamic interactions; see SI for further details. As shown in Fig. 1b, we qualitatively capture the transition from seemingly random motion at low density of swimmers, to bacterial turbulence at higher densities of pushers, visible as large-scale fluid vortices and jets. Quantitatively, this is characterised in Fig. 1c as a rapid deviation of the root-mean-square fluid velocity $\langle U^{2}\rangle^{1/2}$ for pushers at $\rho_{\mathrm{body}}\approx 0.02$ (i.e., close to the experimentally observed transition densities Sokolov and Aranson (2012); Gachelin et al. (2014)), from the expected $\rho^{1/2}$ behaviour at low density to a state with much larger velocity fluctuations; for pullers, correlations lead to $\langle U^{2}\rangle^{1/2}$ increasing slower than $\rho^{1/2}$. Figure 1c also highlights the strong system-size dependence of these results: velocity fluctuations change appreciably when going from a box length of $L=25$ ($N\simeq 10^{3}-10^{4}$, comparable to previous particle-resolved studies Saintillan and Shelley (2011); Lushi and Peskin (2013); Krishnamurthy and Subramanian (2015)) to $L=200$ ($N\simeq 10^{6}-10^{7}$), indicating the need for very large-scale simulations (at least $L=100$ and $N\simeq 10^{5}$) to describe the collective motion in these systems. We postpone further characterisation of the turbulent state to future studies, and in the following we focus on velocity fluctuations in the pre-transitional region ($\rho_{\mathrm{body}}\leq 0.02$) and the buildup of correlations that lead to bacterial turbulence.

In Fig. 2, we make the observation that, in this density regime, interactions between swimmers are dominated by their mutual rotation, while the effect of advection is secondary. There we also show results obtained for so-called shakers – particles that apply forces to the fluid but do not swim. These behave very similarly to swimmers, not only qualitatively Hatwalne et al. (2004) but also quantitatively, indicating that the effect of swimming is subdominant. We now proceed to develop the kinetic theory, where the latter observation will prove to be mathematically convenient.

Kinetic theory. The central theoretical tool developed in this Letter is a kinetic theory that describes the suspension at a coarse-grained level. We assume that the swimmer density $\rho$ is sufficiently low that the single-swimmer Fourier-space velocity field $\mathbf{u}_{\mathbf{k}}(\mathbf{p})$ can be described by that of a regularised point dipole Cortez et al. (2005):

[TABLE]

Here, $k=|\mathbf{k}|$, $A(x)=x^{2}K_{2}(x)/2$ with $K_{2}$ the modified Bessel function of the second kind, and $\epsilon$ is a regularisation parameter. The starting point of any kinetic theory is standard Balescu (1997) and is briefly summarised in the following. The dynamics of the system is described by the master equation Gardiner (2004) for the $N$-body probability density function (PDF) $f_{N}(\mathbf{z}_{1},...,\mathbf{z}_{N},t)$, with $\mathbf{z}_{i}=(\mathbf{r}_{i},\mathbf{p}_{i})$. We then introduce the reduced PDFs $f_{s}=N!/(N-s)!\int d\mathbf{z}_{s+1}...d\mathbf{z}_{N}f_{N}$ and derive a BBGKY hierarchy for the $f_{s}$ SI . This hierarchy is conveniently written in terms of connected correlations $g_{s}$ defined by $f_{2}(\mathbf{z}_{1},\mathbf{z}_{2},t)=f_{1}(\mathbf{z}_{1},t)f_{1}(\mathbf{z}_{2},t)+g_{2}(\mathbf{z}_{1},\mathbf{z}_{2},t)$ and similarly for higher orders. To close the exact BBGKY hierarchy in a controlled way, we follow the approximation scheme used in equilibrium Balescu (1997); Nicholson (1983); Heyvaerts (2010); Campa et al. (2009); Bouchet et al. (2010) and non-equilibrium Nardini et al. (2012a, b) systems with long-range interactions, such as plasmas and self-gravitating systems. We consider the limit of a large number of swimmers, $N\gg 1$, at a fixed density $\rho$; our approach is perturbative in the small parameter $1/N$, not in $\rho$, and, for this reason, we get accurate predictions even close to the onset of bacterial turbulence. In this limit, and in the absence of hidden divergences 222This hypothesis is checked a posteriori showing that the final results keep the order of magnitude assumed. , we can demonstrate that $f_{1}\sim\mathcal{O}(1)$ and $g_{s}\sim\mathcal{O}(1/N^{s-1})$. At leading order, $g_{s}=0$ for $s\geq 2$, we obtain the mean-field approximation already analysed in the literature for swimmers Saintillan and Shelley (2008); Subramanian and Koch (2009); Hohenegger and Shelley (2010); Saintillan and Shelley (2013); Krishnamurthy and Subramanian (2015) and suspensions of passive rods Doi and Edwards (1988). In this Letter we go beyond the mean-field approximation, retaining $g_{2}$ and discarding $g_{s}$ for $s\geq 3$ SI ; henceforth, we use $f\equiv f_{1}$ and $g\equiv g_{2}$.

The perturbative analysis summarized above shows that $g$ solves

[TABLE]

where 333We assume here $f=n/4\pi$, which simplifies significantly both the expressions for $C$ and for $L_{f}$ and corresponds to a suspension which is isotropic and homogeneous on average. Their expressions for general $f$ is given in SI .

[TABLE]

and

[TABLE]

where we used the notation $\partial^{\alpha}\equiv\mathbb{P}^{\alpha\beta}\frac{\partial}{\partial p^{\beta}}$. Here, $U_{\mathrm{mf}}^{\alpha}[h](\mathbf{r})=\int d\mathbf{r}_{1}d\mathbf{p}_{1}\,\,u^{\alpha}(\mathbf{r}-\mathbf{r}_{1},\mathbf{p}_{1})h(\mathbf{r}_{1},\mathbf{p}_{1},t)$, and $L_{f}[h]$ is the mean-field operator $V$ linearised close to $f$ and acting on the function $h$, where $V=\nabla^{\alpha}(\dot{r}_{\mathrm{mf}}^{\alpha}f)+\partial^{\alpha}(\dot{p}_{\mathrm{mf}}^{\alpha}f)+\lambda f-(\lambda/4\pi)\int d\mathbf{p}f$ and $(\dot{\mathbf{r}}_{\mathrm{mf}},\dot{\mathbf{p}}_{\mathrm{mf}})$ are given by Eqs. (1) with $U^{\alpha}(\mathbf{r}_{i})$ replaced by $U_{\mathrm{mf}}^{\alpha}[f](\mathbf{r})$; $L_{f}^{(i)}$ acts on $\mathbf{z}_{i}$.

Motivated by our numerical observation that the nature of the transition and the properties of the suspension can be understood in the absence of self-propulsion, we consider only the case of $v_{s}=0$. We stress that, while being significantly more complex, all the results presented below can also be obtained for $v_{s}>0$, as will be shown in a forthcoming publication Nardini et al. .

We now introduce an Itō white noise $\eta$ with covariance $\mathbb{E}[\eta\,\eta]=C\,\delta(t-t^{\prime})$, where $\mathbb{E}$ denotes the average over $\eta$, and $C$ is given by Eq. (The role of correlations in the collective behaviour of microswimmer suspensions). Formally, the connected correlator $g$ can be written as $g(\mathbf{z}_{1},\mathbf{z}_{2},t)=\mathbb{E}[\delta f(\mathbf{z}_{1},t)\,\delta f(\mathbf{z}_{2},t)]$, where $\delta f$ solves

[TABLE]

Because $g$ is the covariance of density fluctuations close to $f$, these fluctuations are described by the random field $\delta f$. Being small, they are given by a linear stochastic process, Eq. (6), although the variance of the noise $C$ is non-trivial and could not have been guessed a priori. It is also remarkable that, as a result of the coarse graining procedure, the density fluctuations are described by a stochastic process even when the underlying microscopic dynamics are deterministic ($\lambda=0$). Equations (3) and (6) can be solved exactly to yield $\delta f=e^{-tL_{f}}\delta f(t=0)+\int_{0}^{t}e^{-(t-s)L_{f}}\eta(s)ds$ and $g=\int_{0}^{t}e^{-sL_{f}^{(1)}}e^{-sL_{f}^{(2)}}C\,ds$ 444We assume here $g(t=0)=0$. A different initial condition for $g$ only amounts to a transient, exponentially decaying contribution which does not affect the stationary properties.. As $\widetilde{\delta f}=e^{-tL_{f}}\widetilde{\delta f}(t=0)$ solves $\partial_{t}\widetilde{\delta f}+L_{f}[\widetilde{\delta f}]=0$, to compute $g$ we only need the solution $\widetilde{\delta f}$ of the above deterministic dynamics, with appropriate initial conditions set by $C$. It turns out that $\widetilde{\delta f}$ can be found exactly for a generic initial condition SI . For the fluctuations of the fluid velocity in Fourier-Laplace space $\widetilde{\delta U}^{\alpha}_{\mathbf{k}}(\omega)=\int d^{2}\mathbf{p}^{\prime}\,u^{\alpha}_{\mathbf{k}}(\mathbf{p}^{\prime})\,\widetilde{\delta f}_{\mathbf{k}}(\mathbf{p}^{\prime},\omega)$, we then obtain the following closed expression

[TABLE]

where, for the regularised dipolar field in Eq. (2), we have $C_{0}(k,\omega)=1-(\kappa n/5)A(k\epsilon)/(-i\omega+\lambda)$. Equation (7) is valid under the assumption that the dynamical state described by $f$ is linearly stable, which corresponds to the zeros $\omega_{*}=\omega_{R}+i\omega_{I}$ of $C_{0}(k,\omega)$ having negative imaginary parts in the Laplace domain, $\omega_{I}<0$. These zeros are given by $\omega_{R}=0$ and $\omega_{I}(k\epsilon)=-\lambda+\kappa nA(k\epsilon)/5$, which implies that a suspension of pullers ($\kappa<0$) is always stable, while pusher suspensions ($\kappa>0$) are stable only if $n<5\lambda/\kappa$, in agreement with earlier results Saintillan and Shelley (2008); Subramanian and Koch (2009); Hohenegger and Shelley (2010); Saintillan and Shelley (2013).

We note the emergence of the characteristic time-scale $\omega_{I}^{-1}$, which describes the typical time for a small fluctuation of the fluid velocity to relax. At $\rho=0$ it reduces to $\omega_{I}^{-1}=\lambda^{-1}$, while upon increasing the density, it decreases for pullers and increases for pushers. This suggests qualitative differences between the statistical properties of suspensions of pushers and pullers even below the onset of bacterial turbulence.

Fluid velocity variance. We now compare the kinetic theory developed above with the results of LB simulations by computing the variance $\langle U^{2}\rangle$ of the fluid velocity: $\langle U^{2}\rangle=\langle U^{2}\rangle_{\mathrm{mf}}+\langle U^{2}\rangle_{\mathrm{corr}}$, where $\langle U^{2}\rangle_{\mathrm{mf}}=n/(4\pi(2\pi)^{3})\int d\mathbf{k}d\mathbf{p}\,u_{\mathbf{k}}^{\alpha}(\mathbf{p})u_{-\mathbf{k}}^{\alpha}(\mathbf{p})$ is the mean-field contribution and $\langle U^{2}\rangle_{\mathrm{corr}}$ contains corrections induced by swimmer-swimmer correlations. Using the formal solution for $g$ and Eqs. (The role of correlations in the collective behaviour of microswimmer suspensions) and (7) we obtain $\langle U^{2}\rangle_{\mathrm{mf}}=v_{0}^{2}\frac{\rho\kappa_{n}^{2}}{15\pi^{2}}\int_{0}^{\infty}dkA^{2}(\epsilon_{n}k)$ and $\langle U^{2}\rangle_{\mathrm{corr}}=-v_{0}^{2}\frac{\rho^{2}\kappa_{n}^{3}}{75\pi^{2}}\int_{0}^{\infty}dk\frac{A^{2}(\epsilon_{n}k)}{\omega_{n}(\epsilon_{n}k)}$, which is well approximated by SI

[TABLE]

where the upper sign corresponds to pushers and the bottom one to pullers. Here, we have introduced the non-dimensional units $\kappa_{n}=\kappa/(l^{2}v_{0})$, $\epsilon_{n}=\epsilon/l$, $\omega_{n}=\omega_{I}l/v_{0}$, and $\lambda_{n}=\lambda l/v_{0}$, using $l$ and $l/v_{0}$ as the respective length- and time-scales. The density $\rho_{c}\equiv 5\lambda_{n}/|\kappa_{n}|$ corresponds to the onset of collective motion for pushers, and acts as a characteristic density scale for pullers. Fig. 2 shows excellent agreement between this prediction and the LB data, even close to the onset of turbulence, emphasising the non-perturbative nature of our approach. Equation (8) and the data in Fig. 2 can also be used to assess the relative importance of correlations: for $\rho\lesssim 0.1\rho_{c}$, pushers and pullers behave equivalently and follow the mean-field prediction, while above it, correlations have to be taken into account to obtain correct quantitative predictions. We furthermore note that, using the body volume $(4\pi/3)a^{3}$ with $a\approx 0.3$ SI , we obtain $\rho_{c}^{\mathrm{body}}\approx 0.023$, in good agreement with the experimentally observed transitional volume fractions of $2\%$ in E. coli Gachelin et al. (2014) and B. subtilis Sokolov and Aranson (2012).

Enhanced diffusivity. We now consider the advection of a passive tracer with the dynamics $\dot{\mathbf{x}}=\mathbf{U}(\mathbf{x},t)$ immersed in a suspension of shakers. Its long-time motion is diffusive Thiffeault (2015), and we use the kinetic theory developed here to calculate its effective diffusion constant $D_{h}$. Denoting by $\mathbf{x}_{t}$ the tracer position at time $t$, we have $\langle|\mathbf{x}_{t}-\mathbf{x}_{0}|^{2}\rangle=6D_{h}t$, with $D_{h}$ related to the fluid velocity autocorrelation by Kubo (1966)

[TABLE]

Here, $\mathcal{C}(\mathbf{x}_{s},s,\mathbf{x}_{0},0)=\mathbb{E}\left[\delta U^{\alpha}(\mathbf{x}_{s},s)\,\delta U^{\alpha}(\mathbf{x}_{0},0)\right]$ and $\delta\mathbf{U}$ is the fluctuating fluid velocity obtained from the solution of Eq. (6). The computation of $D_{h}$ is greatly simplified by iteratively inserting the solution of the tracer dynamics into Eq. (9), recalling that $\delta\mathbf{U}$ is small, and thus Taylor expanding around $\mathbf{x}_{0}$ (see Bouchet and Dauxois (2005a, b)). At leading order, we obtain $D_{h}\simeq\frac{1}{3}\int_{0}^{\infty}ds\,\mathcal{C}(\mathbf{x}_{0},s,\mathbf{x}_{0},0)$, where corrections due to higher-order correlations have been discarded in agreement with the approximations made in the kinetic approach. Given the solution of the linear dynamics, the diffusivity reads $D_{h}/(lv_{0})=-\frac{\rho\kappa_{n}^{2}}{45\pi^{2}}\int_{0}^{\infty}dkA^{2}(\epsilon_{n}k)/\omega_{n}(k\epsilon_{n})+\frac{\rho^{2}\kappa_{n}^{3}}{225\pi^{2}}\int_{0}^{\infty}dkA^{2}(\epsilon_{n}k)/\omega_{n}^{2}(k\epsilon_{n})$, and is approximated by SI

[TABLE]

As in Eq. (8), the top (bottom) sign corresponds to pushers (pullers). In the low density limit, $D_{h}\to D_{h}^{\mathrm{free}}+\mathcal{O}(\rho^{2})$, where $D_{h}^{\mathrm{free}}/(lv_{0})=\langle U^{2}\rangle_{\mathrm{mf}}/3\lambda_{n}v_{0}^{2}=7\rho\kappa_{n}^{2}/(2048\epsilon_{n}\lambda_{n})$ is the diffusivity obtained from discarding interactions between shakers. An interesting observation is that $D_{h}^{\mathrm{free}}$ is, contrary to the variance of the fluid velocity, sensitive to the presence of self-propulsion: for shakers, $D_{h}^{\mathrm{free}}$ diverges as $\epsilon_{n}\to 0$, while an analogous computation for swimmers gives $D_{h}^{\mathrm{free}}/lv_{0}=\rho\kappa_{n}^{2}/48\pi$, in agreement with earlier theoretical predictions Lin et al. (2011); Miño et al. (2011); Jepson et al. (2013); Pushkin and Yeomans (2013); Kasyap et al. (2014); Morozov and Marenduzzo (2014). We have confirmed this difference of $D_{h}^{\mathrm{free}}$ between shakers and swimmers through LB simulations at low density.

In Fig. 3, we show that the enhanced diffusivity of tracers measured in LB simulations is perfectly described by the kinetic theory even close to the onset of turbulence. We further observe that $D_{h}$ deviates from $D_{h}^{\mathrm{free}}$ even for small densities of shakers, again highlighting the importance of correlations. Moreover, Eq. (10) correctly predicts how $D_{h}$ depends differently on $\rho$ for pushers and pullers – an effect that has only briefly been discussed in the literature Underhill et al. (2008); Ishikawa et al. (2010); Kurtuldu et al. (2011); Krishnamurthy and Subramanian (2015); Saintillan and Shelley (2011).

An approach previously used by several authors Lin et al. (2011); Miño et al. (2011); Jepson et al. (2013); Pushkin and Yeomans (2013); Kasyap et al. (2014); Morozov and Marenduzzo (2014) to predict $D_{h}$ considers tracer displacements due to scattering from a single swimmer. This leads to $D_{h}\simeq D_{h}^{\mathrm{free}}\sim\rho\kappa_{n}^{2}$, and reflects the pusher-puller symmetry upon time-reversal in the Stokes equation. The extension of this argument to scattering by any finite number of swimmers presents a conceptual problem: the tracer displacement due to a scattering event by a collection of pushers can also be obtained in a suspension of pullers whose initial positions are set equal to the final positions of the pushers, with their orientations reversed. This argument thus suggests that $D_{h}$ should be equal for pushers and pullers at all densities, at odds with the numerical data in Fig. 3. The caveat in this argument is that it assumes a uniform sampling of initial conditions for the swimmers, while correlations between them will make some configurations more probable. When taken into account properly, these correlations break the pusher-puller symmetry even at moderate densities, as in Eq. (10).

Conclusions. In this Letter, we have presented results from a novel kinetic theory and unprecedently large particle-resolved simulations of microswimmer suspensions, describing quantitatively the fluctuations and correlations that arise at intermediate swimmer density. We have numerically shown that the collective motion in swimmer suspensions is the result of their rotational dynamics in the flow created by other swimmers, while self-propulsion and spatial correlations play subdominant roles. We calculated the fluid velocity fluctuations and the enhanced diffusivity of tracer particles, and found significant deviations from the mean-field predictions even at moderate swimmer densities. We demonstrated that swimmer-swimmer correlations are responsible for these deviations and should thus be taken into account well below the onset of bacterial turbulence. Understanding such correlations is a prerequisite for a deeper understanding of the turbulent state itself, in particular with regards to the presence or absence of a finite, characteristic length-scale in collective motion of microswimmers Dombrowski et al. (2004); Sokolov et al. (2007); Gachelin et al. (2014); Creppy et al. (2015); Sokolov and Aranson (2012); Dunkel et al. (2013); Wensink et al. (2012).

Acknowledgements. Discussions with Mike Cates and Joost de Graaf are kindly acknowledged. JS is financed by a grant from the Swedish Research Council (2015-05449), CN by EPSRC grant EP/J007404, RWN by Intel through EPCC’s Intel Parallel Computing Centre, and AM by EPSRC through grant number EP/I004262/1. CN acknowledges the hospitality provided by DAMTP, University of Cambridge while most of this work was being done and the support of an Aide Investissements d’Avenir du LabEx PALM (ANR-10-LABX-0039-PALM). Research outputs generated through the EPSRC grant EP/I004262/1 can be found at http://dx.doi.org/10.7488/ds/1703.

Author contributions. JS and CN contributed equally to this work.