1. Introduction
Understanding the mass transport of colloidal and molecular species in non-equilibrium environments is crucial for various processes, ranging from active intracellular transport (Mogre, Brown & Koslover Reference Mogre, Brown and Koslover2020) to the dispersion of nutrients in world oceans (Katija Reference Katija2012). Apart from their practical importance, the transport properties of tracer particles in generic ‘active baths’ have attracted much interest from a statistical physics perspective, where they can be viewed as a minimal example of particles driven by external, non-equilibrium noise (Argun et al. Reference Argun, Moradi, Pinçe, Bagci, Imparato and Volpe2016; Park et al. Reference Park, Paneru, Kwon, Granick and Pak2020). Beyond the level of tracer particles driven by generic non-equilibrium noise, the archetypical example of a tracer particle in an active bath is a collection of point-like tracers being advected by a set of microswimmers such as bacteria or algae (Lauga & Powers Reference Lauga and Powers2009). When swimming through a viscous fluid, these swimmers create long-range flow fields that advect the tracers, leading to tracer dynamics that is ballistic at short times and diffusive over time scales longer than the autocorrelation time of the local flow field (Lin, Thiffeault & Childress Reference Lin, Thiffeault and Childress2011). Realisations of this system have been studied extensively both experimentally, typically in suspensions of E. coli bacteria (Wu & Libchaber Reference Wu and Libchaber2000; Kim & Breuer Reference Kim and Breuer2004; Drescher et al. Reference Drescher, Dunkel, Cisneros, Ganguly and Goldstein2011; Jepson et al. Reference Jepson, Martinez, Schwarz-Linek, Morozov and Poon2013; Koumakis et al. Reference Koumakis, Lepore, Maggi and Di Leonardo2013; Miño et al. Reference Miño, Dunstan, Rousselet, Clément and Soto2013, Reference Miño, Mallouk, Darnige, Hoyos, Dauchet, Dunstan, Soto, Wang, Rousselet and Clément2011; Patteson et al. Reference Patteson, Gopinath, Purohit and Arratia2016; Peng et al. Reference Peng, Lai, Tai, Zhang, Xu and Cheng2016; Semeraro, Devos & Narayanan Reference Semeraro, Devos and Narayanan2018) or Chlamydomonas algae (Leptos et al. Reference Leptos, Guasto, Gollub, Pesci and Goldstein2009; Yang et al. Reference Yang, Peng, Liu, Tang, Xu and Cheng2016; Ortlieb et al. Reference Ortlieb, Rafaï, Peyla, Wagner and John2019; von Rüling, Kolley & Eremin Reference von Rüling, Kolley and Eremin2021), and theoretically, with microswimmers typically being modelled either as force dipoles acting on the surrounding fluid (Pushkin, Shum & Yeomans Reference Pushkin, Shum and Yeomans2013; Pushkin & Yeomans Reference Pushkin and Yeomans2013; Morozov & Marenduzzo Reference Morozov and Marenduzzo2014; Nordanger, Morozov & Stenhammar Reference Nordanger, Morozov and Stenhammar2022), as spherical ‘squirmers’ with an imposed slip velocity along their body (Thiffeault & Childress Reference Thiffeault and Childress2010; Lin et al. Reference Lin, Thiffeault and Childress2011; Thiffeault Reference Thiffeault2015), or as needle-shaped ‘slender swimmers’ with imposed stresses along their body lengths (Saintillan & Shelley Reference Saintillan and Shelley2012; Krishnamurthy & Subramanian Reference Krishnamurthy and Subramanian2015). While the details of these three microswimmer models differ, the results regarding enhanced tracer diffusion are largely generic and consistent with experimental results, which have shown the swimmer-induced, hydrodynamic diffusivity 
$D_A$ to scale linearly with microswimmer density 
$n$ in the dilute limit where swimmer–swimmer correlations can be neglected (Thiffeault & Childress Reference Thiffeault and Childress2010; Lin et al. Reference Lin, Thiffeault and Childress2011; Miño et al. Reference Miño, Dunstan, Rousselet, Clément and Soto2013). In this limit, a fruitful way of calculating 
$D_A$ is to consider the net displacement due to binary swimmer–tracer scattering events (Pushkin & Yeomans Reference Pushkin and Yeomans2013; Morozov & Marenduzzo Reference Morozov and Marenduzzo2014); two examples of resulting (deterministic) tracer trajectories for scattering events are shown in figure 1. For a non-tumbling swimmer, starting and ending at 
$x = \pm \infty$, the resulting loop is closed, leading to a vanishing tracer net displacement 
$\varDelta$ (figure 1a). For tumbling swimmers with a finite persistence length, the trajectory is, however, punctuated midway through the tracer loop, leading to significantly larger values of 
$\varDelta$ (figure 1b). The resulting value of 
$D_A$ due to a large set of such scattering events can then be obtained by summing over all possible sets of scattering parameters.
Figure 1. Tracer trajectories for infinite and finite swimmer paths. (a) Typical trajectory for a non-diffusing tracer advected by a non-tumbling, point-dipole swimmer following an effectively infinite, straight path. (b) Corresponding trajectory terminated due to a tumbling event. Note that, per (1.1), the effective tracer diffusion 
$D_A$ is independent of tumbling rate 
$\lambda$ for sufficiently high swimming speeds when averaged over all possible swimmer–tracer configurations, even though the net displacement 
$\varDelta$ is much larger for the tumbling swimmer. Tracer trajectories were obtained through direct numerical integration using a non-regularised dipolar flow field, as described by Morozov & Marenduzzo (Reference Morozov and Marenduzzo2014).
In spite of the dependence of 
$\varDelta$ on the microswimmer tumbling rate 
$\lambda$ for single scattering events such as that in figure 1, Pushkin & Yeomans (Reference Pushkin and Yeomans2013) showed that, in the limit of large swimming speeds 
$v_s$, 
$D_A$ is in fact independent of 
$\lambda$ when summed over all possible swimmer trajectories. This result was later generalised by Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020) to arbitrary swimming speeds, leading to the following approximate expression for 
$D_A$:
\begin{equation} D_A \approx \frac{7\kappa^2 n}{2048 \lambda \varepsilon + 336{\rm \pi} v_s}. \end{equation}
Here, 
$\kappa$ is the magnitude of the microswimmer dipole (in units of volume over time), and 
$\varepsilon$ is a characteristic linear size of the microswimmer, which we take to be equal to the short-range regularisation length of the dipolar flow field; unless stated otherwise, in the following we will use 
$\kappa$ and 
$\varepsilon$ to non-dimensionalise the numerical data.
Rather than using the scattering approach outlined above, (1.1) was derived by formulating a kinetic theory for the spatiotemporal correlations of the disturbance velocity field 
${\boldsymbol U}$ created by the swimmers. Due to the linearity of Stokes flow, we can write 
${\boldsymbol U}$ as a superposition of the individual swimmer flow fields 
${\boldsymbol u}_s$:
\begin{equation} {\boldsymbol U}({\boldsymbol r},t) = \sum_{i=1}^{N} {\boldsymbol u}_s({\boldsymbol r};{\boldsymbol r}_i,{\boldsymbol p}_i), \end{equation}
where 
${\boldsymbol r}_i$ and 
${\boldsymbol p}_i$ are, respectively, the position and orientation of swimmer 
$i$. Knowing the statistical properties of 
${\boldsymbol U}$, 
$D_A$ can be calculated via the Green–Kubo relation
\begin{align} D_A = \frac{1}{3} \int_{0}^{\infty}\langle \dot{{\boldsymbol r}}_T(t) \boldsymbol{\cdot} \dot{{\boldsymbol r}}_T(0)\rangle \, {\rm d} t = \frac{1}{3} \int_{0}^{\infty}\langle {\boldsymbol U}({\boldsymbol r}_T,t) \boldsymbol{\cdot} {\boldsymbol U}({\boldsymbol r}_T,0) \rangle \, {\rm d} t \equiv \frac{1}{3} \int_{0}^{\infty}C_T(t) \, {\rm d} t, \end{align}
where, in the second equality, we have assumed point-like tracers advected by the disturbance flow, so that 
$\dot {{\boldsymbol r}}_T = {\boldsymbol U}({\boldsymbol r}_T)$, and the third equality defines the velocity autocorrelation function 
$C_T(t)$ in the co-moving tracer frame. While yielding results for dilute suspensions identical to those for the scattering approach discussed above, kinetic theories are more readily extended to accommodate the effect of swimmer–swimmer correlations due to the mutual advection and reorientation of swimmers (Škultéty et al. Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020). Importantly, these interactions break the symmetry between rear-actuated (pusher) swimmers, such as most bacteria, and front-actuated (puller) ones, such as Chlamydomonas, leading to a super-linear scaling of 
$D_A$ with 
$n$ for pushers, and a corresponding sub-linear scaling for pullers (Stenhammar et al. Reference Stenhammar, Nardini, Nash, Marenduzzo and Morozov2017).
Equation (1.1) shows two qualitatively different regimes for high and low 
$v_s$: for 
$v_s \!\rightarrow\! 0$ – the so-called shaker limit – the dominant mechanism controlling the decorrelation of 
$C_T(t)$ is tumbling of the swimmer. For fast swimmers, with 
$v_s \gg \lambda \varepsilon$, the decorrelation of the fluid velocity is instead dominated by the swimmer self-propulsion and is thus independent of 
$\lambda$, and (1.1) reduces to the expression derived by Pushkin & Yeomans (Reference Pushkin and Yeomans2013). A third, somewhat less explored, mechanism affecting 
$C_T(t)$ is Brownian translational diffusion of the tracer: even though Brownian diffusion does not affect the statistics of the flow field 
${\boldsymbol U}({\boldsymbol r},t)$ as measured in the lab frame, the positional noise of the tracer particle will cause it to cross the streamlines of the disturbance flow, thus perturbing its trajectory compared to the athermal case shown in figure 1, leading to a lower 
$D_A$. Since the Brownian diffusion constant 
$D_0$ depends inversely on the tracer radius per the Stokes–Einstein relation, the magnitude of this effect is expected to be significant primarily for small tracer particles, and it has been hypothesised to explain the non-monotonic size dependence of enhanced tracer diffusion seen in experiments with colloids in E. coli suspensions, where Patteson et al. (Reference Patteson, Gopinath, Purohit and Arratia2016) observed a maximum in 
$D_A$ for a tracer radius of approximately 5 
$\mathrm {\mu }$m. Furthermore, the effect of Brownian motion on enhanced tracer diffusion was analysed theoretically by Kasyap, Koch & Wu (Reference Kasyap, Koch and Wu2014) for a model of slender-body swimmers, showing that 
$D_A$ is a non-monotonic function of 
$D_0$, with 
$D_A$ first showing a small increase for intermediate 
$D_0$, before falling below the athermal value as 
$D_0$ grows larger. In a more recent study of finite-size spherical tracers in microswimmer suspensions, Dyer & Ball (Reference Dyer and Ball2021) analysed numerically the combined effect of thermal fluctuations and near-field flows on the size-dependent tracer dynamics, finding non-monotonic behaviour similar to that observed experimentally. In this paper, we will revisit the problem of the interplay between Brownian and hydrodynamic diffusion for the simple case of point-like tracers immersed in a dilute suspension of microswimmers described via a regularised dipolar flow field. Using kinetic theory and large-scale lattice Boltzmann (LB) simulations of E. coli-like suspensions, we show that the effect of Brownian diffusion on active diffusion is practically relevant only whenever 
$v_s < \lambda \varepsilon$, which corresponds to extremely slow (or frequently tumbling) swimmers. For biologically relevant values of 
$v_s$ and 
$\lambda$, swimming is instead the dominant decorrelation mechanism, so that 
$D_A$ becomes independent of both 
$D_0$ and 
$\lambda$. In contrast to previous studies, our results thus indicate that the effect of Brownian motion on the enhanced diffusion is in fact negligible for most microswimmer realisations, and thus that the experimentally observed non-monotonic size dependence of 
$D_A$ on tracer size reported by Patteson et al. (Reference Patteson, Gopinath, Purohit and Arratia2016) has other explanations.
2. Model and method
We consider a collection of 
$N$ non-interacting microswimmers at number density 
$n = N/V$ moving through a three-dimensional viscous fluid of viscosity 
$\mu$. Each microswimmer is composed of two equal and opposite point forces of magnitude 
$F$ separated by a length 
$\ell$, and swims with a constant speed 
$v_s$. The resulting reduced hydrodynamic dipole strength is 
$\kappa = F\ell /\mu$. Furthermore, the swimming direction 
${\boldsymbol p}_i$ relaxes through Poisson-distributed random tumbles with uncorrelated directions, occurring with average frequency 
$\lambda$.
The position 
${\boldsymbol r}_T$ of a point-like tracer obeys the equation of motion
\begin{equation} \dot{{\boldsymbol r}}_T = {\boldsymbol U}({\boldsymbol r}_T) + \sqrt{2D_0}\,\boldsymbol{\eta}, \end{equation}
where 
$\boldsymbol {\eta }$ is a unit-variance white noise, 
$\delta$-correlated in time, and 
$D_0$ is the Brownian diffusion constant. Thus, the effect of Brownian motion is fully contained in the tracer dynamics, while we assume the effects of thermal fluctuations on the fluid and on the pairwise swimmer–tracer dynamics to be subdominant. The fluid disturbance velocity 
${\boldsymbol U}({\boldsymbol r}_T)$ due to the presence of all swimmers is given by (1.2) and can, in principle, be summed up explicitly on each time step. However, to avoid the (prohibitively costly) pairwise summation over all swimmers and tracers, we instead solve numerically for the flow field using an efficient point-force implementation of the LB method described previously (Nash, Adhikari & Cates Reference Nash, Adhikari and Cates2008; Bárdfalvy et al. Reference Bárdfalvy, Anjum, Nardini, Morozov and Stenhammar2020). Our simulations comprise a system with periodic boundaries and size 
$100^3$ lattice units. In LB units, set by the LB lattice spacing 
${\rm \Delta} l$ and time step 
${\rm \Delta} t$, the microswimmer density was kept constant at 
$n = 0.01$ unless stated otherwise, corresponding to 
$N = 10^4$ microswimmers. Furthermore, 
$N_T = 5 \times 10^4$ point tracers were included for statistical averaging. In addition to the LB simulations, in § 3 we will furthermore extend the kinetic theory developed previously by Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020) to the case of microswimmers undergoing Brownian diffusion, and show that in the limit of non-interacting swimmers, the derived expression for 
$D_A$ is equivalent to that for Brownian tracers in a suspension of non-Brownian microswimmers.
Inserting (2.1) into the Green–Kubo relation (1.3) yields
\begin{equation} D_{Tot} = D_0 + \frac{1}{3} \int_{0}^{\infty} C_T(t) \, {\rm d} t = D_0 + D_A, \end{equation}
where 
$C_T$ was defined in (1.3). Thus, to obtain 
$D_A$, we evaluate numerically the time correlation of the disturbance velocity measured in the co-moving tracer frame. Since the tracer position 
${\boldsymbol r}_T$ depends on 
$D_0$, 
$C_T(t)$, and thus 
$D_A$, will depend implicitly on 
$D_0$.
To characterise the system, we will use three dimensionless quantities. First, we define the Péclet number, which measures the relative importance of active and thermal forces, as
\begin{equation} {Pe} \equiv \frac{D_A(D_0=0)}{D_0}, \end{equation}
where 
$D_A(D_0=0)$ is the active diffusivity of a tracer immersed in an equivalent microswimmer suspension but in the absence of Brownian tracer motion. It should be noted that our definition of 
$Pe$ is qualitatively different from that of Kasyap et al. (Reference Kasyap, Koch and Wu2014), who, instead of 
$D_A(D_0=0)$, use the swimming speed 
$v_s$ to characterise the active forces. We, however, argue that the activity experienced by the tracers depends on the magnitude of the velocity fields generated by the swimmers, and is thus dependent on 
$\kappa$ and encoded in 
$D_A(D_0 = 0)$, while 
$v_s$ is instead a measure of the swimmers’ self-propulsion. For experimental realisations of microswimmers, 
$\kappa$ and 
$v_s$ are directly proportional to each other; however, the specific relation between them will nevertheless be specific to each type (or species) of swimmer, and decoupling them enables us to study separately the effects of self-propulsion and fluid forcing, as we demonstrate further below.
Second, we measure the change in active diffusion due to Brownian motion through the quantity
\begin{equation} \xi \equiv \frac{D_A(D_0)}{D_A(D_0 = 0)}. \end{equation}
In the limit 
${Pe} \rightarrow \infty$, where Brownian motion becomes negligible, we thus expect that 
$\xi \rightarrow 1$ as 
$D_A$ approaches its non-Brownian value. Finally, in accordance with Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020), we account for the effect of microswimmer self-propulsion using the reduced swimmer persistence length 
$L$, defined by
\begin{equation} L \equiv \frac{v_s}{\varepsilon \lambda}. \end{equation}3. Kinetic theory
In this section, we will outline the main steps in the derivation of 
$D_A$ for a suspension of Brownian tracer particles immersed in a dilute microswimmer suspension, whose dynamics are governed by (2.1). Just as in our previous works (Stenhammar et al. Reference Stenhammar, Nardini, Nash, Marenduzzo and Morozov2017; Škultéty et al. Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020), we describe the flow field measured at 
${\boldsymbol r}$ due to a swimmer with position 
${\boldsymbol r}_i$ and orientation 
${\boldsymbol p}_i$ by a regularised dipolar flow field 
${\boldsymbol u}_s({\boldsymbol r})$:
\begin{equation} {\boldsymbol u}_s({\boldsymbol r};{\boldsymbol r}_i,{\boldsymbol p}_i) = \frac{\kappa}{8{\rm \pi}} \left[ 3\,\frac{({\boldsymbol p}_i \boldsymbol{\cdot} {\boldsymbol r}')^2 {\boldsymbol r}' + \varepsilon^2 ({\boldsymbol p}_i \boldsymbol{\cdot} {\boldsymbol r}'){\boldsymbol p}_i }{(r'^2 + \varepsilon^2)^{5/2}} - \frac{{\boldsymbol r}'}{(r'^2 + \varepsilon^2)^{3/2}} \right], \end{equation}
where 
${\boldsymbol r}' = {\boldsymbol r}-{\boldsymbol r}_i$, 
$r' = |{\boldsymbol r}'|$, and 
$\varepsilon$ is the regularisation length. Our starting point is the derivation of Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020), where we formulated and solved a kinetic theory describing the fluctuations of the velocity field 
${\boldsymbol U}({\boldsymbol r},t)$ due to a superposition of single-swimmer flow fields. In the limit of non-interacting swimmers, which is the case that we consider here, the temporal correlations of the steady-state velocity field 
${\boldsymbol U}$ measured in the lab frame, 
$C_U(t) \equiv \langle {\boldsymbol U}({\boldsymbol r},t) \boldsymbol {\cdot } {\boldsymbol U}({\boldsymbol r},0) \rangle$, are given by
\begin{equation} C_U(t) = \frac{\kappa^2 n}{15 {\rm \pi}^2 \varepsilon} \int_0^{\infty} A^2(\zeta) \, {\rm e}^{-\tau}\, \mathcal{F}(L\zeta\tau) \, {\rm d} \zeta, \end{equation}where
\begin{equation} A(x) = \tfrac{1}{2} x^2\,K_2(x), \end{equation}
with 
$K_2$ being the modified Bessel function of the second kind, and
\begin{equation} \mathcal{F}(x) = 15\,\frac{(5x^2-12)\sin x - x(x^2-12)\cos x}{x^5}, \end{equation}
defined such that 
$\mathcal {F}(0) = 1$. Furthermore, we used the dimensionless variables 
${L = v_s/(\lambda \varepsilon )}$, 
$\tau = t\lambda$ and 
$\zeta = k\varepsilon$, where 
$k = |{\boldsymbol k}|$ is the wavevector magnitude. Equation (3.2) can be expressed equivalently in closed form in terms of elliptic integrals, as given by (72) of Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020). Equation (3.2) contains decorrelation of the flow field due to two separate mechanisms: exponential decay of 
$C_U(t)$ due to tumbling, and a more complex, oscillatory behaviour due to swimming, encoded in the function 
$\mathcal {F}$.
To obtain the hydrodynamic diffusivity of a passive tracer, Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020) used a stationary tracer approximation, implying that the tracer advection by the swimmer flow field is negligible compared to the self-propulsion of the swimmer. This implies that 
${\boldsymbol r}_T$ remains effectively constant over the time it takes for 
${\boldsymbol U}$ to relax, such that 
$C_T(t) = \langle {\boldsymbol U}({\boldsymbol r}_T[t],t) \boldsymbol {\cdot } {\boldsymbol U}({\boldsymbol r}_T[t=0],0) \rangle \approx \langle$ 
${\boldsymbol U}({\boldsymbol r}_T[t=0],t) \boldsymbol {\cdot } {\boldsymbol U}({\boldsymbol r}_T[t=0],0) \rangle = C_U(t)$, where, in the last equality, we have made the additional assumption that the tracers are distributed homogeneously in space so that they sample an unweighted spatial average of the flow field. Thus if we can replace the correlation function 
$C_T$ in the co-moving tracer frame with that in the stationary lab frame, 
$C_U$, then we can insert (3.2) into the Green–Kubo relation (1.3) and integrate over time to yield the following expression for 
$D_A$, identical to (85) of Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020):
\begin{equation} D_A = \frac{\kappa^2 n}{45 {\rm \pi}^2 \lambda \varepsilon} \int_0^{\infty} A^2(\zeta)\, \mathcal{G} ( L\zeta ) \, {\rm d} \zeta, \end{equation}where
\begin{equation} \mathcal{G}(x) = \frac{5}{2}\,\frac{3x + 2x^3 -3(1+x^2)\arctan x}{x^5}, \end{equation}
defined such that 
$\mathcal {G}(0) = 1$. Furthermore, by matching the asymptotic behaviours for 
$L\rightarrow 0$ and 
$L \rightarrow \infty$, (3.5) can be approximated by the simple expression given in (1.1)
The generalisation of (3.5) to the case of Brownian tracers might seem straightforward, but unfortunately is not. Since the approximation 
$C_T(t) = C_U(t)$ amounts to the tracer remaining effectively stationary during a swimmer–tracer scattering event, this approximation will, by construction, not capture any effects on 
$D_A$ coming from Brownian diffusion across streamlines. This can be realised easily by noticing that 
$C_U(t)$ is solely a property of the swimmer suspension, and will be strictly unaffected by the tracer dynamics; thus 
$D_A$ in (3.5) remains unaffected by the inclusion of tracer diffusion. Instead of going beyond the stationary tracer approximation, we circumvent this problem by noticing that in a suspension of non-interacting microswimmers, the dynamics of a Brownian tracer will be statistically identical to that of a non-Brownian tracer in a suspension of Brownian swimmers with the same translational diffusivity 
$D_0$. This is because the single-swimmer flow field in (3.1) depends solely on the separation vector 
${\boldsymbol r}_T - {\boldsymbol r}_i$, implying that the flow field experienced by a diffusing tracer (i.e. noise acting on 
${\boldsymbol r}_T$) is identical to that experienced by a non-Brownian tracer sampling the flow field from a swimmer with the same noise instead applied to 
${\boldsymbol r}_i$. In a non-interacting microswimmer suspension, this equivalence is exact as long as the noise has zero mean and identical spectral properties, and we verify it numerically in figure 3. However, it breaks down as soon as swimmer–swimmer correlations become significant, since swimmer diffusion will then affect the magnitude of such correlations, which tracer diffusion will not.
Thus, as outlined in Appendix A, we instead calculate 
$C_U(t)$ for the case of a suspension of diffusing swimmers, described by the dynamics
\begin{equation} \dot{{\boldsymbol r}}_i = v_s {\boldsymbol p}_i + \sqrt{2D_0}\,\boldsymbol{\eta},\end{equation}including the same tumbling dynamics as before. This yields the following generalisation of (3.2):
\begin{equation} C_U(t;D_0) = \frac{\kappa^2 n}{15 {\rm \pi}^2 \varepsilon} \int_0^{\infty} A^2(\zeta) \, {\rm e}^{-(1+\tilde{D}\zeta^2) \tau}\, \mathcal{F}(L\zeta\tau) \, {\rm d} \zeta, \end{equation}
where we have additionally defined the non-dimensional diffusivity 
$\tilde {D} = D_0/(\lambda \varepsilon ^2)$. Since the effect of Brownian diffusion is now fully incorporated into the properties of 
$C_U$, we again use the stationary-tracer approximation and insert this expression into the Green–Kubo relation (1.3), leading to the following expression for 
$D_A$ in the presence of Brownian diffusion:
\begin{equation} D_A(D_0) = \frac{\kappa^2 n}{45 {\rm \pi}^2 \lambda \varepsilon} \int_0^{\infty} \frac{A^2(\zeta)}{1 + \tilde{D}\zeta^2}\,\mathcal{G}\left( \frac{L\zeta}{1 + \tilde{D}\zeta^2} \right) {\rm d} \zeta. \end{equation}In § 4, we evaluate (3.8) and (3.9) numerically, and compare the results with direct numerical simulations of microswimmer suspensions.
4. Results and discussion
In figure 2, we begin by verifying (3.5) for the hydrodynamic diffusion coefficient 
$D_A$ in the limit 
$D_0 = 0$. In figure 2(a), we demonstrate the 
$\lambda ^{-1}$ dependence of 
$D_A$ in the shaker limit 
$v_s = 0$, while figure 2(b) shows its more complex dependence on 
$v_s$ for constant 
$\lambda$. Apart from numerically verifying the kinetic theory expression (3.5), these results illustrate how 
$D_A$ decreases abruptly due to the temporal decorrelation of the flow field 
${\boldsymbol U}({\boldsymbol r}_T,t)$ induced respectively by tumbling and swimmer self-propulsion. The slight deviation between the theoretical curve and simulation results at small 
$v_s$ in figure 2(b) is likely due to the specific form for the short-range regularisation, which becomes important as 
$v_s \rightarrow 0$. In our derivation of (3.5), we use the regularised flow field (3.1), based on the regularisation first introduced by Cortez, Fauci & Medovikov (Reference Cortez, Fauci and Medovikov2005). In the LB simulations, we instead use a numerical interpolation scheme based on a regularisation of the 
$\delta$ function (Peskin Reference Peskin2002) acting separately on the two point forces that make up each microswimmer. Unlike the expression in (3.1), this numerical regularisation does not allow a direct mapping (or adjustment) of the regularisation length 
$\varepsilon$. We thus do not expect perfect agreement between kinetic theory and simulation in the low-
$v_s$ regime where the short-range regularisation becomes important, and therefore treat 
$\varepsilon$ as a fitting parameter when comparing data from LB simulations with kinetic theory predictions. However, we find that the fitted value of 
$\varepsilon$ varies only slightly (
$\varepsilon \in [1.9\,{\rm \Delta} l, 2.5\,{\rm \Delta} l]$) for the values of 
$v_s$ used throughout this work, in good accordance with the fact that the regularised 
$\delta$ function is interpolated over a support of 
$2\,{\rm \Delta} l$ in each Cartesian direction; for a more in-depth discussion of the effect of the interpolation scheme on the tracer dynamics, see de Graaf & Stenhammar (Reference de Graaf and Stenhammar2017).
Figure 2. Effective tracer diffusion in the absence of Brownian motion. Panel (a) shows the 
$\lambda ^{-1}$ dependence of 
${D_A(D_0=0)}$ for shakers with 
$v_s = 0$, and panel (b) its 
$v_s$ dependence at constant 
$\lambda = 10^{-4}$. Symbols denote simulation results, and solid lines show results from (3.5) using 
$\varepsilon$ as a fitting parameter. Error bars represent one standard deviation as obtained from averaging over four separate runs with different initial conditions. The results are non-dimensionalised in terms of 
$\kappa$ and 
$\varepsilon$.
In figure 3, we study the additional effect of varying the Brownian diffusion coefficient 
$D_0$, as encoded in (3.9). From the data in figure 3(a), it is clear that for large enough 
$D_0$, the active diffusivity 
$D_A$ decreases compared to its non-Brownian value. To enable an easier analysis of the effect of varying swimming speed, in figure 3(b) we present the same data instead plotted as a function of the reduced variables 
$\xi$ and 
$Pe$. For shakers with 
$L = 0$ (blue curve in figure 3b), 
$D_A$ is reduced compared to its non-Brownian value (
$\xi < 1$) as soon as 
${Pe} < 1$, reaching a value as low as 
$\xi = 0.2$ for 
${Pe} \approx 10^{-3}$. However, for finite values of 
$v_s$, this effect on 
$D_A$ occurs for gradually lower values of 
$Pe$; for the fastest swimmers considered here, with 
$L = 4.0$, no significant reduction of 
$D_A$ is observed even for 
$Pe$ as low as 
$10^{-3}$. Instead, we observe a small but significant increase in the active diffusion compared to its non-Brownian value, in accordance with what was observed previously for slender swimmers by Kasyap et al. (Reference Kasyap, Koch and Wu2014); we discuss this effect further below. Crucially, a reduced persistence length of 
$L = 4$ nevertheless corresponds to relatively slow swimming from a biological perspective: According to the approximate calculation in Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020), 
$L$ for E. coli bacteria lies somewhere in the range between 5 and 20, indicating that the effect of Brownian motion on active diffusion is likely negligible in suspensions of swimming bacteria due to their fast self-propulsion. In our LB simulations, studying values higher than 
$L \approx 4$ is challenging, as these large swimming speeds both require very large systems to avoid significant finite-size effects and yield artefacts due to the effect of finite Reynolds number in the swimmer–tracer scattering dynamics (de Graaf & Stenhammar Reference de Graaf and Stenhammar2017). We nevertheless studied 
$\xi ({Pe})$ numerically using (3.9) for larger values of 
$L$, verifying that both the peak and the subsequent decrease in 
$\xi$ continue to move to even lower values of 
$Pe$ as 
$L$ is increased.
Figure 3. Brownian motion suppresses active diffusion for slow swimming speeds. (a) Values of 
$D_A$ measured from LB simulations (diamonds) and calculated from (3.9) (solid lines), both expressed in LB units. (b) The same data but expressed in the dimensionless quantities 
$\xi$ and 
$Pe$. For very slow swimmers with 
$L \lesssim 1$, 
$D_A$ is reduced compared to the non-Brownian value (
$\xi = 1$) whenever 
${Pe} < 1$, while for faster swimmers, significantly lower values of 
$Pe$ are necessary to affect 
$D_A$. The circles for 
$L = 0$ and 2.4 correspond to the hydrodynamic diffusion of non-Brownian tracers measured in a suspension of Brownian swimmers of the same 
$D_0$, verifying the statistical equivalence between tracer and swimmer diffusion in the non-interacting limit. Error bars represent one standard deviation as obtained from averaging over four separate runs with different initial conditions. The results in (a) are non-dimensionalised using 
$\kappa$ and 
$\varepsilon$.
To understand further the mechanism behind the reduction in 
$D_A$ with 
$D_0$, we consider the two autocorrelation functions 
$C_T(t)$ and 
$C_U(t)$, which measure the fluid autocorrelation in the co-moving tracer frame and in the lab frame, respectively. Figure 4 shows these correlation functions for 
$L = 0$ and 
$L = 2.4$, with figures 4(a,b) corresponding to LB results for 
$C_T$, and figures 4(c,d) to kinetic theory results from (3.8) for 
$C_U$. First, we notice that the two sets of curves are very similar, implying that the stationary tracer approximation 
$C_T \approx C_U$ is indeed accurate. Second, we notice that the decay of the correlation function is significantly faster for swimmers than for shakers, again illustrating that self-propulsion acts an efficient decorrelation mechanism for 
${\boldsymbol U}$. The effect of finite 
$D_0$ for shakers (figures 4a,c) is simply to decrease the relaxation time of the exponential decay, in accordance with the 
$L = 0$ limit of (3.8). For swimmers, the situation is, however, more complex: For short times, the flow field decays faster with decreasing 
$Pe$, while the long-time tail of 
$C_T$ and 
$C_U$ instead becomes somewhat more extended with decreasing 
$Pe$. For fast enough swimmers, the latter effect leads to the local maximum at 
$\xi > 1$ for intermediate 
$Pe$ observed in figure 3(b) for the two highest 
$L$. Finally, we note that the equal-time fluid velocity variance 
$\langle U^2({\boldsymbol r}_T) \rangle$, corresponding to the 
$t = 0$ values of 
$C_T$ and 
$C_U$, is independent of 
$Pe$. This means that, regardless of the ratio between diffusive and active motion, the tracer particles sample the overall flow field homogeneously. This fact is non-trivial, since dry active particles that move autonomously on a solid substrate are known to preferentially sample regions where they move slowly (Stenhammar et al. Reference Stenhammar, Wittkowski, Marenduzzo and Cates2016). Our results thus highlight that this generic mechanism is absent for inertialess point tracers advected by an incompressible fluid. However, for microswimmer systems where entrainment due to tracers being captured by the near-field flows of passing swimmers (Jeanneret et al. Reference Jeanneret, Pushkin, Kantsler and Polin2016), we would expect the fluid flow sampled by tracers to be significantly different from the average flow field in the system.
Figure 4. Brownian motion decorrelates tracer trajectories. (a,b) Fluid velocity autocorrelation 
$C_T(t)$ in the co-moving tracer frame measured from LB simulations for (a) shakers (
$L=0$), and (b) swimmers with 
$L = 2.4$ at indicated values of 
$Pe$. The dotted line shows the correlation function 
$C_U(t)$ of the fluid velocity in the lab frame, demonstrating that the stationary-tracer approximation 
$C_T(t) \approx C_U(t)$ is excellent in the absence of Brownian tracer diffusion (
$Pe \rightarrow \infty$). (c,d) Corresponding lab-frame correlation functions 
$C_U(t)$ obtained from kinetic theory (3.9) for a suspension of diffusing swimmers, as described in § 3. Throughout, 
$C_T$ and 
$C_U$ are non-dimensionalised using 
$\kappa$ and 
$\varepsilon$.
In figure 5, we study the dependence of the suppression of active diffusion on the microswimmer density 
$n$ in the shaker limit 
$L=0$. At first sight, this dependence might appear trivial, since 
$D_A$ is well known to be linearly dependent on 
$n$ (Jepson et al. Reference Jepson, Martinez, Schwarz-Linek, Morozov and Poon2013; Miño et al. Reference Miño, Dunstan, Rousselet, Clément and Soto2013) in the limit of non-interacting swimmers, a fact that is unaffected by Brownian motion as shown by (3.9). Since 
$\xi$ measures the ratio between the Brownian and non-Brownian values of 
$D_A$, one would naively expect 
$\xi$ to be independent of 
$n$. However, since the Péclet number itself, as defined in (2.3), increases with 
$n$ for constant 
$D_0$, the relative effect of Brownian motion on 
$D_A$ is in fact a complex function of 
$n$ even for non-interacting swimmers, as shown in figure 5. More specifically, according to figure 5(b), the suppression of active diffusion becomes more significant with increasing microswimmer density. In physical units, the highest concentration considered here (
${n\varepsilon ^3 \approx 0.08}$) corresponds approximately to a bacterial concentration of 
$10^9$ ml
$^{-1}$, which is somewhat higher than the highest E. coli concentration considered by Jepson et al. (Reference Jepson, Martinez, Schwarz-Linek, Morozov and Poon2013), but still within the range of concentrations where swimmer–swimmer correlations are reasonably small (Stenhammar et al. Reference Stenhammar, Nardini, Nash, Marenduzzo and Morozov2017). In summary, our results thus show that for Brownian diffusion to have any measurable influence on hydrodynamic diffusion, it is necessary to create a system with a relatively high density of very slow microswimmers; as we discuss in § 5, this set of parameters is likely not achievable for suspensions of biological microswimmers.
Figure 5. Reduction of 
$D_A$ varies with microswimmer density. (a) Plots of 
$\xi$ as a function of 
$Pe$ for shakers (
$L=0$) at various densities 
$n$, as indicated. (b) Plot of 
$\xi$ as a function of 
$n$, at a fixed value of 
${Pe} = 0.013$. Simulation data are given by the symbols, with error bars obtained from averaging over four separate runs with different initial conditions, while solid lines are computed from (3.9).
5. Summary and conclusions
In this study, we have demonstrated a number of theoretical and computational results regarding the effect of Brownian diffusion on the swimmer-induced hydrodynamic diffusion of tracer particles in a suspension of dipolar microswimmers. Our key finding is that the effect of Brownian diffusivity 
$D_0$ on the activity-induced, hydrodynamic diffusivity 
$D_A$ is significant only when the Péclet number as defined by (2.3) is below unity, meaning that Brownian diffusivity needs to dominate hydrodynamic diffusivity. However, the necessary requirement 
${Pe} < 1$ is sufficient only in the shaker limit 
$v_s \rightarrow 0$: for swimmers with persistence lengths larger than the organism size, significantly lower values of 
$Pe$ are required to perturb the tracer trajectories sufficiently to affect 
$D_A$. This conclusion is analogous to the independence of 
$D_A$ from the tumbling rate 
$\lambda$ for large 
$v_s$ illustrated in (1.1). Whenever 
$v_s$ is large, the decorrelation of 
${\boldsymbol U}$ by swimming will dominate the decorrelation due to tumbling and translational diffusion, and the dependence on 
$\lambda$ and 
$D_0$ will thus vanish in the limit 
$v_s \rightarrow \infty$. While this effect is expected, what is perhaps surprising is the rather moderate values of 
$L = v_s / (\varepsilon \lambda )$ necessary to render the coupling between 
$D_0$ and 
$D_A$ negligible, as illustrated in figure 3(b). To put these values into perspective, we use the conservative estimate 
$L = 5$ for E. coli. By virtue of figure 3(b), we require that 
${Pe} \leq 10^{-3}$ for Brownian motion to have an effect of 
${\sim }5$ % on 
$D_A$. Using as an example the minimum value 
$D_A \approx 10^{-2}\ \mathrm {\mu }{\rm m}^2\ {\rm s}^{-1}$ measured by Jepson et al. (Reference Jepson, Martinez, Schwarz-Linek, Morozov and Poon2013) in a three-dimensional E. coli suspension, this requirement thus implies that 
$D_0 \geq 10\ \mathrm {\mu }{\rm m}^2\ {\rm s}^{-1}$, which by virtue of the Stokes–Einstein relation corresponds to a tracer radius 
$R_0 \sim 20$ nm. While significantly smaller than the particles used in typical measurements on colloidal tracers (Leptos et al. Reference Leptos, Guasto, Gollub, Pesci and Goldstein2009; Miño et al. Reference Miño, Dunstan, Rousselet, Clément and Soto2013; Patteson et al. Reference Patteson, Gopinath, Purohit and Arratia2016), this diffusion coefficient is close to the value of 
$D_0$ measured for dextran in E. coli suspensions by Kim & Breuer (Reference Kim and Breuer2004). It is also fully feasible to realise such low Péclet numbers for micron-sized spheres by instead decreasing the bacterial density to very low values; however, measuring the correction to 
$D_A$ for 
${Pe} \sim 10^{-3}$ represents a major difficulty, since it amounts to measuring a 
$\sim$5 % deviation of an effective diffusivity that is itself a thousand times smaller than the Brownian diffusion. It would thus require an extremely accurate determination of 
$D_0$, which then needs to be subtracted from the total measured diffusion constant to determine 
$D_A$. Obtaining this accuracy in a colloidal suspension would be very challenging due to particle polydispersity, interactions with boundaries, temperature gradients, and other system-specific complications. We thus conclude that, for typical 
$L$ values relevant for biological microswimmers, our results imply that the effect of Brownian motion on 
$D_A$ is likely negligible for all practical purposes. To observe experimentally the reduction in 
$D_A$, one would instead need to study a system of dipolar shakers, which stir up the surrounding fluid without self-propelling. While this is a somewhat exotic type of system, it could potentially be realised by anchoring molecular motors or biological microswimmers to a surface. In a biological setting, the shaker limit furthermore resembles previously developed models of enzymes anchored to lipid bilayers that induce dipolar flows through cyclical conformation changes (Hosaka, Komura & Mikhailov Reference Hosaka, Komura and Mikhailov2020).
Our results are qualitatively consistent with the previous theoretical results obtained by Kasyap et al. (Reference Kasyap, Koch and Wu2014), including the non-monotonic behaviour of 
$D_A$ with respect to 
$D_0$ at high swimming speed. Their results, however, differ in three important ways. First, they consider a ‘slender swimmer’ model, where the fluid is forced via a stress applied along a line representing the swimmer body, rather than by two point forces. While this model should lead to a dipolar flow in the far field, the near-field differences are significant, as illustrated in figure 8 of Kasyap et al. (Reference Kasyap, Koch and Wu2014). Furthermore, they consider only the fast-swimming limit where, according to (1.1), 
$D_A$ for dipolar microswimmers is independent of 
$\lambda$. Finally, they parametrise their model in a qualitatively different way to us: as discussed above, they adopt a definition of 
${Pe}$ based on 
$v_s$ rather than on the tracer diffusivity as in (2.3). In this description, 
$v_s$ is furthermore directly coupled to 
$\kappa$, so that changing 
$Pe$ simultaneously changes the activity of the bath (via 
$\kappa$) and the swimming speed 
$v_s$. These differences make it difficult to compare directly with our results, as we consider the effects of fluid advection and self-propulsion separately via 
$Pe$ and 
$L$, respectively.
Thus, even though Brownian motion is unlikely to provide a significant dependence of 
$D_A$ on 
$R_0$ for tracers immersed in suspensions of biological microswimmers, there are several other mechanisms that need to be studied to explain the non-monotonic dependence observed experimentally (Patteson et al. Reference Patteson, Gopinath, Purohit and Arratia2016) and computationally (Dyer & Ball Reference Dyer and Ball2021). First, the effect of tracer entrainment by the near field flows of the swimmer is strongly dependent on the size ratio between the swimmer and the tracer (Jeanneret et al. Reference Jeanneret, Pushkin, Kantsler and Polin2016), although we expect this term to be small for micron-sized tracer particles in E. coli suspensions. Second, the finite size of the tracer will change the equation of motion (2.1) into the Faxén equation that takes into account the nonlinearity of the flow field (Kim & Karrila Reference Kim and Karrila1991), an effect that was included implicitly in the wavelet Monte Carlo simulations by Dyer & Ball (Reference Dyer and Ball2021), and, together with tracer entrainment, is a significant explanation of their observed 
$R_0$ dependence of 
$D_A$. Finally, non-hydrodynamic interactions such as direct collisions, electrostatic interactions and artefacts due to container walls are likely to depend in a non-trivial manner on the tracer size for each system in question. We thus conclude by noting that significant further experimental and theoretical work is necessary to disentangle the system-specific properties from the generic properties of tracer dynamics in microswimmer suspensions.
Acknowledgements
Discussions with C. Nardini are kindly acknowledged. The computations were enabled by resources provided by LUNARC. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.
Funding
This work was financed through the Knut and Alice Wallenberg Foundation (project grant KAW 2014.0052). J.S. acknowledges financial support from the Swedish Research Council (project no. 2019-03718).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Kinetic theory for 
$D_A(D_0)$
In this appendix, we demonstrate the main steps in the derivation of (3.8). The analysis follows closely a similar derivation presented in Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020), albeit with two major differences. First, in addition to that work, we include the effect of microswimmer Brownian diffusivity, as discussed in § 3. Second, we consider the case of non-interacting microswimmers, which simplifies the analysis significantly. Due to the similarity with the derivation in Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020), we here present the key steps of the derivation and refer the interested reader to that paper for technical details.
The quantity of interest is the fluid velocity autocorrelation function 
$C_U(t)$, formally defined as
\begin{equation} C_U(t) = \lim_{t'\rightarrow\infty} \frac{1}{V}\int {\rm d} {\boldsymbol r}\, \overline{U^{\alpha}( {\boldsymbol r}, t' )\,U^{\alpha}( {\boldsymbol r}, t+t')}, \end{equation}
where 
$V$ is the volume of the system, and the overbar denotes an average over the stochastic history, i.e. the history of tumble events, and the long-time limit ensures independence of the initial state of the system. Here and in the following, the superscript indices denote Cartesian components of vectors. For a given state of the system, the instantaneous fluid velocity 
$U^{\alpha }$ at a position 
${\boldsymbol r}$ and time 
$t$ is obtained as a superposition of individual velocity fields generated by the microswimmers
\begin{equation} U^{\alpha}( {\boldsymbol r}, t ) = \sum_{i=1}^N u_s^\alpha({{\boldsymbol r}}; {{\boldsymbol r}}_i(t),{{\boldsymbol p}}_i(t)). \end{equation}
Here, 
${\boldsymbol r}_i(t)$ gives the instantaneous position of particle 
$i$, while the unit vector 
${\boldsymbol p}_i(t)$ gives its instantaneous orientation; the index 
$i= 1, \dots, N$ enumerates the particles, where 
$N$ is the total number of microswimmers. In the following, we assume 
${\boldsymbol u}_s$ to be given by the regularised hydrodynamic dipole (3.1).
Time evolution of the suspension comprises spatial motion of the microswimmers according to (3.7), and their random re-orientation with rate 
$\lambda$, as discussed in § 1. Alternatively, this dynamics can be described by the master equation
\begin{equation} \partial_t F_N + \sum_{i=1}^N ( v_s p_i^\alpha\,\nabla_i^\alpha - D_0\,\nabla_i^2 ) F_N ={-} N \lambda F_N + \frac{\lambda}{4{\rm \pi}} \sum_{i=1}^N \int \mathrm{d} {\boldsymbol p}_i\,F_N,\end{equation}
where 
$F_N=F_N({{\boldsymbol r}}_1,\ldots,{{\boldsymbol r}}_N,{{\boldsymbol p}}_1,\ldots,{{\boldsymbol p}}_N,t)$ is the 
$N$-particle probability distribution function. Here, 
${\boldsymbol \nabla }_i$ denotes spatial derivatives with respect to the coordinates of particle 
$i$. As shown in Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020), the same dynamics can be conveniently encoded by the equation
\begin{equation} \partial_t h + v_s p^\alpha\,\nabla^\alpha h - D_0\,\nabla^2 h + \lambda h - \frac{\lambda}{4{\rm \pi}} \int {\rm d} {\boldsymbol p}\, h = \chi({\boldsymbol r},{\boldsymbol p}, t),\end{equation}
where 
$h=h({\boldsymbol r},{\boldsymbol p}, t)$ is an auxiliary field related to the Klimontovich correlation function (Klimontovich & Silin Reference Klimontovich and Silin1962), and 
$\chi$ is a noise term with the properties
\begin{align} \langle \chi({\boldsymbol r},{\boldsymbol p},t) \rangle &= 0,\end{align}
\begin{align} \langle \chi({\boldsymbol r},{\boldsymbol p}, t)\,\chi({\boldsymbol r}',{\boldsymbol p}', t')\rangle &= \frac{n}{4{\rm \pi}}\,\delta(t-t') \left[ 2 \lambda \delta({\boldsymbol r} - {\boldsymbol r}') \left( \delta({\boldsymbol p} - {\boldsymbol p}' ) - \frac{1}{4{\rm \pi}} \right) \right. \nonumber\\ &\quad- \left.D_0 \delta({\boldsymbol p} - {\boldsymbol p}' ) ( \nabla^2 + \nabla'^{2}) \delta({\boldsymbol r} - {\boldsymbol r}' ) \vphantom{\frac{1}{4{\rm \pi}}} \right]. \end{align}The advantage of this representation lies in its direct relation to the phase-space density correlation function (Klimontovich & Silin Reference Klimontovich and Silin1962), which allows us to express the fluid velocity autocorrelation function as
\begin{align} C_U(t) &= \lim_{t'\rightarrow\infty} \frac{1}{V}\int {\rm d} {\boldsymbol r} \int {\rm d} {\boldsymbol r}'\, {\rm d} {\boldsymbol r}'' \, {\rm d} {\boldsymbol p}'\, {\rm d} {\boldsymbol p}''\,u_s^\alpha({\boldsymbol r}; {\boldsymbol r}',{\boldsymbol p}')\, u_s^\alpha({\boldsymbol r}; {\boldsymbol r}'',{\boldsymbol p}'') \nonumber\\ & \quad \times \langle h({\boldsymbol r}',{\boldsymbol p}', t)\,h({\boldsymbol r}'',{\boldsymbol p}'', t+t') \rangle_\chi, \end{align}
where the angle brackets denote the average of the (fictitious) noise 
$\chi$ (see Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020) for details). The linear equation (A4) is solved by introducing the Fourier–Laplace transform 
$\hat {h}$ of the auxiliary field 
$h$:
\begin{equation} \hat{h}({\boldsymbol k}, {\boldsymbol p}, s) = \int_0^\infty {\rm d} t \, {\rm e}^{{-}s t} \int {\rm d} {\boldsymbol r} \,{\rm e}^{{\rm i} {\boldsymbol k}\boldsymbol{\cdot}{\boldsymbol r}}\, h({\boldsymbol r}, {\boldsymbol p}, t), \end{equation}which yields
\begin{equation} \hat{h}({\boldsymbol k}, {\boldsymbol p}, s) = \frac{\hat{\chi}({\boldsymbol k}, {\boldsymbol p}, s)}{\sigma({\boldsymbol k}, {\boldsymbol p}, s)} + \frac{\lambda}{4{\rm \pi}\,\sigma({\boldsymbol k}, {\boldsymbol p}, s)}\,\frac{\displaystyle \int {\rm d} {\boldsymbol p}'\,\frac{\hat{\chi}({\boldsymbol k}, {\boldsymbol p}', s)}{\sigma({\boldsymbol k}, {\boldsymbol p}', s)}}{\displaystyle 1-\frac{\lambda}{4{\rm \pi}} \int \frac{{\rm d} {\boldsymbol p}' }{\sigma({\boldsymbol k}, {\boldsymbol p}', s)}}. \end{equation}
Here, 
$\sigma ({\boldsymbol k},{\boldsymbol p}, s) = s + \lambda + D_0 k^2 + i v_s ({\boldsymbol k}\boldsymbol {\cdot }{\boldsymbol p})$, 
$\hat {\chi }({\boldsymbol k}, {\boldsymbol p},s)$ is the Fourier–Laplace transform of the noise, and we have dropped the initial condition 
$\hat {h}({\boldsymbol k}, {\boldsymbol p}, t=0)$, which does not contribute in the large-
$t$ limit.
As demonstrated by Škultéty et al. (Reference Škultéty, Nardini, Stenhammar, Marenduzzo and Morozov2020), only the first term in (A9) contributes to 
$C_U(t)$. Performing the Fourier–Laplace transform in (A7) and combining it with (A9) yields
\begin{align} C_U(t) &= \frac{ n \kappa^2}{16{\rm \pi}^4 } \lim_{t'\rightarrow\infty} {\mathcal{L}}^{{-}1}_{s_1,t'} {\mathcal{L}}^{{-}1}_{s_2,t'+t} \int {\rm d} {\boldsymbol k}\, \frac{A^2(k\varepsilon)}{k^4}\,( \lambda + D_0 k^2) \nonumber\\ &\quad \times \int {\rm d} {\boldsymbol p}\, (\boldsymbol k \boldsymbol{\cdot} {\boldsymbol p})^2 \left[1- \frac{(\boldsymbol k \boldsymbol{\cdot} {\boldsymbol p})^2}{k^2}\right] \frac{1}{s_1+s_2} \nonumber\\ & \quad \times \frac{1}{\lambda + s_1 + D_0 k^2 + i v_s (\boldsymbol k\boldsymbol{\cdot} {\boldsymbol p})}\, \frac{1}{\lambda + s_2 + D_0 k^2 - i v_s (\boldsymbol k\boldsymbol{\cdot} {\boldsymbol p})}. \end{align}
Performing the inverse Laplace transforms, denoted symbolically by 
${\mathcal {L}}^{-1}$, and integrating over 
$\boldsymbol p$ and the orientation of 
$\boldsymbol k$, we finally arrive at (3.8) in the main text.
Open access
