Enhanced velocity fluctuations in interacting swimmer suspensions

TL;DR

This paper investigates how interactions among micro-swimmers affect velocity fluctuations, revealing divergence in variance and the impact of orientation decorrelation, with implications for passive tracer dynamics.

Contribution

It introduces a theoretical framework for understanding velocity fluctuations in interacting swimmer suspensions, including run-and-tumble particles, and explains experimental observations.

Findings

Velocity variance diverges logarithmically with system size in interacting suspensions.

Inclusion of orientation decorrelation mechanisms arrests the divergence.

Passive tracer dynamics exhibit a crossover between ballistic and diffusive regimes.

Abstract

A dilute non-interacting suspension of micro-swimmers exhibits a finite velocity variance and short-ranged correlations that decay over a swimmer length. For a suspension of interacting straight swimmers, however, pair-interactions leads to a non-decaying velocity covariance, and a variance that diverges logarithmically with system size. The divergence is arrested on inclusion of orientation decorrelation mechanisms. Results for suspensions of run-and-tumble particles (RTPs) are presented, where the underlying straight-swimmer divergence leads to a broad cross-over between the ballistic and diffusive regimes, of immersed passive tracers, in the limit of long run lengths. Our analysis explains long-standing experimental observations of a volume-fraction dependent crossover time for passive tracer dynamics.

Equations

$${\mathbf{u}}_i(\mathbf{x},t)={\int }_{-L/2}^{L/2}f(s)\text{mathsfbi}G(\mathbf{x}-{\mathbf{x}}_i-s{\mathbf{p}}_i)\cdot {\mathbf{p}}_i\mathbf{d}s,$$

$$\mathbf{u}(\mathbf{x},t)\equiv \sum _{i=1}^N{\mathbf{u}}_i(\mathbf{x},t)=\sum _{i=1}^N{\int }_{-L/2}^{L/2}f(s)\text{mathsfbi}G(\mathbf{x}-{\mathbf{x}}_i-s{\mathbf{p}}_i)\cdot {\mathbf{p}}_i\mathbf{d}s,$$

$$⟨\mathbf{u}(\mathbf{x},t)\cdot \mathbf{u}({\mathbf{x}}^{\prime },t)⟩$$

$$⟨\mathbf{u}(\mathbf{x})\cdot \mathbf{u}({\mathbf{x}}^{\prime })⟩$$

$$⟨\hat{\mathbf{u}}(\mathbf{k})\cdot \hat{\mathbf{u}}({\mathbf{k}}^{\prime })⟩$$

$$⟨\hat{\mathbf{u}}(\mathbf{k})\cdot \hat{\mathbf{u}}({\mathbf{k}}^{\prime })⟩=\frac{n{L}^3}{16{\pi }^7(\ln \kappa {)}^2{k}^4}\mathbf{\delta }(\mathbf{k}+{\mathbf{k}}^{\prime })\int \mathrm{d}{\mathbf{p}}_1\left[1-\frac{{\left(\mathbf{k}\cdot {\mathbf{p}}_1\right)}^2}{k^2}\right]\frac{{\sin }^4\left(\frac{\pi \mathbf{k}\cdot {\mathbf{p}}_1}{2}\right)}{(\mathbf{k}\cdot {\mathbf{p}}_1{)}^2}.$$

$$⟨\mathbf{u}(\mathbf{x})\cdot \mathbf{u}(\mathbf{x}+\mathbf{r})⟩{|}_{uncorr}=\frac{n{L}^3}{4{\pi }^6(\ln \kappa {)}^{2}r}{\int }_0^{\infty }\mathrm{d}k\frac{\sin (2\pi kr)}{k^5}{\int }_{-1}^1\mathrm{d}\mu \left(\frac{1-{\mu }^2}{{\mu }^2}\right){\sin }^4\left(\frac{\pi k\mu }{2}\right).$$

$$⟨\hat{\mathbf{u}}(\mathbf{k})\cdot \hat{\mathbf{u}}({\mathbf{k}}^{\prime })⟩{|}_{corr}$$

$$⟨\hat{\mathbf{u}}(\mathbf{k})\cdot \hat{\mathbf{u}}({\mathbf{k}}^{\prime })⟩{|}_{corr}=\frac{\mathbf{\delta }(\mathbf{k}+{\mathbf{k}}^{\prime })}{4{\pi }^6(\ln \kappa {)}^3{k}^4}$$

$${\text{bnabla}}_{\mathbf{r}}\cdot \{[(U{\mathbf{p}}_2+{\mathbf{u}}_1)-(U{\mathbf{p}}_1+{\mathbf{u}}_2)]{\Omega }_2\}{\text{bnabla}}_{{\mathbf{p}}_1}\cdot ({\dot{\mathbf{p}}}_{12}{\Omega }_2)+{\text{bnabla}}_{{\mathbf{p}}_2}\cdot ({\dot{\mathbf{p}}}_{21}{\Omega }_2)=0.$$

$$({\mathbf{p}}_2-{\mathbf{p}}_1)\cdot {\nabla }_{\mathbf{r}}{\Omega }_2^{(1)}=-\frac{(n{L}^3{)}^2}{(4\pi {)}^2}\left[{\nabla }_{{\mathbf{p}}_1}\cdot {\dot{\mathbf{p}}}_{12}+{\nabla }_{{\mathbf{p}}_2}\cdot {\dot{\mathbf{p}}}_{21}\right],$$

$$2\pi i\mathbf{k}\cdot ({\mathbf{p}}_2-{\mathbf{p}}_1){\hat{\Omega }}_2^{(1)}=-\frac{(n{L}^3{)}^2}{(4\pi {)}^2}\left[{\nabla }_{{\mathbf{p}}_1}\cdot {\widehat{\dot{\mathbf{p}}}}_{12}+{\nabla }_{{\mathbf{p}}_2}\cdot {\widehat{\dot{\mathbf{p}}}}_{21}\right].$$

$${\hat{\Omega }}_2^{(1)}=c\delta \left[\mathbf{k}\cdot ({\mathbf{p}}_2-{\mathbf{p}}_1)\right]-\frac{(n{L}^3{)}^2}{(4\pi {)}^{2}2\pi i\mathbf{k}\cdot ({\mathbf{p}}_2-{\mathbf{p}}_1)}\left[{\nabla }_{{\mathbf{p}}_1}\cdot {\widehat{\dot{\mathbf{p}}}}_{12}+{\nabla }_{{\mathbf{p}}_2}\cdot {\widehat{\dot{\mathbf{p}}}}_{21}\right],$$

$$\underset{z\to -\infty }{\lim }{\Omega }_2^{(1)}\equiv \underset{z\to -\infty }{\lim }\int \mathrm{d}\mathbf{k}\exp \left[2\pi i\mathbf{k}\cdot \mathbf{x}\right]{\hat{\Omega }}_2^{(1)}=0.$$

$$\underset{z\to -\infty }{\lim }$$

$$c=\frac{(n{L}^3{)}^2}{32{\pi }^2|{\mathbf{p}}_2-{\mathbf{p}}_1|}\underset{k_z\to 0}{\lim }\left[{\nabla }_{{\mathbf{p}}_1}\cdot {\widehat{\dot{\mathbf{p}}}}_{12}+{\nabla }_{{\mathbf{p}}_2}\cdot {\widehat{\dot{\mathbf{p}}}}_{21}\right].$$

$$c$$

$${\hat{\Omega }}_2^{(1)}$$

$$⟨\mathbf{u}(\mathbf{x})\cdot \mathbf{u}(\mathbf{x}+\mathbf{r})⟩{|}_{corr}=\frac{3(n{L}^3{)}^2}{128{\pi }^6(\ln \kappa {)}^{3}r}{\int }_0^{\infty }\mathrm{d}k\frac{\sin (2\pi kr)}{k^6}$$

$$⟨\mathbf{u}(\mathbf{x})\cdot \mathbf{u}(\mathbf{x})⟩{|}_{corr}=\frac{3(n{L}^3{)}^2}{64{\pi }^5(\ln \kappa {)}^3}{\int }_0^{\infty }\mathrm{d}k\frac{1}{k^5}{\int }_{-1}^1\mathrm{d}\mu (1-{\mu }^2{)}^2\sin (\pi k\mu )\frac{1}{{\left(\frac{\pi }{2}k\mu \right)}^3}{\sin }^6\left(\frac{\pi }{2}k\mu \right).$$

$$\left(\frac{U\tau }{L}\right)({\mathbf{p}}_2-{\mathbf{p}}_1)\cdot {\text{bnabla}}_{\mathbf{r}}{\Omega }_2^{(1)}$$

$$\left(\frac{U\tau }{L}\right)2\pi i\mathbf{k}\cdot ({\mathbf{p}}_2-{\mathbf{p}}_1){\hat{\Omega }}_2^{(1)}$$

$$\left(\frac{U\tau }{L}\right)2\pi i\mathbf{k}\cdot ({\mathbf{p}}_2-{\mathbf{p}}_1){\hat{\Omega }}_2^{(1)}+2{\hat{\Omega }}_2^{(1)}=-\frac{(n{L}^3{)}^2}{(4\pi {)}^2}\left(\frac{U\tau }{L}\right)\left[{\text{bnabla}}_{{\mathbf{p}}_1}\cdot {\widehat{\dot{\mathbf{p}}}}_{12}+{\text{bnabla}}_{{\mathbf{p}}_2}\cdot {\widehat{\dot{\mathbf{p}}}}_{21}\right].$$

$${\hat{\Omega }}_2^{(1)}$$

$${\hat{\Omega }}_2^{(1)}{|}_{st}$$

$$⟨\hat{\mathbf{u}}(\mathbf{k})\cdot \hat{\mathbf{u}}({\mathbf{k}}^{\prime })⟩{|}_{corr}$$

$$⟨\mathbf{u}(\mathbf{x})\cdot \mathbf{u}(\mathbf{x}+\mathbf{r})⟩{|}_{corr}$$

$$⟨\mathbf{u}(\mathbf{x})\cdot \mathbf{u}(\mathbf{x})⟩{|}_{corr}$$

Full text

1 Introduction

Suspensions of rear-actuated swimming microorganisms (pushers), such as bacteria E. coli and B. subtilis, exhibit a state of large-scale coherent motion that arises, in part, due to long-ranged hydrodynamic interactions (Simha & Ramaswamy, 2002; Toner et al., 2005; Saintillan & Shelley, 2007, 2008; Underhill et al., 2008; Subramanian & Koch, 2009; Ramaswamy, 2010; Subramanian & Nott, 2011; Koch & Subramanian, 2011; Underhill & Graham, 2011; Marchetti et al., 2013; Krishnamurthy & Subramanian, 2015; Clement et al., 2016). Experiments (Dombrowski et al., 2004; Gachelin et al., 2014) and simulations (Saintillan & Shelley, 2007; Wensink et al., 2012; Krishnamurthy & Subramanian, 2015) have shown that the transition to collective dynamics occurs beyond a threshold concentration, leading to ‘bacterial turbulence’ (Wensink et al., 2012); very recent experiments, in fact, point to the dominant role of hydrodynamic interactions in this regard (Colin et al., 2019). Collective motion has important consequences for both transport and rheology, with experiments and mean-field theories having shown a reduction in viscosity leading to apparent superfluidity and unexpected shear-banding behavior (Sokolov & Aranson, 2009; Saintillan, 2010; Lopez et al., 2015; Nambiar et al., 2017; Saintillan, 2018, 2010; Nambiar et al., 2018; Laxminarsimharao et al., 2018; Guo et al., 2018). This transition has often been inferred from an anomalous enhancement in the diffusivities of passive tracer particles (Wu & Libchaber, 2000; Chen et al., 2007b; Underhill et al., 2008; Lin et al., 2011; Jepson et al., 2013; Kasyap et al., 2014; Thiffeault, 2015; Krishnamurthy & Subramanian, 2015; Stenhammar et al., 2017; Bárdfalvy et al., 2019).

A large body of theoretical work studying swimmer suspensions relies on phenomenological and mean-field models (Saintillan & Shelley, 2008; Underhill et al., 2008; Subramanian & Koch, 2009; Simha & Ramaswamy, 2002; Aranson et al., 2007; Wensink et al., 2012; Dunkel et al., 2013a, b; Marchetti et al., 2013; Heidenreich et al., 2016). For a dilute swimmer suspension, mean-field theory has shown that it is the long-ranged hydrodynamic interactions (Underhill et al., 2008; Saintillan & Shelley, 2008; Subramanian & Koch, 2009) between pushers and the mutually reinforcing orientation and velocity fluctuations that are responsible for the said transition to collective motion. In this paper, we go beyond this limiting mean-field assumption, and demonstrate the crucial role played by hydrodynamically induced swimmer-correlations. The interplay of swimming with long-ranged hydrodynamic interactions is shown, for the first time, to lead to a larger velocity variance, and thence, larger tracer diffusivities, in suspensions of run-and-tumble pushers vis-a-vis pullers (the swimming mechanism for pullers is front-actuated, as is typically the case for algae; see Kessler (1986); Saintillan & Shelley (2007); Lauga & Powers (2009); Saintillan (2010); Koch & Subramanian (2011); Guasto et al. (2012); Krishnamurthy & Subramanian (2015); Goldstein (2015)). Interactions also lead to a logarithmic divergence of the velocity variance in the straight-swimmer limit, that has been observed in earlier simulations (Underhill & Graham, 2011). The prediction of larger velocity fluctuations in pusher suspensions is also consistent with more recent simulations which have revealed a bifurcation of the fluid velocity variance in pusher and puller suspensions beyond a threshold (Krishnamurthy & Subramanian, 2015; Stenhammar et al., 2017).

In a homogeneous Stokesian suspension of sedimenting particles with a random microstructure, the velocity variance is predicted to diverge linearly with system size (Caflisch & Luke, 1985; Hinch, 1988; Koch & Shaqfeh, 1991; Nicolai et al., 1995; Ladd, 1996; Segre et al., 1997; Ladd, 1997; Levine et al., 1998; Ramaswamy, 2001; Guazzelli & Hinch, 2011; Goldfriend et al., 2017). The resolution of this divergence has been a long-standing theoretical challenge (Hinch, 1988; Ramaswamy, 2001; Guazzelli & Hinch, 2011). It arises due to the long-ranged O($1/r$) disturbance velocity fields of the individual particles acting as point-forces (monopoles), as shown in figure 1a. The fluid velocity variance at a given point $\boldsymbol{x}$ may be estimated as $\langle\boldsymbol{u}(\boldsymbol{x})\cdot\boldsymbol{u}(\boldsymbol{x})\rangle\equiv\int nL^{3}U_{s}^{2}(1/r)^{2}\mathrm{d}\boldsymbol{r}\sim(nL^{3})U_{s}^{2}L_{box}/L$. Here, $U_{s}$ is the mean sedimenting speed, $L_{box}$ the system size and $nL^{3}$ the hydrodynamic volume fraction, with $n$ being the particle number density, and $L$ a characteristic particle size. Unlike passive particles, microswimmers, both pushers and pullers, act as force-dipoles in the far-field, leading to a disturbance velocity field that decays more rapidly as O($1/r^{2}$). An argument along lines similar to that for the passive Stokesian suspension above gives $\langle\boldsymbol{u}(\boldsymbol{x})\cdot\boldsymbol{u}(\boldsymbol{x})\rangle\equiv\int nL^{3}U^{2}(1/r)^{4}\mathrm{d}\boldsymbol{r}\sim(nL^{3})U^{2}$, $U$ being the swimming speed, implying that the velocity variance remains finite in the limit $nL^{3}\ll 1$ with correlations between swimmers being neglected at leading order (the divergence at small $r$ implied in the scaling integral above is regularized once the finite swimmer size is accounted for); the covariance $\langle\boldsymbol{u}(\boldsymbol{x})\cdot\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle\sim(nL^{3})\int U^{2}/(r-r^{\prime})^{4}\mathrm{d}\boldsymbol{r}^{\prime}$ at this order exhibits an O($1/r$) decay for $r\gg L$ (Underhill et al., 2008; Underhill & Graham, 2011). Introducing pair-level correlations in straight-swimmer suspensions, as shown in figure 1b, leads to a pair-orientation probability density that decays as O($1/r^{2}$) in the far-field. In this scenario, where a given swimmer interacts pair-wise with a cloud of swimmers surrounding it, one has the fluid velocity variance due to ‘each cloud’ scaling as $[\boldsymbol{u}(\boldsymbol{x})\cdot\boldsymbol{u}(\boldsymbol{x})]_{cloud}\equiv nL^{3}\int U^{2}r^{\prime 2}\mathrm{d}r^{\prime}/(r^{4}r^{\prime 2})\sim nL^{3}U^{2}/r^{3}$, which when integrated over all such correlation clouds gives $\langle\boldsymbol{u}(\boldsymbol{x})\cdot\boldsymbol{u}(\boldsymbol{x})\rangle\equiv(nL^{3})^{2}\ln(L_{box}/L)$; a logarithmically divergent variance. For passive suspensions, early theoretical analyses have predicted a Debye-like screening of the long-ranged hydrodynamic interactions at distances of O$(nL^{2})^{-1}$ (Koch & Shaqfeh, 1991), leading to a finite variance. In contrast, the divergence in experiments is cut off due to the development of a container-scale stratification (Luke, 2000; Ramaswamy, 2001; Guazzelli & Hinch, 2011). Clearly, the fact that correlations in active suspensions (of straight-swimmers) act to yield a divergent variance is in sharp contrast to the nature of velocity fluctuations in passive suspensions, and highlights the novel consequences of activity (swimming).

The manuscript is broadly organized into three sections. In §2, we examine the fluid velocity covariance. Following the formulation of the general expression, in §2.1, we determine the covariance for suspensions of non-interacting swimmers, in which case the covariance is independent of the mechanism by which the swimmer orientation decorrelates. Next, in §2.2, we examine the covariance for the more involved case of a dilute suspension of interacting swimmers. Here, it is necessary to discuss the case of straight-swimmers and run-and-tumble particles (RTPs) separately, and this is done in §2.2.1 and §2.2.2, respectively. The calculation of the correlated contribution, at O$(nL^{3})^{2}$, requires the steady state pair probability density in position-orientation (phase) space, for both straight-swimmers and RTPs, which in turn requires obtaining an expression for the rotation rate experienced by a slender swimmer, due to the disturbance velocity field generated by another; the latter is derived in appendix A. For straight swimmers, the phase space probability density is directly obtained in terms of generalized functions. For RTPs, the pair-probability density is again obtained in closed form, now as a function of the dimensionless parameter $U\tau/L$ which measures the run length in units of the swimmer size. Straight-swimmers correspond to the (singular) limit $U\tau/L\to\infty.$ The resulting expression for the correlated component of the variance is shown to scale linearly with $U\tau/L$ in the limit $U\tau/L\ll 1$, while increasing logarithmically in the limit $U\tau/L\gg 1$, corresponding to persistent swimmers. For $U\tau/L\ll 1$, the covariance in suspensions of RTPs transitions directly from the variance plateau for $r\ll 1$ to an O($1/r$) far-field decay for $r\gg 1$. However, for large $U\tau/L$, we observe a weak decay of the covariance in the interval $1\ll r\ll U\tau/L$ that is intermediate between the variance plateau ($r\ll 1$) and the aforementioned O($1/r$) far-field decay regime ($r\gg U\tau/L$); the covariance remains non-decaying in the singular straight-swimmer limit. For RTPs in the limit $U\tau/L\gg 1$, we present an alternate derivation of both the variance and the covariance using a matched asymptotic expansions approach in §2.2.3, which is shown to compare well with the exact result for $U\tau/L\geq$ O(1). Importantly, this approach shows that the aforementioned weak intermediate decay of the covariance for $U\tau/L\gg 1$ is a logarithmic one. Next, in §3, we study the diffusivity of immersed passive tracers. Here too, we first consider suspensions of non-interacting swimmers in §3.1, and then interacting swimmers in §3.2. Finally, in §4, we present concluding remarks and a course for future work. We also briefly discuss the nature of pair-orientation correlations in Appendix B.

2 The fluid velocity covariance

We begin with a discrete formulation applicable to a suspension of $N$ slender swimmers, each having a length $L$ and swimming with speed $U$ directed along its axis. The configuration of the swimmer suspension is characterized by the positions ($\boldsymbol{x}_{i}$) and orientations ($\boldsymbol{p}_{i}$) of the $N$ swimmers, with the swimmer number density field being defined as $c(\boldsymbol{x},t)=\sum_{i=1}^{N}\boldsymbol{\delta}(\boldsymbol{x}-\boldsymbol{x}_{i})$. The disturbance velocity and pressure fields $\boldsymbol{u}_{i}$ and $P_{i}$, due to the $i^{th}$ swimmer, at leading logarithmic order, satisfy the Stokes equations forced by a line distribution of Stokeslets along the swimmer axis, and the equation of continuity: {subeqnarray} -\bnablaP_i + η∇^2 ui(x) &= ∫-L/2^L/2 f(s) p_i δ(x -x_i - sp_i) ds,

\bnabla⋅u_i = 0.

Here, $\eta$ is the viscosity of the suspending fluid and $\boldsymbol{\delta}(\boldsymbol{z})$ represents the Dirac-delta function. Assuming the swimmers to be fore-aft symmetric and force-free, the linear force density along the swimmer axis (the axial coordinate being $s$) may be expressed as $f(s)=\eta Usgn(s)/\ln\kappa$. For non-fore-aft symmetric swimmers, there is an O(1) change in the force density, and thence the disturbance fluid velocity it generates, however, none of the principal conclusions detailed below change on account of swimmer asymmetry. The velocity disturbance satisfying (2) may be written in terms of the Oseen-Burgers tensor $\mathsfbi{G}(\boldsymbol{x})=1/(8\pi\eta r)[\mathsfbi{I}+\boldsymbol{xx}/r^{2}]$ as:

[TABLE]

and therefore, the suspension velocity field is expressible as:

[TABLE]

where the time dependence comes from the evolving swimmer positions and orientations. Now, the fluid velocity covariance in a swimmer suspension is defined as:

[TABLE]

where the angular brackets denote an ensemble average ($\langle\boldsymbol{\cdot}\rangle$) over the configurations of swimmers $i,j$. Since the focus in this section is on the single-time covariance, we will neglect mention of the time dependence from hereon. Within the continuum framework ($N\gg 1$), the above average is expressible in terms of the configurational (position-orientation) probability density function. To see this, one may split the double summation in (3) into two distinct contributions; one containing only the diagonal terms, with $i=j$, which involves the singlet probability density defined as $\Omega_{1}=\langle\sum_{i=1}^{N}\delta(\boldsymbol{x}_{1}-\boldsymbol{x}_{i})\delta(\boldsymbol{p}_{1}-\boldsymbol{p}_{i})\rangle$; and a second contribution containing the off-diagonal terms $i\neq j$ that involves the pair probability density defined as $\Omega_{2}$ = $\langle\sum_{i,j=1;i\neq j}^{N}\boldsymbol{\delta}(\boldsymbol{x}_{1}-\boldsymbol{x}_{i})\boldsymbol{\delta}(\boldsymbol{p}_{1}-\boldsymbol{p}_{i})\boldsymbol{\delta}(\boldsymbol{x}_{2}-\boldsymbol{x}_{j})\boldsymbol{\delta}(\boldsymbol{p}_{2}-\boldsymbol{p}_{j})\rangle$. The expression (3) above may be written in terms of these probability densities as:

[TABLE]

Since we examine a spatially homogeneous isotropic suspension, one has $\Omega_{1}=nL^{3}/(4\pi)$, and the pair probability density function $\Omega_{2}$ above only depends on the relative separation $\boldsymbol{r}=\boldsymbol{x}_{2}-\boldsymbol{x}_{1}$ of the pair of swimmers under consideration. The assumption of a statistically homogeneous and isotropic suspension of swimmers may only remain valid over a time smaller than O$(nUL^{2})^{-1}$. Evidence from simulations of spherical squirmers (Alarcón & Pagonabarraga, 2013; Evans et al., 2011; Oyama et al., 2017; Alarcón et al., 2017) does suggest that hydrodynamic interactions, even in a dilute setting, can induce global orientational order over times long compared to the aforementioned scale. For slender swimmers, this assumption is not a limiting one, and the isotropic state appears to be stable to orientation perturbations at leading logarithmic order. We therefore proceed with the above expression for $\Omega_{1}$.

It is convenient to solve for the pair-probability density in Fourier space, and towards this end, the Fourier transformed covariance is expressible as:

[TABLE]

Here, the Fourier transformed quantities are defined as $\hat{f}(\boldsymbol{k})=\int\exp(-2\pi i\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x})f(\boldsymbol{x})\mathrm{d}\boldsymbol{x}$, and the factor $\delta(\boldsymbol{k}+\boldsymbol{k}^{\prime})$ represents the translational invariance due to the assumed spatial homogeneity. Using the aforementioned form of $f(s)$, and the Fourier transformed Oseen-Burgers tensor, $\hat{\mathsfbi{G}}(\boldsymbol{k})=1/(4\pi^{2}\eta k^{2})(\mathsfbi{I}-\hat{\boldsymbol{k}}\hat{\boldsymbol{k}})$, we get:

[TABLE]

as the Fourier transformed velocity covariance in a suspension of slender swimmers. Note that (6) has been rendered non-dimensional by using $L$ and $U$, respectively, as length and velocity scales.

In what follows, we characterize the variance and covariance first for a suspension of non-interacting swimmers in §2.1, and then, for suspensions of interacting straight and run-and-tumble swimmers, respectively, in §2.2.1 and §2.2.2. For the latter case, we also present an alternate derivation of the covariance using a matched asymptotic expansions approach in §2.2.3, applicable to RTP suspensions with large swimmer run lengths.

2.1 Suspensions of non-interacting swimmers

In the absence of interactions, only the term involving the single swimmer probability density is relevant, so that (6) takes the simplified form:

[TABLE]

An inverse Fourier transform of (7) yields:

[TABLE]

Owing to isotropy, the covariance only depends on the scalar distance $r$. The variance is obtained on setting $r=0$ in (8), which gives $\langle\boldsymbol{u}(\boldsymbol{x})\cdot\boldsymbol{u}(\boldsymbol{x})\rangle|_{uncorr}=nL^{3}/[96\pi(\ln\kappa)^{2}]$. In the far-field ($r\gg 1$), we recover the limiting form of the covariance for point dipoles, given by: $\langle\boldsymbol{u}(\boldsymbol{x})\cdot\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{uncorr}=nL^{3}/[120\pi(\ln\kappa)^{2}r]$ (Underhill & Graham, 2011). The aforementioned variance and far-field covariance have been plotted alongside (8), evaluated numerically, in figure 2. Note that the covariance, as given by (8), depends only on the single-swimmer statistics at a given instant, and therefore, remains the same both for straight-swimmers and run-and-tumble swimmers provided the tumbles in the latter case are assumed to be instantaneous (so that any disturbance field generated during a tumble is neglected).

2.2 Suspensions of interacting swimmers

Herein, we consider the covariance in a suspension of hydrodynamically interacting slender swimmers to O$(nL^{3})^{2}$, with the contribution at this order requiring the analysis of pairwise interactions between slender swimmers. Therefore, the additional correlated contribution involving the pair-probability density in (6) needs to be calculated. This correlated contribution is given by:

[TABLE]

where $\hat{\Omega}_{2}$ is to be determined.

Considering slender swimmers, with an aspect ratio $\kappa\gg 1$, leads to logarithmically weak interactions on length scales of O($L$) (scales that contribute dominantly to the velocity variance, as may be verified posteriori), and this allows one to expand $\Omega_{2}$ as a series in $(\ln\kappa)^{-1}$: $\Omega_{2}=\Omega_{2}^{(0)}+1/(\ln\kappa)\Omega_{2}^{(1)}+\ldots$; here, $\Omega_{2}^{(0)}$ is the product of two singlet probability densities, $(nL^{3})^{2}/(4\pi)^{2}$, and denotes the uncorrelated base state at leading order. $\Omega_{2}^{(1)}\equiv\Omega_{2}^{(1)}(\boldsymbol{r},\boldsymbol{p}_{1},\boldsymbol{p}_{2})$ gives the correlated probability, at O($1/\ln\kappa$), of finding two swimmers with orientations $\boldsymbol{p}_{1}$ and $\boldsymbol{p}_{2}$, separated by $\boldsymbol{r}$. Typical bacteria such as E. coli and B. subtilis are quite slender (Berg, 2004; Gachelin et al., 2014), and this assumption reasonably reproduces experimentally-measured disturbance velocity fields to distances of O($L$) (Kasyap et al., 2014). Note that the contribution due to $\Omega_{2}^{(0)}$ in (9) is identically zero. This is as expected since the velocity variance in an uncorrelated swimmer suspension must be proportional to $nL^{3}$, with no constraints of diluteness involved.

Therefore, to leading logarithmic order, one only needs to determine $\Omega_{2}^{(1)}(\boldsymbol{p}_{1},\boldsymbol{p}_{2};\boldsymbol{k})$ to characterize pair correlations in the swimmer suspension, and thence, determine the correlated contribution to the covariance in (9), which is now expressible as:

[TABLE]

At this order, there is a crucial difference in the pair-correlations that develop in suspensions of straight-swimmers and RTPs, and therefore, the two cases are treated separately.

2.2.1 Suspensions of interacting straight-swimmers

For straight-swimmers, $\Omega_{N}$ satisfies the Liouville equation, and evolves due to swimming, and due to convection and rotation of each swimmer by the disturbance velocity fields due to the remaining swimmers. In the dilute limit, integrating over the degrees of freedom of the remaining $N-2$ swimmers, while neglecting three-swimmer interactions, we obtain the equation for the pair probability density $\Omega_{2}$, at steady state, as:

[TABLE]

The terms within braces in (11) denote the convection of $\Omega_{2}$ by the relative velocity of the swimmer pair that includes contributions due to both swimming ($U\boldsymbol{p}_{1}$, $U\boldsymbol{p}_{2}$) and the disturbance velocity fields ($\boldsymbol{u}_{2}$, $\boldsymbol{u}_{1}$). The third and fourth terms denote rotation of the swimmer orientations due to the disturbance velocity fields, with $\dot{\boldsymbol{p}}_{ij}$ denoting the rotation of swimmer $i$ by the disturbance velocity field due to swimmer $j$.

Using the aforementioned series expansion of $\Omega_{2}$ and neglecting the convection by the O$(\ln\kappa)^{-1}$ disturbance velocity fields in (11), at leading logarithmic order, pair-correlations develop along straight-swimming trajectories. Thus, at O($1/\ln\kappa$), for hydrodynamically interacting slender swimmers, $\Omega_{2}^{(1)}$ satisfies:

[TABLE]

on applying the non-dimensionalization mentioned below (6). In Fourier space, (12) takes the form:

[TABLE]

The solution for $\hat{\Omega}_{2}^{(1)}$ from (13) is expressible as:

[TABLE]

where the first term on the right side, involving the Dirac delta function, is the homogeneous solution, and the expressions for $\hat{\dot{\boldsymbol{p}}}_{12}$ and $\hat{\dot{\boldsymbol{p}}}_{21}$ that appear in the particular solution are given in appendix A. The constant $c$ is determined from the fact that the swimmers are uncorrelated prior to interaction; that is, $\Omega_{2}^{(1)}\to 0$ for swimmers when they are infinitely separated in the upstream direction. Choosing a Cartesian coordinate system with the $z$ axis along the relative swimming velocity vector, that is $\hat{\boldsymbol{z}}=(\boldsymbol{p}_{2}-\boldsymbol{p}_{1})/|\boldsymbol{p}_{2}-\boldsymbol{p}_{1})|$, implies that $z\to-\infty$ denotes upstream infinity. Therefore, the aforementioned constraint of uncorrelated swimmers infinitely far upstream may be given as:

[TABLE]

One can represent $\boldsymbol{x}=(\boldsymbol{x}^{\perp},z)\equiv(\boldsymbol{x}\boldsymbol{\cdot}[\mathsfbi{I}-\hat{\boldsymbol{z}}\hat{\boldsymbol{z}}],z)$, and since (15) is valid for arbitrary $\boldsymbol{k}^{\perp}$, one need only consider the integral over $k_{z}$, whence one obtains:

[TABLE]

While the first term within brackets in (LABEL:omega21d) may be evaluated readily, for $z\to-\infty$, the dominant contribution to the second term arises when $k_{z}\rightarrow 0$ such that $k_{z}z$ remains finite. This results in the following expression for $c$:

[TABLE]

We now choose $\boldsymbol{x}^{\perp}=(x,y)$ such that $\hat{\boldsymbol{x}}=(\boldsymbol{p}_{2}\wedge\boldsymbol{p}_{1})/|\boldsymbol{p}_{2}\wedge\boldsymbol{p}_{1}|$ and $\hat{\boldsymbol{y}}=(\boldsymbol{p}_{2}+\boldsymbol{p}_{1})/|\boldsymbol{p}_{2}+\boldsymbol{p}_{1}|$, are the unit vectors in the plane orthogonal to $\hat{\boldsymbol{z}}$. For $k_{z}\to 0$, we have, $\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_{1}=\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_{2}$ with each of them being given by $k_{y}|\boldsymbol{p}_{2}+\boldsymbol{p}_{1}|/2$. Using these relations in (17), and after some algebra, $c$ is given by:

[TABLE]

From (68) and (18) we finally get the O($1/\ln\kappa$) correction to the Fourier transformed pair-probability density, in (14), as:

[TABLE]

where the first term on the right side of (LABEL:eq:omega21Full) is the homogeneous solution, with $c$ being replaced by (18); the integral involving the argument of $PV[\cdot]$ needs to be interpreted as a Cauchy principal value integral. The analysis of the homogeneous solution above is crucial, since it is this term alone that contributes to the covariance.

On using (LABEL:eq:omega21Full) in (10), and on applying the inverse Fourier transform, we obtain the O$(nL^{3})^{2}$ fluid velocity covariance in an interacting straight-swimmer suspension to be:

[TABLE]

Again, the variance is obtained by setting $r=0$ in (20), which gives:

[TABLE]

Notice that the integrand for $k$ in (21) scales as $1/k$ as $k\to 0$, implying a logarithmic divergence. Thus, the O$(nL^{3})^{2}$ contribution to the covariance is not well defined for any $r$, and in particular, results in a logarithmically divergent variance on setting $r=0$. The physical arguments for this logarithmic divergence may now be laid out in a little more detail than in the introduction. As mentioned therein, the divergence arises from the slowO($1/r^{2}$) decay of the pair-probability density in the far-field. This can be readily inferred from (LABEL:eq:omega21Full), where power counting in $k$, as $k\to 0$ suggests that $\hat{\Omega}_{2}^{(1)}\sim$ O($1/k$) (which in real space would imply $\Omega_{2}^{(1)}\sim$ O($1/r^{2}$) in the far-field). Alternatively, one may also observe this from the physical space form, given by (12), wherein integrating along the trajectory of relative swimming, that is: $\Omega_{2}^{(1)}=-(nL^{3})^{2}/[(4\pi)^{2}U|\boldsymbol{p}_{2}-\boldsymbol{p}_{1}|]\int_{-\infty}^{z}\mathrm{d}z^{\prime}(\bnabla_{\boldsymbol{p}_{1}}\boldsymbol{\cdot}\dot{\boldsymbol{p}}_{12}+\bnabla_{\boldsymbol{p}_{2}}\boldsymbol{\cdot}\dot{\boldsymbol{p}}_{21})$, with $\dot{\boldsymbol{p}}_{ij}\sim$ O($1/r^{3}$), the velocity gradient scaling in the far-field, leads to the O($1/r^{2}$) decay. In this scenario, where a given swimmer interacts pair-wise with a cloud of swimmers surrounding it (see figure 1b), one has the fluid velocity variance due to a “correlated cloud”, of radius $r$, scaling as $[\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x})]_{cloud}\equiv nL^{3}\int U^{2}r^{\prime 2}\mathrm{d}r^{\prime}/(r^{4}r^{\prime 2})\sim nL^{3}U^{2}/r^{3}$, which when integrated over all such correlation clouds in position space, gives $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x})\rangle|_{corr}\equiv(nL^{3})^{2}\ln(L_{box}/L)$; a logarithmically divergent variance.

In practice, the logarithmic divergence above will be cut off at an appropriate screening length $L_{screen}$, such that the total variance, including the uncorrelated, given by (8), and the correlated contribution above, takes the form $c_{1}(nL^{3})+(nL^{3})^{2}\{c_{2}\ln[L_{screen}/L]+c_{3}\}$, with the constants $c_{i}$’s being functions of swimmer aspect ratio. The screening length depends on the particular scenario. In box-size-limited simulations, the largest admissible wavelength is set by the computational domain, so $L_{screen}=k_{min}^{-1}\equiv L_{box}$. In figure 3, we highlight this box-size-dependent divergence of the correlated contribution as a function of $k_{min}$, with $k_{min}^{-1}$ replacing 0 as the lower limit of the Fourier integral in (21), and thereby, enforcing a long wavelength cut-off. For suspensions of interacting straight-swimmer, the box-size limitation arises provided the swimmer mean free path, which is O$(nL^{2})^{-1}$, is larger than $L_{box}$. For $L_{box}\gg(nL^{2})^{-1}$, these suspensions become linearly unstable to orientation perturbations, and one expects a transition to collective motion (Saintillan & Shelley, 2008; Subramanian & Koch, 2009); in simulations, this is characterized by a box size dependent sharp increase in the fluid velocity variance and tracer diffusivities (Saintillan & Shelley, 2011; Krishnamurthy & Subramanian, 2015; Stenhammar et al., 2017). As inherent in our homogeneous suspension assumption, our analysis therefore caters to the former scenario (upto $L_{box}\sim(nL^{2})^{-1}$ or for time scales $\sim$ O$(nUL^{2})^{-1}$ for larger box sizes). Simulations of interacting regularized force dipole swimmers carried out by Underhill & Graham (2011), at $nL^{3}\sim 10^{-2}$ (and in the regime $L_{box}<(nL^{2})^{-1}$), do point to a logarithmic divergence of the variance, as may be inferred from the O($1/r^{2}$) far-field decay of the spatial orientation correlations, and thence, the pair-probability density (see figure 7 therein). Interestingly, the authors observe a logarithmically diverging variance in simulations, even at higher volume fractions ($nL^{3}\sim 10^{-1}$), for a range of box sizes $L_{box}\geq(nL^{2})^{-1}$. For slender swimmers, owing to the weak pair correlations (at least by a factor of $\ln\kappa$), the orientation decorrelation is expected to occur over a logarithmically large number of pair interactions, implying a larger effective mean free path. Thus, our analysis may remain valid for suspension of slender swimmers at even higher volume fractions than for the regularized point dipoles considered in earlier simulations (Underhill et al., 2008; Underhill & Graham, 2011; Stenhammar et al., 2017).

For real bacteria, intrinsic decorrelation mechanisms such as rotary diffusion or run-and-tumble dynamics, lead to $L_{screen}$ = $U/D_{r}$ or $U\tau$, where $D_{r}$ is the rotary diffusivity and $\tau$ the mean run duration. Note that the rotary diffusion above could have an entirely hydrodynamic origin. For slender straight-swimmer suspensions not limited by box size, one expects the logarithmic divergence to nevertheless be cut off at $L_{screen}\sim$ O$(U/D_{r}^{h})$, where $D_{r}^{h}\sim$ O$(nUL^{2}/(\ln\kappa)^{2})$ is a hydrodynamic rotary diffusivity arising from logarithmically weak pairwise interactions between slender swimmers each of which lead to an O($1/\ln\kappa$) angular displacement. This diffusivity has been calculated earlier (Subramanian & Koch, 2009), and accounting for this finite $D_{r}^{h}$ leads to a suspension velocity variance of the form: $c_{1}(nL^{3})+(nL^{3})^{2}\{c_{2}\ln[(\ln\kappa)^{2}/(nL^{3})]+c_{3}\}$. As will be seen in the next subsection, the covariance integral, (20), in the strict straight-swimmer limit ($U\tau/L=\infty$) is non-decaying. However, the straight-swimmer limit is better interpreted as a limiting case of RTPs, in which case one finds an intermediate asymptotic regime where the covariance exhibits a weak logarithmic decay with separation.

2.2.2 Suspensions of interacting RTPs

The kinetic equation for the pair-probability density for a suspension of RTPs includes additional terms compared to (12) that describe the tumble dynamics of the individual bacteria, and is given by:

[TABLE]

The bracketed terms on the left-hand side of (22) correspond to the swimmers undergoing tumbles in accordance with Poisson statistics (Subramanian & Koch, 2009; Othmer et al., 1988), with there being no correlation between pre- and post-tumble orientations (random tumbles); $\tau$ is the mean run duration (that is, the mean interval between successive random tumbles). In the absence of additional forcing, tumbling causes the single-swimmer distribution to relax to isotropy on a time scale of O($\tau$). Here, $U\tau/L$ is a non-dimensional mean run length; thus, $U\tau/L\rightarrow\infty$ in the limit of straight-swimmers. Again, on Fourier transforming, one gets:

[TABLE]

One can obtain an analytical solution for $\hat{\Omega}_{2}^{(1)}(\boldsymbol{p}_{1},\boldsymbol{p}_{2};\boldsymbol{k})$ in (23), by determining the Green’s function of the linear operator on the left-hand side. This requires one to first solve the time dependent operator, and then obtain the steady state pair probability density by taking the long time limit. The solution procedure involves representing the Green’s function as a superposition of the eigenfunctions. The singular nature of the linear operator implies that both the discrete and (singular) continuous spectrum need to be examined. For the sake of brevity, we will be reporting this detailed Green’s function analysis elsewhere (Garg et al., 2020).

For purposes of the covariance calculation, it is worth noting that the symmetry of the orientation dilatation forcing terms on the right-hand side of (23) ensures that one does not require information of the full Green’s function. A representation of $\hat{\Omega}_{2}^{(1)}$ as a convolution of the Green’s function with the right-hand side of (23) (as briefly mentioned above) implies that only those eigenfunctions, that the orientation dilatation terms project onto, contribute. To see this, we first choose a coordinate system with its polar axis along $\hat{\boldsymbol{k}}$, such that $\mu_{i}\equiv\cos\theta_{i}=\hat{\boldsymbol{k}}\boldsymbol{\cdot}\boldsymbol{p}_{i}$, with $\theta_{i}$ being the polar angle, and $\phi_{i}$ representing the azimuthal angle measured in a plane perpendicular to $\boldsymbol{k}$. In this coordinate system, the orientation dilatation terms are proportional to $\cos\phi_{i}$, and hence, vanish when integrated in orientation space. Thus, one need only solve for a reduced pair probability density, which accounts for eigen modes proportional to $\cos\phi_{i}$. Now, eigen modes proportional to $\cos(n\phi_{i})$ or $\sin(m\phi_{i})$ remain unaffected by the integral terms for $m,n\neq 0$, since the latter vanish for all such cases; physically, the orientation dilatation terms lead to pure-orientation correlations without concentration (number density) perturbations. Thus, one may set $\int\hat{\Omega}_{2}^{(1)}\mathrm{d}\boldsymbol{p}_{1}=\int\hat{\Omega}_{2}^{(1)}\mathrm{d}\boldsymbol{p}_{2}=0$, and the governing equation for $\hat{\Omega}_{2}^{(1)}$ reduces to:

[TABLE]

Solving (24), yields the following exact expression:

[TABLE]

The pair-correlation function for RTPs above does not involve any concentration fluctuations, as may seen by integrating (LABEL:eq:omega21RTP5) over orientation space. From (LABEL:eq:omega21RTP5), one can also infer that for large $U\tau/L$ (the straight-swimmer limit), $\hat{\Omega}_{2}^{(1)}\sim$ O($1/k$) as $k\to 0$, implying that $\Omega_{2}^{(1)}\sim$ O($1/r^{2}$) in the far-field; however, for rapid tumblers ($U\tau/L\ll 1$), $\Omega_{2}^{(1)}$ exhibits a more rapid far-field decay than an algebraic one. The polar and nematic correlations between swimmer orientations may be derived using the above expression for $\hat{\Omega}_{2}^{(1)}$ and are discussed in Appendix B.

One may recover the straight-swimmer limiting case, obtained in §2.2.1, from (LABEL:eq:omega21RTP5), in the limit of large $U\tau/L$. To see this, we first rewrite (LABEL:eq:omega21RTP5) as:

[TABLE]

where $\vartheta=L/(\pi kU\tau)$. In the limit $U\tau/L\gg 1$, $\vartheta\ll 1$, implying that $i\vartheta\to 0$. Interpreting the bracketed term on the right-hand side of (LABEL:eq:omega21RTP5Recovery1), in this limit, needs careful consideration. Naively setting $\vartheta=0$ only leads to the particular solution corresponding to the straight-swimmer pair probability density, as given in (LABEL:eq:omega21Full). To obtain the correct answer in the straight-swimmer limiting form of $\hat{\Omega}_{2}^{(1)}$, one needs to apply the Plemelj-Sokhotski formula (Gakhov, 1966; Fokas & Ablowitz, 2012), which may written as $\lim_{\epsilon\to 0}1/(x+i\epsilon)=-\delta(x)+PV[1/x]$, and this yields

[TABLE]

One may now readily note that the term involving $\boldsymbol{\delta}(\hat{\boldsymbol{k}}\boldsymbol{\cdot}(\boldsymbol{p}_{2}-\boldsymbol{p}_{1}))$ in (LABEL:eq:omega21RTP5Recovery2) constitutes the homogeneous solution, whereas, the term involving $PV[\cdot]$ is the particular solution of the straight-swimmer pair probability density given by (LABEL:eq:omega21Full).

Now, using the coordinate system described above, the Fourier transformed covariance from (10) is expressible as:

[TABLE]

On inverse Fourier transforming (28) one obtains the correlated contribution to the fluid velocity covariance for interacting RTPs to be:

[TABLE]

where $j_{0}(z)=\sin z/z$ is the spherical Bessel’s function of the first kind (Gradshteyn & Ryzhik, 1980). The correlated variance is obtained by setting $r=0$ in (LABEL:covarRTP), and is expressible as:

[TABLE]

In figure 4, we plot the correlated variance given by (LABEL:varRTP_corr). In accordance with earlier scaling arguments, the correlated variance in suspensions of interacting RTPs takes the form, $(nL^{3})^{2}[c_{1}\ln[U\tau/L]+c_{2}]$ for large $U\tau/L$. In contrast, for $U\tau/L\ll 1$, owing to the more rapid decay of $\Omega_{2}^{(1)}$, the correlated variance scales as, $c_{1}(U\tau/L)(nL^{3})^{2}$ (see inset of figure 4). Next, in figure 5, the correlated covariance is plotted for pusher-type suspensions. For rapid tumblers ($U\tau/L\ll 1$), the covariance directly transitions from the initial variance plateau to an O($1/r$) decay for $r\gg 1$. In this limit, one can obtain a simplified form of the correlated covariance from (LABEL:covarRTP), on setting $i(U\tau/L)k(\mu_{2}-\mu_{1})$ to zero, which gives:

[TABLE]

Setting $r=0$ in (31), readily yields the correlated variance to be of the form $c_{1}(U\tau/L)(nL^{3})^{2}$, consistent with the scaling arguments above. In contrast, for $U\tau/L\gg 1$, there emerges a weak intermediate regime for $1\ll r\ll U\tau/L$, which delays the onset of the eventual O($1/r$) far-field decay (this weak scaling of the covariance in the intermediate region is a logarithmic one, scaling as O($\ln[r/(U\tau)]$), and will be derived in §2.2.3). One may simplify (LABEL:covarRTP) to obtain the limiting form of this latter far-field dipole asymptote, which may be expressed as: $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{corr}=3(nL^{3})^{2}(U\tau/L)(28800\pi r)^{-1}$. The inset in figure 5 plots the correlated covariance as a function of $r(U\tau/L)^{-1}$, highlighting the far-field collapse.

One may now write down the expression for the total variance in a suspension of interacting RTPs by adding the variance expression for suspensions of non-interacting swimmers given below (8) (see §2.1), and the correlated contribution to the variance given by (LABEL:varRTP_corr) above. This is expressible as:

[TABLE]

Expectedly, the total variance scales as $c_{1}(nL^{3})+(nL^{3})^{2}[c_{2}\ln[U\tau/L]+c_{3}]$ for large $U\tau/L$, and as $c_{1}(nL^{3})+c_{2}(U\tau/L)(nL^{3})^{2}$ for $U\tau/L\ll 1$. An important implication of including correlations is the scaling of the variance with the dipole strength. The O($nL^{3}$) contribution of the variance has a dimensional scaling of $U^{2}$, implying that it scales as the square of the dipole strength $D$ (measured in units of $\eta UL^{2}$). However, the correlated contribution to the variance is proportional to $U^{3}$, and thence, to $D^{3}$. Thus, the variance for a suspension of pairwise interacting RTPs depends on the sign of $D$ ($D<0$ for pushers and $>0$ for pullers), and therefore, on the swimming mechanism. Consistent with recent simulations (Krishnamurthy & Subramanian, 2015; Stenhammar et al., 2017), at O$(nL^{3})^{2}$, (LABEL:varRTP_Total) predicts enhanced fluctuations in pusher suspensions as shown in figure 6.

The total covariance in a suspension of interacting pusher-type RTPs can be similarly expressed by combining (8) and (LABEL:covarRTP), which gives:

[TABLE]

Figure 7 shows a plot of the total covariance as given by (33). Note that the chosen $nL^{3}$ is 5 (and, therefore, a little greater than unity); there is evidence, from earlier calculations for suspensions of slender fibers, that the dilute theory, with pair-interactions accounted for, remains valid for $nL^{3}$ up until around 7 (Mackaplow & Shaqfeh, 1996). In figure 7, for large $U\tau/L$, while the variance plateau grows slowly with increase in $U\tau/L$, owing to the $\ln U\tau/L$ scaling of the correlated component, the far-field O($1/r$) scaling grows faster on account of the linear scaling with $U\tau/L$ (see the expression for the far-field asymptote below (31)). In addition, for large $U\tau/L$, there is an intermediate regime of shallow decay for separations $1\ll r\ll U\tau/L$, owing to the O($\ln[r/(U\tau)]$) scaling of the correlated covariance, which delays the onset of the eventual O($1/r$) far-field decay of the covariance. Therefore, as $U\tau/L\to\infty$, one asymptotes to a non-decaying covariance, as expected for the singular straight-swimmer limit discussed earlier.

2.2.3 Suspensions of interacting RTPs - large $U\tau/L$ covariance using matched asymptotic expansions

Until now, we have established the detailed expressions for the covariance in suspensions of hydrodynamically interacting straight-swimmers in §2.2.1, and in suspensions of interacting RTPs for arbitrary $U\tau/L$ in §2.2.2. While the limiting form of the covariance in the rapid tumbling limit ($U\tau/L\ll 1$) is readily inferred from the general expression, and is stated in (31), the form of the covariance in the straight-swimmer limit is more subtle. Owing to a separation of scales between the swimmer size and the run length, one may, however, use the method of matched asymptotic expansions to obtain the covariance in the limit $U\tau/L\gg 1$. In a reference frame moving with one of the swimmers (the test swimmer), the physical domain may be divide into two regions as follows:

An inner region $r\sim$ O(1) $\ll U\tau/L$, where the pair of swimmers move relative to each other predominantly in straight trajectories. In this region, one may approximate the correlations induced by the swimmer disturbance velocity fields as those along straight-swimmer trajectories. The expressions pertaining to this region have already been obtained as part of the straight-swimmer analysis given in §2.2.1. 2. 2.

An outer region $r\gg U\tau/L$, where orientation decorrelation due to random tumbles need to be taken into account. However, given the large separations, the pair of swimmers may be treated as point force-dipoles. Neglect of the finite swimmer size allows for a considerable simplification of the expressions involved.

A schematic highlighting the above demarcation of the physical domain, for purposes of the matched asymptotic analysis, is shown in figure 8. One expects that in the inner region the covariance has a non-decaying character on account of the strong correlations between interacting straight-swimmers. In the outer region, there is a transition to the eventual O($1/r$) decay for $r\gg U\tau/L$, characteristic of interacting RTP’s. The covariance plotted in figure 5 does conform to these scalings. However, in the interval $1\ll r\ll U\tau/L$, corresponding to the region of overlap, there is only a weak decay with separation. The analysis presented in this section enables us to determine the functional form of the covariance in this intermediate region, where the swimmers may be approximated as straight-swimming point force-dipoles.

Now, the point force-dipole, used for modeling the swimmers in the intermediate and outer regions above, has velocity disturbance and pressure fields that satisfy the following equations: {subeqnarray} -\bnablaP_i + η∇^2 u_i(x) &= D p_ipi**:\bnabla**xδ(x-x_i),

\bnabla**⋅u**_i = 0.

where $D=1/(4\ln\kappa)$ is the non-dimensional dipole strength and $\boldsymbol{p}_{i}$ denotes the orientation of the $i^{th}$ dipole swimmer. The Fourier transformed velocity covariance for interacting point force-dipoles is then expressible as:

[TABLE]

where $\hat{\Omega}_{2\,dp}^{(1)}$ denotes the Fourier transformed pair-probability density, and the subscript $dp$ denotes dipoles. The use of the Fourier transformed point-dipole fields obtained from (8), along with the procedure outlined in §2.2.2, leads to the following simplified expression for $\hat{\Omega}_{2\,dp}^{(1)}$ in the outer region:

[TABLE]

Based on the aforementioned arguments, and with the intent of arriving at an asymptotic expression for large $U\tau/L$, one may partition the original covariance integral in (10) into inner and outer contributions in the following manner:

[TABLE]

where we have introduced a parameter $\epsilon\in(0,1)$ (Subramanian & Brady, 2006). In (36), the subscripts $SS$ and $dpT$ denotes slender straight-swimmers and tumbling force-dipole swimmers; here, $\zeta_{SS}(\boldsymbol{k},\boldsymbol{p}_{1},\boldsymbol{p}_{2})$ and $\zeta_{dpT}(\boldsymbol{k}_{1},\boldsymbol{p}_{1},\boldsymbol{p}_{2})$, respectively, represent the factors multiplying $\boldsymbol{\delta}(\boldsymbol{k}+\boldsymbol{k}^{\prime})$ in (10) and (34). Note that, in terms of the wavenumber $k$, the inner and outer regions above correspond to $kL$ and $kU\tau$, respectively, being O(1), and this differing non-dimensionalzation is accounted for in (36). Thus, in the second integral, the scalings are the same as before, with $r$ being the radial distance scaled by $L$; on the other hand, in the first integral $r_{1}=rL/(U\tau)$.

On naively applying the limit $U\tau/L\to\infty$ in (36), one encounters divergent integrals. This is because $\zeta_{dpT}$ diverges for $k_{1}\to\infty$ (point dipoles exhibit a divergence as $r_{1}\to 0$), while $\zeta_{SS}$ diverges for $k\to 0$ (the logarithmic divergence of straight-swimmers as $r\to\infty$). Subtracting the straight-swimming point force dipole contribution (corresponding to the contribution of the overlap region described above) from the respective integrands can, however, account for the said divergences. Therefore, one can recast (36) as:

[TABLE]

such that the divergent contributions are now contained in the final two terms on the right hand side of (37). Here, the subscript $dpS$ denotes straight-swimming point force-dipoles. Now, the added integral terms are divergent both as $k\to 0$ and as $k\to\infty$, and the Heaviside functions $H(x)$ included in (37) ensure that these added terms solely account for the divergences already present (in the limit $U\tau/L\to\infty$) in (37). Note that both $\zeta_{dpS}(\boldsymbol{k},\boldsymbol{p}_{1},\boldsymbol{p}_{2})$ and $\zeta_{dpT}(\boldsymbol{k},\boldsymbol{p}_{1},\boldsymbol{p}_{2})$ have a common form as $U\tau/L\to\infty$. To see this one may use the Plemlj-Sokhotski formula on the pair-probability density appearing in $\zeta_{dpT}$, as explained in §2.2.2 for slender swimmers. Alternatively, one could also determine the pair-probability in $\zeta_{dpS}$ following the method outlined in §2.2.1, which gives:

[TABLE]

The ‘difference integrals’ in (37) that involve the Heaviside functions are now convergent, and one may therefore set $(U\tau/L)^{\epsilon}=\infty$ and $(U\tau/L)^{(\epsilon-1)}=0$ in (37) to get:

[TABLE]

We designate the first two terms on the right side as $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{corr}^{dpT-dpS}$, the next two as $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{corr}^{SS-dpS}$ and the last two $U\tau/L$-dependent terms as $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{corr}^{Div}$, and evaluate them separately. Note that terms containing $\epsilon$ are now restricted to $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{corr}^{Div}$. To obtain a leading order $\epsilon$-independent form from (39), we first evaluate $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{corr}^{Div}$. This will also yield natural choices of $L_{o}$ and $U_{p}$ that are most convenient for the analysis.

After some algebra the $U\tau/L$-dependent terms on the right-hand side of (39) are expressible as:

[TABLE]

where $Ci$ represents the Cosine integral (Gradshteyn & Ryzhik, 1980). We now use the additive method of constructing a uniformly valid approximation of (40) at leading order in $\epsilon$ (Van Dyke, 1964; Subramanian & Brady, 2006). This requires approximate forms of (40) valid, respectively, in the inner, outer and matching regions.

In the inner region, where $r_{1}\ll 1$ and $r\sim$ O(1) one gets:

[TABLE]

where we have used the approximation $Ci(x)\approx j_{0}(x)$, for large $x$. In the outer region, one has $r_{1}\sim$ O(1) and $r\gg 1$, and at leading order (40) simplifies to:

[TABLE]

In the matching region, with $r_{1}\ll 1$ and $r\gg 1$, one gets

[TABLE]

Here, $\gamma_{E}$ is the Euler-Mascheroni constant. From (41)-(43), we find that $L_{o}=U_{p}=1$ is a convenient choice for the analysis, using which, the leading order $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{corr}^{Div}$ is expressible as:

[TABLE]

With $L_{o}$ and $U_{p}$ fixed, we now evaluate the remaining terms in (39). Following standard procedure presented in §2.2.1-§2.2.2, and additionally using (34), (35) and (38), $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{corr}^{dpT-dpS}$ can be simplified as:

[TABLE]

Similarly, using (20) for slender straight-swimmers, and (34) and (38) for straight-swimming point force-dipole swimmers, $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{corr}^{SS-dpS}$ can be written as:

[TABLE]

Now, using (44)-(46) in (39), we finally get the asymptotic form of $\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle$ in the limit $U\tau/L\gg 1$ to be:

[TABLE]

It is worth reiterating that while (LABEL:covarRTP) is the exact expression for the velocity covariance valid for arbitrary $U\tau/L$, (47) is valid only for large $U\tau/L$. In figure 9, we compare the correlated contributions to the covariance calculated from (LABEL:covarRTP) and (47). The matched asymptotic expansion approach yields excellent agreement with the exact variance, derived in §2.2.2 for the three largest values of $U\tau/L$, and remains a good approximation down to $U\tau/L=7.5$. Up until $U\tau/L\sim$ O(1), the covariance smoothly transitions to an O($1/r$) far-field scaling following the variance plateau. However, for $U\tau/L\gg 1$, the covariance exhibits a weak intermediate logarithmic decay for $1\ll r\ll U\tau/L$, following which, it again decays as O($1/r$) in the far-field. In this regime, $\langle\boldsymbol{u}(\boldsymbol{x})\cdot\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle|_{corr}^{Div}$ is given by (43), and with $r_{1}\equiv r*L/(U\tau)$, the weak decay arises from the O($\ln[rL/(U\tau)]$) scaling of the O$(nL^{3})^{2}$ term in the covariance. The aforementioned weak intermediate logarithmic decay of the covariance is compared to the exact solution, for $U\tau/L=750$, in figure 10.

Therefore, for $U\tau/L\gg 1$, the behavior of the covariance may be summarized as follows:

For $r\ll 1$, it exhibits a finite variance plateau, with the variance itself scaling as $\ln U\tau/L$. 2. 2.

For $1\ll r\ll U\tau/L$, it exhibits a weak decay of O($\ln[rL/(U\tau)]$). 3. 3.

For $r\gg U\tau/L$, corresponding to the far-field, it decays as O($1/r$).

3 Tracer diffusivity

The velocity fluctuations analyzed in §2 convect a passive tracer immersed in the suspension. The dynamics of the passive-tracer is governed by,

[TABLE]

where $\boldsymbol{u}(\boldsymbol{x})$ is the suspension velocity field defined in (2). During a tracer-slender swimmer interaction, the convection velocity of the tracer is O$(U/\ln\kappa)$ while the relative tracer-swimmer velocity is O$(U)$. Thus, the tracer displacement during an interaction event is negligible and the tracer statistics may be calculated within the Eulerian approximation (Kasyap et al., 2014). The expression for the mean-squared displacement (MSD) of the tracer, at time $t$, is then given in terms of the time-dependent fluid velocity correlation function as,

[TABLE]

where the angular-brackets again denote an ensemble average (Balakrishnan, 2008; Zwanzig, 2001). In terms of the Fourier-transformed variables, the tracer mean-squared displacement is finally given as,

[TABLE]

Following the discussion in §2, we analyze the tracer MSD till O($(nL^{3})^{2}$); thus it may be written as,

[TABLE]

where $\langle x^{2}(t)\rangle|_{uncorr}\sim$ O($nL^{3}$) is the tracer MSD in a non-interacting suspension, while $\langle x^{2}(t)\rangle|_{corr}\sim$ O($(nL^{3})^{2}$) denotes the correlated contribution.

3.1 Suspensions of non-interacting swimmers

For non-interacting swimmers, the Fourier-transformed velocity correlation function is given by

[TABLE]

where $\Omega_{1}=n/4\pi$ is the steady-state singlet probability density and $\hat{G}_{1}(\boldsymbol{p}_{1},t_{1}|\boldsymbol{p}_{1}^{\prime},t_{2};\boldsymbol{k})$ is the Fourier transformed transition probability of finding a swimmer with orientation $\boldsymbol{p}_{1}$ at time $t_{1}$ starting with a swimmer with orientation $\boldsymbol{p}_{1}^{\prime}$ at $t_{2}$. The governing equation for $\hat{G}_{1}$ is then the same as the equation governing the relaxation of the singlet probability density and is given by,

[TABLE]

The orientation-space eigenmodes of the operator governing the transition probability density, $\hat{G}_{1}$ in (53), that contribute to the velocity-correlation function have no concentration fluctuations (see the discussion in §2.2.2). Thus, for these modes, the inverse-tumbling term in (53) vanishes and the expression for $\hat{G}_{1}$ is then easily given as,

[TABLE]

We emphasize that this expression for $\hat{G}_{1}$ does not capture the relaxation of the concentration modes. The complete expression for $\hat{G}_{1}$ is needed to characterize the time-dependent evolution of a initially localized population of run-and-tumble swimmers and will be reported separately (Garg et al., 2020). Using the above result, and upon simplifying, we obtain the expression for the MSD of a passive tracer as,

[TABLE]

At times much shorter than the velocity correlation time, the tracer shows ballistic motion and the MSD simplifies to,

[TABLE]

The magnitude of the MSD in the ballistic regime is hence set by the fluid velocity variance discussed in section §2. On the other hand at times much longer than the velocity correlation time, the tracer shows diffusive motion with $\langle x^{2}(t)\rangle|_{uncorr}=D_{t}|_{uncorr}t$, where the effective diffusivity is given by,

[TABLE]

The diffusivity may also be derived directly using the Green-Kubo formula (Balakrishnan, 2008) and has also been obtained by Kasyap et al. (2014) using a different approach. On the other hand, to the best of our knowledge, the complete expression for the MSD in (LABEL:eq:msd_sing) has not been given in the literature even for non-interacting swimmers. As discussed earlier in §2.1, the integrals in (LABEL:eq:msd_sing) may be numerically evaluated with the choice of a $\boldsymbol{k}$-aligned spherical coordinate system.

3.2 Suspensions of interacting swimmers

For interacting swimmers, the Fourier-transformed velocity correlation function at O$((nL^{3})^{2})$ is given by

[TABLE]

where $\hat{\Omega}_{2}^{(1)}$ is the correlated pair probability density derived in §2.2.2 and $\hat{G}_{2}$ gives the Fourier transformed transition probability of finding a pair of swimmers with orientations $\boldsymbol{p}_{1}$ and $\boldsymbol{p}_{2}$ at time $t_{1}$ starting with a pair of swimmers with orientations $\boldsymbol{p}_{1}^{\prime}$ and $\boldsymbol{p}_{2}^{\prime}$ at time $t_{2}$. The governing equation for $\hat{G}_{2}$ is then same as the equation governing the relaxation of the pair probability density and is given by,

[TABLE]

The inverse tumbling terms in (59) again don’t contribute (see §2.2.2), and the expression for $\hat{G}_{2}$ relevant to evaluating the velocity correlation function is given by,

[TABLE]

Using this result and simplifying, we obtain the final expression for the MSD of a passive tracer at O$((nL^{3})^{2})$ as

[TABLE]

At times much shorter than the velocity correlation time, the tracer shows ballistic motion and the MSD simplifies to,

[TABLE]

The magnitude of the MSD in the ballistic regime is again set by the fluid velocity variance, with the relevant $U\tau/L$ scalings discussed in section §2. On the other hand, for times much longer than the velocity correlation time, the tracer shows diffusive motion with $\langle x^{2}(t)\rangle|_{corr}=D_{t}|_{corr}t$ for any finite swimmer run-length ($U\tau/L$). The effective diffusivity is given by,

[TABLE]

The integrals in (LABEL:eq:msd_pair) and (63) are again evaluated numerically (see §2.1 for the coordinate system used). The total mean-squared displacement of the tracer, to O$(nL^{3})^{2}$, is obtained by combining the expressions in (LABEL:eq:msd_sing) and (LABEL:eq:msd_pair), and is plotted in figure 11, as a function of $nL^{3}$, for a suspension of pushers.

Surprisingly, figure 11 shows that the time taken to transition from the ballistic to the diffusive regime increases with increasing volume-fraction of the swimmers. This surprising behavior can be explained by noting the differing time scales for the decay of the velocity correlations at O$(nL^{3})$ and O$(nL^{3})^{2}$. At O$(nL^{3})$, this time-scale ($t_{c}$) is set by the swimmer-tracer interaction time. For rapid tumblers ($U\tau/L\ll 1$), the interaction is cut off by the decorrelation of the swimmer orientation, implying $t_{c}\sim$ O($\tau$). On the other hand for straight swimmers ($U\tau/L\gg 1$), the distance that the tracer is convected asymptotes to a finite value for tracer-swimmer separations greater than O$(L)$, and thus $t_{c}\sim$ O($L/U$). At O$(nL^{3})^{2}$, the decay of the velocity correlations, irrespective of $U\tau/L$, only occurs due to orientation decorrelation of the swimmers on the time scale $\tau$. As a result for $nL^{3}$ of order unity, the transition time between the ballistic and diffusive regimes diverges as $U\tau/L$. The resulting broad cross-over gives the impression of an $nL^{3}$-dependent anomalous exponent in the interval $L/U\ll t\ll\tau$ (see the curve for $nL^{3}=2.5$ in figure 11). For rapid tumblers $(U\tau/L\ll 1)$, the transition time is O$(\tau)$ regardless of $nL^{3}$, and there is no intermediate anomalous scaling.

The inset of figure 11 shows the O$(nL^{3})^{2}$ correlated contribution to the tracer diffusivity ($D_{t}|_{corr}$). In the straight-swimmer limit, $D_{t}|_{corr}$ diverges linearly in $U\tau/L$. This scaling may be rationalized by starting from $D_{t}\sim U_{t}^{2}t_{c}$ where $U_{t}$ is the scale of the velocity convecting the tracer and $t_{c}$ is the time scale of decay of the velocity correlations. $U_{t}$ is O$(U/\ln\kappa)$ regardless of the swimmer volume-fraction $(nL^{3})$, while the $t_{c}$ scaling is discussed above. The resulting scaling for tracer diffusivity in an interacting suspension, to O$(nL^{3})^{2}$, is therefore given by,

[TABLE]

and

[TABLE]

and confirms the scalings shown in figure 11 (inset) at O$(nL^{3})^{2}$. The tracer diffusivity variation with volume fraction also shows a pusher-puller bifurcation similar to the velocity variance as shown in figure 12.

Several experiments probing bacterial suspension dynamics with passive tracers have reported an intermediate super-diffusive regime and a volume fraction dependent cross-over time (Wu & Libchaber, 2000; Kim & Breuer, 2004; Patteson et al., 2016; Valeriani et al., 2011; Peng et al., 2016; Chen et al., 2007a; Argun et al., 2016). To the best of our knowledge, this is the first theoretical explanation for these observations. Our analysis thus shows that in a bath of interacting persistent swimmers ($U\tau/L\gg 1$), the temporal and spatial correlations of the (effective) noise experienced by the passive tracer increase to O$(\tau)$ and O$(U\tau)$ respectively. The simplistic approximation for the active noise often found in the literature is thus seen to be insufficient (Maggi et al., 2014; Argun et al., 2016; Kanazawa et al., 2019). Further, when confined by a harmonic potential with a trap radius of O$(U\tau)$, the probability distribution of the tracer displacements is thus expected to become pronouncedly non-Gaussian with increasing $nL^{3}$. This has indeed been seen in experiments (Krishnamurthy et al., 2016; Argun et al., 2016).

4 Conclusion

In this paper, we have shown that long-ranged pair-correlations lead to divergent velocity fluctuations and tracer diffusivities in straight-swimmer suspensions. For RTPs, this divergence is replaced by a logarithmic increase of the velocity variance, and a linear increase of the tracer diffusivity, with the run length for $U\tau/L\gg 1$. Correlated fluctuations are also crucially dependent on the swimming mechanism (pusher vis-a-vis puller), being larger for pusher suspensions. Thus, in contrast to recent work emphasizing the role of anisotropic tracers (Peng et al., 2016; Yang et al., 2016), spherical tracers already discriminate between pusher and puller suspensions due to correlated fluctuations.

While our detailed analysis is limited to slender swimmers, the scaling arguments presented only involve far-field hydrodynamics. It is thus expected that they would be applicable to suspensions of general anisotropic swimmers (Toppaladoddi & Balmforth, 2014; Kyoya et al., 2015) provided far-field interactions dominate over near-field ones. This assertion is supported by the evidence of logarithmic scaling for the velocity variance in simulations of regularized point dipoles by Underhill & Graham (2011). The spherical squirmer model often used in studies of active suspensions may therefore represent an exceptional scenario and the results obtained need to be extrapolated with care (Ishikawa et al., 2006; Ishikawa & Pedley, 2007, 2008; Mehandia & Nott, 2008; Alarcón & Pagonabarraga, 2013; Delmotte et al., 2018).

In this paper, we have examined correlations due to direct pairwise hydrodynamic interactions between swimmers, while not accounting for those induced (indirectly) by long-wavelength fluctuations in the suspension velocity and orientation fields. It is well-known that these fluctuations lead to an instability of a pusher suspension, in turn implying diverging contributions of the ‘collective effects’ to statistical correlations at the threshold (Subramanian & Koch, 2009; Stenhammar et al., 2017; Qian et al., 2017; Bárdfalvy et al., 2019). This divergence is distinct from the one described here. Nevertheless, the framework needed for this calculation requires the $U\tau/L$-dependent pair-probability, determined here, as an input, and will be reported in future-work. Our work thus lays a rigorous foundation for understanding hydrodynamic interactions between anisotropic swimmers, and their role in bacterial turbulence (Wensink et al., 2012).

Appendix A Rotation rates of slender swimmers

Herein, we derive the expression for the rates of rotation $\dot{\boldsymbol{p}}_{ij}$ for a pair of interacting slender swimmers. One can define the rate of rotation of swimmer $1$ due to swimmer $2$ as $\dot{\boldsymbol{p}}_{12}=\boldsymbol{\omega}_{1}\wedge\boldsymbol{p}_{1}$, with $\boldsymbol{\omega}_{1}$ representing its angular velocity. The velocity at any point $s$ on a rigid slender body may be written as $\boldsymbol{u}_{1}^{s}(s\boldsymbol{p}_{1})=\boldsymbol{u}_{1}^{T}+\boldsymbol{\omega}_{1}\wedge\boldsymbol{p}_{1}s$, where $\boldsymbol{u}_{1}^{T}$ is the translational velocity of the geometric center ($s=0$). From viscous slender body theory (Batchelor, 1970; Kim & Karrila, 1991), one has the relation: $\boldsymbol{u}_{1}^{s}(s\boldsymbol{p}_{1})-\boldsymbol{u}^{\infty}(s\boldsymbol{p}_{1})=(\ln\kappa/\eta)f(s)\boldsymbol{p}_{1}\boldsymbol{\cdot}(\mathsfbi{I}+\boldsymbol{p}_{1}\boldsymbol{p}_{1})$ to leading logarithmic order. Here, $\boldsymbol{u}^{\infty}$ is the ambient field which in the present problem is the velocity field generated by the swimmer $2$, that is, $\boldsymbol{u}^{\infty}\equiv\boldsymbol{u}_{2}$; the force density $f(s)=\eta Usgn(s)/\ln\kappa$. To determine, $\boldsymbol{u}_{2}$, one solves the Stokes equations (2) at any location ($-\boldsymbol{r}+s\boldsymbol{p}_{1}$) along swimmer 1; that is:

[TABLE]

Again, it is convenient to solve (66) in Fourier space, and the Fourier transformed ambient field along the length of swimmer $1$, $\hat{u}_{2}(\boldsymbol{k})$, can be written as:

[TABLE]

on applying the non-dimensionalizations stated below (6). The force-free condition ensures that the swimmer propagates with the sum of the swimming and disturbance velocities averaged over its length. However, given that the latter is O($1/\ln\kappa$) smaller, it may be neglected. Applying the torque-free condition, $\int_{-1}^{1}s\boldsymbol{p}_{1}\wedge[\boldsymbol{u}_{1}^{s}(s\boldsymbol{p}_{1})-\boldsymbol{u}_{2}(s\boldsymbol{p}_{1})]\mathrm{d}s=0$, results in the following expression for $\hat{\dot{\boldsymbol{p}}}_{12}$:

[TABLE]

where $j_{1}(x)=\sin(z)/z^{2}-\cos(z)/z$ is the spherical Bessel’s function of the first kind (Gradshteyn & Ryzhik, 1980). The Fourier transformed rate of rotation of the second swimmer $\hat{\dot{\boldsymbol{p}}}_{21}$ is given by interchanging the indices $1$ and $2$ in (68).

Appendix B Orientation Correlations

Herein, we interpret the orientation correlations that develop in a suspension of pairwise interacting RTPs, by characterizing the nematic and polar correlations. In a suspension of interacting swimmers, the nematic correlation, at the leading order, between two swimmers at a distance $r$ is given by,

[TABLE]

whereas, the polar correlation, at the leading order, is expressible as,

[TABLE]

where we have used the normalisation condition

[TABLE]

Using the result for $\hat{\Omega}^{(1)}_{2}$ from (LABEL:eq:omega21RTP5), (69) and (70) simplify to,

[TABLE]

and

[TABLE]

respectively. These expressions are numerically evaluated using the coordinate system described in detail in §2.2. Figures 13 and 14 show the spatial decay of the nematic and polar correlations, respectively, for a range of $U\tau/L$ in a suspension of pushers. The orientation correlations are qualitatively similar to those obtained in previous numerical computations (Saintillan & Shelley, 2007; Bárdfalvy et al., 2019).

For small $U\tau/L$, the orientational correlations are constant for a swimmer length ($r\sim 1$). However, for large distances ($r\gg 1$), each swimmer sees the other swimmer as a rapidly tumbling dipole and hence the correlations decay to zero faster than any power law. On the other hand, in the limit of straight-swimmers ($U\tau/L\gg 1$), there is an intermediate region ($1\ll r\ll U\tau/L$) where the correlations decay as a power law ($\mathcal{P},\mathcal{N}\sim 1/r^{2}$). For distances $r\gg U\tau/L$, the correlations are again screened due to tumbling. Further, we see that the orientational correlations, for any $r$, asymptote to a constant value for straight swimmers ($U\tau/L\gg 1$) unlike the velocity variance. For rapid tumblers, on the other hand, as expected the nematic correlation, $\mathcal{N}(r)/(U\tau/L)\sim$ O$(1)$, dominates over the polar correlation, $\mathcal{P}(r)/(U\tau/L)\sim$ O$((U\tau/L)^{2})$.