Floating active carpets drive transport and aggregation in aquatic ecosystems

Abstract

Communities of swimming microorganisms often thrive near liquid–air interfaces. We study how such ‘active carpets’ shape their aquatic environment by driving biogenic transport in the water column beneath them. The hydrodynamic stirring that active carpets generate leads to diffusive upward fluxes of nutrients from deeper water layers, and downward fluxes of oxygen and carbon. Combining analytical theory and simulations, we examine the biogenic transport by studying fundamental metrics, including the single and pair diffusivity, the first passage time for particle pair encounters and the rate of particle aggregation. Our findings reveal that the hydrodynamic fluctuations driven by active carpets have a region of influence that reaches orders of magnitude further in distance than the size of the organisms. These non-equilibrium fluctuations lead to a strongly enhanced diffusion of particles, which is anisotropic and space dependent. Fluctuations also facilitate encounters of particle pairs, which we quantify by analysing their velocity pair correlation functions as a function of distance between the particles. We found that the size of the particles plays a crucial role in their encounter rates, with larger particles situated near the active carpet being more favourable for aggregation. Overall, this research broadens our comprehension of aquatic systems out of equilibrium and how biologically driven fluctuations contribute to the transport of fundamental elements in biogeochemical cycles.

1. Introduction

Swimming microorganisms that accumulate into dense floating films are widespread at interfaces of aquatic environments (Durham & Stocker 2012; Zhan et al. 2014; Mathijssen et al. 2016; Sengupta, Carrara & Stocker 2017; Vaccari et al. 2017; Desai & Ardekani 2020; Sepúlveda Steiner, Bouffard & Wüest 2021). The emergence of these ‘active carpets’ (ACs) (Mathijssen et al. 2018a; Guzmán-Lastra, Löwen & Mathijssen 2021) results from various factors, including hydrodynamic interactions, mechanical responses, taxis and the optimisation of metabolic activities, to name a few (Berke et al. 2008; Durham & Stocker 2012; Sommer et al. 2017; Ahmadzadegan et al. 2019; Desai & Ardekani 2020). Recently, Mathijssen et al. (2018a) 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. 2013; Mathijssen et al. 2018a; Kanale et al. 2022), and generate anisotropic material transport in the water column (Guzmán-Lastra et al. 2021). 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. 2011; Omar et al. 2018; Dani, Yeganeh & Maldarelli 2022; Gokhale et al. 2022; Madden et al. 2022). 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. 2016; Mathijssen, Jeanneret & Polin 2018b; Pellicciotta et al. 2020; Jin et al. 2021; Bárdfalvy et al. 2024; Škultéty et al. 2024), 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. 2011, 2013; Pushkin & Yeomans 2013; Morozov & Marenduzzo 2014; de Graaf & Stenhammar 2017).

Although particles disperse due to active diffusion, their final concentration profile is restricted because of fluid confinement (Morozov & Marenduzzo 2014; Hamada, Cueto-Felgueroso & de Anna 2020). In the case of ACs settling near no-slip boundaries, we find non-Boltzmannian concentration distributions (Guzmán-Lastra et al. 2021). 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 2008; Takatori & Brady 2015; Ortlieb et al. 2019).

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 1992; Hill & Pedley 2005; Dabiri 2010; Wang & Ardekani 2015; Simoncelli, Thackeray & Wain 2017; Sommer et al. 2017; Javadi et al. 2020; Noto & Ulloa 2023). Yet, biomixing owing to active turbulence (Alert, Casademunt & Joanny 2022) in ACs remain unexplored.

A rigorous quantification of biogenic mixing requires the analysis of tracer pair dispersion, rather than single-particle diffusivity (Belan & Kardar 2019). 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 2009; Camassa et al. 2019; Font-Muñoz et al. 2019; Maheshwari et al. 2019; Arguedas-Leiva et al. 2022; Cruz & Neuer 2022). Recent studies have deepened our insights into the dynamic clustering and aggregation of colloids (Angelani et al. 2011; Jia et al. 2019; Dani et al. 2022). Notably, these studies have revealed the role of colloid size in governing pairwise interactions with microswimmers (Gokhale et al. 2022; Kushwaha et al. 2023). Factors such as tumbling rate and chirality have been found to influence short-range encounters, diffusion rates and mixing (Belan & Kardar 2019; Madden et al. 2022). 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 1998; Gonzalez, Aponte-Rivera & Zia 2021; Dani et al. 2022).

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 2020). 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. (2011). 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. 2016)

(2.1)\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. 2016). 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 1983).

Figure 1. (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

(2.2) \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

(2.3)\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

(2.4a)$$\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}$$

(2.4b)$$\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:

(2.5a)$$\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}$$

(2.5b)$$\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. (2011), 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

(2.6)\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. 2021), 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 2014; Wioland, Lushi & Goldstein 2016; Bianchi, Saglimbeni & Di Leonardo 2017; Ahmadzadegan et al. 2019). 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

(3.1)\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:

(3.2)\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

(3.3)\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

(3.4)\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.

Figure 2. (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

(3.5)\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}$

(3.6)\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. 2023). 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:

(3.7)\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

(3.8)\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

(3.9)\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

(3.10)\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$

(3.11)\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 (2019), 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. 2017). This enhancement, observed to be at least an order of magnitude, was documented in ‘Lago di Cadagno’, Switzerland (Sepúlveda Steiner, Bouffard & Wüest 2019). Similar diffusion enhancement has been measured in laboratory experiments driven by E. coli (Singh et al. 2021). 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. 2021; Sepúlveda Steiner et al. 2023).

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

(3.12)\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

(3.13)\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

(3.14)\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.

Figure 3. (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. 1998; Font-Muñoz et al. 2019). 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.