1. Introduction
Swimming microorganisms that accumulate into dense floating films are widespread at interfaces of aquatic environments (Durham & Stocker Reference Durham and Stocker2012; Zhan et al. Reference Zhan, Sardina, Lushi and Brandt2014; Mathijssen et al. Reference Mathijssen, Doostmohammadi, Yeomans and Shendruk2016; Sengupta, Carrara & Stocker Reference Sengupta, Carrara and Stocker2017; Vaccari et al. Reference Vaccari, Molaei, Niepa, Lee, Leheny and Stebe2017; Desai & Ardekani Reference Desai and Ardekani2020; Sepúlveda Steiner, Bouffard & Wüest Reference Sepúlveda Steiner, Bouffard and Wüest2021). The emergence of these ‘active carpets’ (ACs) (Mathijssen et al. Reference Mathijssen, Guzmán-Lastra, Kaiser and Löwen2018a; Guzmán-Lastra, Löwen & Mathijssen Reference Guzmán-Lastra, Löwen and Mathijssen2021) results from various factors, including hydrodynamic interactions, mechanical responses, taxis and the optimisation of metabolic activities, to name a few (Berke et al. Reference Berke, Turner, Berg and Lauga2008; Durham & Stocker Reference Durham and Stocker2012; Sommer et al. Reference Sommer2017; Ahmadzadegan et al. Reference Ahmadzadegan, Wang, Vlachos and Ardekani2019; Desai & Ardekani Reference Desai and Ardekani2020). Recently, Mathijssen et al. (Reference Mathijssen, Guzmán-Lastra, Kaiser and Löwen2018a) introduced a framework to quantify how ACs generate hydrodynamic fluctuations in their surrounding fluid. Depending on the swimming strategy of the microbes, ACs can either attract or repel fluid in the presence of orientational and/or density gradients (Lambert et al. Reference Lambert, Picano, Breugem and Brandt2013; Mathijssen et al. Reference Mathijssen, Guzmán-Lastra, Kaiser and Löwen2018a; Kanale et al. Reference Kanale, Ling, Guo, Fürthauer and Kanso2022), and generate anisotropic material transport in the water column (Guzmán-Lastra et al. Reference Guzmán-Lastra, Löwen and Mathijssen2021). Previous studies have shown that ACs can facilitate nutrient transport and enhance diffusion orders of magnitude larger than thermal fluctuations. However, it remains unexplored how such ACs can mobilise particles to and from liquid–air interfaces and impact the clustering dynamics of suspended matter within aquatic systems.
The active material transport induced by living and artificial microswimmers stands at the forefront of fluid dynamics research (Angelani et al. Reference Angelani, Maggi, Bernardini, Rizzo and Di Leonardo2011; Omar et al. Reference Omar, Wu, Wang and Brady2018; Dani, Yeganeh & Maldarelli Reference Dani, Yeganeh and Maldarelli2022; Gokhale et al. Reference Gokhale, Li, Solon, Gore and Fakhri2022; Madden et al. Reference Madden, Wang, Simmchen and Luijten2022). It encompasses diverse phenomena, starting with active diffusion in which tracer particles are entrained due to hydrodynamic flows induced by swimming microorganisms (Jeanneret et al. Reference Jeanneret, Pushkin, Kantsler and Polin2016; Mathijssen, Jeanneret & Polin Reference Mathijssen, Jeanneret and Polin2018b; Pellicciotta et al. Reference Pellicciotta, Hamilton, Kotar, Faucourt, Delgehyr, Spassky and Cicuta2020; Jin et al. Reference Jin, Chen, Maass and Mathijssen2021; Bárdfalvy et al. Reference Bárdfalvy, Škultéty, Nardini, Morozov and Stenhammar2024; Škultéty et al. Reference Škultéty, Bárdfalvy, Stenhammar, Nardini and Morozov2024), generating an effective enhanced diffusion over the summation of these encounters. This active diffusion has been measured experimentally for various microorganisms and via numerical simulations, predicting a fundamental relationship between the active flux, the persistence of swimmers and the relative distance between microswimmers (Miño et al. Reference Miño, Mallouk, Darnige, Hoyos, Dauchet, Dunstan, Soto, Wang, Rousselet and Clement2011, Reference Miño, Dunstan, Rousselet, Clément and Soto2013; Pushkin & Yeomans Reference Pushkin and Yeomans2013; Morozov & Marenduzzo Reference Morozov and Marenduzzo2014; de Graaf & Stenhammar Reference de Graaf and Stenhammar2017).
Although particles disperse due to active diffusion, their final concentration profile is restricted because of fluid confinement (Morozov & Marenduzzo Reference Morozov and Marenduzzo2014; Hamada, Cueto-Felgueroso & de Anna Reference Hamada, Cueto-Felgueroso and de Anna2020). In the case of ACs settling near no-slip boundaries, we find non-Boltzmannian concentration distributions (Guzmán-Lastra et al. Reference Guzmán-Lastra, Löwen and Mathijssen2021). This counterintuitive behaviour suggests the existence of an active temperature gradient, which has not been measured before in the context of biogenic hydrodynamic diffusion (Loi, Mossa & Cugliandolo Reference Loi, Mossa and Cugliandolo2008; Takatori & Brady Reference Takatori and Brady2015; Ortlieb et al. Reference Ortlieb, Rafaï, Peyla, Wagner and John2019).
When a microorganism swims, it induces local flows that, when superimposed over an entire colony of microswimmers, may generate coherent flows and large-amplitude disturbances. Examples of the impact of collective motions on scalar transport at intermediate Reynolds numbers, 
$Re\sim O(1)$, include the biogenic mixing induced by zooplankton diel vertical migration and bioconvection observed in natural aquatic systems and laboratory experiments (Pedley & Kessler Reference Pedley and Kessler1992; Hill & Pedley Reference Hill and Pedley2005; Dabiri Reference Dabiri2010; Wang & Ardekani Reference Wang and Ardekani2015; Simoncelli, Thackeray & Wain Reference Simoncelli, Thackeray and Wain2017; Sommer et al. Reference Sommer2017; Javadi et al. Reference Javadi, Arrieta, Tuval and Polin2020; Noto & Ulloa Reference Noto and Ulloa2023). Yet, biomixing owing to active turbulence (Alert, Casademunt & Joanny Reference Alert, Casademunt and Joanny2022) in ACs remain unexplored.
A rigorous quantification of biogenic mixing requires the analysis of tracer pair dispersion, rather than single-particle diffusivity (Belan & Kardar Reference Belan and Kardar2019). Pair dispersion provides information not only about the mixing capacity of a system but also about particle aggregation. The latter is a crucial mechanism for sustaining aquatic life and a fundamental skill of biological cells (Burd & Jackson Reference Burd and Jackson2009; Camassa et al. Reference Camassa, Harris, Hunt, Kilic and McLaughlin2019; Font-Muñoz et al. Reference Font-Muñoz, Jeanneret, Arrieta, Anglès, Jordi, Tuval and Basterretxea2019; Maheshwari et al. Reference Maheshwari, Sunol, Gonzalez, Endy and Zia2019; Arguedas-Leiva et al. Reference Arguedas-Leiva, Słomka, Lalescu, Stocker and Wilczek2022; Cruz & Neuer Reference Cruz and Neuer2022). Recent studies have deepened our insights into the dynamic clustering and aggregation of colloids (Angelani et al. Reference Angelani, Maggi, Bernardini, Rizzo and Di Leonardo2011; Jia et al. Reference Jia, Huang, Lan, Ren, Jiang and Lee2019; Dani et al. Reference Dani, Yeganeh and Maldarelli2022). Notably, these studies have revealed the role of colloid size in governing pairwise interactions with microswimmers (Gokhale et al. Reference Gokhale, Li, Solon, Gore and Fakhri2022; Kushwaha et al. Reference Kushwaha, Semwal, Maity, Mishra and Chikkadi2023). Factors such as tumbling rate and chirality have been found to influence short-range encounters, diffusion rates and mixing (Belan & Kardar Reference Belan and Kardar2019; Madden et al. Reference Madden, Wang, Simmchen and Luijten2022). However, our current understanding of the hydrodynamic interactions induced by microswimmers living near fluid–fluid interfaces on mixing and colloid clustering remains limited (Wang, Wexler & Zhou Reference Wang, Wexler and Zhou1998; Gonzalez, Aponte-Rivera & Zia Reference Gonzalez, Aponte-Rivera and Zia2021; Dani et al. Reference Dani, Yeganeh and Maldarelli2022).
In this paper, we build upon the singularity method to delve into the dynamics of ACs living near fluid–air interfaces and their impact on the tracer dynamics. Our objectives are twofold. Firstly, we seek to expand the theory of ACs by providing analytical solutions for the velocity fluctuations induced by a colony of microswimmers inhabiting fluid–air interfaces. Secondly, we investigate suspended tracer dynamics, quantifying universal metrics, including the single and pair diffusivities, and particle aggregation driven by ACs. We start by presenting an overview of the theoretical framework for ACs and the numerical methods applied to analyse both ACs and the tracer dynamics. Finally, we report our findings and draw conclusions.
2. Active carpets
We consider microswimmers in a three-dimensional (3-D) semi-infinite fluid that self-organise into an AC, as illustrated in figure 1(a). The fluid–air interface is located at 
$z=0$, which satisfies non-penetrative and free-slip boundary conditions. The AC is situated just beneath the interface, at a distance 
$z=h$, where 
$h$ is the typical swimmer size. We use a Cartesian coordinate system where the positive 
$z$ direction points downward from the liquid–air interface so that 
$z$ indicates depth. To analyse the flow generated by each swimmer, we utilise a multipole expansion of the fundamental solution of the Stokes equations together with the method of images (see e.g. Lauga Reference Lauga2020). For a swimmer located at 
$\boldsymbol {r}_s$, an image swimmer is located at the mirror position 
$\boldsymbol {r}_s^*={\mathsf {M}}\boldsymbol {\cdot }\boldsymbol {r}_s$, where 
${\mathsf {M}}$ represents the mirror matrix 
${\mathsf {M}}=\text {diag}(1,1,-1)$. Microswimmers, such as motile bacteria, can be considered as force dipoles in the far-field approximation, as shown by Drescher et al. (Reference Drescher, Dunkel, Cisneros, Ganguly and Goldstein2011). Defining the orientation of the swimmer 
$\boldsymbol {\hat {p}}$ as the unit vector along the axis that connects the flagella with the head of the swimmer, the fluid flow at a position 
$\boldsymbol {r}$ caused by a swimmer located at position 
$\boldsymbol {r}_s$ is (Mathijssen et al. Reference Mathijssen, Doostmohammadi, Yeomans and Shendruk2016)
\begin{equation} \boldsymbol{u}_s(\boldsymbol{r},\boldsymbol{r}_s,\boldsymbol{\hat{p}})={-}\kappa\hat{\boldsymbol{p}}\boldsymbol{\cdot}\boldsymbol{\nabla}_{s} \boldsymbol{v}-\kappa\hat{\boldsymbol{p}}\boldsymbol{\cdot}\boldsymbol{\nabla}^*_{s}\boldsymbol{v}^*, \end{equation}
where 
$\boldsymbol {v} = {\mathsf {G}}(\boldsymbol {r}-\boldsymbol {r}_s)\boldsymbol {\cdot }\hat {\boldsymbol {p}}$ is the Stokeslet flow and 
$\boldsymbol {v}^{*} = {\mathsf {M}}{\mathsf {G}}(\boldsymbol {r}-\boldsymbol {r}_s^*)\boldsymbol {\cdot } \hat {\boldsymbol {p}}$ is its image, with 
${\mathsf {G}}(\boldsymbol {x})=({1}/{|\boldsymbol {x}|})({\mathsf {I}}+\boldsymbol {x} \boldsymbol {x}/|\boldsymbol {x}|^2)$ the Oseen tensor. The coefficient 
$\kappa$ is the dipole strength that can be expressed in terms of the exerted force by the swimmer in the fluid, its characteristic length and the fluid viscosity (see e.g. Mathijssen et al. Reference Mathijssen, Doostmohammadi, Yeomans and Shendruk2016). In general, the dipole strength sign depends on the swimming strategy and the swimmer body geometry: 
$\kappa > 0$ represents a pusher, whereas 
$\kappa < 0$ represents a puller (Happel & Brenner Reference Happel and Brenner1983).
(a) An AC near the liquid–air interface generates flows that stir particles suspended below. (b) Strength of hydrodynamic fluctuations as a function of distance from the AC. Markers represent simulation results and lines represent the theoretical prediction from (2.5b). (c) Single-particle diffusivity in the 
$z$ direction. Circles are results from simulated mean-squared displacements, the squares are results from simulating the variance following (3.3) and the solid orange line represents the theoretical prediction from (2.5b) and (3.3).
We consider an AC composed of point dipoles, each of them located at positions 
$\boldsymbol {r}_s(t)=(x_s(t),y_s(t),h)$ and oriented along 
$\boldsymbol {\hat {p}}(t)=(p_x(t),p_y(t),0)$, forming a flat monolayer of microswimmers, as sketched in figure 1(a). At any time 
$t$, microswimmers exert a flow field given by (2.1) so the total flow field generated by an AC is given by 
$\boldsymbol {u}(\boldsymbol {r},t)=\sum _s\boldsymbol {u}_s(\boldsymbol {r},\boldsymbol {r}_s(t),\hat {\boldsymbol {p}}(t))$ with
\begin{align}
\boldsymbol{u}_s(\boldsymbol{r},\boldsymbol{r}_s,\boldsymbol{\hat{p}})
&= \kappa \left[(x-x_s)
\left(\left(\frac{-1}{r_{-}^{3/2}}+\frac{-1}{r_+^{3/2}}\right)+
3(\,\hat{\boldsymbol p}\boldsymbol{\cdot} (\boldsymbol
r-\boldsymbol r_s))^2 \left(\frac{1}{{r_{-}^{5/2}}}+
\frac{1}{r_+^{5/2}}\right)\right)\hat{x} \right.\nonumber\\
&\quad +(y-y_s)\left(
\left(\frac{-1}{r_{-}^{3/2}}+\frac{-1}{r_+^{3/2}} \right)+
3(\,\hat{\boldsymbol p}\boldsymbol{\cdot} (\boldsymbol
r-\boldsymbol r_s))^2 \left(\frac{1}{{r_{-}^{5/2}}}+
\frac{1}{r_+^{5/2}}\right)\right)\hat{y} \nonumber\\ &\quad
\left.+\left(\frac{(h-z)}{r_{-}^{3/2}}+\frac{(h+z)}{r_+^{3/2}}
+ 3(\,\hat{\boldsymbol p}\boldsymbol{\cdot} (\boldsymbol
r-\boldsymbol r_s))^2\left(\frac{h+z}{r_+^{5/2}}-
\frac{h-z}{r_{-}^{5/2}} \right)\right)\hat{z}\right],
\end{align}
where 
$\hat {\boldsymbol p}\boldsymbol {\cdot } (\boldsymbol r-\boldsymbol r_s)=p_x(x-x_s)+p_y(y-y_s)$, and 
$r_{\pm } = (x-x_s)^{2} + (y-y_s)^{2} + (z \pm h)^{2}$.
We consider the case where the induced flow is probed by fluid parcels that are far from the free-slip surface compared with the position of the swimmers. Therefore, we make the substitution 
$z_s=\varepsilon z$ with 
$\varepsilon \ll 1$, and perform a Taylor expansion on the flow exerted by a single swimmer, obtaining
\begin{equation} \boldsymbol{u}_s(\boldsymbol{r},\boldsymbol{r}_s,\boldsymbol{\hat{p}})=\boldsymbol{u}_{s,\varepsilon=0}+ \left.\frac{\partial \boldsymbol{u}_s}{\partial \varepsilon}\right\rvert_0\varepsilon+O(\varepsilon^2). \end{equation} Evaluating the velocity field at 
$\boldsymbol {r}=(0,0,z)$, with 
$\boldsymbol {r}_s=(\rho _s\cos \theta _s,\rho _s\sin \theta _s,h)$ and 
$\boldsymbol {\hat p}=(\cos \phi _s, \sin \phi _s,0)$, the approximated far-field velocity field is
$$\begin{gather} u_{s,\rho_s} = \frac{\kappa\rho_{s} [2z^{2} - \rho_{s}^2(1 + 3\cos(2(\theta_{s} - \phi_{s})))]}{(z^{2} + \rho_{s}^2)^{5/2}}(\cos \theta_{s},\sin \theta_{s}), \end{gather}$$
$$\begin{gather}u_{s,z} = \frac{z\kappa[{-}2z^{2} + \rho_{s}(1+3\cos(2(\theta_{s} - \phi_{s})))]}{(z^{2} + \rho_{s}^2)^{5/2}}. \end{gather}$$
Our AC is characterised by a uniform number density equal to 
$n$ organisms per unit area, with uniformly distributed orientations and positions, so the AC probability density function is 
$f(\boldsymbol {r}_s,\boldsymbol {\hat {p}})=n/2{\rm \pi}$. The ensemble-averaged flow field induced by an AC is 
$\langle \boldsymbol {u} \rangle = \int _{0}^{\infty } \int _{-{\rm \pi} }^{{\rm \pi} } \int _{-{\rm \pi} }^{{\rm \pi} }f\boldsymbol {u}_s(\boldsymbol {r},\boldsymbol {r}_s,\hat {\boldsymbol {p}}_s)\,{\rm d}\phi _s\,{\rm d}\theta _s \,\rho _s{\rm d}\rho _s =0.$ Thus, the mean flow at any instant is equal to zero for an infinite uniform AC, as expected by symmetry. However, the variances of the hydrodynamic fluctuations are not. We can derive them as follows:
$$\begin{gather} \langle u_i u_j\rangle= \int_{0}^{\infty} \int_{-{\rm \pi}}^{\rm \pi} \int_{-{\rm \pi}}^{\rm \pi} \frac{n}{2{\rm \pi}} u_{s,i} u_{s,j}\, {\rm d}\phi_s\,{\rm d}\theta_s \,\rho_s{\rm d}\rho_s, \end{gather}$$
$$\begin{gather}\langle u_x^2\rangle= \langle u_y^2\rangle=\frac{11{\rm \pi} n\kappa^2}{16z^2},\quad \langle u_z^2\rangle=\frac{9{\rm \pi} n\kappa^2}{8z^2}. \end{gather}$$
Thus, the intensity of the hydrodynamic fluctuations is proportional to the surface density of microswimmers, 
$n$. Moreover, it scales with the square of the individual flow strength, 
$\kappa ^2$, and is the same for swimmers of the puller and pusher type. Finally, the fluctuations scale with the inverse of the distance squared, 
$1/z^2$, so they are stronger near the AC and weaker farther away.
To test our analytical theory, we perform numerical simulations where the ACs are implemented explicitly as a large set of discrete, non-interacting point dipoles. They are randomly distributed on a square horizontal surface with a side length of 
$L=Nh$, where 
$N$ is a large integer that varies from 
$10^3$ to 
$10^4$. The surface density of the microswimmers within the carpet is fixed to 
$n=0.1/h^2$, so the total number of organisms in the carpet is given by 
$N_s= n L^2 = 0.1 N^{2}$.
To evaluate the variances of the hydrodynamic fluctuations shown in (2.5b), we compute the average 
$\langle {\cdot } \rangle$ over a large ensemble of independent AC configurations. In each AC configuration, 
$N_s$ microswimmers are given new random positions and random orientations within the 
$x$–
$y$ plane where the colony lives.
In what follows, we choose parameters that are typical for Escherichia coli bacteria. That is, 
$h$ is taken equal to its body width, 
$h=1\,\mathrm {\mu }{\rm m}$, which corresponds to bacteria swimming just beneath the air–water interface, and the dipole strength 
$\kappa =p/(8{\rm \pi} \eta )$ is taken equal to the value measured by Drescher et al. (Reference Drescher, Dunkel, Cisneros, Ganguly and Goldstein2011), i.e. 
$\kappa = 30\,\mathrm {\mu } {\rm m}^3\,{\rm s}^{-1}$. Moreover, to simplify notation in the remainder of the paper, we will use ‘natural units’, where units of length correspond to micrometres and units of time correspond to seconds. With this, 
$h=1$ and 
$\kappa =30$, and all subsequent numerical values appearing in the paper must be completed with micrometres and seconds, according to their dimensions.
We then consider 
$N_t$ neutrally buoyant and non-Brownian particles to gather statistics for the stochastic tracer dynamics. The equation of motion of the 
$\alpha$th tracer is controlled by the flows induced by the AC
\begin{equation} \frac{{\rm d}\boldsymbol{r}_\alpha}{{\rm d}t}= \sum_s^{N_s}\boldsymbol{u}_s(\boldsymbol{r}_\alpha,\boldsymbol{r}_s,{\boldsymbol{\hat{p}}}), \end{equation}
where 
$\boldsymbol {u}_s(\boldsymbol {r}_\alpha,\boldsymbol {r}_s,\boldsymbol {\hat {p}})$ is the velocity field produced by an individual microswimmer at a position 
$\boldsymbol {r}_s$ and orientation 
$\boldsymbol {\hat {p}}$, acting on a tracer placed at 
$\boldsymbol {r}_\alpha$, using either the full expression from (2.1) or the far-field approximation from (2.4).
3. Results
3.1. Hydrodynamic fluctuations
We first focus on examining the velocity variance. Numerically, we produce a large ensemble of AC configurations and evaluate each variance component for a given initial tracer depth 
$z_0$. The simulation box has a size 
$L=10^3$, in which 
$N_s=10^5$ swimmers were placed. Figure 1(b) compares the analytical solution for the variance found in (2.5b) (black and dashed lines) with full numerical simulations (markers). Notice that the horizontal fluctuations, 
$\langle u_x^2\rangle$ and 
$\langle u_y^2\rangle$, are approximately 60 % larger than the vertical fluctuations, 
$\langle u_z^2\rangle$, showing a significant anisotropic behaviour. The results show that the far-field theory is less accurate close to the AC, 
$z_0\lesssim 10$, but it offers a good approximation at long distances from the AC.
These results show that the hydrodynamic fluctuations are space dependent in the vertical direction, decaying as 
$\langle u_i^2\rangle \propto 1/z_0^2$ through the fluid column. Compared with a solid no-slip boundary, where 
$\langle u_i^2\rangle \propto 1/z_0^4$ (Guzmán-Lastra et al. Reference Guzmán-Lastra, Löwen and Mathijssen2021), these fluctuations near a liquid–air interface are much longer in range, giving rise to a larger region of influence.
In Appendix A, we also consider more detailed models of ACs that more closely resemble experimental observations (Li & Ardekani Reference Li and Ardekani2014; Wioland, Lushi & Goldstein Reference Wioland, Lushi and Goldstein2016; Bianchi, Saglimbeni & Di Leonardo Reference Bianchi, Saglimbeni and Di Leonardo2017; Ahmadzadegan et al. Reference Ahmadzadegan, Wang, Vlachos and Ardekani2019). Instead of modelling the AC as a perfectly planar thin 2-D sheet with swimmers that are aligned parallel to the surface, we also consider the cases of a thick AC, a thin AC of swimmers with non-zero orientation angles and the two combined cases. In all cases, we observe that the variances do not change significantly compared with the initial case, so henceforth we continue with that AC description.
3.2. Single-particle diffusivity
Secondly, we examine the diffusion of particles that is caused by the hydrodynamic fluctuations generated by the AC. To explicitly obtain their diffusivity, we consider the single-particle mean-squared displacement (MSD), given by
\begin{equation} \langle \Delta r_{i}\Delta r_{j} \rangle= \sum_{t' = 0}^{\tilde{n}\tau} \sum_{t'' = 0}^{\tilde{n}\tau} \langle u_{i}(\boldsymbol{r}_{0}, t')u_{j}(\boldsymbol{r}_{0}, t'') \rangle \Delta t' \Delta t'', \end{equation}
where 
$\tilde {n}$ is the total number of time steps in the simulations and 
$\tau$ is selected such that it matches the characteristic time of swimmer reorientation. We consider the case where the flow fields are uncorrelated between two consecutive time steps. Each time step corresponds to a different independent snapshot of the AC with microswimmers that have randomly sampled positions and orientations. Then we have the relation 
$\langle u_{i}(\boldsymbol {r}_{0}, t')u_{j}(\boldsymbol {r}_{0}, t'') \rangle = \langle u_{i}(\boldsymbol {r}_{0})u_{j}(\boldsymbol {r}_{0}) \rangle \delta _{t't''}$, leading to the following expression for the MSD:
\begin{equation} \langle \Delta r_{i}\Delta r_{j} \rangle = \langle u_{i}(\boldsymbol{r}_{0})u_{j}(\boldsymbol{r}_{0})\rangle \sum_{t' = 0}^{\tilde{n}\tau}\sum_{t'' = 0}^{\tilde{n}\tau} \delta_{t't''}\Delta t'\Delta t'' = \langle u_iu_j \rangle\tilde{n}\tau^{2} = \langle u_iu_j \rangle\tau t_{f}, \end{equation}
where 
$t_{f} = \tilde {n}\tau$ denotes the final integration time. Therefore, we have 
$\langle \Delta r_{i}^{2}\rangle = \langle u_{i}^{2} \rangle \tau t_{f} \equiv 2D_{i}^{s}t_{f}$, so the single-particle diffusivity is
\begin{equation} D_{i}^{s} \equiv \tfrac{1}{2}\langle u_{i}^{2} \rangle \tau. \end{equation}
To verify this result numerically, we place tracers initially at a depth 
$z_0$ and we compute their MSD for 
$\tilde {n}=10^2$ time steps. Figure 1(c) shows the single-particle diffusivity 
$D_{i}^{s}$ for 
$i=z$ as a function of 
$z_{0}$ obtained from the simulations and from the theory. The simulated observable 
$\langle u_z^{2} \rangle \tau /2$ is shown as green squares, and the diffusivity, 
$D_{z}^{s}$, obtained from the simulated MSDs, is displayed as blue circles. These markers collapse onto a single curve, which is compared with the theory from (3.3) shown by the orange line. Indeed, the results exhibit strong agreement between theory and simulations. As predicted, the diffusivity decays with an exponent of 
$-2$ for large 
$z_0$. This implies that particle stirring significantly intensifies near the AC, potentially enhancing particle encounters, which we discuss next.
3.3. Velocity pair correlation
Thirdly, we look at how pairs of suspended particles respond to fluid distortions caused by the AC. Although the generated flow is time uncorrelated, the spatial structure is more complex: particles that are initially located close to each other follow adjacent streamlines that are almost parallel to each other, so their motions remain correlated at short distances. By contrast, with larger initial separations, the particles follow different streamlines and separate further from each other. This relative motion is quantified by the velocity pair correlation function
\begin{equation} g_{i,p}(d_0, z_0) = \frac{\langle u_{i}(\boldsymbol{r}_1) u_{i} (\boldsymbol{r}_2)\rangle}{\sqrt{\langle u_{i}^{2}(\boldsymbol{r}_1)\rangle \langle u_{i}^{2}(\boldsymbol{r}_2) \rangle}}, \end{equation}
where the initial positions of the two particles are 
$\boldsymbol {r}_1 = (d_0/2\sqrt {2}, d_0/2\sqrt {2}, z_0)$ and 
$\boldsymbol {r}_2 = (-d_0/2\sqrt {2}, -d_0/2\sqrt {2}, z_0)$, which are located at the same distance from the carpet 
$z_0$ and separated laterally by 
$d_0$. Notice that the velocity pair correlation is symmetric in the 
$x$ and 
$y$ directions.
Numerically, we sample over 
$10^{4}$ ensembles, each with 
$N_s=10^7$ swimmers, to measure the velocity using (2.1) at the positions 
$\boldsymbol {r}_1$ and 
$\boldsymbol {r}_2$ given above, with 
$z_0$ in the range 
$[2,20]$ and 
$d_0$ in the range [
$10^{-3}$, 
$10^3$]. The resulting velocity pair correlations, in the vertical and horizontal orientations, are shown as circles in figure 2(a,c). As predicted, we find that pairs become progressively more uncorrelated with distance 
$d_0$. Moreover, the separation required for the pairs to become uncorrelated grows with 
$z_0$. Notably, when the velocity pair correlations are plotted in terms of the scaled distance 
$d_0/z_0$, there is a collapse for all depths onto a single curve, as shown in the insets of figure 2(a,c). Globally, velocities decorrelate at different rates, depending on the velocity direction, with particles decorrelating at distances larger than 
$d_0\gtrsim 7 z_0$ in the horizontal direction and 
$d_0\gtrsim 4 z_0$ in the vertical direction.
(a) Velocity pair correlation in the planar coordinate, 
$g_{\rho }$, as a function of horizontal particle separation. Coloured circles represent simulation results for different depths, ranging from 
$z_0 = 2$ (blue) to 
$z_0 = 20$ (light green), each separated by 
$\Delta z_0=2$. Lines represent the best fit to the model from (3.5). Inset: these curves collapse when plotted as 
$d_0/z_0$. The yellow line is the best fit from (3.5), and the pink dashed line is the semi-analytical result obtained by numerically integrating 
$g_{\rho }(d_0, z_0)$. (b) Pair diffusivity in the planar coordinate, 
$D^{p}_{\rho }$, as a function of horizontal particle separation, normalised by the single-particle diffusivity. Markers represent simulation values for different depths, 
$z_0$. The black dashed line is the theoretical asymptotic value from 
$D_{i}^{p\infty }$. (c,d) Show the pair correlation and the pair diffusion in the vertical coordinate, 
$g_z$ and 
$D^p_{z}$, respectively.
First, we propose a heuristic model for the collapsed velocity correlation of the form
\begin{equation} g_{i,p}^{model}(\xi) = \exp(-\xi)[A_1 \sin(\xi) + \cos(\xi)] + \exp(-\xi/A_2)[A_3 + A_4\xi], \end{equation}
where 
$\xi = d_0/z_0$. Using least squares to fit this model to our simulation data, we find that the optimal fitting parameters are in the horizontal direction 
$A_1 = 0.097$, 
$A_2 = 1.3$, 
$A_3 = -0.014$, 
$A_4 = 0.89$ and in the vertical direction 
$A_1 = -0.18$, 
$A_2 = 0.52$, 
$A_3 = -0.09$, 
$A_4 = 1.35$. The resulting model is shown as coloured lines for the different 
$z_0$ values in figure 2(a,c), and also in the inset.
Analytical progress can be made by using the far-field approximation, where 
$g_{i,p}$ is obtained by evaluating the velocity from (2.4) and the variances from (2.5) at the positions 
$\boldsymbol {r}_{1}$ and 
$\boldsymbol {r}_{2}$
\begin{equation} g_{i,p}(d_0, z_0) = \frac{1}{\langle u_i^{2} \rangle} \int_{0}^{2{\rm \pi}}\int_{0}^{2{\rm \pi}}\int_{0}^{\infty} u_{i}(\boldsymbol{r}_1) u_{i}(\boldsymbol{r}_2) \frac{n}{2{\rm \pi}} \rho_{s}\,{\rm d}\rho_{s}\,{\rm d}\theta_{s}\,{\rm d}\phi_{s}. \end{equation}
It is challenging to find a closed-form expression for the above integral, but it can be evaluated numerically; the semi-analytical result is shown in the inset of figure 2(a). The theoretical result, which is in strong agreement with our simulations, underscores the crucial role of the vertical distance to the AC in our study: it not only determines the strength of the hydrodynamic fluctuations but also the relative motion of the flow structures, with strong correlation at larger distances for flows parallel to the AC and fast decaying correlations in the vertical direction, which is also evident with the heuristic model in (3.5) with the parameter 
$A_2$ representing the rate at which particles decorrelate with distance, where 
$A_2^{\rho }=1.3$ and 
$A_2^{z}=0.52$. Particles at greater depths tend to move together synchronously, which directly impacts their pair diffusivity, as we show next.
3.4. Pair diffusion
We now explore how an AC drives can mix two particles up, as a function of their distance. This is quantified by the ‘pair diffusivity’. Numerically, we sample the carpet as described previously. Two tracers are put initially at the positions 
$\boldsymbol {r}_1$ and 
$\boldsymbol {r}_2$, after which they follow the equation of motion (2.6). The integration time step is 
$\Delta t = 10^{-3}$. Averages are performed over 
$300$ pair trajectories. The pair diffusivity 
$D_{i}^p$ is obtained by measuring the pair's distance squared, 
$\Delta \rho ^2(t)=(|\boldsymbol {r}_1(t)-\boldsymbol {r}_2(t)|-d_0)^2$. After the simulation ends, each 
$\Delta \rho ^2$ curve is fitted with a power-law function of the form 
$4D_{i}^p t^{b}$. Results show that all curves follow a diffusive regime (
$b \sim 1$).
Figure 2(b,d) displays the resulting pair diffusivity 
$D_{i}^p$, normalised by the single-particle diffusivity 
$D^{s}_{i}$, as a function of the initial horizontal particle separation 
$d_0$. For small separations, the particles follow almost the same streamlines, so their relative motion and pair diffusivity are small. Indeed, it is difficult to mix particles that are close to each other, as we observe in the kitchen when stirring macaron batter (Mathijssen et al. Reference Mathijssen, Lisicki, Prakash and Mossige2023). At large particle separations, it gets easier to mix them, so the pair diffusivity increases with 
$d_0$ until it plateaus.
To predict this analytically, we define 
$\boldsymbol {r}_{12}= \boldsymbol {r}_{1}- \boldsymbol {r}_{2}$. Similar to the single-particle MSD calculation, the pair MSD is then computed as follows:
\begin{equation} \langle \Delta r_{12i} \Delta r_{12j} \rangle = \langle \Delta r_{1i} \Delta r_{1j} \rangle + \langle \Delta r_{2i}\Delta r_{2j}\rangle - \langle \Delta r_{1i}\Delta r_{2j} \rangle - \langle \Delta r_{2i}\Delta r_{1j} \rangle, \end{equation}
where 
$\Delta r_{\alpha i}$ denotes the displacement of the particle 
$\alpha$ along coordinate 
$i$ and analogously for 
$\Delta r_{12 i}$. On average, terms 1 and 2 on the right-hand side of (3.7) are the same, and terms 3 and 4 are also equal to each other. Therefore, 
$\langle \Delta r_{12i} \Delta r_{12j}\rangle = 2\langle \Delta r_{1i}\Delta r_{1j}\rangle - 2\langle \Delta r_{1i}\Delta r_{2j} \rangle$. Considering 
$i = j$, the first term is 
$\langle \Delta r_{i}\Delta r_{i} \rangle = 2 D_{i}^{s} t$, with 
$D_i^{s}$ the single diffusivity in the 
$i$-coordinate. The second term can be estimated in a similar way as in the single MSD calculation
\begin{equation} \langle \Delta r_{1i}\Delta r_{2i} \rangle = \sum_{t' = 0}^{\tilde{n}\tau}\sum_{t'' = 0}^{\tilde{n} \tau} \langle u_{1i}u_{2i} \rangle \Delta t'\Delta t''\ ({\rm with}\ \tau=\Delta t). \end{equation}
The function to be integrated is the complete spatio-temporal pair correlation velocity, 
$\langle u_{1i}u_{2i} \rangle (d_0, z_0, t) = \langle u_i^{2} \rangle g_{i,p}(d_0, z_0) \delta _{t't''}$, where 
$g_{i,p}$ is the spatial pair velocity correlation defined in (3.4). With this, we obtain
\begin{equation} \langle \Delta r_{1i}\Delta r_{2i} \rangle = \langle u_{i}^{2} \rangle g_{i,p} \sum_{t'=0}^{\tilde{n}\tau}\sum_{t''=0}^{\tilde{n}\tau} \delta_{t't''} \Delta t' \Delta t'' = \langle u_{i}^{2} \rangle g_{i,p} \tau t_{f}. \end{equation}
Thus, 
$\langle (\Delta r_{12i})(\Delta r_{12j})\rangle = 2(2D_{i}^{s} - \langle u_{i}^{2} \rangle \tau g_{i,p})t_{f}$, where we obtain the pair diffusivity
\begin{equation} D_{i}^{p} \equiv 2D_{i}^{s} - \langle u_{i}^{2} \rangle \tau g_{i,p} = 2D_{i}^{s}(1 - g_{i,p}). \end{equation}
This expression depends upon the variance and single diffusivity relation in (3.3). Since the function 
$g_{i,p}(d_0, z_0)$ decays exponentially with the distance 
$d_0$, the asymptotic pair diffusivity is readily predicted for large 
$d_0$
\begin{equation} \lim_{d_0 \to \infty} D_{i}^{p} \equiv D_{i}^{p\infty} = 2D_{i}^{s}. \end{equation}
Figure 2(b,d) shows this theoretical result for the pair diffusivity (3.10), using the 
$g_{i,p}^{model}$ in (3.5). This result agrees well with the simulations, as shown for various distances 
$z_0$.
Interestingly, the observed diffusive behaviour in our study bears a resemblance to turbulent systems, as seen in Belan & Kardar (Reference Belan and Kardar2019), where a pair of tracers within a 3-D bath of microswimmers also exhibited an asymptotic diffusivity twice that of their self-diffusion, reaching numerical values of 
$D^p_{\infty }\approx 2\,\mathrm {\mu } {\rm m}^2\,{\rm s}^{-1}$, which is of the same order of magnitude as in our case when we choose 
$\tau =\tau _R=10$ as a typical memory time for E. coli bacteria and 
$n = 0.01, \kappa = 30, z_0 = 10, d_0 = 10$, we then get 
$D_{\rho }^p\approx 1.39\,\mathrm {\mu } {\rm m}^2\,{\rm s}^{-1}$, with horizontal and vertical self-diffusivities 
$D^s_{\rho } \approx 1.94\,\mathrm {\mu }\text {m}^2\,\text {s}^{-1}$, 
$D^s_z \approx 3.18\,\mathrm {\mu }\text {m}^2\,\text {s}^{-1}$, respectively, which are one order of magnitude larger than thermal diffusion for a passive spherical micron-sized particle 
$D_{T} \approx 0.22\,\mathrm {\mu }\text {m}^2\,\text {s}^{-1}$.
To put these numbers in perspective, in natural aquatic environments, bioconvection propelled by Chromatium Okenni (with a body length of 1–4 
$\mathrm {\mu }$m) has been demonstrated to notably boost the actual diffusion of heat within their ecological niche (Sommer et al. Reference Sommer2017). This enhancement, observed to be at least an order of magnitude, was documented in ‘Lago di Cadagno’, Switzerland (Sepúlveda Steiner, Bouffard & Wüest Reference Sepúlveda Steiner, Bouffard and Wüest2019). Similar diffusion enhancement has been measured in laboratory experiments driven by E. coli (Singh et al. Reference Singh, Patteson, Maldonado, Purohit and Arratia2021). Remarkably, the thermal diffusion augmentation resulting from the collective swimming of microorganisms is not significantly different from the enhancement in vertical heat diffusion around the thermocline induced by wind-driven shear flows in strongly stratified lakes, as evidenced in Lac Léman (Fernández Castro et al. Reference Fernández Castro, Bouffard, Troy, Ulloa, Piccolroaz, Sepúlveda Steiner, Chmiel, Serra Moncadas, Lavanchy and Wüest2021; Sepúlveda Steiner et al. Reference Sepúlveda Steiner, Forrest, McInerney, Fernández Castro, Lavanchy, Wüest and Bouffard2023).
Our results reveal that the distance parameter 
$z_0$ significantly affects diffusivity. The mixing becomes more intense and converges to the asymptotic value when particles are closer to the AC. This phenomenon can be attributed to amplified hydrodynamic fluctuations near the carpet, which disrupt spatial correlations in fluid flows. The parameter 
$d_0$ is also crucial; closer particles exhibit lower diffusion from each other compared with those at greater distances, resulting in synchronised movement and more similar trajectories, with vertical pair diffusion reaching the asymptotic value faster than horizontal pair diffusion as expected from the pair velocity correlation behaviour. Consequently, the pair diffusivity might offer fundamental insights into aggregation phenomena, which we study next.
3.5. Particle aggregation
We investigate whether fluctuations caused by ACs can initiate the aggregation of spherical particles of finite radius 
$r_{0}$. To model this, we consider a short-ranged pairwise sticky force between the particles. We consider the Morse model, given by
\begin{equation} \boldsymbol{F}_{M}(\boldsymbol r) = 2UW[1 - \exp({-W(r-r_{eq})})]\exp({-W(r - r_{eq})})\hat{\boldsymbol{r}}, \end{equation}
where 
$r=|\boldsymbol {r}|$ is the distance between the tracers. Here, 
$r_{eq}$ is the pair equilibrium distance, while 
$U$ and 
$W$ denote the depth and width of the potential well, respectively. The force is scaled with the Stokes mobility, so it has units of velocity. By fixing 
$r_{eq}$ to twice the radius 
$r_{0}$, 
$U$ to 
$10^{-9}$ and 
$W$ to 30, we set the potential minima along 
$r$, allowing for the precise adjustment of equilibrium distances for tracers with different radii 
$r_{0}$.
Initially, at 
$t=0$, the particles are placed in a square lattice of 
$N_t = 19\times 19 = 361$ particles located at the same depth 
$z_0$. The distance between their centres is equal to 
$d_{0} = \tilde {d}_0 + 2r_{0}$, where 
$\tilde {d}_0$ is their border-to-border distance. We always set 
$\tilde {d}_0 = 0.8$ and vary 
$d_0$ or equivalently 
$r_{0}$. For 
$t>0$, the particles start moving in the 
$x-y$ plane at a fixed depth 
$z_0$, obeying the equation of motion (2.6), with the addition of the inter-particle sticking forces that allow them to collide and aggregate with each other
\begin{equation} \frac{{\rm d}\boldsymbol{r}_\alpha}{{\rm d}t} = \sum_s^{N_{s}} \boldsymbol{u}_s(\boldsymbol{r}_\alpha,\boldsymbol{r}_s,\boldsymbol{\hat{p}})+ \sum_{\alpha{\neq} \beta}^{N_{t}} \boldsymbol{F}_M^{\alpha\beta} (\boldsymbol{r}_\alpha, \boldsymbol{r}_\beta, t). \end{equation}We investigate the particle aggregation by examining the histogram of first passage times, defined as the time taken for particles to collide.
Results are shown in the inset of figure 3(a); the red dashed line shows the mean first passage time, 
$\tau _{FP}$. This observable can also be estimated analytically by equating the pair MSD in (3.10) to the border-to-border distance between particles, 
$\tilde {d_0}^{2} \approx \langle \Delta r_{1} \Delta r_{2} \rangle = 4 D_{\rho }^p\tau _{FP}$, so the mean first passage time is
\begin{equation} \tau_{FP} \approx \frac{\tilde{d_0}^{2}}{4 D_{\rho}^p}. \end{equation}
The theoretical prediction for 
$\tau _{FP}$ agrees well with the numerical simulations, as shown for 
$z_0 = 3$ by the solid green line in figure 3(a). Therefore, 
$\tau _{FP}$ provides a reliable estimate for the average collision time between two particles due to the hydrodynamic stirring induced by the AC.
(a) Mean first passage time of particle collisions as a function of the horizontal particle separation, 
$d_0$, for particles initially located at 
$z_0 = 3$. Squares represent simulations and the line is the prediction 
$\tau _{FP} \approx \tilde {d_0}^{2}/4 D_{\rho }^p.$ Inset: histogram of first passage times. The vertical red dashed line is the mean first passage time. (b) Clearance time against particle radius for three depths. Symbols show simulation results, and solid lines show 
$\tau _c\sim 0.7\tau _{FP}$. (c) Average accumulated number of collisions over time. Markers show averaged simulation results for different particle radii, ranging from 
$a = 3$ (triangles) to 
$a = 10$ (squares). The background colours represent different depths, ranging from 
$z_0 = 3$ (light blue) to 
$z_0 = 9$ (pink).
In addition, we calculate the average ‘clearance time’ denoted as 
$\tau _c$. This is the time taken for half of the particles in the suspension to collide, representing a fundamental time scale indicator for aggregation and cluster formation (Wang et al. Reference Wang, Wexler and Zhou1998; Font-Muñoz et al. Reference Font-Muñoz, Jeanneret, Arrieta, Anglès, Jordi, Tuval and Basterretxea2019). The graph shown in figure 3(b) illustrates the clearance time as a function of particles radius 
$r_{0}$ and depth 
$z_{0}$, where we assume that particles stick together after a collision. Our results show that, as the size of the particle 
$r_{0}$ increases, the clearance time decreases. Moreover, 
$\tau _{FP}$ and 
$\tau _{c}$ are proportional to each other since they both arise from the same physical process. Empirically, we observe that 
$\tau _c \sim 0.7 \tau _{FP}$ for every 
$z_0$, as shown by the lines in figure 3(b), which agrees well with the simulations.
Last, we measure the time evolution of the average accumulated number of collisions 
$\langle N_{cols}\rangle$ for this configuration. Figure 3(c) shows 
$\langle N_{cols}\rangle$ over time for each depth 
$z_0$. The results show a close correlation between 
$\langle N_{cols}\rangle$ and 
$\tau _c$. Collisions become more frequent as we move closer to the AC and increase with larger particle radius. The latter suggests that the pair correlation function 
$g_{\rho }$ in (3.4) governs both the collisions and the measured aggregation time scale. As the particle size increases, the correlation among their centres decreases, making it easier for them to follow different trajectories and ultimately leading to more collisions. Likewise, the intense fluctuations near ACs disrupt this correlation, increasing the likelihood of collisions.
4. Concluding remarks
Here, we have established analytical connections between the flows generated by ACs and the aggregation dynamics of suspended particles. For this, we have investigated broadly used fundamental metrics that characterised the suspended particle dynamics: the single and pair diffusivities and particle aggregation. Our key findings include: (i) an analytical solution for hydrodynamic fluctuations produced by ACs near fluid–air interfaces; (ii) the emergence of space-dependent and anisotropic diffusion, which decreases quadratically with distance; (iii) the role of hydrodynamic fluctuations in facilitating pairing encounters, where particle aggregation is favoured for large compared with small particles; (iv) the mean first passage time between collisions decreases as particles move farther apart; (v) in the close vicinity of ACs, intense hydrodynamic stirring accelerates clearance times and particle aggregation processes.
Although our study focused on ACs, the employed methodology can be applied to other swimmer configurations. In particular (3.10), which relates the pair diffusivity to the velocity pair correlation, is a valuable tool to study dynamical aggregation driven by active fluctuations. This research highlights the pivotal role of biologically driven flows in the transport and spatial organisation of particles in aquatic systems, serving as a noteworthy example of an out-of-equilibrium system that remains analytically tractable.
Acknowledgements
We are grateful for the helpful feedback given H. Löwen and F. Barros. We also thank the anonymous reviewers for their constructive feedback.
Funding
G.A., R.S. and F.G.-L. have received support from the ANID – Millennium Science Initiative Program – NCN19 170, Chile. F.G.-L. was supported by Fondecyt Iniciación No. 11220683. H.N.U. and A.J.T.M.M. were supported by start-up grants from the University of Pennsylvania. A.J.T.M.M. acknowledges funding from the United States Department of Agriculture (USDA-NIFA AFRI grants 2020-67017-30776 and 2020-67015-32330), the Charles E. Kaufman Foundation (Early Investigator Research Award KA2022-129523) and the University of Pennsylvania (University Research Foundation Grant and Klein Family Social Justice Award). This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Thick and tilted active carpet
In nature and laboratory experiments, swimming microorganisms form spatial distributions that are more complex than those of a single monolayer. They can create ‘thick’ films where the number density of microswimmers decays with depth. Also, it has been observed that microswimmers do not always orient parallel to air–liquid interfaces; observations also show that, in some cases, their average angle orientation is slightly out of the plane, which we denote as ‘tilted’ swimmers (Li & Ardekani Reference Li and Ardekani2014; Wioland et al. Reference Wioland, Lushi and Goldstein2016; Bianchi et al. Reference Bianchi, Saglimbeni and Di Leonardo2017; Ahmadzadegan et al. Reference Ahmadzadegan, Wang, Vlachos and Ardekani2019). We extended our results to three different cases: a thick carpet, a tilted carpet and a thick tilted carpet. Hence, we performed simulations to quantify the changes in the hydrodynamic fluctuations.
For the thick carpet, we consider swimmers that do not form a planar 2-D sheet, but they are distributed in three dimensions according to a vertical exponential function
\begin{equation} f(z_s) = 2\exp({-}2(z_s-h)), \end{equation}
where 
$z_s$ is the swimmer's vertical position. Using this distribution, in practice, swimmers are found up to approximately 
$z_s\sim 5h$. In this scenario, we defined the ‘average carpet position’ at 
$\bar {z} = 2.5h$. With this, we computed the variances for different depths in the same way as in the manuscript. Figure 4 displays these results in purple symbols. We notice that the thick carpet variances conserve the thin AC behaviour with vertical fluctuations larger than planar fluctuations decaying with the same power law. We also observe, quantitatively, some changes close to the AC. In contrast, for far distances, the values are alike.
Comparison between variances measured when the AC is thick and have planar swimmers (purple markers), when the AC is flat and have tilted swimmers (light-green markers), when the AC is thick and have tilted swimmers (orange markers) and when the AC is flat and have planar swimmers (hollow markers, the manuscript case).
For the tilted case, we consider microswimmers with vertical out-of-plane angles from 
$-10^{\circ }$ to 
$10^{\circ }$, distributed according to a Gaussian with zero mean and a standard deviation of 
$5^{\circ }$. Figure 4 displays the results in green symbols. Coloured and hollow markers match perfectly well; therefore, the 
$z_0$ dependency holds. Only a slight decrease in magnitude is observed in the tilted case. This change could be understood as a geometric factor (a little less than the unity) multiplying the planar AC used in this manuscript. In addition to the first case, we also performed numerical simulations using a thick carpet containing tilted swimmers (both features together) to ensure completeness. In figure 4 with orange markers, it can be seen that this case is almost the same as the thick carpet with planar swimmers case. Hence, we can conclude that a thick carpet affects hydrodynamic fluctuations more than a carpet with tilted swimmers. In all cases, however, the effect is small and does not change the order of magnitude of the variance.
Therefore, such variants in the AC properties do not obscure the most fundamental results and concluding remarks presented in the manuscript. Other AC distributions can be used, but any density profile that is concentrated near the interface will produce similar results to that of a Dirac delta distribution, with subdominant corrections. Moreover, they do not break the anisotropic behaviour of fluctuations, highlighting the universality of our results.
Open access
