3.1. Suspensions of non-interacting swimmers

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

(3.7)\begin{equation} \langle\hat{\boldsymbol{u}}(\boldsymbol{k})\boldsymbol{\cdot} \hat{\boldsymbol{u}}(\boldsymbol{k}^\prime)\rangle\vert_{uncorr} = \frac{nL^3}{16{\rm \pi}^7 (\ln\kappa)^2 k^4}\boldsymbol{\delta}(\boldsymbol{k}+\boldsymbol{k}^\prime)\int \mathrm{d}\boldsymbol{p}_1 \left[1 - \frac{\left(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)^2}{k^2}\right]\frac{\sin^4\left(\dfrac{{\rm \pi} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1}{2}\right)}{(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1)^2}. \end{equation}

An inverse Fourier transform of (3.7) yields

(3.8)\begin{align} \langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle\vert_{uncorr} = \frac{nL^3}{4{\rm \pi}^6 (\ln\kappa)^2 r}\int_0^\infty \mathrm{d}k \frac{\sin(2{\rm \pi} kr)}{k^5}\int_{-1}^1\mathrm{d}\mu\left(\frac{1-\mu^2}{\mu^2}\right)\sin^4\left(\frac{{\rm \pi} k \mu}{2}\right). \end{align}

Owing to isotropy, the covariance only depends on the scalar distance $r$. The variance is obtained on setting $r = 0$ in (3.8), which gives $\langle \boldsymbol {u}(\boldsymbol {x})\boldsymbol {\cdot }\boldsymbol {u}(\boldsymbol {x})\rangle \vert _{uncorr} = nL^3/[96{\rm \pi} (\ln \kappa )^2]$. One can simplify (3.8) in the far-field ($r\gg 1$) by adopting $\tilde {k} = k r$ as the new variable (in the outer integral) instead of $k$, and approximating $\sin [{\rm \pi} \tilde {k} \mu /(2r)] \approx {\rm \pi}\tilde {k} \mu /(2r)$ (that is, seeking a limiting form of the expression for $\tilde {k}\sim O(1)$, but $r\gg 1$). This yields the far-field limiting form of the covariance which, expectedly, is similar to that obtained earlier for point dipoles (Underhill & Graham Reference Underhill and Graham2011; Zaid et al. Reference Zaid, Dunkel and Yeomans2011), and is given by $\langle \boldsymbol {u}(\boldsymbol {x})\boldsymbol {\cdot }\boldsymbol {u}(\boldsymbol {x}+ \boldsymbol {r})\rangle \vert _{uncorr} = nL^3/[480{\rm \pi} (\ln \kappa )^2 r]$. Recognizing that the (non-dimensional) dipole strength for slender swimmers is $D = 1/(4\ln \kappa )$, one in fact recovers the exact covariance for point force-dipoles, $\langle \boldsymbol {u}(\boldsymbol {x})\boldsymbol {\cdot }\boldsymbol {u} (\boldsymbol {x}+\boldsymbol {r})\rangle \vert _{uncorr} = nL^3 D^2/[30{\rm \pi} r]$. The aforementioned variance and far-field covariance have been plotted alongside (3.8), evaluated numerically, in figure 2. Note that the covariance, as given by (3.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).

Figure 2.

The ${O}(nL^3)$ velocity covariance in a dilute non-interacting suspension of slender swimmers. The covariance transitions from an ${O}(nL^3)/(\ln \kappa )^2$ variance plateau to a far-field ${O}(1/r)$ asymptote. The swimmer aspect ratio is chosen to be $\kappa = 8$.

3.2. Suspensions of interacting straight-swimmers

Herein, we consider the covariance in a suspension of hydrodynamically interacting slender straight-swimmers to ${O}(nL^3)^2$. The contribution at this order involves the correlations that develop due to pairwise interactions between slender straight-swimmers, and these are characterized by the Fourier-transformed pair-probability density $\hat {\varOmega }_2(\boldsymbol {p}_1, \boldsymbol {p}_2; \boldsymbol {k})$ in (3.6). The correlated component of the Fourier transformed velocity variance is given by

(3.9)\begin{align} \langle\hat{\boldsymbol{u}}(\boldsymbol{k})\boldsymbol{\cdot} \hat{\boldsymbol{u}}(\boldsymbol{k}^{\prime})\rangle\vert_{corr} &= \frac{\delta(\boldsymbol{k}+\boldsymbol{k}^\prime)}{4{\rm \pi}^6(\ln\kappa)^2 k^4}\left[ \int\mathrm{d}\boldsymbol{p}_1\,\mathrm{d}\boldsymbol{p}_2\left[ \left( {\boldsymbol{\mathsf{I}}} - \frac{\boldsymbol{k}\boldsymbol{k}}{k^2}\right)\right]\boldsymbol{:}\boldsymbol{p}_1\boldsymbol{p}_2 \frac{\sin^2\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)}{ \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1}\right.\nonumber\\ & \quad\times\!\left.\frac{\sin^2\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_2\right)}{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2}\, \hat{\varOmega}_2(\boldsymbol{p}_1, \boldsymbol{p}_2; \boldsymbol{k})\right], \end{align}

where $\hat {\varOmega }_2$ is to be determined.

As already discussed in § 2, the disturbance field of a slender swimmer is only $O(\ln \kappa )^{-1}$ on length scales of $O(L)$, and this leads to logarithmically weak interactions between swimmers on scales of order the swimmer length. The slender swimmer assumption (Berg Reference Berg2004; Gachelin et al. Reference Gachelin, Rousselet, Lindner and Clement2014) has, in fact, reasonably reproduced experimentally measured disturbance velocity fields to distances of ${O}(L)$ (Kasyap et al. Reference Kasyap, Koch and Wu2014). This allows one to expand $\varOmega _2$ as a series in $(\ln \kappa )^{-1}$ about the base uncorrelated state. One writes $\varOmega _2 = \varOmega _2^{(0)} + 1/(\ln \kappa ) \varOmega _2^{(1)} + \ldots\,$, where $\varOmega _2^{(0)}$ is the product of two singlet probability densities, $(nL^3)^2/(4{\rm \pi} )^2$, and denotes the uncorrelated base state at leading order; $\varOmega _2^{(1)} \equiv \varOmega _2^{(1)}(\boldsymbol {r}, \boldsymbol {p}_1, \boldsymbol {p}_2)$ gives the correlated probability, at ${O}(\ln \kappa )^{-1}$, of finding two swimmers with orientations $\boldsymbol {p}_1$ and $\boldsymbol {p}_2$, separated by $\boldsymbol {r}$. Note that the contribution due to $\varOmega _2^{(0)}$ in (3.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 $\hat {\varOmega }_2^{(1)}(\boldsymbol {p}_1, \boldsymbol {p}_2; \boldsymbol {k})$, the Fourier transform of the $O(\ln \kappa )^{-1}$ pair-probability density above, in order to characterize pair-correlations in the swimmer suspension and, thence, determine the correlated contribution to the covariance in (3.9). The leading-order contribution being zero, (3.9) is now expressible as

(3.10)\begin{align} \langle\hat{\boldsymbol{u}}(\boldsymbol{k})\boldsymbol{\cdot}\hat{\boldsymbol{u}}(\boldsymbol{k}^\prime) \rangle\vert_{corr} &= \frac{\boldsymbol{\delta}(\boldsymbol{k}+\boldsymbol{k}^\prime)}{4{\rm \pi}^6(\ln\kappa)^3 k^4} \left[\int\mathrm{d}\boldsymbol{p}_1\,\mathrm{d}\boldsymbol{p}_2 \frac{1}{\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_1} \frac{1}{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2} \sin^2 \left(\frac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right) \sin^2\left(\frac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2\right)\right. \nonumber\\ &\quad \times\!\left.\left({\boldsymbol{\mathsf{I}}} - \frac{\boldsymbol{kk}}{k^2}\right)\boldsymbol{:}\boldsymbol{p}_1\boldsymbol{p}_2 \hat{\varOmega}^{(1)}_2(\boldsymbol{p}_1, \boldsymbol{p}_2; \boldsymbol{k}) \right ]. \end{align}

Now, in a macroscopic suspension of $N$ ($\gg 1$) interacting slender straight-swimmers, $\varOmega _N$ satisfies the Liouville equation, that 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, one obtains the equation for the pair-probability density $\varOmega _2$, at steady-state, as

(3.11)\begin{equation} \boldsymbol{\nabla}_{\boldsymbol{r}}\boldsymbol{\cdot}\{[(U\boldsymbol{p}_2 + \bar{\boldsymbol{u}}_{21}) - (U\boldsymbol{p}_1+\bar{\boldsymbol{u}}_{12})]\varOmega_2\} + \boldsymbol{\nabla}_{\boldsymbol{p}_1}\boldsymbol{\cdot}(\dot{\boldsymbol{p}}_{12} \varOmega_2) + \boldsymbol{\nabla}_{\boldsymbol{p}_2}\boldsymbol{\cdot}(\dot{\boldsymbol{p}}_{21} \varOmega_2) = 0. \end{equation}

The terms within braces in (3.11) denote the convection of $\varOmega _2$ by the relative velocity of the swimmer pair that includes contributions due to swimming ($U\boldsymbol {p}_i$) and due to the disturbance velocity fields $\bar {\boldsymbol {u}}_{ij}$ (of swimmer $j$ at the location of swimmer $i$); the overbar here denotes an average over the length of swimmer $i$, a result derivable from viscous slender-body theory (Batchelor Reference Batchelor1970). These disturbance fields are $O(\ln \kappa )^{-1}$ since we are considering separations of $O(L)$ and larger, and are thus asymptotically small compared with the swimming-induced convection. The third and fourth terms denote rotation of the swimmer orientations due to the disturbance velocity fields, with $\dot {\boldsymbol {p}}_{ij}$ again denoting the rotation of swimmer $i$ by the disturbance velocity field due to swimmer $j$; the expression for the $\dot {\boldsymbol {p}}_{ij}$, in Fourier space, are given in appendix A.

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

(3.12)\begin{equation} (\boldsymbol{p}_2-\boldsymbol{p}_1)\boldsymbol{\cdot}\nabla_{\boldsymbol{r}}\varOmega_2^{(1)} = - \frac{(nL^3)^2}{(4{\rm \pi})^2} \left[\boldsymbol{\nabla}_{\boldsymbol{p}_{1}}\boldsymbol{\cdot}\dot{\boldsymbol{p}}_{12} + \boldsymbol{\nabla}_{\boldsymbol{p}_2}\boldsymbol{\cdot}\dot{\boldsymbol{p}}_{21}\right ], \end{equation}

in non-dimensional terms. Fourier transforming (3.12) leads to

(3.13)\begin{equation} 2{\rm \pi} \textrm{i}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_2-\boldsymbol{p}_1)\hat{\varOmega}_2^{(1)} = - \frac{(nL^3)^2}{(4{\rm \pi})^2} \left[\boldsymbol{\nabla}_{\boldsymbol{p}_1}\boldsymbol{\cdot}\hat{\dot{\boldsymbol{p}}}_{12} + \boldsymbol{\nabla}_{\boldsymbol{p}_2}\boldsymbol{\cdot}\hat{\dot{\boldsymbol{p}}}_{21}\right], \end{equation}

where the convection in physical space reduces to a multiplicative $\boldsymbol {k}$ factor in Fourier space. The general solution of (3.13) may be written as

(3.14)\begin{equation} \hat{\varOmega}_{2}^{(1)} = C\delta\left[\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_2-\boldsymbol{p}_1)\right] - \frac{(nL^3)^2}{(4{\rm \pi})^2 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_2-\boldsymbol{p}_1)}\left[ \boldsymbol{\nabla}_{\boldsymbol{p}_1}\boldsymbol{\cdot}\hat{\dot{\boldsymbol{p}}}_{12} + \boldsymbol{\nabla}_{\boldsymbol{p}_2}\boldsymbol{\cdot}\hat{\dot{\boldsymbol{p}}}_{21}\right], \end{equation}

where the first term on the right-hand side, involving the Dirac delta function, is the homogeneous solution (using the identity $z\delta (z) = 0$; see Gelfand & Shilov (Reference Gelfand and Shilov1964)), and the expressions for $\hat {\dot {\boldsymbol {p}}}_{12}$ and $\hat {\dot {\boldsymbol {p}}}_{21}$ that appear in the second term, 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, $\varOmega _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 the aforementioned upstream infinity. Therefore, the constraint of uncorrelated swimmers infinitely far upstream may be given as

(3.15)\begin{equation} \lim_{z\to-\infty}\varOmega_2^{(1)} \equiv\lim_{z\to-\infty}\int \mathrm{d}\boldsymbol{k}\exp\left[2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{x}\right]\hat{\varOmega}_2^{(1)} = 0. \end{equation}

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

(3.16)\begin{align} &\lim_{z\to-\infty}\int \mathrm{d}k_z\exp\left[2{\rm \pi} \textrm{i} k_z z\right]\nonumber\\ &\quad \times\!\left[C \delta[k_z |\boldsymbol{p}_2-\boldsymbol{p}_1|] - \frac{(nL^3)^2}{(4{\rm \pi})^2 2{\rm \pi} \textrm{i} k_z|\boldsymbol{p}_2-\boldsymbol{p}_1|}\left[\boldsymbol{\nabla}_{\boldsymbol{p}_1}\boldsymbol{\cdot} \hat{\dot{\boldsymbol{p}}}_{12} + \boldsymbol{\nabla}_{\boldsymbol{p}_2}\boldsymbol{\cdot}\hat{\dot{ \boldsymbol{p}}}_{21} \right]\right] = 0. \end{align}

While the first term within brackets in (3.16) 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$:

(3.17)\begin{equation} C = -\frac{(nL^3)^2}{32{\rm \pi}^2} \lim_{k_z\to 0}\left[\boldsymbol{\nabla}_{\boldsymbol{p}_1}\boldsymbol{\cdot}\hat{\dot{\boldsymbol{p}}}_{12} + \boldsymbol{\nabla}_{\boldsymbol{p}_2}\boldsymbol{\cdot}\hat{\dot{\boldsymbol{p}}}_{21} \right]. \end{equation}

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 (3.17), and after some algebra, $C$ is given by

(3.18)\begin{align} C &= \frac{3 (nL^3)^2}{8 {\rm \pi}^5}\left[\frac{|\boldsymbol{p}_2-\boldsymbol{p}_1|^2}{[\boldsymbol{k}\wedge(\boldsymbol{p}_2- \boldsymbol{p}_1)]\boldsymbol{\cdot}[\boldsymbol{k}\wedge(\boldsymbol{p}_2-\boldsymbol{p}_1)]}\right] \frac{1}{\left(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_1+\boldsymbol{p}_2)\right)} \sin^2\left(\frac{\rm \pi}{4}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_1+\boldsymbol{p}_2)\right) \nonumber\\ &\quad \times\sin\left(\frac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_1+\boldsymbol{p}_2)\right)\left[\boldsymbol{p}_2 \boldsymbol{\cdot}\boldsymbol{p}_1 - \frac{\left(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_1+\boldsymbol{p}_2) \right)^2|\boldsymbol{p}_2-\boldsymbol{p}_1|^2}{4\left(\boldsymbol{k} \wedge(\boldsymbol{p}_2-\boldsymbol{p}_1)\right)\boldsymbol{\cdot}\left(\boldsymbol{k} \wedge(\boldsymbol{p}_2-\boldsymbol{p}_1)\right)}\right]. \end{align}

Using (3.18) in (3.14), the Fourier-transformed pair-probability density at, $O(\ln \kappa )^{-1}$, may be expressed as

(3.19)\begin{align} \hat{\varOmega}_2^{(1)} &= \frac{3 (nL^3)^2}{8{\rm \pi}^5}\left[\frac{|\boldsymbol{p}_2-\boldsymbol{p}_1|^2}{\left(\boldsymbol{k}\wedge(\boldsymbol{p}_2-\boldsymbol{p}_1)\right)\boldsymbol{\cdot}\left(\boldsymbol{k}\wedge(\boldsymbol{p}_2-\boldsymbol{p}_1)\right)}\right] \frac{1}{\left(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_1+\boldsymbol{p}_2)\right)} \sin^2\left(\frac{\rm \pi}{4}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_1+\boldsymbol{p}_2)\right) \nonumber\\ &\quad \times\sin\left(\frac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_1+\boldsymbol{p}_2)\right) \left[\boldsymbol{p}_2\boldsymbol{\cdot}\boldsymbol{p}_1 - \frac{\left(\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_1+\boldsymbol{p}_2)\right)^2|\boldsymbol{p}_2- \boldsymbol{p}_1|^2}{4\left(\boldsymbol{k}\wedge(\boldsymbol{p}_2-\boldsymbol{p}_1)\right)\boldsymbol{\cdot} \left(\boldsymbol{k}\wedge(\boldsymbol{p}_2-\boldsymbol{p}_1)\right)}\right] \boldsymbol{\delta}\left[\boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_2-\boldsymbol{p}_1)\right] \nonumber\\ &\quad - \frac{3 \textrm{i} (nL^3)^2}{32{\rm \pi}^6 k^3}PV \left[\frac{1}{\hat{\boldsymbol{k}}\boldsymbol{\cdot}(\boldsymbol{p}_2-\boldsymbol{p}_1)}\right] \left(\boldsymbol{\mathsf{I}}-\hat{\boldsymbol{k}}\hat{\boldsymbol{k}}\right) \boldsymbol{:}\boldsymbol{p}_1\boldsymbol{p}_2 \left[\frac{1}{\left(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)} \sin^2\left(\frac{\rm \pi}{2} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)\sin\left({\rm \pi} \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_2\right)\right. \nonumber\\ &\quad +\!\left.\frac{1}{\left(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2\right)} \sin^2\left(\frac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2\right) \sin\left({\rm \pi} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right) \right] , \end{align}

where both terms in (3.19) are generalized functions. Thus, the integral involving the argument of $PV[\cdot ]$ needs to be interpreted in the sense of a Cauchy principal value. The analysis of the homogeneous solution above is crucial, since it is this term alone that contributes to the covariance.

One can now use (3.19) in (3.10), and express the Fourier transformed correlated covariance in a $\boldsymbol {k}$-aligned polar coordinate system. Denoting the polar angles of the two swimmers as $\theta _1$ and $\theta _2$, and their relative azimuthal angle as $\phi$, one has $\boldsymbol {p}_1 = \cos \theta _1 \hat {\boldsymbol {k}} + \sin \theta _1\cos \phi \hat {\boldsymbol {k}}_\perp + \sin \theta _1\sin \phi \hat {\boldsymbol {k}}_\perp ^\prime$ and $\boldsymbol {p}_2 = \cos \theta _2 \hat {\boldsymbol {k}} + \sin \theta _2 \hat {\boldsymbol {k}}_\perp$. The correlated component of the covariance in Fourier space may then be written as

(3.20)\begin{align} &\langle\hat{\boldsymbol{u}}(\boldsymbol{k})\boldsymbol{\cdot}\hat{\boldsymbol{u}}(\boldsymbol{k}^\prime)\rangle\vert_{corr} \nonumber\\ &\quad =\frac{3(nL^3)^2\boldsymbol{\delta}(\boldsymbol{k}+\boldsymbol{k}^\prime)}{256{\rm \pi}^6(\ln\kappa)^3 k^7} \left [ \int_{-1}^{1}\mathrm{d}\mu_2(1-\mu_2^2)^2 \sin({\rm \pi} k \mu_2)\sin^3\left(\frac{\rm \pi}{2}k\mu_2\right) \left(j_0\left(\frac{\rm \pi}{2}k\mu_2\right)\right)^3\right. \nonumber\\ &\qquad+\!\left.\textrm{i}\int_{-1}^1\mathrm{d}\mu_1\int_{-1}^1 \mathrm{d}\mu_2 \frac{(1-\mu_1^2)^2(1-\mu_2^2)^2}{\mu_2-\mu_1} j_0\left(\frac{\rm \pi}{2}k\mu_1\right) j_0 \left(\frac{\rm \pi}{2}k\mu_2\right)\sin^2\left(\frac{\rm \pi}{2}k\mu_1\right)\right.\nonumber\\ &\qquad \times\!\left. \sin^2\left(\frac{\rm \pi}{2}k\mu_2\right)\!\!\vphantom{\int_{-1}^{1}}\left\{\vphantom{\int_{-1}^{1}} j_0\left(\frac{\rm \pi}{2}k\mu_1\right) \cos\left(\frac{\rm \pi}{2}k\mu_2\right) + j_0\left(\frac{\rm \pi}{2}k\mu_2\right) \cos\left(\frac{\rm \pi}{2}k\mu_1\right)\right\} \right ], \end{align}

where $\mu _i \equiv \cos \theta _i$ and $j_0(z) = \sin z/z$ is the spherical Bessel function of the first kind (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik1980). Here, we have performed the simpler $\phi$ integration for both the homogeneous and the particular solutions, and the delta function in the homogeneous solution of $\hat {\varOmega }_2^{(1)}$ has enabled us to evaluate the $\mu _1$ integral (the first term on the right-hand side of (3.20)). The contribution from the particular solution of $\hat {\varOmega }_2^{(1)}$, represented by the second term on the right-hand side of (3.20) vanishes upon integration. To see this, one may transform from ($\mu _1$, $\mu _2$) to ($\mu _+$, $\mu _-$), where $\mu _+ = \mu _1+\mu _2$ and $\mu _- = \mu _2-\mu _1$. In terms of the new integration variables, the particular solution integrand is odd in $\mu _-$ and therefore vanishes on applying the integration limits, implying that the ${O}(nL^3)^2$ contribution to the covariance comes only from the homogeneous solution as stated above. On inverse Fourier transforming (3.20), one obtains the $O(nL^3)^2$ component of the covariance in physical space, which is expressible as

(3.21)\begin{align} \langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle\vert_{corr} &= \frac{3(nL^3)^2}{128{\rm \pi}^6(\ln\kappa)^3 r}\int_0^\infty \mathrm{d}k \frac{\sin(2{\rm \pi} kr)}{k^6} \int_{-1}^1 \mathrm{d}\mu (1-\mu^2)^2 \sin({\rm \pi} k \mu)\nonumber\\ &\quad \times\sin^3\left(\frac{\rm \pi}{2}k\mu\right)\left(j_0\left(\frac{\rm \pi}{2}k\mu_2\right)\right)^3. \end{align}

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

(3.22)\begin{align} \langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x})\rangle\vert_{corr} = \frac{3(nL^3)^2}{64{\rm \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({\rm \pi} k \mu) \frac{1}{\left(\dfrac{\rm \pi}{2}k\mu\right)^3}\sin^6\left(\frac{\rm \pi}{2}k\mu\right). \end{align}

The integrand for $k$ in (3.22) scales as $1/k$ as $k\to 0$, implying that the integral itself is logarithmically divergent. 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$, in agreement with the scaling arguments in § 2. The physical arguments for this logarithmic divergence may now be laid out in a little more of an equation-specific manner than in § 2. As mentioned therein, the divergence arises from the slow ${O}(1/r^2)$ decay of the pair-probability density in the far-field. One may observe this from the physical space form of $\varOmega _2^{(1)}$, given by (3.12), wherein integrating along the trajectory of relative swimming, that is $\varOmega _{2}^{(1)} = -(nL^3)^2/[(4{\rm \pi} )^2|\boldsymbol {p}_2-\boldsymbol {p}_1|]\int _{-\infty }^z \mathrm {d}z^\prime (\boldsymbol {\nabla }_{\boldsymbol {p}_1}\boldsymbol {\cdot }\dot {\boldsymbol {p}}_{12} + \boldsymbol {\nabla }_{\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. This can be readily inferred from (3.19) as well, where power counting in $k$, as $k\to 0$ suggests that $\hat {\varOmega }_2^{(1)}\sim {O}(1/k)$ (which in real space would imply $\varOmega _2^{(1)}\sim {O}(1/r^2)$ in the far-field); here, we have used the scaling $\delta (k) \sim O(1/k)$ to arrive at the aforementioned argument. In this scenario, where a given swimmer interacts pairwise 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 contribution given by (3.8), and the correlated one given above, takes the form $c_1 (nL^3) + (nL^3)^2\{c_2 \ln [L_{screen}/L] + c_3\}$, with the constants $c_i$ being functions of swimmer aspect ratio. The screening length depends on the particular physical 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}$ replacing 0 as the lower limit of the Fourier integral in (3.22) and, thereby, enforcing a long wavelength cutoff. For suspensions of interacting straight-swimmers, the box-size limitation arises provided the swimmer mean free path, which is ${O}(nL^2)^{-1}$, is larger than $L_{box}$. Note that, for $L_{box} \gg (nL^2)^{-1}$, such suspensions are always unstable to sufficiently long wavelength orientation and velocity perturbations regardless of the swimmer volume fraction (Saintillan & Shelley Reference Saintillan and Shelley2008; Subramanian & Koch Reference Subramanian and Koch2009); this may be seen from the threshold criterion for instability (collective motion) given in the introduction, but with $\tau$ being set to infinity for straight-swimmers. One expects a transition to collective motion to be accompanied by a sharp box-size-dependent increase in the fluid velocity variance and tracer diffusivities (Saintillan & Shelley Reference Saintillan and Shelley2011; Krishnamurthy & Subramanian Reference Krishnamurthy and Subramanian2015; Stenhammar et al. Reference Stenhammar, Nardini, Nash, Marenduzzo and Morozov2017). The longest wavelength orientation (and velocity) modes that control the onset of collective motion and, thence, the amplitude of the velocity variance, have a nematic symmetry. The assumption of a homogeneous isotropic state in the analysis above therefore concerns the long-ranged correlations that develop either prior to the onset of collective motion (implying box sizes $L_{box} \sim O(nL^2)^{-1}$), or on length scales smaller than those corresponding to the aforementioned long wavelength (nematic) modes after onset. Simulations of interacting regularized force dipole swimmers carried out by Underhill & Graham (Reference Underhill and Graham2011), 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, Underhill & Graham (Reference Underhill and Graham2011)). 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}\geqslant (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. Reference Underhill, Hernandez-Ortiz and Graham2008; Underhill & Graham Reference Underhill and Graham2011; Stenhammar et al. Reference Stenhammar, Nardini, Nash, Marenduzzo and Morozov2017; Bárdfalvy et al. Reference Bárdfalvy, Nordanger, Nardini, Morozov and Stenhammar2019).

Figure 3.

The ${O}(nL^3)^2$ correlated straight-swimmer variance plotted against the infrared cutoff ($k_{min}$) used to evaluate (3.22). As expected, it exhibits a logarithmic divergence with decreasing $k_{min}$. The swimmer aspect ratio is chosen to be $\kappa = 8$.

For real bacteria, intrinsic decorrelation mechanisms such as rotary diffusion and/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. Recall that in § 2, we presented scaling arguments for the latter case. It is worth pointing out here 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 Reference Subramanian and Koch2009), 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/(n L^3)] + c_3\}$; although, for reasons discussed above, this estimate will be significantly modified by collective motion (authors' unpublished observations). In the next subsection, we analyse a suspension of interacting RTPs in detail, deriving the above asymptotic forms with the precise expressions for the constants that appear in them. The covariance integral, (3.21), in the strict straight-swimmer limit ($U\tau /L = \infty$) is ill-posed, and therefore, the straight-swimmer limit is better interpreted as a limiting case of persistent RTPs analysed below.

3.3. Suspensions of interacting RTPs

The covariance expression in a suspension of interacting slender RTPs remains the same as derived above, and is given by (3.10). However, the kinetic equation governing the pair-probability density must include additional terms compared to (3.12) that describe the tumble dynamics of the individual bacteria. Therefore, in this case, $\hat {\varOmega }_2^{(1)}$ satisfies

(3.23)\begin{align} &\left(\frac{U\tau}{L}\right)(\boldsymbol{p}_2-\boldsymbol{p}_1)\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{r}} \varOmega_2^{(1)} + \left(\varOmega_2^{(1)} -\frac{1}{4{\rm \pi}}\int\mathrm{d}\boldsymbol{p}_1 \varOmega_2^{(1)}\right) + \left(\varOmega_2^{(1)} -\frac{1}{4{\rm \pi}}\int\mathrm{d}\boldsymbol{p}_2 \varOmega_2^{(1)}\right)\nonumber\\ &\quad = - \frac{(nL^3)^2}{(4{\rm \pi})^2} \left(\frac{U\tau}{L}\right) \left[\boldsymbol{\nabla}_{\boldsymbol{p}_1} \boldsymbol{\cdot}\dot{\boldsymbol{p}}_{12} + \boldsymbol{\nabla}_{\boldsymbol{p}_2}\boldsymbol{\cdot} \dot{\boldsymbol{p}}_{21}\right ]. \end{align}

The bracketed terms on the left-hand side of (3.23) correspond to the individual swimmers undergoing tumbles in accordance with Poisson statistics (Subramanian & Koch Reference Subramanian and Koch2009; Othmer, Dunbar & Alt Reference Othmer, Dunbar and Alt1988), 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$ corresponds to the limiting case of straight-swimmers analysed in the previous subsection. Again, on Fourier transforming, one gets

(3.24)\begin{align} &\left(\frac{U\tau}{L}\right)2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_2-\boldsymbol{p}_1)\hat{\varOmega}_2^{(1)} + \left(2 \hat{\varOmega}_2^{(1)} -\frac{1}{4{\rm \pi}}\int\mathrm{d}\boldsymbol{p}_1 \hat{\varOmega}_2^{(1)} -\frac{1}{4{\rm \pi}}\int\mathrm{d}\boldsymbol{p}_2 \hat{\varOmega}_2^{(1)}\right) \nonumber\\ &\quad = - \frac{(nL^3)^2}{(4{\rm \pi})^2}\left(\frac{U\tau}{L}\right) \left[\boldsymbol{\nabla}_{\boldsymbol{p}_1}\boldsymbol{\cdot}\hat{\dot{\boldsymbol{p}}}_{12} +\boldsymbol{\nabla}_{\boldsymbol{p}_2}\boldsymbol{\cdot}\hat{\dot{\boldsymbol{p}}}_{21}\right]. \end{align}

One may, in fact, obtain an analytical solution for $\hat {\varOmega }_2^{(1)}(\boldsymbol {p}_1, \boldsymbol {p}_2; \boldsymbol {k})$ in (3.24), for arbitrary $U\tau /L$, by determining 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 the spectral representation of the Green's function in terms of an eigenfunction expansion. Further, the singular nature of the linear operator implies that both the discrete and (singular) continuous spectra need to be accounted for in this expansion. For the sake of brevity, we will be reporting this detailed Green's function analysis elsewhere (authors' unpublished observations).

For purposes of the present calculation, one notes that the symmetry of the orientation dilatation forcing terms on the right-hand side of (3.24) ensures that one does not require information contained in the full Green's function. A representation of $\hat {\varOmega }_2^{(1)}$ as a convolution of the Green's function with the right-hand side of (3.24) implies that only those eigenfunctions in the aforementioned spectral expansion, that the orientation dilatation terms project onto, contribute. To see this, we again choose a spherical polar coordinate system with its polar axis along $\hat {\boldsymbol {k}}$, with $\theta _1$ and $\theta _2$ representing the polar angles, and $\phi$ being the relative azimuthal angle between the two swimmers, measured in a plane perpendicular to $\boldsymbol {k}$ (similar to that described in § 3.2). In this coordinate system, the orientation dilatation terms are proportional to $\cos \phi$ and, hence, vanish when integrated in orientation space. Thus, one only needs to solve for the eigenmodes proportional to $\cos \phi$ in the full expansion for the pair-probability density. Now, eigenmodes proportional to $\cos (n \phi )$ or $\sin (m\phi )$ remain unaffected by the integral terms in (3.24) 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 {\varOmega }_2^{(1)}\,\mathrm {d}\boldsymbol {p}_1 = \int \hat {\varOmega }_2^{(1)}\,\mathrm {d}\boldsymbol {p}_2 = 0$ in (3.24), and the governing equation for $\hat {\varOmega }_2^{(1)}$ (for purposes of determining the covariance) reduces to the following much simpler algebraic form:

(3.25)\begin{align} \left(\frac{U\tau}{L}\right)2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_2-\boldsymbol{p}_1)\hat{\varOmega}_2^{(1)} + 2 \hat{\varOmega}_2^{(1)} = - \frac{(nL^3)^2}{(4{\rm \pi})^2}\left(\frac{U\tau}{L}\right) \left[\boldsymbol{\nabla}_{\boldsymbol{p}_1}\boldsymbol{\cdot}\hat{\dot{\boldsymbol{p}}}_{12} +\boldsymbol{\nabla}_{\boldsymbol{p}_2}\boldsymbol{\cdot}\hat{\dot{\boldsymbol{p}}}_{21}\right]. \end{align}

Solving (3.25) is straightforward, and yields the following exact expression:

(3.26)\begin{align} \hat{\varOmega}_2^{(1)} &= \frac{3 (nL^3)^2}{32{\rm \pi}^5 k^2}\left(\frac{U\tau}{L}\right)\left(\frac{1}{{\rm \pi} \textrm{i} (U\tau/L)\boldsymbol{k}\boldsymbol{\cdot} (\boldsymbol{p}_2-\boldsymbol{p}_1) + 1}\right) \left({\boldsymbol{\mathsf{I}}}- \hat{\boldsymbol{k}}\hat{\boldsymbol{k}}\right)\boldsymbol{:}\boldsymbol{p}_2\boldsymbol{p}_1 \nonumber\\ &\quad\times\! \left[\frac{1}{\left(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)} \sin^2\left(\frac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right) \sin\left({\rm \pi} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2\right) + \frac{1}{\left(\boldsymbol{k}\boldsymbol{ \cdot}\boldsymbol{p}_2\right)} \sin^2\left(\frac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_2\right)\sin\left({\rm \pi} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right) \right] . \end{align}

The pair-correlation function for RTPs above does not involve any concentration fluctuations, as may be seen by integrating (3.26) over orientation space ($\boldsymbol {p}_1$ and $\boldsymbol {p}_2$). This implies that the pair distribution function, $g(r) = (nL^3)^{-2}\int \varOmega _2\,\mathrm {d}\boldsymbol {p}_1 \,\mathrm {d}\boldsymbol {p}_2$ is just unity. Therefore, pair-correlations between slender swimmers involve only the orientation degrees of freedom, with the swimmer concentration (positional) fluctuations still predicted to be similar to that in the absence of interactions (and, therefore, Poissonian in character); this is consistent with the simulation results of Underhill & Graham (Reference Underhill and Graham2011). In fact, recent simulations show that the pair distribution function remains close to unity even after the onset of collective motion (Krishnamurthy & Subramanian Reference Krishnamurthy and Subramanian2015; Bárdfalvy et al. Reference Bárdfalvy, Nordanger, Nardini, Morozov and Stenhammar2019). From (3.26), one can also infer that for large $U\tau /L$, $\hat {\varOmega }_2^{(1)}\sim {O}(1/k)$ in the range $U\tau /L\gg 1/k\gg 1$ and, therefore, $\varOmega _2^{(1)}$ has the character of an ${O}(1/r^2)$ decay for $1\ll r\ll U\tau /L$. Recall the discussion in § 3.2, where this algebraic decay of $\varOmega _2^{(1)}$ led to a logarithmically divergent variance in the limit $U\tau /L \to \infty$ (straight-swimmers); in appendix B, we recover the straight-swimmer form of the pair-probability density from (3.26). For any large but finite $U\tau /L$, however, $\varOmega _2^{(1)}$ exhibits faster-than-algebraic decay for $r \gg U\tau /L$, on account of decorrelations induced by tumbling. In the other limit of rapid tumblers, corresponding to $U\tau /L \ll 1$, $\varOmega _2^{(1)}$ directly transitions to a faster-than-algebraic decay for $r\gg 1$.

Now, using (3.26) in (3.10), and on inverse Fourier transforming, one obtains the correlated contribution to the fluid velocity covariance in a suspension of interacting RTPs to be

(3.27)\begin{align} \langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x} + \boldsymbol{r})\rangle\vert_{corr} &= \frac{3(nL^3)^2}{128{\rm \pi}^6(\ln\kappa)^3 r}\left(\frac{U\tau}{L}\right) \int_0^\infty \mathrm{d}k \frac{\sin(2{\rm \pi} k r)}{k^5}\int_{-1}^1\,\mathrm{d}\mu_1\int_{-1}^1\mathrm{d}\mu_2 (1-\!\mu_1^2)\nonumber\\ &\quad \times(1-\mu_2^2)\left(\frac{1}{1 + {\rm \pi}\textrm{i}(U\tau/L)k(\mu_2-\mu_1)}\right) \sin^2\left( \frac{\rm \pi}{2}k\mu_1\right)\sin^2\left(\frac{\rm \pi}{2}k\mu_2\right)\nonumber\\ &\quad \times j_0\left(\frac{\rm \pi}{2}k\mu_1\right) j_0\left(\frac{\rm \pi}{2}k\mu_2\right) \left[\cos\left(\frac{\rm \pi}{2}k\mu_2\right)j_0\left(\frac{\rm \pi}{2}k\mu_1\right) \right.\nonumber\\ &\quad+\!\left.\cos\left(\frac{\rm \pi}{2}k\mu_1\right)j_0\left(\frac{\rm \pi}{2}k\mu_2\right)\right], \end{align}

where $\mu _i\equiv \cos \theta _i$. The correlated variance is obtained by setting $r=0$ in (3.27), and is expressible as

(3.28)\begin{align} \langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x})\rangle\vert_{corr} &= \frac{3(nL^3)^2}{64{\rm \pi}^5(\ln\kappa)^3}\left(\frac{U\tau}{L}\right) \int_0^\infty \mathrm{d}k \frac{1}{k^4}\int_{-1}^1\mathrm{d}\mu_1\int_{-1}^1\mathrm{d}\mu_2 (1-\mu_1^2)(1-\mu_2^2) \nonumber\\ &\quad\times \!\left(\frac{1}{1 \!+\! {\rm \pi}\textrm{i}(U\tau/L)k(\mu_2-\mu_1)}\right) \sin^2\left(\frac{\rm \pi}{2}k\mu_1\right) \sin^2\left(\frac{\rm \pi}{2}k\mu_2\right) j_0\left(\frac{\rm \pi}{2}k\mu_1\right) \nonumber\\ &\quad \times j_0\left(\frac{\rm \pi}{2}k\mu_2\right)\left[\cos\left(\frac{\rm \pi}{2}k\mu_2\right)j_0 \left(\frac{\rm \pi}{2}k\mu_1\right)+ \cos\left(\frac{\rm \pi}{2}k\mu_1\right)j_0\left(\frac{\rm \pi}{2}k\mu_2\right)\right]. \end{align}

In figure 4, we plot the correlated variance given by (3.28). For $U\tau /L \ll 1$, owing to the aforementioned rapid decay of $\varOmega _2^{(1)}$, the correlated variance scales as $c_1 (U\tau /L)(nL^3)^2$ (see inset of figure 4). This estimate arises from neglecting the factor ${\rm \pi} \textrm {i}(U\tau /L)k(\mu _2-\mu _1)$ in the denominator of the right-hand side in (3.28). For large $U\tau /L$, in accordance with earlier scaling arguments, the correlated variance in suspensions of interacting RTPs takes the form $c_2 (nL^3)^2/(\ln \kappa )^3 \ln (U\tau /L)$. As stated above, this is a consequence of the $O(1/k)$ scaling of $\varOmega _2^{(1)}$ for persistent RTPs, which leads to a logarithmic divergence in (3.28) that is eventually cut off for $k \sim O(U\tau /L)^{-1}$.

Figure 4.

The $O(nL^3)^2$ correlated variance, as given by (3.28), plotted as a function of $U\tau /L$. The small (dotted) and large (dash-dotted) $U\tau /L$ asymptotes are shown; the inset highlights the linear scaling of the correlated variance with $U\tau /L$ for $U\tau /L\ll 1$.

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)$ far-field decay for $r\gg 1$. The nature of this far-field decay remains the same irrespective of $U\tau /L$. For $U\tau /L\ll 1$, the far-field decay occurs when $r\gg 1$, whereas, for $U\tau /L\gg 1$, this only occurs when $r\gg U\tau /L$. To obtain this far-field dipole asymptote, one may rescale the outer integral variable to be $\tilde {k} = kr$, and seek the asymptotic form of (3.27) for $\tilde {k}\sim O(1)$ but $r\gg 1$. Using the approximations $\sin [{\rm \pi} \tilde {k}\mu _i/(2r)]\approx {\rm \pi}\tilde {k}\mu _i/(2r)$, $j_0[{\rm \pi} \tilde {k}\mu _i/(2r)] \approx 1$ and neglecting the term involving $(\mu _2-\mu _1)$, one finally obtains the covariance as $\langle \boldsymbol {u}(\boldsymbol {x})\boldsymbol {\cdot }\boldsymbol {u}(\boldsymbol {x} + \boldsymbol {r})\rangle \vert _{corr} = 3 (nL^3)^2 (U\tau /L)(28\,800 {\rm \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 of the far-field asymptotes for $r\gg U\tau /L$. Both from figure 5 and the inset therein, one notes, for $U\tau /L \gg 1$, the emergence of an intermediate asymptotic regime following the variance plateau, characterized by a weak decay in the range $1\ll r\ll U\tau /L$. As is shown in detail in appendix C, the weak decay of the covariance in the intermediate region is a logarithmic one, scaling as ${O}(\ln [r/(U\tau )])$.

Figure 5.

The O $(nL^3)^2$ component of the (normalized) fluid velocity covariance as a function of $r$. In the inset the covariance curves collapse onto the far-field dipole asymptote (represented by the $\square$ symbols) stated below (3.27), for $r\gg U\tau /L$, when the abscissa is scaled with $(U\tau /L)^{-1}$.

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 (3.8) (see § 3.1), and the correlated contribution to the variance given by (3.28) above. This is expressible as

(3.29)\begin{align} &\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x})\rangle\nonumber\\ &\quad= \frac{(nL^3)}{96{\rm \pi}(\ln\kappa)^2} \pm \frac{3(nL^3)^2}{64{\rm \pi}^5(\ln\kappa)^3}\left(\frac{U\tau}{L}\right) \int_0^\infty \frac{\mathrm{d}k}{k^4} \int_{-1}^1\mathrm{d}\mu_1\int_{-1}^1\mathrm{d}\mu_2 (1-\mu_1^2)(1-\mu_2^2) \nonumber\\ &\quad\times\!\left(\frac{1}{1 + {\rm \pi}\textrm{i}(U\tau/L)k(\mu_2-\mu_1)}\right) \sin^2\left(\frac{\rm \pi}{2}k\mu_1\right)\sin^2 \left(\frac{\rm \pi}{2}k\mu_2\right)\nonumber\\ &\quad \times j_0\left(\frac{\rm \pi}{2}k\mu_1\right) j_0\left(\frac{\rm \pi}{2}k\mu_2\right) \left[\cos\left(\frac{\rm \pi}{2}k\mu_2\right)j_0\left(\frac{\rm \pi}{2}k\mu_1\right) + \cos\left(\frac{\rm \pi}{2}k\mu_1\right)j_0\left(\frac{\rm \pi}{2}k\mu_2\right)\right], \end{align}

where the positive and negative signs pertain to suspensions of pushers and pullers, respectively. Based on the preceding arguments, in a suspension of pusher-type swimmers, the total variance scales as $c_1(nL^3)/(\ln \kappa )^2 + (nL^3)^2/(\ln \kappa )^3[c_2 \ln [U\tau /L] + c_3]$ for large $U\tau /L$, and as $c_1(nL^3)/(\ln \kappa )^2 +c_2 (U\tau /L)(nL^3)^2/(\ln \kappa )^3$ 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 involves a dimensional scaling of $(U/\ln \kappa )^2$, implying that it scales as the square of the dipole strength $D$ (measured in units of $\eta U L^2$). However, the correlated contribution to the variance is proportional to $(U/\ln \kappa )^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 Reference Krishnamurthy and Subramanian2015; Stenhammar et al. Reference Stenhammar, Nardini, Nash, Marenduzzo and Morozov2017), at ${O}(nL^3)^2$, (3.29) predicts enhanced fluctuations in suspensions of pusher-type swimmers. In figure 6, this enhancement is approximately 30 % for an $nL^3$ of approximately 4, and for $U\tau /L = 125$. As mentioned in § 2, the enhancement is much greater for tracer diffusivities on account of the linear scaling in $U\tau /L$ (see § 4.3).

Figure 6.

The total velocity variance plotted as a function of the swimmer volume fraction ($nL^3$) at $U\tau /L = 125$, for both pushers and pullers, along with the uncorrelated contribution (solid line). The swimmer aspect ratio is chosen to be $\kappa = 8$.

The total covariance in a suspension of interacting RTPs can be similarly expressed by combining (3.8) and (3.27), which gives

(3.30)\begin{align} &\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u}(\boldsymbol{x}+\boldsymbol{r})\rangle\nonumber\\ &\quad= \frac{1}{4{\rm \pi}^6 r}\int_0^\infty \mathrm{d}k \frac{\sin(2{\rm \pi} kr)}{k^5}\left\{ \frac{(nL^3)}{(\ln\kappa)^2}\int_{-1}^1\mathrm{d}\mu\left(\frac{1-\mu^2}{\mu^2}\right)\sin^4\left(\frac{{\rm \pi} k \mu}{2}\right)\right.\nonumber\\ &\qquad \pm \frac{3(nL^3)^2}{32(\ln\kappa)^3}\left(\frac{U\tau}{L}\right) \int_{-1}^1 \mathrm{d}\mu_1\int_{-1}^1\mathrm{d}\mu_2 (1-\mu_1^2)(1-\mu_2^2) \nonumber\\ &\qquad\times \!\left(\frac{1}{1 + {\rm \pi}\textrm{i}(U\tau/L)k(\mu_2-\mu_1)}\right) \sin^2\left(\frac{\rm \pi}{2}k\mu_1\right) \sin^2\left(\frac{\rm \pi}{2}k\mu_2\right)\nonumber\\ &\qquad \times \!\left. j_0\left(\frac{\rm \pi}{2}k\mu_1\right) j_0\left(\frac{\rm \pi}{2}k\mu_2\right) \left[\cos\left(\frac{\rm \pi}{2}k\mu_2\right)j_0\left(\frac{\rm \pi}{2}k\mu_1\right) + \cos\left(\frac{\rm \pi}{2}k\mu_1\right)j_0\left(\frac{\rm \pi}{2}k\mu_2\right)\right]\right\}, \end{align}

where the positive sign, corresponding to an increase in the covariance due to correlations, pertains to suspensions of pushers. Figure 7 shows a plot of (3.30) as a function of $r$ for $nL^3 = 5$, for pusher suspensions (the positive sign in (3.30)). Note that the chosen $nL^3$ is a little greater than unity. There is evidence, from earlier calculations for suspensions of non-Brownian slender fibres, that the dilute theory, with pair interactions accounted for, remains valid for $nL^3$ up until around 7 (Mackaplow & Shaqfeh Reference Mackaplow and Shaqfeh1996); the extended regime of validity for the dilute approximation is also evident for Brownian fibre suspensions, where the transition of the single fibre self-diffusivity from the isolated Stokes–Einstein value (valid in the dilute regime) to an $O(nL^3)^{-2}$ scaling arising from steric hindrance, characteristic of the semidilute regime, occurs only at an $nL^3$ of approximately 16 (see chapter 10 in Larson (Reference Larson2013)). In figure 7, for large $U\tau /L$, while the variance plateau only grows slowly with increasing $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 (3.28)). This differential growth results in the total covariance developing a shoulder in the intermediate regime ($1\ll r\ll U\tau /L$).

Figure 7.

The total fluid velocity covariance for a suspension of pusher-type RTPs plotted as a function of $r$ over a range of $U\tau /L\in [0.125, 750]$, and for $nL^3 = 5$. The red (dash-dotted) arrow highlights the direction of increasing $U\tau /L$. The swimmer aspect ratio is chosen to be $\kappa = 8$.

We now analyse the weak decay of the correlated component of the covariance in the intermediate asymptotic regime (corresponding to the shoulder) mentioned in the above discussion. Although the exact solution given by (3.27) is valid for arbitrary $U\tau /L$, the form of the intermediate regime is not obvious. Therefore, for large $U\tau /L$, we derive an alternate simpler approximation which reveals an intermediate asymptotic regime of shallow decay for separations in the range $1\ll r\ll U\tau /L$ (as implied by figures 5 and 7). In this limit, owing to a separation of scales corresponding to those of order the swimmer size and those comparable to the run length, the method of matched asymptotic expansions is applicable. The details of the analysis are quite involved, and are therefore relegated to appendix C. The nature of the different asymptotic regions involved may, however, be readily explained. In a reference frame moving with one of the swimmers (termed the ‘test swimmer’), the physical domain may be divided into the following two regions.

1. (i) 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 that develop along straight-swimmer trajectories. The expressions pertaining to this region have already been obtained as part of the straight-swimmer analysis given in § 3.2.

2. (ii) An outer region $r\sim O(U\tau /L) (\gg 1)$, 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 of the rotation rates ($\dot {\boldsymbol {p}}_{ij}$) involved (see appendix A).

The matching region, which therefore corresponds to straight-swimming point force-dipoles interacting in the range $1\ll r\ll U\tau /L$, is the region where the expressions in both the inner and outer regions must reduce to a common form. A schematic highlighting the aforementioned demarcation of the physical domain, is shown in figure 8 (see also appendix C). 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 \geqslant U\tau /L$, characteristic of interacting RTPs; the covariance plotted in figure 5 does conform to these scalings. The matched asymptotic expansions calculation in appendix C yields the following expression for the correlated component of the covariance:

(3.31)\begin{align} &\langle\boldsymbol{u}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{u} (\boldsymbol{x}+\boldsymbol{r})\rangle_{corr-matched}\nonumber\\ &\quad = \frac{(nL^3)^2D^3}{{\rm \pi}^2 r}\left\{\int_0^1\mathrm{d}k \int_{-1}^1\mathrm{d}\mu\left[\frac{3}{2{\rm \pi}^4}\frac{\sin(2{\rm \pi} k r)}{k^6}(1-\mu^2)^2 \right. \right.\nonumber\\ &\qquad \times \!\left.\sin({\rm \pi} k \mu)\left(j_0\left(\frac{\rm \pi}{2}k\mu\right)\right)^3 \sin^3\left(\frac{\rm \pi}{2}k\mu\right)- \frac{1}{105}\frac{\sin(2{\rm \pi} k r)}{k^2}\right]\nonumber\\ &\qquad+\int_1^\infty \mathrm{d}k\frac{\sin(2{\rm \pi} k r)}{k^6}\int_{-1}^1 \mathrm{d}\mu\frac{3}{2{\rm \pi}^4}(1-\mu^2)^2 \left.\left(j_0\left(\frac{\rm \pi}{2}k\mu\right)\right)^3\sin^3\left( \frac{\rm \pi}{2}k\mu\right)\right\} \nonumber\\ &\qquad + (nL^3)^2 D^3\left\{\frac{3}{16{\rm \pi}^2r_1}\int_0^\infty \frac{1}{k_1}\sin(2{\rm \pi} k_1 r_1)\, \mathrm{d}k_1\int_{-1}^1\mathrm{d}\mu_1\right.\nonumber\\ &\qquad\times \vphantom{\int_0^\infty}\int_{-1}^1\mathrm{d}\mu_2 \frac{\mu_1^2\mu_2^2(1-\mu_1^2)(1-\mu_2^2)}{1 + {\rm \pi}^2\left(\dfrac{U\tau}{L}\right)k_1^2(\mu_2-\mu_1)^2} \left.-\frac{2}{105{\rm \pi}}\left[j_0(2{\rm \pi} r_1) - C_i(2{\rm \pi} r_1)\right]\right\}\nonumber\\ &\qquad + \frac{2(nL^3)^2 D^3}{105{\rm \pi}}\left\{\ln(2{\rm \pi} r) - j_0(2{\rm \pi} r ) + \int_0^{2{\rm \pi} r } \frac{\cos(t) -1}{t}\,\mathrm{d}t \right. \nonumber\\ &\left.\qquad \vphantom{\frac{2(nL^3)^2 D^3}{105{\rm \pi}}}+ j_0(2{\rm \pi} r_1) - C_i(2{\rm \pi} r_1) + \gamma_E \right\}, \end{align}

where $C_i(x)$ represents the cosine integral (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik1980), $r_1 = r L/(U\tau )$ and $k_1$ the corresponding wavenumber; $D = 1/(4\ln \kappa )$ is the non-dimensional dipole strength corresponding to dipole swimmers. There are three sets of braced terms on the right-hand side of (3.31), and we briefly explain these in the following. In applying the matched asymptotic expansions method, one needs to take care of the apparent divergences that arise in a naïve straight-swimmer approximation. In the outer region, the slender straight-swimmer covariance has a non-decaying character, whereas, in the inner region, the far-field dipole approximation diverges, owing to the $O(1/r^2)$ scaling of the dipole velocity field. Therefore, one must ensure that both of the aforementioned divergences are accounted for and isolated in the matching scheme. In (3.31), the first set of braces on the right-hand side corresponds to the slender straight-swimmer contribution with a far-field divergence due to straight-swimming dipoles subtracted off. To avoid the (obvious) divergence of the subtracted dipoles in the inner region $k\to \infty$ ($r\to 0$), we set an upper bound on the $k$-integration of the subtracted component. The second set of braces corresponds to the far-field dipole contribution, now, with the inner-region divergence due to straight-swimming dipoles subtracted off. Thus, in this case, we set a lower bound in the $k_1$ integration of the subtracted quantity to again ensure that the region $k_1\to 0$ ($r_1\to \infty$) is removed. We add back these subtracted contributions in the third set of braces in (3.31). Note that it is the third grouping alone that depends on $U\tau /L$, and exhibits the required logarithmic scaling.

Figure 8.

A schematic showing the three distinct physical regions for $U\tau /L\gg 1$. The inner region corresponds to $r\sim {O}(1)$, where we have slender straight-swimmers, and the correlations develop over straight trajectories. The outer region $r\sim U\tau /L$ corresponds to the run-and-tumble dynamics of point force-dipoles. In the intermediate region of overlap, $1 \ll r \ll U\tau /L$, the inner and outer solutions must assume a common asymptotic form.

In figure 9, we compare the correlated contributions to the covariance calculated from (3.27) and (3.31). The matched asymptotic expansions approach yields excellent agreement with the exact covariance for the three largest values of $U\tau /L$, and remains a good approximation down to $U\tau /L = 7.5$. As stated in the aforementioned discussion, it is the straight-swimming point force-dipole contribution in the third set of braces on the right-hand side of (3.31), that captures the functional form in the intermediate asymptote regime. This term simplifies to $2(nL^3)^2/(105{\rm \pi} )[1 - \gamma _E - \ln (2{\rm \pi} r_1)]$, thereby exhibiting an ${O}(\ln [rL/(U\tau )])$ decay. Thus, the aforementioned weak decay of the $O(nL^3)^2$ component of the covariance, is a logarithmic one, and as $U\tau /L\to \infty$, one asymptotes to a non-decaying covariance in the straight-swimmer limit. The logarithmic functional form alone is compared with the exact solution for the full covariance, for $U\tau /L = 750$ in figure 10; a good agreement is obtained in the intermediate range of separations $1 \ll r\ll U\tau /L$.

Figure 9.

The fluid velocity covariance for arbitrary $U\tau /L$ from (3.27) (straight curves) is plotted alongside the same determined by the method of matched asymptotic expansion (3.31) (dashed curves), for $U\tau /L = 2.5$, 7.5, 25, 125, 250. The swimmer aspect ratio is chosen to be $\kappa = 8$.

Figure 10.

The correlated contribution to the fluid velocity covariance for $U\tau /L = 750$ and $\kappa = 8$ obtained from (3.27) (straight curves), along with the intermediate logarithmic scaling ($\diamond$) in the regime $1\ll r\ll U\tau /L$ obtained from (C 8) in appendix C.

Therefore, for $U\tau /L\gg 1$, the behaviour of the covariance, in an interacting suspension may be summarized as follows:

1. (i) for $r\ll 1$, it exhibits a finite variance plateau, with the variance itself scaling as $\ln U\tau /L$;

2. (ii) for $1\ll r\ll U\tau /L$, it exhibits a weak decay proportional to ${O}(\ln [rL/(U\tau )])$;

3. (iii) for $r\gg U\tau /L$, corresponding to the far-field, it decays as ${O}(1/r)$.

4.1. Suspensions of non-interacting swimmers

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

(4.5)\begin{equation} \langle\hat{\boldsymbol{u}}(\boldsymbol{k}, t_1)\boldsymbol{\cdot} \hat{\boldsymbol{u}}(-\boldsymbol{k},t_2)\rangle \vert_{uncorr} = \int \hat{\boldsymbol{u}}(\boldsymbol{k})\boldsymbol{\cdot}\hat{\boldsymbol{u}}(-\boldsymbol{k}) \hat{G}_1 \left( \boldsymbol{p}_1, t_1 \vert \boldsymbol{p}_1^\prime, t_2; \boldsymbol{k} \right) \varOmega_1 \,\mathrm{d}\boldsymbol{p}_1 \,\mathrm{d}\boldsymbol{p}_1^\prime, \end{equation}

where $\hat {\boldsymbol {u}}(\boldsymbol {k})$ is the Fourier-transformed disturbance velocity field (see (A 2)), $\varOmega _1 = nL^3/4{\rm \pi}$ is the steady-state singlet probability density and $\hat {G}_1(\boldsymbol {p}_1, t_1\vert \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$ given that the swimmer had an orientation $\boldsymbol {p}_{1}^{\prime }$ at $t_2$. The governing equation for $\hat {G}_1$ is

(4.6)\begin{equation} \frac{\mathrm{d} \hat{G}_1}{\mathrm{d} t} + \left(\frac{U\tau}{L}\right)2 {\rm \pi}\textrm{i} \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{p}_1 \hat{G}_1 + \left(\hat{G}_1 -\frac{1}{4{\rm \pi}}\int\mathrm{d}\boldsymbol{p}_1 \hat{G}_1 \right) = \boldsymbol{\delta}(\boldsymbol{p}_1 - \boldsymbol{p}_1^\prime)\delta(t - t^\prime) . \end{equation}

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

(4.7)\begin{equation} \hat{G}_1 \left( \boldsymbol{p}_1, t \vert \boldsymbol{p}_1^\prime, t^{\prime}; \boldsymbol{k} \right) = \delta(\boldsymbol{p}_1 - \boldsymbol{p}_1^\prime) \exp\left(\left(- 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_1 \frac{U \tau}{L} -1 \right) (t - t^\prime) \right). \end{equation}

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$, including the concentration modes, is needed to characterize the time-dependent evolution of an initially localized population of run-and-tumble swimmers, and will be reported separately (authors' unpublished observations). Using the above result, and upon simplifying, we obtain the expression for the MSD of a passive tracer as

(4.8)\begin{align} \langle x^2(t) \rangle\vert_{uncorr} & = \frac{nL^3}{(\ln\kappa)^2 8{\rm \pi}^7} \int\mathrm{d} \boldsymbol{p}_1 \int \frac{\mathrm{d} \boldsymbol{k}}{k^4} \left( \frac{-1 + \exp\left(-\left( 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_1 \dfrac{U \tau}{L} +1 \right) t \right)}{ \left( 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_1 \dfrac{U \tau}{L} +1 \right)^2} \right) \nonumber\\ &\quad\times \!\left( {\boldsymbol{\mathsf{I}}} - \frac{\boldsymbol{k}\boldsymbol{k}}{k^2}\right)\boldsymbol{:} \boldsymbol{p}_1\boldsymbol{p}_1\frac{\sin^4\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_1\right)}{(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1)^2} \nonumber\\ &\quad + \frac{nL^3}{(\ln\kappa)^2 8{\rm \pi}^7} t \int \mathrm{d} \boldsymbol{p}_1 \int \frac{\mathrm{d} \boldsymbol{k}}{k^4} \left( \dfrac{1}{ 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_1 \frac{U \tau}{L} +1} \right)\left( {\boldsymbol{\mathsf{I}}} - \frac{\boldsymbol{k}\boldsymbol{k}}{k^2}\right)\boldsymbol{:}\boldsymbol{p}_1\boldsymbol{p}_1\nonumber\\ &\quad\times \frac{\sin^4\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)}{(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1)^2} . \end{align}

In (4.8), the term of $O(U\tau /L)$ in the argument of the exponential may be neglected in the limit $U\tau /L \ll 1$, and this leads to a (dimensional) correlation time of $O(\tau )$. For $U\tau /L \gg 1$, the integral in (4.8) has a dominant contribution for $k \sim L/(U\tau )$ and hence the dimensional correlation time is $O(L/U)$ (Kasyap et al. Reference Kasyap, Koch and Wu2014). For times much shorter than the velocity correlation time, the tracer shows ballistic motion and the MSD simplifies to

(4.9)\begin{equation} \langle x^2(t) \rangle\vert_{uncorr} = \frac{nL^3}{(\ln\kappa)^2 16{\rm \pi}^7} t^2 \int \mathrm{d} \boldsymbol{p}_1 \int\frac{\mathrm{d} \boldsymbol{k}}{k^4} \frac{\sin^4\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)}{(\boldsymbol{k} \boldsymbol{\cdot}\boldsymbol{p}_1)^2} \left( {\boldsymbol{\mathsf{I}}} - \frac{\boldsymbol{k}\boldsymbol{k}}{k^2}\right)\boldsymbol{:}\boldsymbol{p}_1\boldsymbol{p}_1. \end{equation}

The coefficient of O$(t^2)$ in (4.9) is the fluid velocity variance given in (3.7). On the other hand, for times much longer than the velocity correlation time, the tracer shows diffusive motion with $\langle x^2(t) \rangle \vert _{uncorr} = 6 D_t\vert _{uncorr} t$, where the effective diffusivity is given by

(4.10)\begin{align} D_t\vert_{uncorr} = \frac{nL^3}{(\ln\kappa)^2 48{\rm \pi}^7} \int \mathrm{d} \boldsymbol{p}_1 \int\frac{\mathrm{d} \boldsymbol{k}}{k^4} \left( \frac{1}{ 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_1 \dfrac{U \tau}{L} +1} \right)\left( {\boldsymbol{\mathsf{I}}} - \frac{\boldsymbol{k}\boldsymbol{k}}{k^2}\right)\boldsymbol{:} \boldsymbol{p}_1\boldsymbol{p}_1\dfrac{\sin^4\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1 \right)}{(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1)^2} . \end{align}

The above expression may also be derived using the Green–Kubo formula (Balakrishnan Reference Balakrishnan2008) and has also been obtained by Kasyap et al. (Reference Kasyap, Koch and Wu2014) using a different approach. On the other hand, to the best of our knowledge, the complete expression for the MSD in (4.8) has not been given in the literature even for non-interacting swimmers. As discussed earlier in § 3.1, the integrals in (4.8) may be numerically evaluated with the choice of a $\boldsymbol {k}$-aligned spherical coordinate system. Recall that the tracer diffusivity ($D_t\vert _{uncorr}$) in (4.10) is non-dimensionalised using $U^2 \tau$ as the scale. For $U\tau /L \ll 1$, the integral in (4.10) is $O(1)$, and hence $D_t\vert _{uncorr} \sim nL^3/(\ln \kappa )^2 U^2 \tau$. For $U\tau /L \gg 1$, the dominant contribution results from $k \sim$ $O(L/U\tau )$ and, using this scale, the integral in (4.10) is $O(L/U\tau )$, with the tracer diffusivity, $D_t\vert _{uncorr} \sim nL^3/(\ln \kappa )^2 U L$ (Kasyap et al. Reference Kasyap, Koch and Wu2014). Recall that these scalings were obtained earlier in § 2 using physical arguments.

4.2. Suspensions of interacting swimmers

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

(4.11)\begin{align} \langle\hat{\boldsymbol{u}}(\boldsymbol{k}, t_1)\boldsymbol{\cdot} \hat{\boldsymbol{u}}(-\boldsymbol{k},t_2)\rangle \vert_{corr} &= \int \hat{\boldsymbol{u}}(\boldsymbol{k})\boldsymbol{\cdot}\hat{\boldsymbol{u}}(-\boldsymbol{k}) \hat{G}_2 \left( \boldsymbol{p}_1, \boldsymbol{p}_2, t_1 \vert \boldsymbol{p}_1^\prime, \boldsymbol{p}_2^\prime, t_2; \boldsymbol{k} \right) \hat{\varOmega}_2^{(1)} \left( \boldsymbol{k}, \boldsymbol{p}_1^\prime, \boldsymbol{p}_{2}^{\prime} \right) \nonumber\\ &\qquad \times\mathrm{d}\boldsymbol{p}_1 \,\mathrm{d}\boldsymbol{p}_1^\prime \,\mathrm{d}\boldsymbol{p}_2 \, \mathrm{d}\boldsymbol{p}_2^\prime, \end{align}

where $\hat {\varOmega }_2^{(1)}$ is the pair-probability density derived in § 3.3 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$, given that the pair had orientations $\boldsymbol {p}_{1}^{\prime }$ and $\boldsymbol {p}_{2}^{\prime }$ at time $t_2$. The governing equation for $\hat {G}_2$ is

(4.12)\begin{align} &\frac{\mathrm{d} \hat{G}_2}{\mathrm{d} t}+ 2{\rm \pi} \textrm{i} \left(\frac{U\tau}{L}\right) \boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_2-\boldsymbol{p}_1)\hat{G}_2 + \left(2 \hat{G}_2 -\frac{1}{4{\rm \pi}}\int\mathrm{d}\boldsymbol{p}_1 \hat{G}_2 -\frac{1}{4{\rm \pi}} \int\mathrm{d}\boldsymbol{p}_2 \hat{G}_2 \right) \nonumber\\ &\quad = \boldsymbol{\delta}(\boldsymbol{p}_2 - \boldsymbol{p}_2^\prime)\boldsymbol{\delta}(\boldsymbol{p}_1 - \boldsymbol{p}_1^\prime)\delta(t - t^\prime). \end{align}

The inverse tumbling terms in (4.12) again don't contribute (see § 3.3), and the expression for $\hat {G}_2$ relevant to evaluating the velocity correlation function is given by

(4.13)\begin{align} \hat{G}_2 \left( \boldsymbol{p}_1, \boldsymbol{p}_2, t \vert \boldsymbol{p^\prime}_1, \boldsymbol{p}_{2}^{\prime}, 0; \boldsymbol{k} \right) = \delta(\boldsymbol{p}_1 - \boldsymbol{p}_1^\prime) \delta( \boldsymbol{p}_2-\boldsymbol{p}_{2}^{\prime}) \exp\left(\left(- 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot}(\boldsymbol{p}_2-\boldsymbol{p}_1) \frac{U \tau}{L} -2 \right) t \right). \end{align}

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

(4.14)\begin{align} \langle x^2(t) \rangle\vert_{corr} & = \frac{1}{(\ln\kappa)^3 2{\rm \pi}^6} \int \mathrm{d} \boldsymbol{p}_1 \int\mathrm{d} \boldsymbol{p}_2 \int \frac{\mathrm{d} \boldsymbol{k}}{k^4} \left( \frac{-1 + \exp\left(\left(- 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot} (\boldsymbol{p}_2 - \boldsymbol{p}_1 ) \dfrac{U \tau}{L} -2 \right) t \right)}{ \left( 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot} (\boldsymbol{p}_2 - \boldsymbol{p}_1) \dfrac{U \tau}{L} +2 \right)^2} \right) \nonumber\\ & \quad\times \!\left( {\boldsymbol{\mathsf{I}}} - \frac{\boldsymbol{k}\boldsymbol{k}}{k^2}\right)\boldsymbol{:} \boldsymbol{p}_1\boldsymbol{p}_2\frac{\sin^2\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)} {(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1)} \dfrac{\sin^2\left(\dfrac{\rm \pi}{2}\boldsymbol{k} \boldsymbol{\cdot}\boldsymbol{p}_2\right)}{(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2)} \hat{\varOmega}_2^{(1)} \left( \boldsymbol{p}_1, \boldsymbol{p}_{2}; \boldsymbol{k} \right) \nonumber\\ &\quad + \frac{1}{(\ln\kappa)^3 2{\rm \pi}^6} t \int \mathrm{d} \boldsymbol{p}_1 \int \mathrm{d} \boldsymbol{p}_2 \int \frac{\mathrm{d} \boldsymbol{k}}{k^4} \left( \frac{1}{ 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot} (\boldsymbol{p}_2 - \boldsymbol{p}_1) \dfrac{U \tau}{L} +2} \right) \nonumber\\ & \quad\times \!\left( {\boldsymbol{\mathsf{I}}} - \frac{\boldsymbol{k}\boldsymbol{k}}{k^2}\right)\boldsymbol{:}\boldsymbol{p}_ 1\boldsymbol{p}_2\frac{\sin^2\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)}{(\boldsymbol{k} \boldsymbol{\cdot}\boldsymbol{p}_1)} \frac{\sin^2\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2 \right)}{(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2)} \hat{\varOmega}_2^{(1)} \left( \boldsymbol{p}_1, \boldsymbol{p}_{2}; \boldsymbol{k} \right). \end{align}

The correlation time is now $O(\tau )$ irrespective of $U\tau /L$, since correlations between swimmers do not decay in the absence of tumbling. For times smaller than $O(\tau )$, the MSD simplifies to

(4.15)\begin{align} \langle x^2(t) \rangle_{corr} & = \frac{t^2}{(\ln\kappa)^3 4{\rm \pi}^6} \int \mathrm{d} \boldsymbol{p}_1 \int \mathrm{d} \boldsymbol{p}_2 \int\frac{\mathrm{d} \boldsymbol{k}}{k^4} \left( {\boldsymbol{\mathsf{I}}} - \frac{\boldsymbol{k}\boldsymbol{k}}{k^2}\right)\boldsymbol{:}\boldsymbol{p}_1\boldsymbol{p}_2 \nonumber\\ &\quad\times \frac{\sin^2\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)}{( \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1)} \frac{\sin^2\left(\dfrac{\rm \pi}{2} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2\right)}{(\boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{p}_2)} \hat{\varOmega}_2^{(1)} \left( \boldsymbol{k}, \boldsymbol{p}_1, \boldsymbol{p}_{2} \right) . \end{align}

The coefficient of $t^2$ is now the correlated component of the velocity variance given in (3.28). On the other hand, for times much longer than the velocity correlation time, the tracer shows diffusive motion with $\langle x^2(t) \rangle \vert _{corr} = 6 D_t\vert _{corr} t$ for any finite swimmer run length ($U\tau /L$). The effective diffusivity is given by

(4.16)\begin{align} D_t\vert_{corr} & = \frac{1}{(\ln\kappa)^3 12{\rm \pi}^6} \int \mathrm{d} \boldsymbol{p}_1 \int\mathrm{d} \boldsymbol{p}_2 \int\frac{\mathrm{d} \boldsymbol{k}}{k^4} \left( \frac{1}{ 2{\rm \pi} \textrm{i} \boldsymbol{k}\boldsymbol{\cdot} (\boldsymbol{p}_2 - \boldsymbol{p}_1) \dfrac{U \tau}{L} +2} \right)\left( {\boldsymbol{\mathsf{I}}} - \dfrac{\boldsymbol{k}\boldsymbol{k}}{k^2}\right)\boldsymbol{:}\boldsymbol{p}_1\boldsymbol{p}_2 \nonumber\\ &\quad \times \frac{\sin^2\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_1\right)}{(\boldsymbol{k} \boldsymbol{\cdot}\boldsymbol{p}_1)} \frac{\sin^2\left(\dfrac{\rm \pi}{2}\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2\right)}{(\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}_2)} \hat{\varOmega}_2^{(1)} \left( \boldsymbol{k}, \boldsymbol{p}_1, \boldsymbol{p}_{2} \right) . \end{align}

The integrals in (4.14) and (4.16) are again evaluated numerically (see § 3.1 for the coordinate system used). Recall that the tracer diffusivity ($D_t\vert _{corr}$) in (4.10) is also non-dimensionalised using $U^2 \tau$ as the scale. For $U\tau /L \ll 1$, the integral in (4.16) is seen to be $O(nL^3)^2 (U\tau /L)$, and hence $D_t\vert _{corr} \sim (nL^3)^2/(\ln \kappa )^3 U^2 \tau (U\tau /L)$. However, for $U\tau /L \gg 1$ the integral in (4.16) is O$(nL^3)^2$, and hence $D_t\vert _{corr} \sim (nL^3)^2/(\ln \kappa )^3 U^2 \tau$. The scaling was earlier obtained in § 2 using physical arguments.

4.3. Effective tracer MSD

The total MSD of the tracer, to ${O}(nL^3)^2$, is obtained by combining the expressions in (4.8) and (4.14), and is plotted in figure 11, as a function of $nL^3$, for a suspension of pushers. 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 behaviour 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 persistent 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 crossover 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.

Figure 11.

The evolution of the MSD of the passive tracer, in a suspension of pushers, with time for $U\tau /L =50$ and $\kappa = 8$ for varying $nL^3$. The inset shows the correlated tracer diffusivity ($D_t\vert _{corr}$) at ${O}(nL^3)^2$ as a function of the run length, $U\tau /L$, for $\kappa =8$; the dashed lines show the asymptotic scalings discussed in the text.

The inset of figure 11 shows the ${O}(nL^3)^2$ correlated contribution to the tracer diffusivity ($D_t\vert _{corr}$). In the straight-swimmer limit, $D_t\vert _{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. Here, $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

(4.17)\begin{equation} \frac{nL^3 U^2 \tau}{(\ln \kappa)^2} \left(\widetilde{d_1} + \frac{nL^3}{\ln \kappa} \frac{U\tau}{L} \widetilde{d_2}\right) \quad\mathrm{for}\ U\tau/L \ll 1 \end{equation}

and

(4.18)\begin{equation} \frac{nL^4U}{(\ln \kappa)^2} \left(d_1 + \frac{nL^3}{\ln \kappa} \frac{U\tau}{L} d_2 \right) \quad\mathrm{for}\ U\tau/L \gg 1, \end{equation}

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. Importantly, we see that correlations significantly enhance the tracer diffusivity in a pusher suspension; figure 12 already shows an enhancement of $55.51 \%$ at $nL^3 = 0.1$ for a swimmer with the run length, $U\tau /L = 50$. Thus, the enhancement for the tracer diffusivity due to correlations is much more pronounced than that for the velocity variance discussed earlier in § 3.3.

Figure 12.

The variation of the diffusivity of a passive tracer versus the volume fraction of swimmers ($nL^3$) for a fixed run length, $U\tau /L =50$, and aspect-ratio $\kappa = 8$.

Several experiments probing bacterial suspension dynamics with passive tracers have reported an intermediate super-diffusive regime and a volume fraction dependent crossover time (Wu & Libchaber Reference Wu and Libchaber2000; Kim & Breuer Reference Kim and Breuer2004; Chen et al. Reference Chen, Lau, Hough, Islam, Goulian, Lubensky and Yodh2007; Valeriani et al. Reference Valeriani, Li, Novosel, Arlt and Marenduzzo2011; Argun et al. Reference Argun, Moradi, Pinçe, Bagci, Imparato and Volpe2016; Patteson et al. Reference Patteson, Gopinath, Purohit and Arratia2016; Peng et al. Reference Peng, Lai, Tai, Zhang, Xu and Cheng2016). 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. Reference Maggi, Paoluzzi, Pellicciotta, Lepore, Angelani and Di Leonardo2014; Argun et al. Reference Argun, Moradi, Pinçe, Bagci, Imparato and Volpe2016; Kanazawa et al. Reference Kanazawa, Sano, Cairoli and Baule2019). 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 (Argun et al. Reference Argun, Moradi, Pinçe, Bagci, Imparato and Volpe2016; Krishnamurthy et al. Reference Krishnamurthy, Ghosh, Chatterji, Ganapathy and Sood2016).