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 