---
abstract: 'Active matter exhibits various forms of non-equilibrium states in the absence of external forcing, including macroscopic steady-state currents. Such states are often too complex to be modelled from first principles and our understanding of their physics relies heavily on minimal models. These have mostly been studied in the case of “dry" active matter, where particle dynamics are dominated by friction with their surroundings. Significantly less is known about systems with long-range hydrodynamic interactions that belong to “wet" active matter. Dilute suspensions of motile bacteria, modelled as self-propelled dipolar particles interacting solely through long-ranged hydrodynamic fields, are arguably the most studied example from this class of active systems. Their phenomenology is well-established: at sufficiently high density of bacteria, there appear large-scale vortices and jets comprising many individual organisms, forming a chaotic state commonly known as *bacterial turbulence*. As revealed by computer simulations, below the onset of collective motion, the suspension exhibits very strong correlations between individual microswimmers stemming from the long-ranged nature of dipolar fields. Here we demonstrate that this phenomenology is captured by the minimal model of microswimmers. We develop a kinetic theory that goes beyond the commonly used mean-field assumption, and explicitly takes into account such correlations. Notably, these can be computed exactly within our theory. We calculate the fluid velocity variance, spatial and temporal correlation functions, the fluid velocity spectrum, and the enhanced diffusivity of tracer particles. We find that correlations are suppressed by particle self-propulsion, although the mean-field behaviour is not restored even in the limit of very fast swimming. Our theory is not perturbative and is valid for any value of the micro-swimmer density below the onset of collective motion. This work constitutes a significant methodological advance and allows us to make qualitative and quantitative predictions that can be directly compared to experiments and computer simulations of micro-swimmer suspensions.'
author:
- Viktor Škultéty
- Cesare Nardini
- Joakim Stenhammar
- Davide Marenduzzo
- Alexander Morozov
bibliography:
- 'Refs.bib'
title: Swimming suppresses correlations in dilute suspensions of pusher microorganisms
---

Introduction {#section:introduction}
============

In recent years active systems emerged as a new state of matter with unique properties that are absent from their passive counterparts [@Ramaswamy2010; @Marchetti2013]. Such systems comprise particles that are capable of extracting energy from their environment and using it to exert forces and torques on their surroundings. The resulting self-propulsion and interactions between particles break detailed balance at the microscopic level, often leading to steady states that are not invariant under time reversal and exhibit macroscopic currents [@Cates2012]. Such currents, or collective motion, have been reported in a variety of systems [@Vicsek2012], including Vicsek particles [@Chate2008], mixtures of microtubules and molecular motors [@Sanchez2012], light-activated colloids [@Palacci2013], Quincke rollers [@Bricard2013; @Karani2019], bacterial colonies [@Zhang2010], sperm cells [@Creppy2015], locusts [@Buhl2006], birds, and fish [@Parrish1997]. The omnipresence of collective motion raises the need to classify various active systems according to common features of their phenomenological behaviour. Marchetti *et al.* [@Marchetti2013] recently introduced two broad universality classes for active systems, “dry" and “wet", comprising particles dominated by friction with their surroundings and long-ranged hydrodynamic interactions, respectively. Each class is expected to be defined by a few, relatively simple model systems, and significant effort has been invested into finding such models. For dry active matter, these include Vicsek-like models [@Vicsek2012; @Chate2020], that describe cases where alignment interactions are dominant, and Active Brownian Particles [@Howse2007; @Romanczuk2012] or Run and Tumble particles [@Schnitzer1993], that describe systems dominated by steric forces randomising their self-propulsion direction either smoothly or in a discontinuous manner. In this work, we study dilute suspensions of motile bacteria that, arguably, play the same role for wet active matter [@Koch2011; @Saintillan2013].

Collective motion in bacteria has been extensively studied in dilute [@Soni2003; @Dombrowski2004; @Gachelin2014] and dense [@Mendelson1999; @Wu2006; @Sokolov2007; @Sokolov2009; @Cisneros2011; @Sokolov2012; @Wensink2012a; @Dunkel2013] suspensions. These studies reveal the following sequence of dynamical states. At very low densities, bacterial suspensions appear featureless and disordered [@Wu2006; @Gachelin2014]. At higher, yet still sufficiently low densities, collective motion sets in on the scale of the system. In this state, bacterial motion takes the form of large-scale jets and vortices with typical speeds that are larger than the swimming speeds of individual organisms [@Soni2003; @Dombrowski2004; @Gachelin2014]. At significantly higher densities, there emerges a typical lengthscale of the vortices, which is comparable to about $5-10$ times the bacterial size [@Sokolov2012; @Dunkel2013; @Ryan2013]. Although this sequence of dynamical states has never been simultaneously observed in a single systematic bulk experiment, with the exception of Sokolov *et al.* [@Sokolov2009], the transition scenario is supported by computer simulations of self-propelled particles interacting through various forms of long-ranged hydrodynamic fields and short-ranged steric repulsion [@Hernandez-Ortiz2005; @Saintillan2007; @Wolgemuth2008; @Underhill2008; @Hernandez-Ortiz2009; @Lushi2013; @Lushi2014; @Krishnamurthy2015; @Wioland2016; @Saintillan2012; @Stenhammar2017; @Theillard2017; @Schwarzendahl2018; @Bardfalvy2019; @Theillard2019].

Bulk experiments with *E.coli* [@Gachelin2014] and *B.subtilis* [@Dombrowski2004] show that the transition to collective motion occurs around a volume fraction of bacterial bodies of about $1-2\%$. At such densities, the typical distance between the organisms is about $5-8$ times their body length, collisions are rare, and the far-field hydrodynamic interactions are thought to be dominant [@Koch2011; @Saintillan2013]. The latter are well-described by a “pusher”-like Stokesian dipolar field [@Lauga2009; @Drescher2011], generated when two point forces of equal magnitude and pointing away from each other are applied to a viscous fluid. Self-propelled pusher-like dipolar particles thus form a minimal model for dilute bacterial suspensions.

The transition to collective motion in dilute bacterial suspensions can be understood in terms of a mean-field kinetic theory [@Koch2011; @Saintillan2013] incorporating the minimal ingredients discussed above. Such theory identifies re-orientation of bacteria in the velocity field created by other organisms as the key ingredient leading to a global isotropic-nematic transition. The globally ordered state is, however, linearly unstable through a long-wavelength generic instability [@AditiSimha2002; @Marchetti2013], and there ensue never-settling dynamics as a compromise between the two instabilities. The critical density of bacteria at the onset of collective motion is determined by the strength of their dipolar interactions, their shape, and the way individual organisms change their orientation: either by occasionally re-orienting in a random way (tumble), or by rotational diffusion [@Saintillan2008; @Saintillan2008a; @Subramanian2009; @Hohenegger2010; @Krishnamurthy2015]. Typically, the critical threshold density is significantly lower in the latter case, and going to zero in the absence of a decorrelation mechanism for individual bacterium orientation. The mean-field kinetic theory has also been extended to systems with steric interactions [@Ezhilan2013; @Ryan2013; @Heidenreich2016; @Reinken2018] and to microswimmers suspended in non-Newtonian fluids [@Bozorgi2013; @Bozorgi2014; @Li2016].

Below the onset of collective motion, the mean-field kinetic theory predicts that the suspension is homogeneous and isotropic, as featureless as a suspension of non-interacting microswimmers. These assumptions are widely used when describing rheological properties of very dilute suspensions [@Hatwalne2004; @Chen2007; @Sokolov2009a; @Saintillan2010; @Underhill2011; @Lopez2015; @Alonso-Matilla2016; @Nambiar2017; @Guo2018; @Nambiar2019; @Liu2019; @Saintillan2018] and enhanced diffusivity of tracer particles [@Wu2000; @Kim2004; @Underhill2008; @Leptos2009; @Dunkel2010; @Ishikawa2010; @Childress2010; @Childress2011; @Kurtuldu2011; @Mino2011; @Mino2013; @Jepson2013; @Pushkin2013; @Pushkin2013jfm; @Morozov2014; @Kasyap2014; @Thiffeault2015; @Patteson2016; @Burkholder2017]. However, recent large-scale Lattice-Boltzmann simulations of dipolar swimmers [@Stenhammar2017; @Bardfalvy2019] revealed the presence of very strong correlations below the onset of collective motion. It was shown that various observables deviate from their mean-field values at any density of microswimmers [@Stenhammar2017], with the deviation diverging in the vicinity of the onset. The origin of such strong correlations can be readily attributed to the slow spatial decay of the dipolar velocity field, implying a simultaneous coupling between all microswimmers in the system. While this argument is intuitive enough, its implementation as a theoretical framework presents major technical challenges [@Stenhammar2017; @Qian2017; @Nambiar2019a], and only simplified cases were studied until now, notably by Stenhammar *et al.* [@Stenhammar2017], who considered a suspension of “shakers" – particles that apply forces to the fluid but do not self-propel.

In this work we develop a kinetic theory that goes beyond the mean-field assumption for the general model of dilute microswimmer suspensions described above. Our theory explicitly includes particle self-propulsion, and is valid at any density of microswimmers below the onset of collective motion. This constitutes simultaneously a significant methodological development compared to the work by Stenhammar *et al.* [@Stenhammar2017], and a major advance in our understanding of one of the key models defining “wet" active matter. Our theory allows us to make explicit predictions for observables that can be directly set against experiments and numerical simulations.

The paper is organised as follows. In Section \[section:theory\] we formulate a kinetic theory for a model suspension of pusher-like dipolar microswimmers. We explicitly find the dynamics of fluctuations around the homogeneous and isotropic state that describe the system below the onset of collective motion. Since our theory differs significantly from the previous work [@Stenhammar2017], we present its derivation in detail. We appreciate, however, that some readers might only be interested in the results of our theory without feeling the need to go through the rather technical Section \[section:theory\]. We, therefore, present our results in a stand-alone Section \[section:results\], which can be read without Section \[section:theory\]. There, we calculate the temporal and spatial correlation functions, fluid velocity variance, energy spectra, and the enhanced diffusivity of tracer particles. We conclude in Section \[section:discussion\], while Appendices contain additional derivations for technically oriented readers.

Kinetic theory of strongly interacting suspensions {#section:theory}
==================================================

Microscopic model {#subsection:microscopic}
-----------------

We consider a collection of $N$ microswimmers contained in a volume $V$ at a finite number density $n=N/V$. The microswimmers are suspended in a Newtonian fluid with the viscosity $\mu$. Each microswimmer is described by its instantaneous position ${\bm x}_i$ and orientation ${\bm p}_i$, that we collectively denote by ${\bm z}_i = \left( {\bm x}_i, {\bm p}_i\right)$, where $i=1\dots N$ enumerates the particles. Within our model, the dynamics of the suspension is governed by the following equations of motion $$\begin{aligned}
\dot{x}_{i}^{\alpha} &= v_{s} p_{i}^{\alpha} + \mathcal{U}^{\alpha}\left( {\bm x}_i \right),   
\label{xdot} \\
\dot{p}_{i}^{\alpha} &= \mathbb{P}_{i}^{\alpha\beta} \left(\mathcal{W}^{\beta\gamma}\left( {\bm x}_i \right) + B \mathcal{E}^{\beta\gamma}\left( {\bm x}_i \right) \right) p_{i}^{\gamma},
\label{pdot}\end{aligned}$$ where the dot denotes the time derivative, the superscript indices denote Cartesian components of vectors, and the subscript indices label the particles. Throughout this work, we utilise the Einstein summation convention for the superscript indices, while no summation is assumed over repeated subscript indices.

Equations of motion and incorporate the following physical ingredients. First of all, each swimmer self-propels with the speed $v_s$ in the direction of its orientation. To induce self-propulsion, swimmers generate long-ranged flows in the suspending fluid [@Lauga2009]. The superposition of these flows at the position of the $i$-th swimmer, $\mathcal{U}^{\alpha}\left( {\bm x}_i \right)$, advects that particle in addition to its self-propulsion, see Eq., and re-orients it according to Jeffrey’s equation . The latter describes the dynamics of a passive particle in an external flow [@kimkarrila], with $$\begin{aligned}
&\mathcal{W}^{\beta\gamma}\left( {\bm x}_i \right) = \frac{1}{2} \left( \nabla^{\gamma}  \mathcal{U}^{\beta}\left( {\bm x}_i \right) - \nabla^{\beta}  \mathcal{U}^{\gamma}\left( {\bm x}_i \right) \right), \\
&\mathcal{E}^{\beta\gamma} \left( {\bm x}_i \right)= \frac{1}{2} \left(\nabla^{\gamma}  \mathcal{U}^{\beta}\left( {\bm x}_i \right) + \nabla^{\beta}  \mathcal{U}^{\gamma}\left( {\bm x}_i \right) \right),
\label{eq:Strain}\end{aligned}$$ being the Cartesian components of the vorticity and rate-of-strain tensors, respectively. In Eq., $\mathbb{P}_{i}^{\alpha\beta} = \delta^{\alpha\beta} - p^{\alpha}_{i}p_{i}^{\beta}$, is the projection operator, $\delta^{\alpha\beta}$ denotes the Kronecker delta, $\nabla_i^\alpha = \partial/\partial x_i^\alpha$, and $B = \left( a^2 - 1\right)/\left( a^2 + 1\right)$ is the measure of the swimmer’s nonsphericity [@kimkarrila] based on its aspect ratio $a$. For strongly elongated particles, $B\rightarrow 1$, while for spheres, $B=0$. Finally, each swimmer randomly changes its orientation with a rate $\lambda$, thus mimicking the run-and-tumble motion commonly exhibited by bacteria [@Berg1993]. We note here that we neglect the effects of rotational and translational diffusion on the particle’s dynamics, and random tumbling is thus the only source of stochasticity in our model.

The velocity field generated by a self-propelled particle sufficiently far away from its surface is often well-described by the field produced by a point dipole with the same position and orientation [@Lauga2009; @Drescher2011]. In a dilute suspension of microswimmers, where the particles are sufficiently separated from each other, we can approximate $\mathcal{U}^{\alpha}\left( {\bm x}_i \right)$ by a sum of dipolar contributions $$\begin{aligned}
 \mathcal{U}^{\alpha}\left( {\bm x}_i \right) = \sum_{j\neq i}^{N} u_{d}^{\alpha}({\bm x}_{i}; {\bm z}_{j}),\end{aligned}$$ where $$\begin{aligned}
{\bm u}_{d}({\bm x}_{i}; {\bm z}_{j}) = \frac{\kappa}{8\pi} \left[ 3
\frac{\left( {\bm p}_j\cdot {\bm x}'\right)^2 {\bm x}' + \epsilon^2 \left( {\bm p}_j\cdot {\bm x}'\right) {\bm p}_j }{\left( x'^2 + \epsilon^2 \right)^{5/2}} \right. \nonumber \\
\left. -\frac{{\bm x}'}{\left( x'^2 + \epsilon^2 \right)^{3/2}} \right]
\label{realspacedipole}\end{aligned}$$ is the velocity field generated at ${\bm x}_{i}$ by a hydrodynamic dipole located at ${\bm x}_{j}$ with the orientation ${\bm p}_{j}$. Here, $\kappa = F l/\mu$ is the reduced dipolar strength, where $F$ is the magnitude of the forces applied to the fluid, $l$ is the dipolar length, and $\mu$ is the viscosity of the fluid; ${\bm x}' = {\bm x}_i -{\bm x}_j$, and $x'$ denotes the length of ${\bm x}'$. The dipole consists of two regularised Stokeslets, that were introduced by Cortez *et al.* [@Cortez2005], with $\epsilon$ being the regularisation length of the order of swimmer size. For pushers, $\kappa > 0$.

The main goal of our work is to calculate spatial and temporal correlations of the fluid velocity in microswimmer suspensions described by the model above. Both quantities can be succinctly expressed through a combined correlation function $$\begin{aligned}
&C(R,T) \nonumber \\
&\qquad = \lim_{t\rightarrow\infty} \frac{1}{V}\int d{\bm x} \,\overline{U^{\alpha}\left( {\bm x}, t \right) U^{\alpha}\left( {\bm x} + {\bm R}, t+T \right)},
\label{CRTgeneral}\end{aligned}$$ where $U^{\alpha}\left( {\bm x}, t \right)$ is the fluid velocity at the position ${\bm x}$ at time $t$, and the large-$t$ limit guarantees independence of the initial conditions. The spatial and temporal correlation functions are trivially recovered by setting $T=0$ and $R=0$, respectively. The bar in Eq. denotes the average over the history of tumble events, and reflect the stochastic nature of our model. To calculate this and similar averages, below we formulate a kinetic theory of microswimmer suspensions based on our macroscopic model. Such theories have been extensively studied at the mean-field level [@Liao2007; @Lau2009; @Zaid2011; @Underhill2011; @Belan2019; @Bardfalvy2019]. Here, we go beyond the mean-field approximation and explicitly take into account strong correlations between the swimmers caused by the long-range nature of their hydrodynamic fields, Eq..

Kinetic theory and BBGKY hierarchy
----------------------------------

The starting point of our theory is the $N$-particle probability distribution function $F_N\left({\bm z}_1, {\bm z}_2, \dots, {\bm z}_N, t \right)$ that gives the geometric probability of the system occupying a particular point in the $6N$-dimensional phase space $\left\{ {\bm z}_1, \dots, {\bm z}_N\right\} $ at time $t$. The $N$-particle probability distribution function is symmetric with respect to swapping particle labels, reflecting their indistinguishability, and is normalised $$\begin{aligned}
\int d{\bm z}_1 \dots d{\bm z}_N F_N\left({\bm z}_1, \dots, {\bm z}_N, t \right) = 1.
\label{normalisation}\end{aligned}$$ Its time dynamics is governed by the Master equation [@Balescu1975] $$\begin{aligned}
\partial_{t} F_{N} & + \sum_{i=1}^{N} \Big[ \nabla_{i}^{\alpha}(\dot{x}_{i}^{\alpha} F_{N}) + \partial_{i}^{\alpha}(\dot{p}_{i}^{\alpha} F_{N}) \Big] \nonumber \\
& = - N \lambda F_{N} +  \frac{\lambda}{4\pi} \sum_{i=1}^{N} \int \mathrm{d} {\bm p}_{i} F_{N},
\label{Liouville}\end{aligned}$$ where we introduced $\partial_{i}^{\alpha} = \mathbb{P}_{i}^{\alpha\beta} \partial/\partial p^\beta_i$. The l.h.s. of Eq. describes the probability fluxes to and from a particular point in the phase space due to the deterministic particle dynamics given by Eqs. and , while the r.h.s. gives the changes of the probability due to random tumbling from and into that phase space point [@Subramanian2009; @Koch2011]. Next, we introduce the $s$-particle correlation functions defined as $$\begin{aligned}
&F_{s}\left({\bm z}_1,\dots, {\bm z}_s, t \right) = \frac{N!}{(N-s)!N^{s}} \nonumber \\
&\qquad\qquad \times \int d{\bm z}_{s+1} \dots d{\bm z}_N F_N\left({\bm z}_1, \dots, {\bm z}_N, t \right), 
\label{Fs}\end{aligned}$$ Below, we will only be interested in the first partial correlation functions $F_1$, $F_2$, and $F_3$, that we further express as $$\begin{aligned}
F_2\left( {\bm z}_1, {\bm z}_2, t\right) = F_1\left( {\bm z}_1, t\right) F_1\left( {\bm z}_2, t \right) +  G\left( {\bm z}_1, {\bm z}_2, t\right),
\label{F2}\end{aligned}$$ and $$\begin{aligned}
F_3\left( {\bm z}_1, {\bm z}_2, {\bm z}_3, t\right) = F_1\left( {\bm z}_1, t\right) F_1\left( {\bm z}_2, t \right) F_1\left( {\bm z}_3, t \right) \nonumber \\
+ G\left( {\bm z}_1, {\bm z}_2, t\right) F_1\left( {\bm z}_3, t \right) + G\left( {\bm z}_1, {\bm z}_3, t\right) F_1\left( {\bm z}_2, t \right) \nonumber \\
+ G\left( {\bm z}_2, {\bm z}_3, t\right) F_1\left( {\bm z}_1, t \right) + H\left( {\bm z}_1, {\bm z}_2, {\bm z}_3, t\right),
\label{F3}\end{aligned}$$ where $G$ and $H$ are the irreducible (connected) correlation functions [@Balescu1975]. The time evolution of $F_s$ can be deduced from the Master equation by integrating it over $\left\{{\bm z}_{s+1},\dots,{\bm z}_N\right\}$. Integrating by parts and using Eqs. and , we obtain the following equations for the one- and two-particle irreducible correlation functions

$$\begin{aligned}
& \partial_t F_1({\bm z}, t) + {\mathcal L}[F_1({\bm z}, t)]({\bm z}) \nonumber \\
& \qquad\qquad = -N \nabla^\alpha \int d{\bm z}' G({\bm z},{\bm z}',t) u^\alpha_d ({\bm x}; {\bm z}')  
- N \mathbb{P}^{\alpha\beta} \frac{\partial}{\partial p^\beta} \int d{\bm z}' G({\bm z},{\bm z}',t)
p^\gamma \mathbb{X}^{\alpha\mu\nu\gamma} \nabla^\mu u^\nu_d ({\bm x}; {\bm z}'), 
\label{EqF1}
\\ 
\nonumber \\
& \partial_t G({\bm z}_1,{\bm z}_2, t) + {\mathcal L}[G({\bm z}_1,{\bm z}_2, t)]({\bm z}_1) + {\mathcal L}[G({\bm z}_1,{\bm z}_2, t)]({\bm z}_2) \nonumber \\
& \qquad \qquad + N \nabla_1^\alpha \left[ F_1({\bm z}_1,t) \int d{\bm z}' G({\bm z}_2,{\bm z}',t) u_d^\alpha({\bm x}_1;{\bm z}')\right]
+ N \nabla_2^\alpha \left[ F_1({\bm z}_2,t) \int d{\bm z}' G({\bm z}_1,{\bm z}',t) u_d^\alpha({\bm x}_2;{\bm z}')\right]
\nonumber \\ 
& \qquad \qquad + N \mathbb{P}_1^{\alpha\beta} \frac{\partial}{\partial p_1^\beta}\left[  F_1({\bm z}_1,t) p_1^\gamma \mathbb{X}_1^{\alpha\mu\nu\gamma} \int d{\bm z}' G({\bm z}_2,{\bm z}',t)\nabla_1^\mu u_d^\nu({\bm x}_1;{\bm z}')\right]
\nonumber \\ 
& \qquad \qquad + N \mathbb{P}_2^{\alpha\beta} \frac{\partial}{\partial p_2^\beta}\left[  F_1({\bm z}_2,t) p_2^\gamma \mathbb{X}_2^{\alpha\mu\nu\gamma} \int d{\bm z}' G({\bm z}_1,{\bm z}',t)\nabla_2^\mu u_d^\nu({\bm x}_2;{\bm z}')\right] \nonumber \\
& \qquad \qquad =  - \mathcal{S}_{1,2}^{F} - \mathcal{S}_{2,1}^{F} -  \mathcal{S}_{1,2}^{G} - \mathcal{S}_{2,1}^{G} - \mathcal{S}_{1,2}^{H} - \mathcal{S}_{2,1}^{H},
\label{EqG}\end{aligned}$$

where we have introduced the operator $$\begin{aligned}
&{\mathcal L}[\Phi]({\bm z}) \nonumber \\
& \qquad =
v_s p^\alpha \nabla^\alpha \Phi({\bm z}) 
+ N \nabla^\alpha \big[ \Phi({\bm z}) \mathcal{U}^\alpha_{\text{MF}}({\bm x})\big] 
+ N \mathbb{P}^{\alpha\beta} \frac{\partial}{\partial p^\beta}\big[ \Phi({\bm z}) p^\gamma \mathbb{X}^{\alpha\mu\nu\gamma} \nabla^\mu \mathcal{U}^\nu_{\text{MF}}({\bm x})\big] 
+ \lambda \Phi({\bm z})
- \frac{\lambda}{4\pi} \int d{\bm p}\, \Phi({\bm z}),
\label{operatorL}\end{aligned}$$

acting on the variable $\bm z$ of an arbitrary function $\Phi=\Phi({\bm z}_1,\dots,{\bm z}_N)$, and defined the mean-field velocity field as $$\begin{aligned}
 \mathcal{U}^\alpha_{\text{MF}}({\bm x}) = \int d{\bm z}' F_1({\bm z}',t) u^\alpha_d({\bm x}; {\bm z}').
 \label{UMF}\end{aligned}$$ The rank-4 tensor $$\begin{aligned}
\mathbb{X}_i^{\alpha\mu\nu\gamma} = \mathbb{P}_i^{\alpha\beta} \left[ \frac{B+1}{2}\delta^{\mu\gamma}\delta^{\nu\beta} + \frac{B-1}{2}\delta^{\mu\beta}\delta^{\nu\gamma} \right],\end{aligned}$$ encodes the tensorial structure of Jeffrey’s equation , and the r.h.s. of Eq. is given in terms of $$\begin{aligned}
 &\mathcal{S}_{i,j}^{F} = 
F_1({\bm z}_j,t)\Big\{ \nabla_i^\alpha \left[ F_1({\bm z}_i,t) u_d^\alpha({\bm x}_i;{\bm z}_j)\right] \nonumber \\
&+  \mathbb{P}_i^{\alpha\beta} \frac{\partial}{\partial p_i^\beta}\left[ F_1({\bm z}_i, t) p_i^\gamma \mathbb{X}_i^{\alpha\mu\nu\gamma} \nabla_i^\mu u_d^\nu({\bm x}_i; {\bm z}_j) \right] \Big\},\end{aligned}$$ $$\begin{aligned}
 & \mathcal{S}_{i,j}^{G} = 
\nabla_i^\alpha \left[ G({\bm z}_i,{\bm z}_j,t) u_d^\alpha({\bm x}_i;{\bm z}_j)\right] \nonumber \\
&+  \mathbb{P}_i^{\alpha\beta} \frac{\partial}{\partial p_i^\beta}\left[ G({\bm z}_i,{\bm z}_j,t) p_i^\gamma \mathbb{X}_i^{\alpha\mu\nu\gamma} \nabla_i^\mu u_d^\nu({\bm x}_i; {\bm z}_j) \right], \end{aligned}$$ and $$\begin{aligned}
 & \mathcal{S}_{i,j}^{H} = 
N \int d{\bm z}' \Big\{
\nabla_i^\alpha \left[ H({\bm z}_i,{\bm z}_j,{\bm z}',t) u_d^\alpha({\bm x}_i;{\bm z}')\right] \nonumber \\
&+  \mathbb{P}_i^{\alpha\beta} \frac{\partial}{\partial p_i^\beta}\left[ H({\bm z}_i,{\bm z}_j,{\bm z}',t) p_i^\gamma \mathbb{X}_i^{\alpha\mu\nu\gamma} \nabla_i^\mu u_d^\nu({\bm x}_i; {\bm z}') \right] \Big\}.\end{aligned}$$ Eqs. and are the beginning of a BBGKY hierarchy of equations for partial distribution functions [@Balescu1975]. As such, they do not form a closed system as they also depend on the three-particle irreducible distribution function $H$.

Before discussing our choice of closure for this system of equations, let us briefly review the predictions of the mean-field approximation to Eqs. and , which consists of neglecting all correlation functions beyond $s=1$. The remaining equation determines the mean-field approximation to the one-particle correlation function $$\begin{aligned}
\partial_t F_1^{\text{MF}}({\bm z}, t) + {\mathcal L}[F_1^{\text{MF}}({\bm z}, t)]({\bm z})  = 0,
\label{EqF1MF}\end{aligned}$$ that has been extensively studied before [@Saintillan2008; @Saintillan2008a; @Subramanian2009; @Hohenegger2010; @Koch2011; @Saintillan2013; @Krishnamurthy2015]. One of the solutions of this equation is given by a constant, which is fixed to $F_1^{\text{MF}}({\bm z}, t)=1/(4\pi V)$ by the normalisation condition Eq.. This solution, which is valid at any number density, corresponds to a homogeneous and isotropic suspension of microswimmers. For pushers $(\kappa>0)$, this state loses its stability [@Saintillan2008; @Saintillan2008a; @Subramanian2009; @Hohenegger2010; @Stenhammar2017] at the critical number density of microswimmers $n_{crit} = 5\lambda/(B\kappa)$, while for pullers $(\kappa<0)$, the homogeneous and isotropic state is always linearly stable within the mean-field approximation.

The homogeneous and isotropic mean-field solution implies that $N F_1^{\text{MF}}\sim n \sim O(1)$ is finite in the thermodynamic limit. This, in turn, implies that, to leading order, $G\sim O(N^{-2})$, $H\sim O(N^{-3})$, etc. A more comprehensive discussion of this statement, together with the required rescaling of the correlation functions, system parameters, and time is given elsewhere [@Stenhammar2017].

Building upon these results, here we *assume* that upon approaching the thermodynamic limit, $F_1$ is well-approximated by $F_1^{\text{MF}}$, since the r.h.s. of Eq. is $O(1/N)$ compared its l.h.s. In the homogeneous and isotropic state, the mean-field velocity vanishes $\mathcal{U}^\alpha_{\text{MF}}({\bm x}) = 0$, since the integral in Eq. is then proportional to the total flow rate through a surface surrounding the dipole. The latter is zero due to incompressibility. Fluctuations around the homogeneous and isotropic state are then governed by Eq. with $F_1=1/(4\pi V)$, and $\mathcal{S}_{i,j}^{G}$ and $\mathcal{S}_{i,j}^{H}$ neglected $$\begin{aligned}
& \partial_t G({\bm z}_1,{\bm z}_2, t) + {\mathcal L}_{12} [G] + {\mathcal L}_{21} [G] \nonumber \\
& \qquad\qquad = \frac{3B}{\left(4\pi V\right)^2} \Big\{ p_1^\mu p_1^\nu \nabla_1^\mu u_d^\nu({\bm x}_1;{\bm z}_2) \nonumber \\
& \qquad\qquad\qquad\qquad\qquad+ p_2^\mu p_2^\nu \nabla_2^\mu u_d^\nu({\bm x}_2;{\bm z}_1) \Big\},
\label{EqGfinal}\end{aligned}$$ where $$\begin{aligned}
& {\mathcal L}_{ij} [G] = 
v_s p_i^\alpha \nabla_i^\alpha G({\bm z}_1,{\bm z}_2, t) \nonumber \\
&\qquad\qquad - \frac{3 n B}{4\pi} p_i^\mu p_i^\nu \int d{\bm z}' G({\bm z}_j,{\bm z}', t) \nabla_i^\mu u_d^\nu({\bm x}_i;{\bm z}') \nonumber \\
& \qquad\qquad + \lambda G({\bm z}_1,{\bm z}_2, t)
- \frac{\lambda}{4\pi} \int d{\bm p}_i\, G({\bm z}_1,{\bm z}_2, t).
\label{Lij}\end{aligned}$$ This equation has previously been derived and analysed for the case of shakers ($v_s=0$) [@Stenhammar2017]. Here, we solve it in the general case $v_s>0$.

Phase-space density fluctuations
--------------------------------

While the two-point distribution function $G$, given by Eq., contains statistical information about fluctuations in the system, it is not straightforward to relate it to the spatial and temporal correlation function $C(R,T)$, Eq., that we seek to calculate. To establish this connection, we introduce a method based on the phase space density $$\begin{aligned}
\varphi({\bm z},t) = \sum_{i=1}^N \delta({\bm z} - {\bm z}_i(t)),\end{aligned}$$ pioneered by Klimontovich [@klimontovich1967book]. Here, $\delta({\bm z})$ is the three-dimensional Dirac delta function. The average of the phase space density is related to $F_1$ as can be seen from $$\begin{aligned}
&\overline\varphi({\bm z},t) = \int d{\bm z}_1 \dots d{\bm z}_N  \sum_{i=1}^N \delta({\bm z} - {\bm z}_i) F_N({\bm z}_1, \dots, {\bm z}_N, t) \nonumber \\
& \qquad\qquad = N F_1({\bm z}, t),\end{aligned}$$ where we used Eq.. Fluctuations of the phase space density can formally be defined as $\delta \varphi = \varphi - \overline\varphi$, and their second moment is given by $$\begin{aligned}
& G_K( {\bm z}', {\bm z}'', t) \equiv \overline{\delta \varphi({\bm z}',t) \delta \varphi({\bm z}'',t)} \nonumber \\
& \qquad = N^2 G ( {\bm z}', {\bm z}'', t) + N F_1({\bm z}',t) \delta({\bm z}' - {\bm z}'').
\label{GK}\end{aligned}$$ Below, we refer to $G_K$ as the Klimontovich correlation function. Its utility is evident if one considers the spatial correlation function $C(R)$, defined in Eq. as $$\begin{aligned}
&C(R) = \lim_{t\rightarrow\infty} \frac{1}{V}\int d{\bm x} \overline{U^{\alpha}\left( {\bm x}, t \right) U^{\alpha}\left( {\bm x} + {\bm R}, t \right)}.\end{aligned}$$ The velocity of the fluid at a position $\bm x$ is given by the superposition of the velocity fields generated by all swimmers $$\begin{aligned}
&U^{\alpha}\left( {\bm x}, t \right) = \sum_{i=1}^N u_d^\alpha({\bm x}; {\bm z}_i(t)) \nonumber \\
&\qquad\qquad = \int d{\bm z}' \varphi({\bm z}',t) u_d^\alpha({\bm x}; {\bm z}'). \end{aligned}$$ Separating the phase space density into its average and fluctuations, $\varphi =\overline\varphi +  \delta \varphi$, the spatial correlation function becomes $$\begin{aligned}
&C(R) = \lim_{t\rightarrow\infty} \frac{1}{V}\int d{\bm x} \int d{\bm z}'d{\bm z}''  u_d^\alpha({\bm x}; {\bm z}') u_d^\alpha({\bm x}+{\bm R}; {\bm z}'')  \nonumber \\
& \qquad\qquad\qquad \times  \Bigg[ \left(\frac{n}{4\pi}\right)^2 + G_K( {\bm z}', {\bm z}'', t)\Bigg].\end{aligned}$$ The integral with the constant term vanishes, demonstrating that $G_K$ fully determines the spatial correlation function.

Time evolution of the Klimontovich correlation function can readily be derived from Eqs. and , yielding $$\begin{aligned}
& \partial_t G_K({\bm z}_1,{\bm z}_2, t) + {\mathcal L}_{12} [G_K] + {\mathcal L}_{21} [G_K] \nonumber \\
& \qquad\qquad = 
2 \lambda \frac{n}{4\pi} \delta({\bm x}_1 - {\bm x}_2 ) \left[ \delta({\bm p}_1 - {\bm p}_2 ) - \frac{1}{4\pi} \right],
\label{EqGK}\end{aligned}$$ where $ {\mathcal L}_{ij}$ is defined in Eq., and we used $F_1 = 1/(4\pi V)$ in the homogeneous and isotropic state. To solve Eq., we introduce an auxiliary field $h({\bm z}_1, t)$, that satisfies the following equation $$\begin{aligned}
 \partial_t h({\bm z}_1, t) + {\mathcal L}_{11} [h] = \xi({\bm z}_1, t),
\label{Eqh}\end{aligned}$$ where $\xi$ is a noise term with the following properties $$\begin{aligned}
& \langle \xi({\bm z}_1, t) \rangle = 0, 
\label{noiseaverage}\\
& \langle \xi({\bm z}_1, t) \xi({\bm z}_2, t')\rangle = 2 \lambda \frac{n}{4\pi} \delta(t-t') \delta({\bm x}_1 - {\bm x}_2 ) 
\nonumber \\
&\qquad\qquad \times \left[ \delta({\bm p}_1 - {\bm p}_2 ) - \frac{1}{4\pi} \right].
\label{noisevariance}\end{aligned}$$ Here, the angular brackets denote the average over the realisations of the noise $\xi$, and should not be confused with the ensemble averages that we denoted by bars in the equations above. Eq. allows us to factorise the Klimontovich correlation function as $$\begin{aligned}
G_K({\bm z}_1,{\bm z}_2, t) = \langle h({\bm z}_1, t) h({\bm z}_2, t)\rangle,\end{aligned}$$ which replaces the deterministic Eq. by a significantly simpler stochastic Eq. with a fictitious noise $\xi$ with properly chosen spectral properties. Remarkably, the non-equal time correlations of the phase space density can be expressed through the same auxiliary field $$\begin{aligned}
\overline{\delta \varphi({\bm z}',t') \delta \varphi({\bm z}'',t'')} = \langle h({\bm z}', t') h({\bm z}'', t'')\rangle,\end{aligned}$$ as implied by a seminal work of Klimontovich and Silin [@Silin1962]. This, finally, leads to a direct relationship between the field $h$, which encodes the statistical properties of fluctuations in the suspension, and the combined correlation function $$\begin{aligned}
&C(R,T) = \lim_{t\rightarrow\infty} \frac{1}{V}\int d{\bm x} \int d{\bm z}'d{\bm z}''  \nonumber \\
&  \times u_d^\alpha({\bm x}; {\bm z}') u_d^\alpha({\bm x}+{\bm R}; {\bm z}'')  \langle h({\bm z}', t) h({\bm z}'', t+T) \rangle.
\label{CRTtmp}\end{aligned}$$

Dynamics of the auxiliary field $h$ {#subsection:h}
-----------------------------------

Here, we explicitly find the solution to Eq. together with Eqs. and . Since Eq. is linear in $h$, we introduce the Fourier $$\begin{aligned}
h({\bm z}, t) = \frac{1}{(2\pi)^3} \int d{\bm k} e^{i {\bm k}\cdot{\bm x}} \hat{h}({\bm k}, {\bm p}, t),
\label{FT}\end{aligned}$$ and the Laplace transforms $$\begin{aligned}
\hat{h}({\bm k}, {\bm p}, s) = \int_0^\infty dt e^{-s t} \hat{h}({\bm k}, {\bm p}, t).\end{aligned}$$ We will also require the Fourier transform of the regularised dipolar field, Eq., which is given by $$\begin{aligned}
& u_d^\nu({\bm x}; {\bm z}') = \frac{-i \kappa}{(2\pi)^3}  \int d{\bm k} \, e^{i {\bm k}\cdot ({\bm x}-{\bm x}')}  \nonumber \\
& \qquad\qquad\qquad \times  \frac{A(k \epsilon)}{k} (\hat{\bm k} \cdot {\bm p}') \left( \delta^{\nu\delta} - \hat{k}^\nu \hat{k}^\delta \right) p'^{\delta},
\label{udfourier}\end{aligned}$$ where $\hat{\bm k}={\bm k}/k$, and $k=|{\bm k}|$. The function $A$, defined as $$\begin{aligned}
A(x)=\frac{1}{2} x^2 K_2(x),\end{aligned}$$ with $K_2(x)$ being the modified Bessel function of the second kind, is close to unity for $x<1$, and quickly approaches zero for $x>1$. It will serve as a regularisation of the integrals over $k$, suppressing contributions from lengthscales smaller than the size of individual microswimmers.

Performing the Fourier and Laplace transforms of Eq., we obtain after re-arranging $$\begin{aligned}
&\hat{h}({\bm k}, {\bm p}, s) = \frac{1}{\sigma({\bm k},{\bm p}, s)}
\Bigg[\hat{h}_0({\bm k}, {\bm p}) + \hat{\xi}({\bm k}, {\bm p},s) \nonumber \\
& \qquad + \frac{\lambda}{4\pi} I^{(0)}({\bm k},s)
+ \frac{15 \lambda}{4\pi} \Delta A(k\epsilon) \bigg\{ (\hat{\bm k}\cdot {\bm p}) I^{(1)}({\bm k},{\bm p}, s) \nonumber \\
& \qquad\qquad\qquad\qquad\qquad - (\hat{\bm k}\cdot {\bm p})^2 I^{(2)}({\bm k}, s)
\bigg\} 
\Bigg].
\label{hfullsolution}\end{aligned}$$ Here, $\hat{\xi}({\bm k}, {\bm p},s)$ is the Fourier-Laplace transform of the noise, $\sigma({\bm k},{\bm p}, s) = s + \lambda + i v_s ({\bm k}\cdot{\bm p})$, and we defined $$\begin{aligned}
I^{(0)}({\bm k},s) &= \int d{\bm p } \,\hat{h}({\bm k}, {\bm p}, s), \label{I0} \\
I^{(1)}({\bm k},{\bm p}, s) &= \int d{\bm p}'  (\hat{\bm k}\cdot {\bm p}')({\bm p}\cdot {\bm p}') \hat{h}({\bm k}, {\bm p}', s), \label{I1} \\
I^{(2)}({\bm k},s) &= \int d{\bm p }\, (\hat{\bm k}\cdot {\bm p})^2 \hat{h}({\bm k}, {\bm p}, s). \label{I2}\end{aligned}$$ In Eq., $\hat{h}({\bm k}, {\bm p}, t=0) = \hat{h}_0({\bm k}, {\bm p})$ denotes some arbitrary initial condition; below we demonstrate that the long-time statistical properties of the suspension are insensitive to $\hat{h}_0({\bm k}, {\bm p})$. In Eq., we have also introduced an important dimensionless parameter $\Delta = n/n_{crit}$, where $n_{crit} = 5\lambda/(B\kappa)$ is the mean-field onset of collective motion in pusher suspensions, $\kappa>0$. For pushers, $\Delta$ measures the dimensionless distance from the onset, with $\Delta=1$ corresponding to the instability. Eq. is a linear integral equation for $\hat{h}({\bm k}, {\bm p}, s)$ and its solution is straightforward. Substituting Eq. into Eqs.-, gives $$\begin{aligned}
&I^{(0)}({\bm k},s) \nonumber \\
&\qquad\qquad= \frac{1}{1-\frac{\lambda}{4\pi}f_0} \int d{\bm p} \frac{\hat{h}_0({\bm k}, {\bm p}) + \hat{\xi}({\bm k}, {\bm p},s)}{\sigma({\bm k},{\bm p}, s)},
\label{I0sol}\end{aligned}$$ $$\begin{aligned}
&I^{(2)}({\bm k},s) = \frac{\lambda}{4\pi} f_1 I^{(0)}({\bm k},s) \nonumber \\
& \qquad\qquad + \int d{\bm p} (\hat{\bm k}\cdot {\bm p})^2\frac{\hat{h}_0({\bm k}, {\bm p}) + \hat{\xi}({\bm k}, {\bm p},s)}{\sigma({\bm k},{\bm p}, s)},
\label{I1sol}\end{aligned}$$ $$\begin{aligned}
&I^{(1)}({\bm k},{\bm p}, s) = \frac{1}{1+\frac{15 \lambda}{8\pi} \Delta A(k\epsilon) (f_2 - f_1)} \Bigg[ \nonumber \\
& \quad \int d{\bm p}' (\hat{\bm k}\cdot {\bm p}')({\bm p}\cdot {\bm p}') \frac{\hat{h}_0({\bm k}, {\bm p}') + \hat{\xi}({\bm k}, {\bm p}',s)}{\sigma({\bm k},{\bm p}', s)} \nonumber \\
& \quad + (\hat{\bm k}\cdot {\bm p}) \bigg\{ \frac{\lambda}{4\pi} f_1 I^{(0)}({\bm k},s) \nonumber \\
& \qquad\qquad\qquad + \frac{15 \lambda}{8\pi} \Delta A(k\epsilon) (f_2 - f_1) I^{(2)}({\bm k},s)
\bigg\} \Bigg],
\label{I2sol}\end{aligned}$$ where $$\begin{aligned}
f_n = 2\pi \int_{-1}^{1} dx \frac{x^{2n}}{s+\lambda + i v_s k x}.\end{aligned}$$

Having found the explicit expression for $\hat{h}({\bm k}, {\bm p},s)$, we proceed to calculate the combined correlation function, Eq.. Below, we show that only a small number of terms from Eqs. and - contribute to $C(R,T)$.

$C(R,T)$ in terms of $\hat{h}({\bm k}, {\bm p},s)$ {#subsection:hhat}
--------------------------------------------------

In what follows, it will be convenient to re-write $C(R,T)$ in terms of the Fourier and Laplace transforms of all quantities. Substituting Eq. into Eq., and using the Fourier representation of the regularised dipolar field, Eq., we obtain $$\begin{aligned}
&C(R,T) = \lim_{t\rightarrow\infty} {\mathcal L}^{-1}_{s_1,t} {\mathcal L}^{-1}_{s_2,t+T} \frac{\kappa^2}{(2\pi)^3 V} \nonumber \\
&\qquad \times \int d{\bm k} e^{-i {\bm k}\cdot{\bm R}} \frac{A^2(k\epsilon)}{k^2} 
\int d{\bm p}_1 d{\bm p}_2 
(\hat{\bm k} \cdot {\bm p}_1) (\hat{\bm k} \cdot {\bm p}_2) \nonumber \\
& \qquad\qquad\qquad \times \big( \delta^{\alpha\beta} - \hat{k}^\alpha \hat{k}^\beta \big) p_1^{\beta} 
\big( \delta^{\alpha\gamma} - \hat{k}^\alpha \hat{k}^\gamma \big) p_2^{\gamma} \nonumber \\
&  \qquad\qquad\qquad \times \langle \hat{h}({\bm k}, {\bm p}_1,s_1)\hat{h}(-{\bm k}, {\bm p}_2,s_2) \rangle_{\hat{\xi}},
\label{CRTlaplace}\end{aligned}$$ where ${\mathcal L}^{-1}_{s,t}$ formally denotes the inverse Laplace transform from $s$ to $t$, given by the Bromwich integral [@Doetsch1974]. The angular brackets $\langle\dots\rangle_{\hat{\xi}}$ denote the average with the Fourier-Laplace components of the noise $\xi$, with the following spectral properties $$\begin{aligned}
& \langle \hat{\xi}({\bm k}, {\bm p}, s) \rangle_{\hat{\xi}} = 0, 
\label{noiseaverageFL}\\
& \langle \hat{\xi}({\bm k}, {\bm p}_1, s_1) \hat{\xi}(-{\bm k}, {\bm p}_2, s_2)\rangle_{\hat{\xi}} = 2 \lambda V \frac{n}{4\pi} 
\nonumber \\
&\qquad\qquad\qquad\qquad \times \frac{1}{s_1+s_2} \left[ \delta({\bm p}_1 - {\bm p}_2 ) - \frac{1}{4\pi} \right],
\label{noisevarianceFL}\end{aligned}$$ obtained by applying the Fourier-Laplace transform to Eqs. and . While the average in Eq. can readily be formed using the solution for $\hat h$ found in Section \[subsection:h\], the result is very cumbersome. Before proceeding, we make two observations that greatly reduce the number of terms contributing to Eq..

First, we observe that $$\begin{aligned}
\int d{\bm p} \, \big( \delta^{\alpha\beta} - \hat{k}^\alpha \hat{k}^\beta \big) p^{\beta} f(\hat{\bm k} \cdot {\bm p}) = 0,\end{aligned}$$ where $f$ is an arbitrary function of $\hat{\bm k} \cdot {\bm p}$. This statement is readily demonstrated by representing $\bm p$ in spherical coordinates with $\hat{\bm k}$ selected along the $z$-axis, and performing the angular integrals component-wise. This result has profound implications for the average $\langle \hat{h}({\bm k}, {\bm p}_1,s_1)\hat{h}(-{\bm k}, {\bm p}_2,s_2) \rangle_{\hat{\xi}}$ in Eq.. Every term in $\hat{h}({\bm k}, {\bm p}_1,s_1)$, Eq., that only depends on ${\bm p}_1$ through its dependence on $(\hat{\bm k} \cdot {\bm p}_1)$ does not contribute to $C(R,T)$, as its integral over ${\bm p}_1$ with the corresponding dipolar field in Eq. vanishes. The same applies to $\hat{h}(-{\bm k}, {\bm p}_2,s_2)$.

The second observation is related to the initial condition. All terms that involve $\hat{h}_0({\bm k}, {\bm p})$ only depend on the Laplace frequency $s$ through $1/\sigma({\bm k},{\bm p}, s)$, and their inverse Laplace transform can be readily performed before any other integration. Since the inverse Laplace transform of $1/(s+a)$ is $e^{-a t}$, where $a$ is a complex number, the dominant long-time behaviour of such terms is given by $e^{-\lambda t}$, where we ignored the subdominant oscillatory dependencies. In Eq., we are interested in the $t\rightarrow$ limit, and these terms also do not contribute to $C(R,T)$.

With these observations in mind, Eq. can be significantly simplified to read $$\begin{aligned}
&\hat{h}({\bm k}, {\bm p}, s) \cong \frac{\hat{\xi}({\bm k}, {\bm p},s)}{\sigma({\bm k},{\bm p}, s)} \nonumber \\
&\qquad\qquad + \frac{ (\hat{\bm k}\cdot {\bm p}) }{\sigma({\bm k},{\bm p}, s)}  
\frac{\frac{15 \lambda}{4\pi} \Delta A(k\epsilon) }{1+\frac{15 \lambda}{8\pi} \Delta A(k\epsilon) (f_2 - f_1)} \nonumber \\
& \qquad\qquad\qquad \times \int d{\bm p}' (\hat{\bm k}\cdot {\bm p}')({\bm p}\cdot {\bm p}') \frac{\hat{\xi}({\bm k}, {\bm p}',s)}{\sigma({\bm k},{\bm p}', s)},
\label{hneeded}\end{aligned}$$ where $\cong$ signifies that we only kept the terms that contribute to $C(R,T)$. Now, the average $\langle \hat{h}({\bm k}, {\bm p}_1,s_1)\hat{h}(-{\bm k}, {\bm p}_2,s_2) \rangle_{\hat{\xi}}$ assumes a tractable form that can be used in Eq.. Separating the terms independent of $\Delta$, we obtain $C(R,T) = C_0(R,T) + C_1(R,T)$. Here, $$\begin{aligned}
&C_0(R,T) = \frac{\lambda n \kappa^2}{16\pi^4 } \lim_{t\rightarrow\infty} {\mathcal L}^{-1}_{s_1,t} {\mathcal L}^{-1}_{s_2,t+T}  \int d{\bm k} e^{-i {\bm k}\cdot{\bm R}} \frac{A^2(k\epsilon)}{k^2} \nonumber \\
&\quad \times \int d{\bm p} 
(\hat{\bm k} \cdot {\bm p})^2 \left[1- (\hat{\bm k} \cdot {\bm p})^2\right] \frac{1}{s_1+s_2}\nonumber \\
&  \qquad \times  
\frac{1}{\lambda + s_1 + i v_s k (\hat{\bm k}\cdot {\bm p})} 
\frac{1}{\lambda + s_2 - i v_s k (\hat{\bm k}\cdot {\bm p})},
\label{C0Laplace}\end{aligned}$$ represents correlations in the fluid created by non-interacting swimmers. The double inverse Laplace transform in the equation above can be performed using the method outlined in Appendix \[appendix:Cesare\]. It yields $$\begin{aligned}
 & \lim_{t\rightarrow\infty} {\mathcal L}^{-1}_{s_1,t} {\mathcal L}^{-1}_{s_2,t+T} \frac{1}{s_1+s_2} \frac{1}{\lambda + s_1 + i v_s k (\hat{\bm k}\cdot {\bm p})} \nonumber \\
& \qquad\times \frac{1}{\lambda + s_2 - i v_s k (\hat{\bm k}\cdot {\bm p})} = \frac{e^{-\lambda T + i v_s k T (\hat{\bm k} \cdot {\bm p})}}{2\lambda}.
\label{dilt_example}\end{aligned}$$ Performing the angular integration, we finally obtain $$\begin{aligned}
& C_0(R,T) = \frac{n \kappa^2 e^{-\lambda T}}{\pi^2}  \int_{0}^\infty dk \frac{\sin{k R}}{k R} A^2(k\epsilon) \nonumber \\
&\quad \times
\frac{ y(12-y^2) \cos{y} - (12-5y^2)\sin{y}
}{y^5}
\Bigg\vert_{y=v_s k T}.
\label{C0final}\end{aligned}$$

All other terms in Eq. correspond to additional correlations generated by the hydrodynamic interactions among the swimmers, and, as such, they are dependent on the dimensionless microswimmer density $\Delta$. Performing the angular integration over ${\bm p}_1$ and ${\bm p}_2$, gives $$\begin{aligned}
&C_1(R,T) \nonumber \\
& = \frac{2 \lambda n \kappa^2 }{15 \pi^2} \lim_{t\rightarrow\infty} {\mathcal L}^{-1}_{s_1,t} {\mathcal L}^{-1}_{s_2,t+T} \int_{0}^\infty dk \frac{\sin{k R}}{k R} A^2(k\epsilon) \nonumber \\
&\times \frac{1}{\lambda + s_1}\frac{1}{\lambda + s_2}\frac{1}{s_1 + s_2} \frac{z_1 \psi(z_1) + z_2 \psi(z_2)}{z_1+z_2}
 \nonumber \\
&\qquad \times\Bigg[ \frac{z_1 \psi(z_1)}{\omega - z_1 \psi(z_1)} + \frac{z_2 \psi(z_2)}{\omega - z_2 \psi(z_2)} \nonumber \\
&\qquad\qquad 
+ \frac{z_1 z_2 \psi(z_1)\psi(z_2)}{\left(\omega - z_1 \psi(z_1)\right)\left(\omega - z_2 \psi(z_2)\right)} \Bigg].
\label{C1laplace}\end{aligned}$$ Here, we introduced $\omega = v_s k/(\lambda \Delta A(k\epsilon))$, and the function $\psi(z)$, defined as $$\begin{aligned}
\psi(z) = \frac{5}{2}\frac{3z+2z^3-3(1+z^2)\arctan{z}}{z^5},
\label{psi}\end{aligned}$$ which is related to $f_2-f_1$ used in the previous Section. The variable $z_i=v_s k/(\lambda + s_i)$ allows us to write Eq. in a compact form but hides its complex dependence on the Laplace frequencies $s_1$ and $s_2$. Its inverse Laplace transform is discussed below.

Approximate double inverse Laplace transform {#subsection:approximateDILT}
--------------------------------------------

The integrand of Eq. is not a rational function of $s_1$ and $s_2$, and we were unable to calculate its double inverse Laplace transform exactly. Instead, here we develop a rational approximation to $\psi(z)$ that will allow us to find $C_1(R,T)$ analytically.

First, we observe that if the poles of an analytic function are known, its large-$t$ behaviour is determined by the pole with the smallest negative real part [@Doetsch1974]. Therefore, the presence of the pole at $-\lambda$ in Eq. makes all poles with real parts smaller than $-\lambda$ irrelevant in the large-$t$ limit. This reflects the fact that individual tumbling events are always a source of de-correlation between microswimmers.

Next, we introduce $L=v_s/(\lambda \epsilon)$, which compares the typical runlength of a swimmer to its size. In this work we consider $L=0-25$, ranging from non-swimming (shaker) particles to wild-type *E.coli* bacteria [@Berg1993]. Contributions to the integrand in Eq. with $k \epsilon>1$ are strongly suppressed by the regularising factor $A(k\epsilon)$, and therefore, when approximating $\psi(z)$, the relevant domain is $-\lambda<Re(s)<0$, with $v_s k/\lambda$ not exceeding $L$.

In Appendix \[appendix:psia\] we show that a surprisingly good approximation to $\psi(z)$ on this domain is given by $$\begin{aligned}
\psi_a(z) = \frac{7}{7+3 z^2}.
\label{psia}\end{aligned}$$ The simple structure of this expression allows us to deduce the pole structure of the integrand in Eq.. Indeed, with $\psi(z)$ replaced by $\psi_a(z)$, and factorising $$\begin{aligned}
\frac{1}{\omega - z \psi(z)} = \frac{7+3z^2}{3\omega\left(z-z_{+}\right)\left(z-z_{-}\right)},\end{aligned}$$ where $$\begin{aligned}
z_{\pm} =\frac{7}{6\omega} \left[ 1\pm\sqrt{1-\frac{12}{7}\omega^2} \right],
\label{z012}\end{aligned}$$ the denominators in Eq. can now be written as products of linear polynomials in $s_1$ and $s_2$. It is now straightforward to perform the inverse Laplace transform of this expression using the method outlined in Appendix \[appendix:Cesare\]. Taking the limit of $t\rightarrow\infty$, finally gives $$\begin{aligned}
&C_1(\rho,\tau) \nonumber \\
& = e^{-\tau} \frac{n \kappa^2 }{15 \pi^2 \epsilon} \int_{0}^\infty d\xi \frac{\sin{\xi \rho}}{\xi \rho} A^2(\xi) 
\Bigg[ 
-\cos{\left(\sqrt{\frac{3}{7}}L\xi\tau\right)} \nonumber \\
&+ \frac{e^{\frac{1}{2}A(\xi) \Delta \tau}}{1-A(\xi)\Delta + \frac{3}{7}L^2\xi^2} \Bigg\{ \nonumber \\
&  \frac{2-A(\xi)\Delta + \frac{6}{7}L^2\xi^2}{2-A(\xi)\Delta}\cosh{\left(\frac{1}{2}A(\xi) \Delta \tau  \sqrt{1-\frac{12 L^2\xi^2}{7 A^2(\xi)\Delta^2}}\right)} \nonumber \\
&\qquad +\frac{\sinh{\left(\frac{1}{2}A(\xi) \Delta \tau  \sqrt{1-\frac{12 L^2\xi^2}{7 A^2\left(\xi\right)\Delta^2}}\right)}}{\sqrt{1-\frac{12 L^2\xi^2}{7 A^2\left(\xi\right)\Delta^2}}}
\Bigg\}
\Bigg],
\label{C1final}\end{aligned}$$ where we changed the integration variable to $\xi = k\epsilon$, and introduced the dimensionless parameters $\rho = R/\epsilon$ and $\tau = \lambda T$. In Appendix \[appendix:dilt\], we verify that Eq. provides a good approximation to the long-time behaviour of Eq..

Results {#section:results}
=======

For the benefit of the readers who have skipped Section \[section:theory\], we repeat our main result, which comprises an explicit expression for the combined correlation function $C(R,T)$, defined in Eq.. It describes the steady-state correlations between the fluid velocity at two points in space separated by a distance $R$, and two instances in time separated by a time-interval $T$. Our theory includes full hydrodynamic interactions between microswimmers and is valid at any density up to the onset of collective motion. The result consists of the non-interacting part, $$\begin{aligned}
& C_0(\rho,\tau) = \frac{n \kappa^2 e^{-\tau}}{\pi^2\epsilon}  \int_{0}^\infty d\xi \frac{\sin{\xi \rho}}{\xi \rho} A^2(\xi) \nonumber \\
&\quad \times
\frac{ y(12-y^2) \cos{y} - (12-5y^2)\sin{y}
}{y^5}
\Bigg\vert_{y=L\xi\tau}, 
\label{ResultsC0}\end{aligned}$$ and the interacting correlation function $C_1(\rho,\tau)$, given in Eq.. Here, $\rho = R/\epsilon$, where $\epsilon$ is a lengthscale comparable to the microswimmer size, and $\tau = \lambda T$, where $\lambda$ is the tumbling rate. The parameter $L=v_s/(\lambda \epsilon)$ compares the typical distance covered by a swimming microorganism between two tumble events to its size. Finally, $\Delta = n/n_{crit}$ is the dimensionless number density of the particles, where $n_{crit} = 5\lambda/(B\kappa)$ is the onset of collective motion for pusher-like microswimmers [@Subramanian2009; @Hohenegger2010; @Stenhammar2017]; the parameter $B$, is defined after Eq.. Our theory is valid for $\Delta<1$.

The full expression, $C(\rho,\tau)=C_0(\rho,\tau)+C_1(\rho,\tau)$, given as a definite integral, constitutes the main technical result of our study. We now explicitly work out its predictions for the spatial and temporal correlation functions, and other experimentally accessible observables. When discussing their physical meaning, we are going to vary the dimensionless persistence length $L$, while keeping all the other parameters of the microswimmers fixed. We note that in reality the dipolar strength and shape of a microorganism uniquely determine its swimming speed, and hence $L$. We, however, see varying $L$ as a tool to disentangle the effects of self-propulsion (ability to change one’s position in space) from the strength of the hydrodynamic disturbances it causes. In particular, we will consider two limiting cases: shakers ($L=0$) and fast swimmers ($L\rightarrow\infty$). The former case corresponds to microswimmers that exert dipolar forces on the fluid but do not self-propel, and only change their positions due to being advected by the velocity fields created by other microswimmers [@Stenhammar2017]. The latter case, while obviously non-physical, is a useful tool to assess the effect of fast swimming on various quantities of interest. Finally, we note that the terms representing hydrodynamic interactions in Eq. are proportional to the swimmer’s nonsphericity $B$ that enters Jeffrey equation, Eq.. The limit of non-interacting microswimmers therefore corresponds to setting $B$ to zero, which, in turn, can be achieved by setting $\Delta=0$, while keeping $n$ finite.

Velocity variance
-----------------

![The fluid velocity variance $\langle U^2 \rangle$ normalised by its non-interacting value $\langle U^2 \rangle_0$ for various values of $L$. The dotted line represents the non-interacting case $\langle U^2 \rangle=\langle U^2 \rangle_0$. Note that the $L\rightarrow\infty$ line turns sharply upwards and diverges in the vicinity of $\Delta=1$ in a way that cannot be resolved on the scale of this graph.[]{data-label="Fig_Variance"}](Variance.pdf){width="0.8\linewidth"}

Our first quantity of interest is the fluid velocity variance, $\langle U^2 \rangle \equiv C(\rho=0,\tau=0)$. In the absence of thermal noise, re-arrangements of the microswimmer positions and orientations is the sole source of fluid velocity fluctuations. For this reason, it was used in previous studies as an order parameter to identify the onset of collective motion [@Stenhammar2017; @Bardfalvy2019]. Summing up Eqs. and , and setting $\rho=0$ and $\tau=0$, we obtain $$\begin{aligned}
& \langle U^2 \rangle  = \frac{\kappa^2 n}{15 \pi^2 \epsilon} \int_{0}^\infty d\xi A^2(\xi) \nonumber \\
& \qquad\qquad \times \frac{2-A(\xi)\Delta+\frac{6}{7}L^2 \xi^2}{\left(2-A(\xi)\Delta\right)\left(1-A(\xi)\Delta+\frac{3}{7}L^2 \xi^2\right)}.
\label{resultvarience}\end{aligned}$$ We evaluate this integral numerically and plot the fluid velocity variance normalised by its value in the non-interacting case, $\langle U^2 \rangle(\Delta=0)\equiv \langle U^2 \rangle_0$, given by [@Bardfalvy2019] $$\begin{aligned}
\langle U^2 \rangle_0 = \frac{\kappa^2 n}{15 \pi^2 \epsilon} \int_{0}^\infty d\xi A^2(\xi) = \frac{21\kappa^2 n}{2048\epsilon}.\end{aligned}$$ Note that $\langle U^2 \rangle_0$ corresponds to a superposition of uncorrelated fluctuations in the fluid velocity, which, by virtue of the central limit theorem, is proportional to $n$. Any deviations of $\langle U^2 \rangle$ from that value signify the presence of correlations.

As can be seen in Fig.\[Fig\_Variance\], the fluid velocity fluctuations exhibit significant correlations at any density of the microswimmers, as was recognised previously [@Stenhammar2017]. Starting from its non-interacting value at $\Delta=0$, the variance increases with $\Delta$, until it diverges at the onset of collective motion. The strongest correlations are exhibited by suspensions of shakers, while swimming acts to reduce correlations. For large but finite values of $L$, the variance increases mildly from its non-interacting value, until it rises sharply in a small vicinity of $\Delta=1$, with the size of this region shrinking with $L$. Interestingly, the rise of $\langle U^2 \rangle_0$ for $\Delta<1$ remains finite even in the $L\rightarrow\infty$ limit. In other words, while swimming clearly reduces correlations, it does not remove them entirely, and the suspension is never described by the mean-field theory.

To determine the scaling of the fluid velocity variance as $\Delta\rightarrow 1$, we observe that in that limit the integrand in Eq. is dominated by small values of $\xi$, where $A(\xi) \approx 1 - \xi^2/4$. Using this approximation in Eq. and replacing the upper integration limit by unity, we obtain $$\begin{aligned}
& \langle U^2 \rangle  \sim \frac{\kappa^2 n}{15 \pi \epsilon} \frac{1}{\sqrt{1+\frac{12}{7}L^2}\sqrt{1-\Delta}}, \quad \Delta \rightarrow 1.
\label{varianceasymptotic}\end{aligned}$$ Therefore, our theory predicts that the fluid velocity variance diverges as $(1-\Delta)^{-1/2}$ in the vicinity of the transition to collective motion, irrespectively of $L$.

Spatial correlations
--------------------

![image](Spatial.pdf){width="1.0\linewidth"}

Our next quantity of interest is the equal-time spatial correlation function, $C(\rho, T=0)$, given by $$\begin{aligned}
&C(\rho) = \frac{\kappa^2 n}{15 \pi^2 \epsilon} \int_{0}^\infty d\xi \frac{\sin{\xi \rho}}{\xi \rho} A^2(\xi) \nonumber \\
& \qquad\qquad \times \frac{2-A(\xi)\Delta+\frac{6}{7}L^2 \xi^2}{\left(2-A(\xi)\Delta\right)\left(1-A(\xi)\Delta+\frac{3}{7}L^2 \xi^2\right)}.
\label{spatialcorrelations}\end{aligned}$$ While this integral cannot be evaluated analytically, a good approximation can be obtained by setting $A(\xi)=1$ in the integrand, yielding $$\begin{aligned}
&C(\rho) \approx \frac{\kappa^2 n}{30 \pi \epsilon \left(1-\Delta\right)\rho} \nonumber \\
&\qquad\qquad \times \left[ 1 - \frac{\Delta}{2-\Delta} \exp{\left\{-\sqrt{\frac{7}{3}}\sqrt{1-\Delta}\frac{\rho}{L} \right\}}\right].
\label{spatialcorrelationsanalytic}\end{aligned}$$ For $\Delta=0$, this equation reproduces the result obtained previously for non-interacting swimmers [@Zaid2011; @Underhill2011; @Belan2019; @Bardfalvy2019].

In Fig.\[Fig\_spatial\] we evaluate Eq. numerically and compare it against the analytic approximation, Eq.; $\kappa^2 n_{crit}/(15 \pi^2 \epsilon)$ is chosen as the normalisation factor. For all values of $L$ and $\Delta$, the approximation works well for all but small spatial separations $\rho$, where the spatial correlation function is, essentially, equal to the fluid velocity variance. As with the fluid velocity variance, the strongest correlations are exhibited by suspensions of shakers, $L=0$. In this case, the spatial correlation function changes very slowly at short distances, and decays as $\rho^{-1}$ at large distances. Close to the onset of collective motion, the typical scale $\rho_0$ at which the crossover occurs can be estimated from Eqs. and , by requiring that $C(\rho_0) = \langle U^2 \rangle$. For $L=0$, this yields $\rho_0 \sim (1-\Delta)^{-1/2}$. This is readily verified by the data in Fig.\[Fig\_spatial\]A: As the system approaches the onset of collective motion, the overall strength of the correlations grows, with the region of strong correlations extending to progressively larger scales.

![The spatial correlation function $C(\rho)$ as a function of the distance $\rho$ for $\Delta=0.9$ and various values of $L$. At sufficiently large distances, $C(\rho)$ recovers the shaker behaviour, while at small distances correlations are suppressed by swimming. Note that the $L\rightarrow\infty$ line, serving as the limit beyond which correlations cannot be suppressed, joins the shaker line at $\rho\rightarrow\infty$.[]{data-label="Fig_spatial_delta09"}](SpatialDelta09.pdf){width="0.8\linewidth"}

The effect of swimming on the behaviour of $C(\rho)$ is demonstrated in Figs.\[Fig\_spatial\]B-\[Fig\_spatial\]C. As $L$ increases, the strongly correlated core at moderate separations shrinks, indicating that the steady growth of orientational correlations is reduced by the mixing introduced by swimming. The overall strength of correlations inside the core also decreases with $L$, reflecting the reduction of the fluid velocity variance by swimming. At large distances, $C(\rho)$ recovers the behaviour seen in shakers, with the crossover distance given by $\rho_1 \sim L (1-\Delta)^{-1/2}$, as can be deduced from the exponential in Eq.. This behaviour is further demonstrated in Fig.\[Fig\_spatial\_delta09\], where we plot $C(\rho)$ for $\Delta=0.9$ and various values of $L$. In the limit of fast swimming, $L\rightarrow\infty$, the correlation function deviates modestly from the non-interacting case for almost all values of $\Delta$, exhibiting a quick rise and the divergence associated with the onset of collective motion only in a very small vicinity of $\Delta=1$.

The data in Fig.\[Fig\_spatial\] and Eq. demonstrate that $C(\rho)$ exhibits an algebraic decay for large distances, and a true correlation length can thus not be defined. A phenomenological correlation length $\xi_{corr}$ can nevertheless be defined as a distance over which $C(\rho)$ decreases by certain amount, as has been employed in [@Gachelin2014; @Bardfalvy2019]. Setting $C(\xi_{corr}) = \alpha \langle U^2 \rangle$, with $\alpha<1$, we obtain $$\begin{aligned}
\xi_{corr}\sim (1-\Delta)^{-1/2}, \quad \Delta\rightarrow1, \end{aligned}$$ similar to any other typical distance discussed above.

Fluid velocity spectrum
-----------------------

![image](Spectra.pdf){width="1\linewidth"}

Next, we discuss the fluid velocity energy spectrum $E(k)$ that is closely related to the spatial correlation function $C(\rho)$. Defined as $$\begin{aligned}
E(k) = 4\pi k^2 \overline{\hat{U}^{\alpha}\left( {\bm k} \right) \hat{U}^{\alpha}\left( -{\bm k} \right)},
\label{eq:spectrum}\end{aligned}$$ this quantity is often used in turbulence research to study the cascade of the kinetic energy [@Townsend1980]. Although the kinetic energy is not a useful concept for Stokesian flows, $E(k)$ provides an insight into the relative strength of fluid motion at various scales. The energy spectrum is proportional to the Fourier transform of $C(\rho)$, and, up to a prefactor is given by the integrand of Eq. $$\begin{aligned}
&E(\xi) = \frac{8\pi}{15} \kappa^2 n A^2(\xi) \nonumber \\
&\qquad\qquad \times \frac{2-A(\xi)\Delta+\frac{6}{7}L^2 \xi^2}{\left(2-A(\xi)\Delta\right)\left(1-A(\xi)\Delta+\frac{3}{7}L^2 \xi^2\right)},
\label{energyspectrum}\end{aligned}$$ where, again, $\xi=k\epsilon$. This expression is plotted in Fig.\[Fig\_Spectra\] for various values of $\Delta$ and $L$.

First, we observe that $E(\xi)$ has significant energy content at all large scales, $\xi<1$, that quickly decays to zero at the organism-size scales, $\xi\sim1$, due to the regularising factor $A(\xi)$. This is not caused by some form of energy cascade, but is due to the nature of the dipolar field created by the microswimmers. Indeed, the dipolar velocity field decays in space as $r^{-2}$, while its Fourier transform scales as $k^{-1}$. Together with the definition of $E(k)$, Eq., this implies that $E(k)\sim k^0$ even for a single microswimmer, i.e. the dipolar field has a constant energy content at every scale.

In the presence of interactions, the energy spectrum of shakers ($L=0$) preserves the overall structure described above, while its absolute value increases with $\Delta$ and, eventually, diverges at $\Delta=1$. For swimmers, the increase in the energy content is mostly confined to large scales, while in the limit of fast swimming (not shown), the rise in the energy content on the approach to the onset of collective motion is confined to the largest scales available ($k\rightarrow0$) and starts to be visible only in a very close vicinity of $\Delta=1$.

Temporal correlations {#section:results:temporal}
---------------------

![image](Temporal.pdf){width="1.0\linewidth"}

The temporal correlation function $C(\tau)=C_0(\rho=0,\tau)+C_1(\rho=0,\tau)$ is given by Eqs. and . The corresponding expressions do not simplify significantly in the limit $\rho=0$, and we do not repeat them here. In Fig.\[Fig\_Temporal\] we plot $C(\tau)$ normalised by its value at $\tau=0$, which is given by the fluid velocity variance $\langle U^2\rangle$. As with the other quantities discussed above, the temporal correlation function exhibits a progressively slower decay as $\Delta$ approaches the onset of collective motion, eventually diverging at $\Delta=1$. For swimmers, this is offset by a decay of $C(\tau)$ at short times that becomes more pronounced as $L$ increases. For very large swimming speeds, the temporal correlations differ only marginally from the non-interacting case for most values of $\Delta$, eventually exhibiting a rapid increase and divergence in a very small vicinity of $\Delta=1$.

To understand the behaviour of $C(\tau)$ at long times, we analyse its individual contributions. The integral in the non-interacting part, $C_0(T)$, can be explicitly evaluated giving $$\begin{aligned}
& C_0(\tau) = \frac{n \kappa^2 }{\pi \epsilon} \frac{e^{-\tau}}{8\alpha^4 (4+\alpha^2)^2} 
\Bigg[ \nonumber \\
& \qquad\qquad\qquad 4(24+8\alpha^2+\alpha^4) \mathbb{E}\left( -\frac{\alpha^2}{4}\right) \nonumber \\
& \qquad\qquad\qquad-(4+\alpha^2)(24+5\alpha^2) \mathbb{K}\left( -\frac{\alpha^2}{4}\right)
\Bigg],
\label{}\end{aligned}$$ where $\alpha=L\tau$, and $\mathbb{K}(x)$ and $\mathbb{E}(x)$ are the complete elliptic integrals of the first and second order, respectively. In the limits of small and large $\alpha$ this equation predicts $$\begin{aligned}
C_0(\tau) \sim \frac{n \kappa^2 }{\pi \epsilon} e^{-\tau} \times
\begin{cases}
\frac{21\pi}{2048},\quad &L\tau\rightarrow0, \\ \\
\frac{1}{4 (L\tau)^3},\quad &L\tau\rightarrow\infty.
\end{cases}\end{aligned}$$ At short times, tumbling is the leading source of decorrelation, while at large $\tau$ the non-interacting temporal correlation function $C_0$ decays as $\tau^{-3}e^{-\tau}$, as reported previously [@Belan2019; @Bardfalvy2019]. The crossover time is set by $\alpha=L\tau=v_s T/\epsilon\sim1$, and corresponds to the time interval needed for a microswimmer to swim its own size.

To understand the large-$\tau$ asymptotic behaviour of $C_1(\tau)$, we observe that $$\begin{aligned}
e^{-\tau} \int_0^\infty d\xi A(\xi)^2 \left\{ 
\begin{array}{c}
\sin{\gamma \tau \xi} \\
\cos{\gamma \tau \xi}
\end{array} \right\} 
\underset{\tau\rightarrow\infty}{\sim}
e^{-\tau}  \left\{ 
\begin{array}{c}
\tau^{-1}\\
\tau^{-5}
\end{array} \right\},\end{aligned}$$ where $\gamma$ is a real constant. This result implies that a trigonometric function in the integrand of Eq. generates a contribution to $C_1(\tau)$ that decays on the same timescale as the non-interacting part $C_0(\tau)$, and does not contribute to the slow decay in Fig.\[Fig\_Temporal\]. In turn, this restricts the integration domain to $\xi \in [0,\xi_*]$, with $$\begin{aligned}
\xi_* = \sqrt{\frac{7}{12}}\frac{\Delta}{L},\end{aligned}$$ which ensures that the arguments of the hyperbolic functions in Eq. are real. Introducing $\zeta = \xi/\xi_*$, $C_1(\tau)$ can be approximated as $$\begin{aligned}
&C_1(\tau)
\underset{\tau\rightarrow\infty}{\sim}
  \frac{n \kappa^2 }{15 \pi^2 \epsilon} e^{-\tau\left(1-\frac{1}{2}\Delta\right)} \xi_* \int_{0}^1 d\zeta  \frac{1}{1-\Delta + \frac{1}{4}\Delta^2\zeta^2} \nonumber \\
& \qquad \times \Bigg\{
\frac{2-\Delta +\frac{1}{2}\Delta^2\zeta^2}{2-\Delta}\cosh{\left(\frac{1}{2}\Delta \tau  \sqrt{1-\zeta^2}\right)} \nonumber \\
&\qquad\qquad\qquad + \frac{1}{\sqrt{1-\zeta^2}} \sinh{\left(\frac{1}{2}\Delta \tau  \sqrt{1-\zeta^2}\right)} 
\Bigg\},
$$ where we used $A(\xi<\xi_*)\sim 1$ for not-too-small values of $L$. In the limit of large $\tau$, this can be further approximated by $$\begin{aligned}
&C_1(\tau)
\underset{\tau\rightarrow\infty}{\sim}
  \frac{n \kappa^2 }{15 \pi^2 \epsilon} e^{-\tau\left(1-\Delta\right)} \frac{\xi_*}{1-\Delta} \int_{0}^1 d\zeta  \,
e^{-\frac{1}{4}\tau \Delta \zeta^2} \nonumber \\
& =  \frac{n \kappa^2 }{15 \pi^2 \epsilon}  \sqrt{\frac{7\pi}{12}\frac{\Delta}{\tau}} \frac{1}{L(1-\Delta)} e^{-\tau\left(1-\Delta\right)} \mathrm{erf}\left( \frac{1}{2}\tau\Delta\right),
\label{C1asymptotics}\end{aligned}$$ where $\mathrm{erf}(x)$ denotes the error function. Predictions of Eq. are plotted in Fig.\[Fig\_Temporal\]B and C as dashed lines. We find a good agreement between its prediction and the true decay of $C(\tau)$ as $\tau\rightarrow\infty$.

To extract the typical timescale $\tau_{corr}$ of the fluid velocity fluctuations on the approach to collective motion, we combine Eqs. and to obtain $$\begin{aligned}
\frac{C(\tau)}{\langle U^2\rangle}
\underset{\tau\rightarrow\infty}{\sim}
\frac{e^{-\tau\left(1-\Delta\right)}}{\sqrt{\tau \left(1-\Delta\right)}}, \quad \Delta \rightarrow 1,\end{aligned}$$ which implies $$\begin{aligned}
\tau_{corr} \sim (1-\Delta)^{-1}.\end{aligned}$$

Enhanced diffusivity
--------------------

As the final observable, we consider here the enhanced diffusivity of a passive tracer particle embedded in a suspension of motile microorganisms. The tracer is assumed to be neutrally buoyant and move due to advection by the velocity fields created by the microswimmers. Brownian diffusion of the tracer is significantly weaker than its enhanced counterpart, and is neglected for simplicity. This problem has been extensively studied both experimentally [@Wu2000; @Kim2004; @Leptos2009; @Kurtuldu2011; @Mino2011; @Mino2013; @Jepson2013; @Patteson2016] and theoretically [@Underhill2008; @Dunkel2010; @Ishikawa2010; @Childress2010; @Childress2011; @Pushkin2013; @Pushkin2013jfm; @Morozov2014; @Kasyap2014; @Thiffeault2015; @Burkholder2017] in the dilute regime, where $\Delta\ll 1$, and for arbitrary densities of shakers [@Stenhammar2017]. Here, we consider the case of arbitrary density $\Delta < 1$ and $L$.

The position of the tracer ${\bm a}(T)$ obeys the following equation of motion $$\begin{aligned}
\dot{\bm a}(T) = {\bm U}({\bm a}(T),T),
\label{tracer_position}\end{aligned}$$ which implies that the tracer is point-like and follows the velocity of the fluid at its position. The long-time behaviour of such a tracer is diffusive [@Wu2000; @Leptos2009; @Thiffeault2015], and the associated diffusion coefficient can be extracted in the usual way $$\begin{aligned}
D = \lim_{T\rightarrow\infty} \frac{1}{6T} \overline{{\bm a}(T)\cdot {\bm a}(T)}.\end{aligned}$$ Here, the bar denotes the average over the history of tumble events, and has the same meaning as in Eq.. Solving formally Eq., ${\bm a}(T)={\bm a}(0) + \int_0^T dt' {\bm U}({\bm a}(t'),t')$, the diffusion coefficient can be written as [@Kubo1966] $$\begin{aligned}
&D = \frac{1}{3} \lim_{t\rightarrow\infty}\int_{0}^\infty dT\, \overline{{\bm U}({\bm a}(t+T),t+T) \cdot {\bm U}({\bm a}(t),t)}.\end{aligned}$$ Here, $t$ is sufficiently large so that any influence of the initial conditions has died away. To proceed, we observe that ${\bm U}({\bm a}(t+T),t+T)$ can be iteratively calculated by substituting the formal solution for ${\bm a}(T)$ into its spatial argument, i.e. $$\begin{aligned}
&{\bm U}({\bm a}(t+T),t+T) = {\bm U}({\bm a}(t),t+T)\nonumber \\
& +\nabla  {\bm U}({\bm a}(t),t+T) \cdot \int_{t}^{t+T} dt' {\bm U}({\bm a}(t'),t') + \cdots.
\label{DUapprox}\end{aligned}$$ As was argued by Pushkin and Yeomans [@Pushkin2013], for very dilute suspensions velocity gradients over the typical distance travelled by the tracer particle during the microswimmer runtime are small compared to the velocity of the fluid at any of these positions, and can be neglected. Therefore, we can approximate the diffusion coefficient as $$\begin{aligned}
& D \approx \frac{1}{3} \int_{0}^\infty dT\, \overline{{\bm U}({\bm a}(t),t+T) \cdot {\bm U}({\bm a}(t),t)}\nonumber \\
&\qquad\qquad = \frac{1}{3} \int_{0}^\infty dT\, C(T).
\label{DvsCT}\end{aligned}$$ As we have seen in Section \[section:results:temporal\], as $\Delta$ increases, the correlation time increases from $\lambda^{-1}$ (corresponding to $\tau_{corr}=1$) in the very dilute regime to progressively larger values, eventually diverging as $\Delta\rightarrow 1$, implying that the second, etc. terms in Eq. grow rapidly in this limit. However, the fluid velocity variance, which sets the magnitude of the leading term in Eq. also diverges as $\Delta\rightarrow 1$. Further work is required to assess the validity of the approximation above for all values of $\Delta$. Here, we proceed by using Eq. with the potential caveat that it might not be accurate in the vicinity of $\Delta=1$.

The integral in Eq. can be evaluated explicitly, leading to $D = D_0 + D_1$, where the non-interacting and interacting contributions are given by $$\begin{aligned}
D_0 = \frac{\kappa^2 n}{45 \pi^2 \lambda \epsilon} \int_{0}^\infty d\xi A^2(\xi) \psi\left( \xi L\right),
$$ and $$\begin{aligned}
&D_1 = \frac{\kappa^2 n \Delta}{45 \pi^2 \lambda \epsilon} \int_{0}^\infty d\xi A^3(\xi) \nonumber\\
&\qquad\qquad \times \frac{ 2-A(\xi)\Delta+\frac{6}{7}L^2 \xi^2 }{\left( 1+\frac{3}{7} L^2\xi^2\right)\left( 1-A(\xi)\Delta+\frac{3}{7}L^2 \xi^2\right)^2},\end{aligned}$$ respectively, and $\psi(x)$ is defined in Eq.. At this point, we would like to comment on the shaker limit of these expressions, when they should reduce to the ones obtained by Stenhammar *et al.* [@Stenhammar2017]. Instead, we observe that the expression for $D_1$ reported there erroneously contained $A^2(\xi)$ instead of $A^3(\xi)$ under the integral. We note, however, that since $A(\xi)$ is a regularised representation of a step function, this has almost no bearing on the numerical evaluation of $D_1$ presented in [@Stenhammar2017].

The integral in the non-interacting part $D_0$ cannot be represented in terms of special functions, but its limiting behaviour can readily be obtained. Combining the asymptotic results for $L=0$ and $L\rightarrow\infty$, results in the following approximation $$\begin{aligned}
D_0 \approx \frac{\kappa^2 n}{\lambda \epsilon} \frac{7}{2048 + 336 \pi L}.
\label{D0approx}\end{aligned}$$

To derive an approximate expression for $D_1$, we set $A(\xi)\approx 1$ under the integral sign, to obtain $$\begin{aligned}
D_1 \approx \frac{\kappa^2 n}{90 \pi \lambda \epsilon \sqrt{1+\frac{12}{7}L^2}} \left\{ \frac{2-\Delta}{\left( 1-\Delta\right)^{3/2}} -2\right\}.
\label{D1approx}\end{aligned}$$

![image](Diffusivity.pdf){width="0.9\linewidth"}

In Fig.\[Fig\_Diffusivity\] we compare the numerical evaluation of $D_0$ and $D_1$ against Eq. and . We observe that while the uniform approximation Eq. does not work well for small but finite values of $L\sim1$, all other values of $L$ are well-represented by the approximation. The interacting part of the diffusivity is well-approximated by Eq..

Finally, we remark that Eq. predicts that $$\begin{aligned}
D_1 \sim (1-\Delta)^{-3/2}, \quad \Delta\rightarrow1,\end{aligned}$$ even though this prediction should be treated with caution, as discussed above.

Discussion and Conclusion {#section:discussion}
=========================

In this work, we have presented a kinetic theory for dilute suspensions of pusher-like microswimmers interacting via long-ranged dipolar fields. Our theory goes beyond the mean-field assumption and explicitly includes correlations between microswimmers. We have overcome a significant technical difficulty in including particle self-propulsion that has limited our previous work to the case of shaker microswimmers [@Stenhammar2017]. Our theory makes explicit predictions for various experimentally relevant observables for any density of microswimmers up to the onset of collective motion. All of its parameters can be independently measured in experiments, and its predictions can be directly compared against experimental data.

The results of our theory, presented in Section \[section:results\], reveal that all observables considered deviate from their mean-field values, which can be recovered from our results by setting $\Delta=0$, indicating that the mean-field theory is incorrect at any density below the onset of collective motion. We have also uncovered the following interplay between the strength of correlations between microswimmers and their self-propulsion speed. For all observables considered, the strongest correlations are exhibited by suspensions of shakers, $L=0$. This can be readily seen by observing that, in the absence of self-propulsion, the microswimmer positions only change due to their mutual advection. In dilute suspension, displacements thus accumulated over one correlation time are small compared to the interparticle distances, and, to first approximation, shaker suspensions perform orientational dynamics only. In turn, this implies that they spend maximum amount of time possible adjusting to the orientational fields created by other microswimmers. In contrast, motile microswimmers are aligning in a local velocity field that constantly changes due to their self-propulsion, implying weaker correlations in such suspensions. This effect becomes stronger as $L$ increases.

The degree to which correlations are suppressed by self-propulsion depends on the nature of the observable. Spatial-like observables (the fluid velocity variance, the energy spectrum, and the spatial correlation function) are significantly reduced as $L$ increases, but do not reach their mean-field values even in the limit $L\rightarrow\infty$. For instance, as can be seen from Fig.\[Fig\_Variance\], the fluid velocity variance is significantly larger than its mean-field value at any density $\Delta$, even in the limit of fast swimming. In a similar fashion, as $L\rightarrow\infty$, the spatial correlation function in Fig.\[Fig\_spatial\_delta09\] does not reduce to its mean-field behaviour, which is given by the $\Delta\rightarrow0$ limit in Fig.\[Fig\_spatial\]. Instead, it recovers the strongly correlated shaker-like behaviour at sufficiently large distances. This can be understood by employing the same argument as above. For any value of $L$, there are such separations $\rho$ that the typical distance travelled by a microswimmer during one correlation time of the suspension is small compared to $\rho$. For such separations, the difference between swimmers and shakers vanishes and $C(\rho)$ recovers its shaker-like behaviour.

On the other hand, temporal-like observables (the temporal correlation function and the enhanced diffusivity of tracer particles) are almost completely suppressed as $L\rightarrow\infty$ for $\Delta<1$, though they still diverge in the limit of $\Delta\rightarrow1$. This behaviour mirrors the dependence of their mean-field values on $L$, which vanish in the limit of fast swimming below the onset of collective motion. An intuitive argument for this behaviour has been put forward by Dunkel *et al.* [@Dunkel2010], who demonstrated that the total displacement of a tracer by a single motile particle vanishes as the length of a straight path covered by the swimmer diverges. This is fundamentally related to the time-reversibility of Stokesian flows. The presence of correlations between microswimmers breaks this time-reversibility: although the pathway between two states in phase space is still reversible, the probabilities of finding the suspension in those states are a priori different. Strong swimming introduces effective phase space mixing and recovers equal a priori probabilities for the phase space states. Again, this argument only holds for $\Delta<1$ when $L$ is large, yet finite.

  **Observable**                **Scaling law for $\Delta\rightarrow1$**
  ----------------------------- ------------------------------------------
  Fluid velocity variance       $(1-\Delta)^{-1/2}$
  (Pseudo-)correlation length   $(1-\Delta)^{-1/2}$
  Correlation time              $(1-\Delta)^{-1}$
  Enhanced diffusivity          $(1-\Delta)^{-3/2}$

  : Critical exponents[]{data-label="table:exponents"}

To gain further insight into the nature of the transition to collective motion exhibited by our model, we extracted the scaling behaviour of the observables considered in this work upon the approach to the onset, $\Delta=1$. All of these quantities diverge at $\Delta=1$ and the values of the critical exponents predicted by our theory are summarised in Table \[table:exponents\]. We want to stress that these exponents rely on the approximation introduced in Section \[subsection:approximateDILT\], and while we are confident that it semi-quantitatively captures the spatial and temporal behaviour of the generalised correlation function $C(\rho,\tau)$ for $\Delta<1$, its quality in the close vicinity of $\Delta=1$ is untested. The values presented in Table \[table:exponents\] should thus be seen as a starting point, and more work is needed to make a conclusive statement about the critical exponents predicted by our model. This is especially important in the context of defining universality subclasses of “wet" active matter models.

In this work, we have only considered pusher-like microswimmers below the onset of collective motion. Recent simulations suggest [@Stenhammar2017; @Bardfalvy2019] that suspensions of pullers also exhibit strong correlations, although their effect is opposite to what is observed for pushers. The results presented in this work cannot be used to study this effect, i.e. by replacing $\Delta$ with $-\Delta$ in the relevant expressions. Instead, to extend our theory to pullers, one would have to re-evaluate the long-term behaviour of the approximate double inverse Laplace transform in Section \[subsection:approximateDILT\] for negative values of $\Delta$.

Finally, we note that the scenario emerging from our work is that microswimmer motility suppresses correlations and brings the suspension closer to the mean-field state. This observation could prove to be instrumental in understanding the origins of a surprising observation made by Stenhammer *et al.* [@Stenhammar2017] regarding the nature of the transition to collective motion in bacteria. There, we performed large-scale Lattice Boltzmann simulations of dilute suspensions of pusher-like microswimmers, and observed that all quantities of interest increased sharply around the threshold predicted by the mean-field theory ($\Delta=1$ in the notation adopted here). This increase, however, fell short of a divergence that should be associated with a true transition. While it would be natural to assume that this is the result of finite-size effects always present in computer simulations, Stenhammar *et al.* [@Stenhammar2017] did not see any sharpening of the curves around the supposed transition point upon a systematic increase of the system size. This phenomenon is currently not understood. Either the finite-size effects are decaying very slowly and realistic-size computer simulations cannot reach such large system sizes, or the presence of strong correlations below the onset of collective motion changes the stability properties of the homogeneous and isotropic state, replacing the mean-field transition with a correlation-smoothened cross-over. The results of this work provide a way to test the latter hypothesis: Since swimming suppresses correlations, sharpening of the potential transition/cross-over should be observed in suspensions with progressively increased swimming speed.

Double inverse Laplace transform of an archetypal term {#appendix:Cesare}
======================================================

Here, we show how to calculate the double inverse Laplace transform of Eq.. The derivation of Eq. is similar, though lengthy, and we do not present it here.

We start by observing that the double inverse Laplace transform in Eq., given in terms of two Bromwich integrals [@Doetsch1974], can be written as $$\begin{aligned}
\label{app:B-1}
\lim_{t\to\infty}
\int_{\Gamma_1}\frac{ds_1}{2\pi i} \frac{e^{s_1 t}}{\lambda + s_1 + i v_s k (\hat{\bm k}\cdot {\bm p})}J(s_1),\end{aligned}$$ where $$\begin{aligned}
\label{app:B-2}
J(s_1)
=
\int_{\Gamma_2}\frac{ds_2}{2\pi i} 
\frac{1}{s_1+s_2}
\frac{e^{s_2 (t+T)}}{{\lambda + s_2 - i v_s k (\hat{\bm k}\cdot {\bm p})}}.\end{aligned}$$ By the definition of the inverse Laplace transform [@Doetsch1974], the contours defining the integrals above have to be chosen such that $\Gamma_2$ passes on the right of $-s_1$ and of $-\lambda +iv_s k(\hat{\bm k}\cdot {\bm p})$, while $\Gamma_1$ should pass on the right of all the poles of $J(s_1)$ and of $-\lambda-i v_s k(\hat{\bm k}\cdot {\bm p})$. Observe that the first condition implies that $\Gamma_2$ should be chosen on the right of $-\Gamma_1$.

Next, we observe that, again from the definition of the inverse Laplace transform, $J(s_1)$ is only defined for $\textrm{Re}(s_1)>0$. To proceed, we follow the method often utilised in plasma physics to describe the Landau damping [@Balescu1975]. We perform the analytic continuation of $J(s_1)$ to purely imaginary values of $s_1$ (recall that the analytic continuation of a complex function defined on an open set is the only function $\hat{J}(s_1)$ that is analytic, defined on a larger set, and equals $J(s_1)$ on the original set), and replace $J(s_1)$ with $\hat{J}(s_1)$ in Eq. . Since $\lambda>0$, the difficulty in performing the analytic continuation of $J(s_1)$ lies in the pole at $s_2=-s_1$ of the integrand from Eq.. We, therefore, define $$\begin{aligned}
& \hat{J}(s_1) = \int^* \frac{ds_2}{2\pi i} 
\frac{e^{s_2 (t+T)}}{s_2+\lambda-i v_s k(\hat{\bm k}\cdot {\bm p})} \nonumber \\
& \qquad\quad +\frac{1}{2} \ell(-s_1)
 \frac{e^{s_2 (t+T)}}{s_2+\lambda-i v_s k(\hat{\bm k}\cdot {\bm p})}\Bigg\vert_{s_2=-s_1},\end{aligned}$$ where $$\begin{aligned}
\ell(s_1) = \begin{cases}
0, & \textrm{Re}(s_1)>0, \\
1, & \textrm{Re}(s_1)=0, \\
2, & \textrm{Re}(s_1)<0.
\end{cases}\end{aligned}$$ The meaning of the integral denoted by $\int^* ds_1$ above depends on the sign of $\textrm{Re}(s_1)$: If $\textrm{Re}(s_1)>0$, it is just a standard complex integral over a contour passing on the right of $-s_1$ and of $-\lambda +iv_s k(\hat{\bm k}\cdot {\bm p})$; If $\textrm{Re}(s_1)=0$, $\int^* ds_1$ stands for a principal value integral; Finally, if $\textrm{Re}(s_1)<0$, $\int^* ds_1$ stands for a standard complex integral over a contour passing on the left of $-s_1$ but on the right of $-\lambda +iv_s k(\hat{\bm k}\cdot {\bm p})$. With the definitions above, it is easy to show that $\hat{J}(s_1)$ is holomorphic in an infinitesimal stripe around $s_1\in \mathbb{R}$. Hence, it is the analytic continuation of $J(s_1)$.

Replacing $J(s_1)$ by $\hat{J}(s_1)$ in Eq., we obtain two terms. The first term, containing $\int^*ds_2$, vanishes for $t\to\infty$, since we are now free to choose the integration contours $\Gamma_1$ and $\Gamma_2$ such that $\textrm{Re}(s_1+s_2)<0$. The other term reads $$\begin{aligned}
& \lim_{t\to\infty} \frac{1}{2}\int_{\Gamma_1}\frac{ds_1}{2\pi i} 
\frac{1}{s_1+\lambda+i v_s k(\hat{\bm k}\cdot {\bm p})} \nonumber \\
&\qquad\qquad\qquad\qquad \times \frac{\ell(-s_1)\,e^{-s_1 T}}{-s_1+\lambda-i v_s k(\hat{\bm k}\cdot {\bm p})}.\end{aligned}$$ Closing the contour at $+\infty$, the only pole contributing to the integral is at $s_1=\lambda-i v_s k(\hat{\bm k}\cdot {\bm p})$, and we obtain $$\begin{aligned}
&\lim_{t\to\infty}
\int_{\Gamma_1}\frac{ds_1}{2\pi i} \frac{e^{s_1 t}}{s_1+\lambda+i v_s k(\hat{\bm k}\cdot {\bm p})}J(s_1) \nonumber \\
&\qquad\qquad\qquad\qquad = \frac{e^{-\lambda T + i v_s k T (\hat{\bm k} \cdot {\bm p})}}{2\lambda}.\nonumber\end{aligned}$$ This completes the proof of the equality in Eq..

Approximating $\psi(z)$ {#appendix:psia}
=======================

Here, we develop an approximation to $\psi(z)$ from Eq.. Our goal is to find a rational function with a pole structure that is similar to the original $\psi(z)$. As discussed in Section \[subsection:approximateDILT\], the relevant domain is set by the values of $z$ given by $z = \beta/(1+s/\lambda)$, with $\beta=v_s k/\lambda$ varying from $0$ to $L=v_s/(\lambda \epsilon)=0-25$, and by the real part of $s$ ranging from $-\lambda$ to $0$.

Our starting point are the observations that as $z\rightarrow0$, $\psi(z) \rightarrow 1-3z^2/7$, while for $z\rightarrow\infty$, $\psi(z) \rightarrow 0$. Both asymptotic behaviours can be combined into $\psi_a(z)= 7/(7+3z^2)$. Now we show that this is a surprisingly good approximation to $\psi(z)$, both reproducing its global shape and having a similar pole structure.

![Comparison between $\psi(z)$ (solid lines) and $\psi_a(z)$ (dotted lines) for $z = \beta/(1+s/\lambda)$ for real values of $s$ and various values of $\beta$.[]{data-label="psi_vs_psia"}](Psia.pdf){width="0.8\linewidth"}

In Fig.\[psi\_vs\_psia\] we compare $\psi(z)$ and $\psi_a(z)$ for real values of $s$. We observe a good agreement between the two functions for various values of $\beta$. Similar, semi-quantitative, degree of agreement is observed for larger values of $\beta$ and also for complex values $s$.

To demonstrate that $\psi_a(z)$ also reproduces the pole structure of $\psi(z)$, we consider a typical term from the analysis in Section \[subsection:hhat\] $$\begin{aligned}
\frac{1}{s+\lambda-\lambda\Delta A(k\epsilon)\psi(z)} = \frac{1}{\lambda \Delta A(k\epsilon)}\frac{z}{\omega - z \psi(z)}.
\label{A_original}\end{aligned}$$ We compute its inverse Laplace transform numerically, using the original function $\psi(z)$, and compare the result with the analytic expression, which we obtain by replacing $\psi(z)$ with $\psi_a(z)$ in the expression above. The latter is straightforward: factorising $\omega \left(7+3z^2\right) - 7 z = 3\omega(z-z_{+})(z-z_{-})$, where $z_{\pm}$ are given in Eq., we obtain $$\begin{aligned}
&\frac{1}{\lambda \Delta A(k\epsilon)}\frac{z}{\omega - z \psi_a(z)} \nonumber \\
&\qquad\qquad = \frac{1}{7}\frac{1}{s+\lambda}\frac{7(s+\lambda)^2 + 3(v_s k)^2}{\left(s+\lambda - \frac{v_s k}{z_{+}}\right)\left(s+\lambda - \frac{v_s k}{z_{-}}\right)}.\end{aligned}$$ Performing the inverse Laplace transform of this expression and introducing the dimensionless time $\tau = \lambda t$ yields $$\begin{aligned}
&e^{-\tau} \Bigg[ 1 +\frac{2}{\sqrt{1-\frac{12\beta^2}{7\Delta^2}}} \exp{\left(\frac{\tau \Delta }{2}\right)} \nonumber \\
& \qquad\qquad\qquad\qquad \times \sinh{\left(\frac{\tau \Delta}{2}\sqrt{1-\frac{12\beta^2}{7\Delta^2}}\right)}\Bigg].
\label{A_test_analytic}\end{aligned}$$ Since both $A(k\epsilon)$ and $\Delta$ take values between $0$ and $1$, for the purpose of comparing to its numerical counterpart, we set $A(k\epsilon)=1$ in the expression above, without loss of generality.

![image](ILTcombined.pdf){width="1.0\linewidth"}

The inverse Laplace transform of the original function Eq. written in terms of the same parameters is given by the Bromwich integral $$\begin{aligned}
\frac{1}{2\pi i} \int_{\gamma-i \infty}^{\gamma+i \infty} d{\tilde s} \frac{e^{\tilde s \tau}}{{\tilde s}+1-\Delta \psi(\beta/(\tilde s+1))},
\label{A_test_numerics}\end{aligned}$$ where $\gamma$ is a real number, chosen to be greater than the real part of any singularity of the integrand [@Doetsch1974]. We perform this integral numerically, using the Gaver–Wynn–Rho algorithm as presented by Valko and Abate [@Valko2005]. Valko and Abate provide an explicit Mathematica function GWR [@math1], which we use here. A Mathematica notebook with the details of this calculation can be found here [@SI].

In Fig.\[ILTcombined\]A we compare Eq. against the numerical Laplace transform of Eq. for $\Delta=0.1$ and $\beta=0.1$, $1$, and $10$. We observe a very good agreement, which is not surprising: At small microswimmer densities, the hydrodynamic interactions between particles affect their dynamics only weakly, and correlations decay as $e^{-\tau}$. This regime does not test the quality of our approximation. A more stringent test is provided, on the other hand, in Fig.\[ILTcombined\]B, were we compare the two Laplace transforms for $\Delta=0.9$. For $\beta<1$, we observe a very good agreement even at such high values of $\Delta$ (close to the mean-field transition). This is the most interesting regime, corresponding to large-scale motion in the suspension, and it is encouraging that our approximation shows quantitative agreement with the numerical data. Note that the black line and the black circles, corresponding to $\beta=0.1$, do not follow $e^{-\tau}$, i.e. our approximation is capable of capturing a non-trivial decay rate. At higher values of $\beta$, corresponding to scales comparable to individual microswimmers, the agreement is semi-quantative, but the overall decay is again close to the tumbling-dominated decay $e^{-\tau}$.

In Appendix \[appendix:dilt\], we assess the quality of our approximation, when used in Eq., which is its ultimate purpose.

Double inverse Laplace transform {#appendix:dilt}
================================

In Section \[subsection:approximateDILT\], we performed the double inverse Laplace transform in Eq. analytically by replacing $\psi(z)$ with $\psi_a(z)$, which led to Eq.. Here, we assess the quality of that approximation by performing the double inverse Laplace transform in Eq. numerically. The relevant part of Eq. reads $$\begin{aligned}
&{\mathcal L}^{-1}_{\tilde{s}_1,\tilde{t}} {\mathcal L}^{-1}_{\tilde{s}_2,\tilde{t}+\tau}
 \frac{1}{1 + \tilde{s}_1}\frac{1}{1 + \tilde{s}_2}\frac{1}{\tilde{s}_1 + \tilde{s}_2} \frac{\tilde{z}_1 \psi(\tilde{z}_1) + \tilde{z}_2 \psi(\tilde{z}_2)}{\tilde{z}_1+\tilde{z}_2}
 \nonumber \\
&\qquad \times\Bigg[ \frac{\tilde{z}_1 \psi(\tilde{z}_1)}{\tilde{\omega} - \tilde{z}_1 \psi(\tilde{z}_1)} + \frac{\tilde{z}_2 \psi(\tilde{z}_2)}{\tilde{\omega} - \tilde{z}_2 \psi(\tilde{z}_2)} \nonumber \\
&\qquad\qquad 
+ \frac{\tilde{z}_1 \tilde{z}_2 \psi(\tilde{z}_1)\psi(\tilde{z}_2)}{\left(\tilde{\omega} - \tilde{z}_1 \psi(\tilde{z}_1)\right)\left(\tilde{\omega} - \tilde{z}_2 \psi(\tilde{z}_2)\right)} \Bigg],
\label{DILTeqnum}\end{aligned}$$ where, in anticipation of performing numerical calculations, we introduced the dimensionless times $\tau = \lambda T$ and ${\tilde t} = \lambda t$, Laplace frequencies ${\tilde s}_{1,2} = s_{1,2}/\lambda$, ${\tilde z}_{1,2} = \beta/(1+{\tilde s}_{1,2})$, and $\tilde{\omega} = \beta/\Delta$, where we absorbed $A(k\epsilon)$ into $\Delta$, as in Appendix \[appendix:psia\]. In what follows, we set $\tilde t=20$ to imitate the limit $\tilde{t}\rightarrow\infty$. The calculations are performed in Mathematics using the combined Fixed-Talbot and Gaver–Wynn–Rho algorithm described by Valko and Abate [@Valko2005]. A Mathematica notebook with the details of this calculation can be found here [@SI]. The results are compared to the relevant part of Eq., recast in the same dimensionless variables $$\begin{aligned}
&e^{-\tau} \Bigg[ -\cos{\left(\sqrt{\frac{3}{7}}\beta\tau\right)} + \frac{e^{\frac{1}{2} \Delta \tau}}{1-\Delta + \frac{3}{7}\beta^2} \Bigg\{ \nonumber \\
&  \qquad\qquad \frac{2-\Delta + \frac{6}{7}\beta^2}{2-\Delta}\cosh{\left(\frac{1}{2}\Delta \tau  \sqrt{1-\frac{12 \beta^2}{7 \Delta^2}}\right)} \nonumber \\
&\qquad\qquad\qquad +\frac{\sinh{\left(\frac{1}{2} \Delta \tau  \sqrt{1-\frac{12 \beta^2}{7 \Delta^2}}\right)}}{\sqrt{1-\frac{12 \beta^2}{7 \Delta^2}}}
\Bigg\}
\Bigg].
\label{DILTeqanal}\end{aligned}$$

![image](DILTcombined.pdf){width="1.0\linewidth"}

The results of the numerical double inverse Laplace transform and its analytical counterpart are shown in Fig.\[DILTcombined\]. As in Appendix \[appendix:psia\], we focus on high values of $\Delta$, which provide the most stringent test of our results. For $\beta\le1$, the analytic approximation agrees quite well with the numerical data, capturing not only the decay rate, but also the oscillatory behaviour, as can be seen from the $\beta=1$ case. These calculations required a very high number of terms, $O(100)$, in the combined Fixed-Talbot and Gaver–Wynn–Rho algorithm [@Valko2005]. For $\beta > 1$, we were unable to obtain converged results for the numerical Laplace transform for any viable number of terms in the numerical algorithm. Nevertheless, the results of Appendix \[appendix:psia\], and the degree of agreement exhibited in Fig.\[DILTcombined\] for the physically most relevant case of $\beta<1$ make us confident that Eq. faithfully reproduces the long-time behaviour of Eq..
