## Abstract

Biological activity is often highly concentrated on surfaces, across the scales from molecular motors and ciliary arrays to sessile and motile organisms. These ‘active carpets’ locally inject energy into their surrounding fluid. Whereas Fick’s laws of diffusion are established near equilibrium, it is unclear how to solve non-equilibrium transport driven by such boundary-actuated fluctuations. Here, we derive the enhanced diffusivity of molecules or passive particles as a function of distance from an active carpet. Following Schnitzer’s telegraph model, we then cast these results into generalised Fick’s laws. Two archetypal problems are solved using these laws: First, considering sedimentation towards an active carpet, we find a self-cleaning effect where surface-driven fluctuations can repel particles. Second, considering diffusion from a source to an active sink, say nutrient capture by suspension feeders, we find a large molecular flux compared to thermal diffusion. Hence, our results could elucidate certain non-equilibrium properties of active coating materials and life at interfaces.

## Introduction

Imagine a particle floating in outer space, surrounded by uniformly distributed stars that all exert a gravitational force on it. On average, the total force is zero, by symmetry, but its variance is infinite, as described by the Holtsmark distribution^{1}. Microscopically, a similar situation occurs when molecular motors or swimming cells generate long-ranged flows and fluctuations in ‘living fluids’^{2,3,4,5,6,7,8,9,10,11,12}. These active suspensions operate far from thermal equilibrium where surprising effects emerge that fundamentally impact biological processes, including enhanced diffusion of enzymes^{13,14,15}, viscous and chaotic mixing by motility^{16,17,18,19,20,21,22,23}, viral infections^{24}, bioconvection^{25,26,27}, oxygen redistribution^{28}, nutrient uptake^{29,30,31,32} and communication via hydrodynamic trigger waves^{33}.

Besides active suspensions, biological activity is commonly concentrated on surfaces. Hence, ‘active carpets’ can form that drive systems out of equilibrium by injecting energy from boundaries. These carpets exist across the length scales: Inside cells, cytoskeletal motors generate cytoplasmic streaming while membranes host catalysing enzymes and transport proteins, all producing non-equilibrium fluctuations and active stresses^{34,35,36,37,38,39,40}. Outside cells, ciliary arrays create globally directed flows across entire organs but locally also facilitate mixing^{41,42,43,44,45,46,47}, which together may enhance pathogen clearance^{48}. Cells themselves accumulate on surfaces too, driving flows to replenish nutrients, including biofilms or microbial colonies^{49,50,51,52,53} and sessile suspension feeders^{54,55,56,57,58}. Multicellular organisms also drive feeding currents from walls, including sponges, reef corals^{59,60}, carpets of upside-down *Cassiopea* jellyfish^{61} and other macrobenthos^{62} at larger Reynolds numbers. Beyond biological hotspots, synthetic active carpets have been developed, including artificial cilia^{63,64}, self-propelled droplets and colloids^{9,10,65,66}, engineered bacterial carpets^{67,68}, molecular motility assays^{69}, light-controlled microfluidic flow networks^{70}, hydrogel actuators, liquid crystal elastomers and other responsive materials^{71,72,73}.

In this article, the properties of non-equilibrium diffusion of molecules or passive particles are examined near such active carpets that generate long-ranged flows. We distinguish between various types of carpets that inject energy into the surrounding medium in different ways: We consider actuators that exert a net force on the liquid (like pumping cilia or sessile suspension feeders), a net torque (like mixing rotors), or a net stress (like swimming bacteria). First, we show how the resulting active fluctuations decay with distance from the surface. Next, we derive the space-dependent diffusivities for the different carpet types, verified by hydrodynamic simulations, and the corresponding generalised Fick’s laws are written out. These laws are then used to solve the active analogue of two classic problems: We first examine sedimentation of particles towards an active carpet. Instead of following the Boltzmann distribution, the particles can be repelled by the fluctuations close to the surface, so the surface features an active self-cleaning effect. We then consider the diffusion of molecules (or particles) from a source to an active sink, such as a surface covered with sessile suspension feeders. The nutrient flux (or particle capture rate) is found to scale quadratically with the active forcing, so it can be significantly larger than the flux due to thermal diffusion. Overall, these problems can be described by relatively simple equations, even if the resulting dynamics feature unexpected solutions.

## Results

### Active carpet definition

We consider an active carpet made of actuators that generate flows at low Reynolds numbers near a planar no-slip surface (Fig. 1a–d). The surface is located at *z* = 0. The actuators have position \({\boldsymbol{r}}_{a}=({x}_{a},{y}_{a},h)\) and orientation unit vector \({\boldsymbol{p}}_{a}\). Each actuator *a* drives an individual flow, \(\boldsymbol{u}(\boldsymbol{r},{\boldsymbol{r}}_{a},{\boldsymbol{p}}_{a})\), evaluated at position \(\boldsymbol{r}\). In general, this flow can be written in terms of the Blake tensor \({\mathcal{B}}(\boldsymbol{r},{\boldsymbol{r}}_{a})\)^{74} and a multipole expansion thereof (see ‘Methods: Individual flow fields’).

For example, molecular motors walking or cells crawling along a substrate can entrain the surrounding fluid (Fig. 1a)^{34,35,36,37,39,40}. This can be described by a ‘parallel Stokeslet’, a point force *f*_{∥} oriented parallel to the surface, giving rise to the flow \({\boldsymbol{u}}_{\parallel }={f}_{\parallel }{\mathcal{B}}\cdot {\boldsymbol{p}}_{a}\). At a larger scale, sessile suspension feeders like *Vorticella* cells generate nutrient currents (Fig. 1b)^{54,55}, which can be described by a point force *f*_{⊥} oriented perpendicular towards the surface, \({\boldsymbol{u}}_{\perp }={f}_{\perp }{\mathcal{B}}\cdot {\boldsymbol{e}}_{z}\). Similarly, the flow generated by an *E. coli* bacterium swimming parallel to a wall (Fig. 1c) is described well^{50} by a ‘Stokes dipole’ flow, given by \({\boldsymbol{u}}_{\text{D}}=\kappa ({\boldsymbol{p}}_{a}\cdot {\boldsymbol{\nabla }}_{a})({\mathcal{B}}\cdot {\boldsymbol{p}}_{a})\), with dipole moment *κ*. Finally, a torque-generating actuator can be described by a Stokes rotlet (Fig. 1d), given by \({\boldsymbol{u}}_{\text{R}}=\varrho (({\boldsymbol{e}}_{x}\cdot {\boldsymbol{\nabla }}_{a})({\mathcal{B}}\cdot {\boldsymbol{e}}_{y})-({\boldsymbol{e}}_{y}\cdot {\boldsymbol{\nabla }}_{a})({\mathcal{B}}\cdot {\boldsymbol{e}}_{x}))/2\), with strength *ϱ*. This could represent rotors on a surface, such as bacteria with tethered flagella, or a carpet of nodal cilia^{45} that move around in circles in the *x**y* plane. A more complex example could be (non-nodal) airway cilia, beating almost perfectly in a plane but with some off-plane fluctuations^{48,75}. For such situations, to establish a more realistic description, one can combine terms from this multipole expansion, using Stokeslets, rotlets, dipoles and high-order terms as needed.

To be explicit, we have written out the full expressions of these first multipoles in ‘Methods: Individual flow fields’, and their corresponding flow fields are illustrated in Fig. 1e–h. Also note, throughout this paper we consider non-dimensional quantities in our simulations and figures, but all results can be predicted by analytical expressions, into which dimensional values can be inserted. Hence, we will discuss our results for typical numbers in biology.

Once the individual flows \(\boldsymbol{u}\) are established, the total flow \(\boldsymbol{v}={\sum }_{a}\boldsymbol{u}\) due to all *N*_{a} actuators is probed by a passive tracer particle as a function of its distance from the surface. For a given carpet architecture, which is defined by the probability density \(F({\boldsymbol{r}}_{a},{\boldsymbol{p}}_{a})\) of finding an actuator at position \({\boldsymbol{r}}_{a}\) and orientation \({\boldsymbol{p}}_{a}\), the average total flow evaluated at position \(\boldsymbol{r}\) is \(\langle \boldsymbol{v}(\boldsymbol{r})\rangle =\int {\boldsymbol{u}}F{\mathrm{d}}{\boldsymbol{r}}_{a}{\mathrm{d}}{\boldsymbol{p}}_{a}\), where 〈…〉 denotes averaging over a statistical ensemble of independent active carpet configurations. If there are any spatial or orientational gradients in the distribution *F*, for example due to bacterial clustering, or topological defect patterns, then long-ranged flows can emerge^{52}. However, in this paper we will consider cases where the actuators are uniformly distributed, so *F* is equal to a constant. Then the mean drift cancels out, so \(\langle \boldsymbol{v}\rangle =0\) in the absence of gradients in the carpet architecture.

### Active fluctuations

Even if the mean flow generated by the actuators is equal to zero, its variance at any one time is not. Hence, the flows can lead to ‘active fluctuations’ that push and pull on particles near the carpet. We first determine the strength of these fluctuations numerically. A carpet is simulated by placing *N*_{a} actuators that are randomly distributed with a uniform surface density *n* within the *x**y* plane. We then evaluate the total flow \(\boldsymbol{v}\) evaluated at particle position \({\boldsymbol{r}}_{0}=(0,0,{z}_{0})\). By repeating this for a large ensemble of *N*_{e} independent carpet configurations, we evaluate the distribution of the total flow, \(\,\text{PDF}\,(\boldsymbol{v})\), analogous to the Holtsmark distribution^{1}. Subsequently, the moments of this total flow distribution (the mean, variance, skewness and kurtosis) are found as a function of distance from the surface, for different carpet types (see ‘Methods: Characterising fluctuations: simulation details’).

We first focus on carpets made of Stokeslets oriented parallel to the surface, with uniformly distributed orientations, so *p*_{z} = 0 and *F* = *n*/2*π*. Figure 1i shows the resulting histogram of the total vertical flow, *v*_{z}, for different values of *z*_{0}. This distribution is an even function, by symmetry of the parallel Stokeslet, so it is immediately clear that the mean flow vanishes. However, the width is finite and decays with distance from the surface, so the active fluctuations are stronger closer to the carpet. This variance can be calculated analytically,

for the directions *i*, *j* ∈ *x*, *y*, *z*, as detailed in ‘Methods: Characterising fluctuations: theory details’. For the vertical component, for instance, this integral yields the variance

Physically, this expression represents how deeply the active boundary can influence the passive bulk fluid. This theoretical result is listed in Table 1 (row 4, column 3), and it is compared with simulations in Fig. 1m (magenta line). The decay with distance is \(1/{z}_{0}^{2}\) for all diagonal components, but the off-diagonal components are zero. Interestingly, the variance in the vertical direction (purple stars) is twice as strong as in the horizontal directions (green points, blue triangles), so \(\langle {v}_{x}^{2}\rangle =\langle {v}_{y}^{2}\rangle =\frac{1}{2}\langle {v}_{z}^{2}\rangle\) for parallel Stokeslets.

Moreover, when repeating this algebra for Stokeslets oriented perpendicular to the surface, we find that \(\langle {v}_{z}^{2}\rangle\) is five times stronger than \(\langle {v}_{x}^{2}\rangle\) (Table 1, fourth column) and the decay is now \(1/{z}_{0}^{4}\). This scaling also holds for the dipole flows (fifth column) and the rotlets (sixth column). Note that rotlets only generate flows in the *x**y* plane, so their resulting variance only has horizontal components. This fact could be exploited to tune the relative magnitude of \(\langle {v}_{x}^{2}\rangle\) and \(\langle {v}_{z}^{2}\rangle\). For example, one could use an active carpet made of both Stokeslets and rotlets and vary their relative prefactors *f*_{⊥} and *ϱ*, or vary their relative densities *n*_{⊥} and *n*_{ϱ}. Therefore, as we discuss next, the diffusion anisotropy may become a tunable.

The distributions of the active fluctuations, \(\,\text{PDF}\,(\boldsymbol{v})\), are not purely Gaussian. They feature skewness and kurtosis (Table 1; bottom rows), but they still obey the central limit theorem because the variance is finite for *z*_{0} > 0. Otherwise, Lévy flights must be considered^{12,16,76,77}. Additionally, these higher-order moments also decay with *z*_{0}, so far away from the active carpet the profiles are more Gaussian (Fig. 1i–l; insets).

### Space-dependent diffusivity

To understand how these active fluctuations may lead to particle diffusion, it is important to realise how exactly they vary over time. Most systems in nature are dynamic, such as bacteria (or synthetic self-propelled particles) that swim over a surface, with some source of stochasticity. As the bacteria move and reorient, the flows they produce change dynamically. Therefore, a particle advected by these flows will trace out a path that is eventually diffusive. As opposed to Brownian motion, though, this diffusion process is not the sum of random kicks that are instantaneous. Instead, the active fluctuations are correlated in time. They have a memory based on the history of the active carpet configuration.

Therefore, we first focus on a simple case with only one memory time scale, *τ*. Consider a carpet of sessile suspension feeders like *Vorticella convallaria*, tethered organisms that produce time-varying nutrient currents (Fig. 2a). These are modelled as *N*_{a} perpendicular Stokeslets that are again uniformly distributed in space, and fixed because they are non-motile. Each organism generates a flow with its own time-varying strength *f*_{⊥} = *f*(*t*), which is independent of the other organisms. These *N*_{a} forces fluctuate according to an Ornstein–Uhlenbeck process, \(\frac{{\mathrm{d}}f}{{\mathrm{d}}t}=-f/\tau +\sigma \eta (t)\), with relaxation time *τ* and strength *σ*. These behavioural quantities are set by the cells’ intrinsic properties, such as the internal biochemistry and biophysics. Supplementary Movie 1 depicts these actuator dynamics as well as the evolution of the total flow \(\boldsymbol{v}(t)\) they produce. We then simulate the motion of tracer particles that do not diffuse because of Brownian thermal fluctuations, but because of the fluctuating flows generated by the active carpet (see ‘Methods: Simulating the diffusivity of particles near carpets made of fluctuating forces’).

The resulting velocity correlation function (VCF) and mean-squared displacement (MSD) are shown in Fig. 2b and c, respectively. The tracer dynamics are ballistic at short times, *t* ≪ *τ*, called the Holtsmark regime. However, at long times they are diffusive with a linear relation 〈Δ*x*^{2}〉 = 2*D*_{xx}*t*, and similarly in the other directions. Hence, we determine the components of the diffusion tensor as a function of *z*_{0} (Fig. 2d). Like the flow variances, the vertical component *D*_{zz} is five times stronger than the horizontal components, leading to anisotropic diffusion, and they all decay as \(1/{z}_{0}^{4}\) in all three directions.

This system can be solved analytically when considering the limit of small displacements, \(\langle {{\Delta }}{\boldsymbol{r}}^{2}\rangle \ll {\boldsymbol{r}}_{0}^{2}\), when the noise amplitude is small. This ensures that we determine the local diffusivity, *D*_{ij}(*z*) with small variations in *z*. Using information from the variance of the active fluctuations and their temporal correlations, the motion can then be integrated (see ‘Methods: Derivation of the mean-squared displacement and space-dependent diffusivity’). This gives the MSD,

which captures both the short-term ballistic motion and the diffusivity after long times (Fig. 2c, dashed red lines). Thus, for the vertical Stokeslets for example, we find the space-dependent diffusion,

which is compared with the simulations in Fig. 2d. Because our theoretical approximation is formulated for small amplitudes of active fluctuations, the expression only holds far from the surface, when \(z\gg \root{6} \of {30\pi n{h}^{4}\langle {f}^{2}\rangle {\tau }^{2}}\) (Eq. (54)). Beyond this distance we find a good agreement between the simulations and the theory.

Overall, we find that the diffusion driven by an active carpet can be much stronger than Brownian thermal diffusion. Considering a carpet of *Vorticella* with cell radius *a* ~25 μm, *h* ~150 μm, *n* ~1/(100 μm)^{2}, *τ* ~1 s and 8*π**μ**f*_{⊥} ~500 pN^{55}, the active diffusion can be dominant up to distances of \({z}_{\text{th}}=\root 4 \of {15\pi n{h}^{4}\langle {f}^{2}\rangle \tau /{D}_{\text{th}}} \sim 1\) mm for small nutrient molecules of *D*_{th} ~ 500 μm^{2}/s. That is much larger than the cell size. Moreover, for micron-sized prey of *D*_{th} ~0.5 μm^{2}/s, we find *z*_{th} ~7 mm, orders of magnitude larger than the organism itself. Hence, we expect these results to be highly relevant across the scales, for non-equilibrium transport from molecular to organismic sizes.

### Diffusion due to motile actuators

These concepts can be extended to a more complex system with motile actuators (Fig. 2e–h). We consider Stokeslets that move with a constant velocity *V* along the surface, subject to rotational diffusion *D*_{r}, which exert on the liquid a constant force *f*_{∥} that is aligned with the direction of motion (see ‘Methods: Simulating the diffusivity of particles near carpets made of moving actuators’). Importantly, this system features multiple time scales: The reorientation time, *τ*_{r} ~ 1/*D*_{r}, and the time taken to move underneath a particle, *τ*_{u} ~ *z*_{0}/*V*. The corresponding tracer dynamics are therefore more complicated to solve in general, but for slow actuators the decorrelation of the carpet memory is primarily controlled by the rotational diffusion, since *τ*_{r} ≪ *τ*_{u} for small *V*. Indeed, as shown in Fig. 2f, in that case the VCF simulated for different positions *z*_{0} (green lines) closely follows the actuators’ rotational correlation function (red dashed). Then the same theory as before can be applied to predict the MSD (Fig. 2g) and the diffusivity (Fig. 2h). In particular, we can approximate the space-dependent diffusion for parallel Stokeslets as

and similarly for other actuator types following Table 1. In the remainder of this paper we stay in the small *V* limit, but future work should also address the case of faster actuators.

Again, the active diffusion is significant compared to thermal diffusion. Inserting some typical values into Eq. (5), say *h* ~ 1 μm, *n* ~ 1 μm^{−2}, *D*_{r} ~ 1 s^{−1} and 8*π**μ**f* ~ 1 pN, we have \({z}_{\text{th}}=\sqrt{6\pi n{h}^{2}{f}_{\parallel }^{2}/{D}_{\text{r}}{D}_{\text{th}}} \sim 7.7\,\upmu {\mathrm{m}}\) for molecules of *D*_{th} ~ 500 μm^{2}/s. Similarly, for micron-sized particles with *D*_{th} ~ 0.5 μm^{2}/s we have *z*_{th} ~ 240 μm, so again much larger than *h*. Since *n**h*^{2} remains constant when scaling down in size, the enhanced transport can be significant even for subcellular systems.

### Generalised Fick’s laws

Once the local diffusivity is determined, we can establish the generalised Fick’s laws (or Fokker-Planck equations) that govern the global stochastic transport. For clarity, we first revise the case of constant diffusion \(\tilde{D}\). Fick’s first law relates gradients in particle concentration \(\varphi (\boldsymbol{r},t)\) and an external drift flow \({\boldsymbol{v}}_{\text{d}}\) to the total flux, \({\boldsymbol{J}}={\boldsymbol{v}}_{\text{d}}\varphi -{\tilde{D}}\boldsymbol{\nabla }\varphi\). Fick’s second law describes how the particle concentration evolves in time, \({\partial }_{t}\varphi =-{\boldsymbol{\nabla }}\cdot ({\boldsymbol{v}}_{\text{d}}\varphi )+\tilde{D}{\nabla }^{2}\varphi\). This follows directly from the first law and the continuity equation, \({\partial }_{t}\varphi =-\boldsymbol{\nabla }\cdot \boldsymbol{J}\).

If the diffusivity is not constant, the spatial dependence may either arise from gradients in the mean speed or an inhomogeneous mean free path^{78}. Indeed, as we saw earlier (Eqs. (4) and (5)), the diffusion tensor \({D}_{ij}(z)={{\mathcal{V}}}_{ij}(z)\tau (z)\) can be written in terms of the variance of the active fluctuations \({\mathcal{V}}\) and the memory time *τ*, which in general can both depend on position. As described in ‘Methods: Generalised Fick’s laws’, the diffusion driven by an active carpet can be analysed by constructing a ‘telegraph model’, following Schnitzer^{78}. The particle flux is then described by the first generalised Fick’s law,

and the second generalised Fick’s law still follows from the continuity equation. Interestingly, these laws can be used to describe the active analogue of several classic problems, as we discuss in the next sections.

### Sedimentation towards an active carpet

As a first application, we consider particles sedimenting towards an active carpet (Fig. 3a, ‘Methods: Sedimentation towards an active carpet: simulation details’). The fluctuations are driven by slowly moving parallel Stokeslets, as before. Typical particle trajectories *z*(*t*) are shown for different values of *v*_{g}, the sedimentation velocity (Fig. 3b). As expected, the particles with the smallest *v*_{g} reach the highest *z* positions (green track) but, surprisingly, they hover above the active carpet at some finite height. This phenomenon is somewhat reminiscent of the Leidenfrost effect, where droplets hover above a hot plate. That is, far from the active carpet the particles fall under gravity, but nearby they are repelled by strong active fluctuations (also see Supplementary Movie 2). Because the system is driven out of equilibrium, the sedimentation profile *φ*(*z*) does not follow the Boltzmann distribution with an exponential decay. Instead, it features a maximum value (Fig. 3c) with much less particles close to the surface than expected.

To quantify this ‘self-cleaning’ effect, we solve the generalised Fick’s laws for this system analytically. For active fluctuations that scale algebraically with distance from the carpet, \({\mathcal{V}}(z)=\tilde{{\mathcal{V}}}/{z}^{\alpha }\) as in Table 1, and similarly for the memory time, \(\tau (z)=\tilde{\tau }/{z}^{\beta }\), the vertical component of the generalised flux (Eq. (6)) becomes

where the constant part of the vertical diffusivity is \(\tilde{D}={\tilde{{\mathcal{V}}}}_{zz}\tilde{\tau }\). We require that *J*_{z} = 0 at steady state, which yields the non-Boltzmannian sedimentation profile

where *φ*_{0} is a normalisation factor. To clarify, *v*_{g} has standard units of m/s and \(\tilde{D}\) has units of m^{2+α+β}/s. In the limit of a constant diffusivity (*α* = *β* = 0) we recover the Boltzmann distribution, \(\varphi (z)={\varphi }_{0}{\mathrm{e}}^{-{v}_{\text{g}}z/\tilde{D}}\). A major different with the sedimentation profile near an active carpet is the *z*^{α/2} factor in Eq. (8), where *α* = 2 for parallel Stokeslets. Therefore, the sedimentation profile features a maximum (Fig. 3c), which is located at

These results agree with our simulations (Fig. 3c, d), for different values of the sedimentation velocity.

The self-cleaning effect can be understood by analysing the generalised flux more carefully (see ‘Methods: Self-cleaning effect’). The first two terms are identical to the case of constant diffusion. However, a third term emerges that represents a flux toward regions where the speed is low. Since the active fluctuations decay with distance, this flux term is always directed away from the active carpet. Importantly, the self-cleaning effect also persists when thermal fluctuations are included in the analysis (see ‘Methods: Sedimentation with active and thermal diffusion’).

### Diffusion from a source to an active carpet sink

Next, we consider particles diffusing from a source to an active carpet sink (Fig. 4a, ‘Methods: Diffusion from a source to an active carpet sink: simulation details’), which could represent particle capture by sessile suspension feeders, or substrate molecules being catalysed by a carpet of enzymes (also see Supplementary Movie 3). Every time a simulated particle diffuses across the absorbing boundary at *z*_{sink} = *h*, we place it back at *z*_{source} = *H*. This way, a diffusive flux emerges towards the carpet, which is equivalent to the particle capture rate. Typical trajectories *z*(*t*) show that the particles spend the majority of time close to the top surface (Fig. 4a), which is also apparent in the concentration profiles *φ*(*z*) that are peaked at the top surface (Fig. 4b). Moreover, the mean first-passage time is much larger for distant sources, so the flux decays rapidly with gap size (Fig. 4c), which may be important for biology or applications in confined spaces.

This system can again be solved analytically (see ‘Methods: Diffusion from a source to an active carpet sink: theory details’). The flux is given by Eq. (7), with *v*_{g} = 0 and with fixed particle concentrations at the source and sink. Hence, we must solve the continuity equation ∂_{z}*J*_{z} = 0 subject to the boundary conditions *φ*(*H*) = *φ*_{+} and *φ*(*h*) = 0. This gives the solution

for *α* ≥ 0 and *β* ≥ −1, or a slightly more complex function for other values. The corresponding solution for the flux, equivalent to the particle capture rate, is

In the limit of constant diffusion and *h* ≪ *H*, we recover the profile *φ*(*z*) ≈ *φ*_{+}*z*/*H*, which increases linearly from sink to source, and the flux \({J}_{z}\approx -{\varphi }_{+}\tilde{D}/H\). However, for an active carpet the profile is more concentrated at the source, in agreement with the simulations (Fig. 4b). Indeed, the flux also decays more rapidly than thermal diffusion with the gap size (Fig. 4c). These predictions (dashed lines) agree well with the simulations.

Finally, we also examine the effect of different force strengths *f*_{∥} on the diffusion to the active sink. Interestingly, we find that the flux scales quadratically with the active fluctuations, so stronger Stokeslets lead to a much larger capture rate (Fig. 4d). Inserting the typical values written below Eq. (5) into Eq. (11), with source concentration *φ*_{+} ~ 1 particle/μm and gap size *H* ~ 10 μm, we find an active flux of ∣*J*_{z}∣ ~ 60 s^{−1}. That is more than the thermal flux for micron-sized and molecular particles, ∣*J*_{z}∣ ~ 0.05 and 50 s^{−1}, respectively. Interestingly, the thermal and active fluctuations can co-operate to give an even larger flux of ~128 s^{−1} together (‘Methods: Diffusion from a source to an active carpet sink: comparison with thermal diffusion’). Hence, by enhancing the local diffusivity, the active carpet may increase its nutrient uptake significantly.

### Advective and diffusive transport

Until now, we have considered active carpets that are spatially uniform, where \(F({\boldsymbol{r}}_{a},{\boldsymbol{p}}_{a})\) is constant so the average flow vanishes, \(\langle {\boldsymbol{v}}\rangle =0\). However, any natural carpet is likely to feature some heterogeneity in its force distribution that can drive local advection flows (see Methods: Advective and di usive transport’). This imposes a constraint on the generalised Fick’s laws we discussed so far: If the particles are stuck in local advection currents, they cannot diffuse around freely. Interestingly, this is also related to the occurrence of quenched disorder (‘Methods: Quenched disorder’). The key question is then how strong these advection flows are compared to the actively driven diffusion.

To investigate this, we consider a carpet composed of a square lattice of perpendicular Stokeslets that all fluctuate about a non-zero mean force \(\bar{f}\) with variance Var(*f*). The actuators are located at (*i**λ*, *j**λ*, *h* = 1) where *i*, *j* are integer numbers and *λ* is the lattice spacing. The total flow driven by these actuators is then described by an advective contribution and a diffusive contribution, \(\boldsymbol{v}={\boldsymbol{v}}_{\text{adv}}+{\boldsymbol{v}}_{\text{diff}}\). As shown in Fig. 5, and described in detail in ‘Methods: Advective and diffusive transport’, the active diffusion is stronger than the advection beyond a specific distance from the carpet (Eq. (92)). As mentioned earlier, we also require *z* to be large (Eq. (54)) in order for the time scale of diffusive transport to be slow compared to the time scale of the active fluctuations themselves. Indeed, this is the case in many biological and engineered settings, but care should be taken that these conditions are satisfied.

## Discussion

In summary, active carpets are ubiquitous in nature, from molecular and cellular to organismic length scales. In this paper we described how these active carpets can drive non-equilibrium fluctuations by locally injecting energy into their surrounding fluid. These fluctuations were quantified with a general theoretical framework in terms of the fundamental solutions of the Stokes equation. Hence, we derived the diffusivity as a function of distance from an active carpet, and we found the corresponding generalised Fick’s laws. The predictive power of these laws was demonstrated by solving them for two archetypal problems: sedimentation and diffusion towards an active sink. Of course, one could extend this straightforwardly to other problems like diffusion away from an active source, potentially in combination with sedimentation, or systems with multiple active boundary conditions in one or more dimensions.

This framework can be applied to a much more general class of problems. Until now we have considered non-interacting actuators, where the fluctuating forces and associated time scales are determined by their intrinsic properties. It would be interesting to consider actuators featuring different types of collective behaviours^{2,3,4,5,6,7,8,9,10,11}, say swarming and clustering, so these forcing properties might emerge collectively. This could lead to interesting types of anomalous diffusion^{79}.

Additionally, one could consider active transport and diffusion in complex geometries (see ‘Methods: Extension to more complex geometries’). Immediate extensions could be carpets on the outside of a sphere, such as cilia covering *Volvox carteri*^{30} and microbes covering marine snow^{80,81}, but also carpets on the inside of a spherical cavity, such as motors on the cytoskeleton^{11}. Once these topologies are understood, one could consider more complex curved spaces, like the highly folded active membranes of the endoplasmic reticulum. Furthermore, besides translational diffusion one could derive the rotational diffusion^{82,83}, for spherical but also for ellipsoidal particles and other shapes^{84}. Finally, flows with finite inertia could be modelled to describe carpets with actuators that are large or perform rapid motions^{33,85}.

Besides hydrodynamic actuators, this framework could equally be applied to different kinds of energy injection from surfaces. For example, one could consider active carpets that generate fluctuating electric or magnetic fields, or catalyse chemical reactions. Conversely, there are also surface-driven thermo-osmotic and acoustic fluctuations, or liquid interfaces with a variable surface tension, which can all induce transport and diffusion.

The understanding of active carpets may therefore prove useful in a wide range of applications. Recent advances in nanotechnology enable the design of artificial cilia^{63,64}, origami micromachines^{86,87}, and four-dimensionally printed active materials^{88}. One may therefore think of ‘active coating materials’, for example with self-cleaning properties based on the described effect of boundary-actuated fluctuations repelling particles. One could also consider biomolecular condensates^{89,90}, including phase-separated biopolymers, where chemical reactions are concentrated in one phase but the products are transported to the other phase by mixing actuators at the interface. This idea of interfacial mixing may also be applied to bacteria moving on water–oil interfaces^{91,92} for bioremediation of oil spills^{93,94}.

## Methods

### Individual flow fields

The flow generated by an individual actuator is written in terms of the Blake tensor^{74}, given by

where the flow is evaluated at position \(\boldsymbol{r}\) and the actuator is located at position \({\boldsymbol{r}}_{a}=({x}_{a},{y}_{a},{z}_{a}=h)\) and oriented along \({\boldsymbol{p}}_{a}=({p}_{x},{p}_{y},{p}_{z})\). The indices *i*, *j*, *k* ∈ {*x*, *y*, *z*} denote Cartesian components, the mirror matrix is M_{jk} = diag(1, 1, −1), and the Oseen tensor is

with distance \(\boldsymbol{d}=\boldsymbol{r}-{\boldsymbol{r}}_{a}\). All the derivatives ∇_{a} of the Oseen tensor in Eq. (12) are taken with respect to the actuator position \({\boldsymbol{r}}_{a}\). As described for example in ref. ^{51}, the individual flows can then be written as a multipole expansion of this Blake tensor.

The flow due to a point force *f*_{∥} oriented parallel to the surface is given by

Similarly, the flow due to a point force *f*_{⊥} oriented perpendicular to the surface is

Note, we have scaled the forces as *f* = *f*_{N}/(8*π**μ*), where *μ* ~ 10^{−3} Pa s is the viscosity of water, so *f*_{N} are in Newtons and *f* have units [m^{2}/s]. The Stokes dipole flow with prefactor *κ* in [m^{3}/s] is found by taking the derivative along the \({\boldsymbol{p}}_{a}\) direction,

Finally, the Stokes rotlet describes an actuator that generates a torque about the *z*-axis. Its flow with prefactor *ϱ* in [m^{3}/s] is defined by

These flows fields are quite complex when evaluated algebraically. Therefore, to make progress we approximate the flows as being evaluated far from the surface, such that *z*_{a} = *ϵ**z* with *ϵ* ≪ 1. For each of the individual actuator flows ((14)–(17)) we then perform a Taylor expansion,

and we only keep the leading-order contribution.

To be explicit, for the parallel Stokeslets evaluated at \(\boldsymbol{r}=(0,0,{z}_{0})\), and using \({\boldsymbol{r}}_{\boldsymbol{a}}=({\rho }_{a}\cos {\theta }_{a},{\rho }_{a}\sin {\theta }_{a},h)\) and \({\boldsymbol{p}}_{\boldsymbol{a}}=(\cos {\phi }_{a},\sin {\phi }_{a},0)\), this gives the Cartesian components

Similarly, for the perpendicular Stokeslets with the same position and orientation \({{\boldsymbol{p}}_{\boldsymbol{a}}}=(0,0,1)\) we have

For Stokes dipoles oriented parallel to the surface, with orientation \({\boldsymbol{p}}_{\boldsymbol{a}}=(\cos {\phi }_{a},\sin {\phi }_{a},0)\), we have

Finally, for the rotlet dipoles we have

and so forth for higher-order multipole flows.

### Characterising fluctuations: simulation details

To characterise the strength of the active fluctuations, we simulate the distribution \(\,\text{PDF}\,(\boldsymbol{v})\) of the total flow evaluated at particle position \({\boldsymbol{r}}_{0}=(0,0,{z}_{0})\). The carpet consists of *N*_{a} = 4*n**L*^{2} actuators that are randomly distributed with uniform surface density *n* within *x*, *y* ∈ [ − *L*, *L*], with carpet size *L* and vertical position *z*_{a} = *h*. Once the carpet configuration is sampled with standard random number generators, we evaluate the total flow \(\boldsymbol{v}={\sum }_{a}\boldsymbol{u}\) due to all *N*_{a} actuators. We repeat this simulation to obtain an ensemble of *N*_{e} total flow velocities, each due to an independent carpet realisation. Hence, we evaluate the distribution \(\,\text{PDF}\,(\boldsymbol{v})\) and its moments as a function of distance from the carpet, *z*_{0}. In non-dimensionalised simulation units, we use the parameters *h* = 1, *n* = 0.1, *N*_{e} = 10^{4}, and *L* = 500 so *N*_{a} = 10^{5}. We have verified that *L* is large enough to avoid edge effects. The actuator strengths used for the four types are \({f}_{\parallel }={f}_{\perp }=\frac{3}{4}\) and *κ* = *ϱ* = 30, respectively. These simulations are presented in Fig. 1. Note the results are shown in simulation units in order to keep the description general. However, for any specific application, dimensional quantities (with standard SI units) can be inserted in all the equations throughout the paper.

### Characterising fluctuations: theory details

We consider active carpets made of uniformly distributed actuators, so the average flow vanishes in the absence of gradients, \(\langle {\boldsymbol{v}}\rangle =0\). This may be demonstrated by integrating any of the expressions Eqs. (19)–(30) with a constant carpet distribution, *F* = *n*/2*π*, which yields the average total flow

The variance tensor of the active fluctuations is calculated in the same way, by evaluating the integral

To give an explicit example, the vertical component of the variance for parallel Stokeslets is given by

This result corresponds to Eq. (2) in the main text, and is compared with simulations in Fig. 1m. We repeat this calculation for the different actuator types, and for the different components, *i*, *j*. Note that the off-diagonal components vanish, so the only non-trivial results are *i* = *j*. These results are all listed in Table 1.

Besides variance, the distributions of the active fluctuations also feature skewness

and kurtosis,

For all these quantities, we find that the only non-zero results are diagonal, where *i* = *j* = *k* = *l*. These results are listed in the bottom rows of Table 1 for all actuator types, and plotted in the insets of Fig. 1i–l.

### Simulating the diffusivity of particles near carpets made of fluctuating forces

We first consider the motion of a tracer particle above a carpet of fluctuating forces, composed of *N*_{a} perpendicular Stokeslets (Eq. (22)) that have a fixed position on the wall. The point forces have vertical position *h* = 1 and horizontally they are randomly distributed over the surface, with uniform density *n* = 0.1 per unit area, and within a simulation box of dimensions *x*, *y* ∈ [−*L*, *L*] with *L* = 500 so *N*_{a} = 4*n**L* = 10^{5}. We have verified that the simulation box is large enough to avoid edge effects by testing that the results are independent of large *L*. The tracer is initially located at (0, 0, *z*_{0}). The forces do not interact with each other. Each perpendicular Stokeslet force *f*_{⊥} = *f*_{a}(*t*) evolves dynamically in time following its own independent Ornstein–Uhlenbeck process.

This Ornstein–Uhlenbeck process is defined as

where *τ* is a relaxation time, *σ* is a constant that relates to the force strength, and *η* is Gaussian white noise with 〈*η*(*t*)〉 = 0 and 〈*η*(*t*_{1})*η*(*t*_{2})〉 = *δ*(*t*_{1} − *t*_{2}). The solution of this stochastic differential equation is

This corresponds to an overdamped relaxation process driven by fluctuations, which admits a Gaussian stationary distribution

with mean 〈*f*〉 = 0 and force magnitude \(\langle {f}^{2}\rangle =\frac{{\sigma }^{2}\tau }{2}\). Importantly, the temporal correlation function at steady state for this Ornstein–Uhlenbeck process is also known exactly,

In our simulations, we use an average force magnitude of 〈*f*^{ 2}〉 = (3/4)^{2} and *τ* = 1/10 so that *σ*^{2} = 2〈*f*^{ 2}〉/*τ*, and the initial force *f*(*t* = 0) for each actuator is drawn randomly from the stationary distribution (Eq. (41)). Besides these active fluctuations, we do not include additional Brownian thermal fluctuations here.

The tracer’s equation of motion is then given by

which is integrated numerically, along with the *N*_{a} Ornstein–Uhlenbeck processes for all the actuators, using a fourth-order Runge–Kutta scheme. This gives one tracer trajectory. We repeat this for an ensemble of *N*_{e} = 100 independent trajectories, each with an independent carpet realisation composed of actuators that are distributed at different random but fixed positions, and with independent Ornstein–Uhlenbeck processes. We then repeat all this for each starting position *z*_{0}.

By averaging over the ensemble of *N*_{e} trajectories we compute the VCFs, \(C(| t^{\prime} -t^{\prime\prime} | )=\langle {v}_{i}(\boldsymbol{r},t^{\prime} ){v}_{j}(\boldsymbol{r},t^{\prime\prime} )\rangle\), shown in Fig. 2b for different values of *z*_{0} ∈ [2, 20]. Similarly, we determine the ensemble-averaged mean-square displacements, 〈Δ*r*_{i}Δ*r*_{j}〉, shown in Fig. 2c. From the latter data we also determine the components of the diffusion tensor, by statistical regression of the expression 2*D*_{ij}*t* = 〈Δ*r*_{i}Δ*r*_{j}〉 in the regime *t* ≫ *τ*, the diffusive limit, which is shown in Fig. 2d.

### Derivation of the MSD and space-dependent diffusivity

Analytical progress is made by considering small tracer displacements,

Then their equation of motion reduces to \(\frac{{\mathrm{d}}\boldsymbol{r}}{{\mathrm{d}}t}\approx \boldsymbol{v}({\boldsymbol{r}}_{0},t)\), plus higher-order terms that can be expanded as a power series in 1/*z*_{0}^{51}. Consequently, the MSD is given by

where we identify the VCF of the total flow at position \({\boldsymbol{r}}_{0}\), that is

The relationship between this ensemble-averaged flow correlation and the Ornstein–Uhlenbeck force correlations is given by

where we used Eq. (42). This expression is shown as a red dashed line in Fig. 2b. Hence, using the variance 〈*v*^{2}〉 = 〈*v*(0)*v*(0)〉 of the active fluctuations (Eq. (33)), for example, for vertical displacements due to the perpendicular Stokeslets, the MSD simplifies to

and similarly for other actuator types. After performing the final integrals, we obtain the MSD

and similarly for the other directions 〈Δ*r*_{i}Δ*r*_{j}〉. This MSD transition from ballistic to diffusive motion is compared with simulations in Fig. 2c (dashed line). The corresponding anisotropic diffusivity is then given by

These expressions are shown in Fig. 2d. The same analysis can also be applied to all other actuator types using the variances listed in Table 1.

Importantly, we should come back to the initial assumption of small displacements (Eq. (44)) and analyse its consequences. Rewriting this condition as *z*^{2} ≫ 2*D*_{zz}*t* = 30*π**n**h*^{4}〈*f*^{ 2}〉*τ**t*/*z*^{4} using Eq. (53), we can rearrange this for a temporal condition, *t*/*τ* ≪ *z*^{6}/(30*π**n**h*^{4}〈*f*^{ 2}〉*τ*^{2}). We also require that *t* ≫ *τ* for the particle motion to transition from ballistic to diffusive motion (see Eq. (51)). In other words, the theory is only expected to hold when the time scale of diffusive transport is slow compared to the time scale of the active fluctuations themselves. This is true far from the surface, where

and similarly for other actuator types. Inserting the values used for simulations in Fig. 2a–d, being *n* = 0.1, *h* = 1, 〈*f*^{ 2}〉 = (3/4)^{2}, *τ* = 0.1, we find *z* ≫ 0.61. This condition ensures that we determine the local diffusivity, *D*(*z*) with small variations in *z*. Once this local diffusivity is determined from the MSDs within this limit, the global stochastic dynamics of particles leaving this local area (with large variations in *z*) can be solved using the space-dependent Fick’s laws, as discussed in ‘Generalised Fick’s laws’.

### Simulating the diffusivity of particles near carpets made of moving actuators

Here, we consider the motion of a tracer particle above a carpet of *N*_{a} parallel Stokeslets (Eq. (19)) that each move with velocity \(\boldsymbol{V}=V{\boldsymbol{p}}_{a}\) along their director, \({\boldsymbol{p}}_{a}\), with a constant speed *V*. They generate a flow with force \(\boldsymbol{f}={f}_{\parallel }{\boldsymbol{p}}_{a}\). The Stokeslets move along the surface, so *p*_{z} = 0, with an orientation angle *ϕ*_{a} = arctan(*p*_{y}/*p*_{x}) that is subject to rotational diffusion. That is, \(\frac{{\mathrm{d}}{\phi }_{a}}{{\mathrm{d}}t}=\sqrt{2{D}_{\text{r}}}\eta (t)\) where *η* is white Gaussian noise with 〈*η*(*t*)〉 = 0 and 〈*η*(*t*_{1})*η*(*t*_{2})〉 = *δ*(*t*_{1} − *t*_{2}). The parallel Stokeslets are all independent of each other. We use the parameters *V* = 1, *D*_{r} = 10, and *f*_{∥} = 3/4, such that *τ*_{r} ≪ *τ*_{u} for all *z*_{0} ∈ [2, 20]. As before, the actuators are distributed randomly with uniform density *n* = 0.1 in a simulation box, *x*, *y* ∈ [−*L*, *L*], with *L* = 500 and *h* = 1. Periodic boundary conditions are imposed so that the surface density *n* = 0.1 of the actuators remains constant. Again, besides these active fluctuations, we do not include additional Brownian thermal fluctuations here.

These *N*_{a} + 1 equations of motion, for the tracer and the actuators, are then integrated as before over a time period *t* ∈ [0, 10] for each tracer trajectory. We repeat this simulation for an ensemble of *N*_{e} = 100 independent trajectories, each with an independent carpet realisation with actuators that are initially distributed at different random positions and orientations. Again, we repeat all this for ten different initial positions *z*_{0} ∈ [2, 20]. The ensemble-averaged results are shown in Fig. 2e–h.

### Generalised Fick’s laws

For clarity, we first revise the case of constant diffusion \(\tilde{D}\) in one spatial dimension, *z*. Then, Fick’s first law relates gradients in concentration *φ*(*z*, *t*) and an external drift *v*_{d} to the total flux,

Together with the continuity equation, ∂_{t}*φ* = −∂_{z}*J*_{z}, this gives Fick’s second law,

which is also known as the Fokker–Planck equation. This is equivalent to motion described by the following Langevin equation, \(\frac{{\mathrm{d}}z}{{\mathrm{d}}t}={v}_{\text{d}}+\sqrt{2\tilde{D}}\eta (t)\), where *η* is Gaussian white noise as defined below Eq. (39).

Rather than being constant in space, the diffusivity of tracers near an active carpet depends continuously on position. Such stochastic processes can often be described by an effective Smoluchowski equation^{78}, rather than standard Langevin methods which make no reference to individual collisions. Here, we follow this approach as a foundation for the generalised Fick’s laws that describe the diffusion of a tracer particle as a function of distance from an active carpet. Since we only have gradients in the vertical direction, we use the short-hand notations *D*(*z*) = *D*_{zz} for the vertical diffusivity and \(v(z)=\sqrt{\langle {v}_{z}^{2}\rangle }=\sqrt{{{\mathcal{V}}}_{zz}}\) for the mean vertical speed.

With this spatial dependence, the question arises whether the second term in Eq. (55) should be interpreted as *D*∂_{z}*φ* or ∂_{z}(*D**φ*), or something in between. This question is not well posed, because it depends on the physical processes in question. In other words, generalising this expression for macroscopic quantities requires partial knowledge of the microscopic mechanism for diffusion. In particular, information is needed about the spatial dependence of the memory time *τ*(*z*), and of the mean vertical speed, *v*(*z*). The diffusivity can then be written as the combination of these ingredients, *D*(*z*) = *v*^{2}(*z*)*τ*(*z*), as shown in Eq. (53). Indeed, either a larger average speed or a longer time between reorientations can give rise to a larger diffusivity. Note that for Ornstein–Uhlenbeck forcing the time scale *τ* does not depend on position, but this need not be true in general. For example, for rapidly moving actuators the smallest decorrelation time is *τ* ~ *z*/*V*. On the other hand, for slowly moving actuators the rotational diffusion constant *D*_{r} sets the smallest memory time.

Using this information about the microscopic interactions, a ‘telegraph model’^{78} can be constructed that describes the space-dependent diffusion process. Inspired by this model, we consider two particle populations, with population densities *φ*_{up}(*z*, *t*) and *φ*_{down}(*z*, *t*), respectively, that either move up or down along *z* with mean speed *v*(*z*) due to the active carpet fluctuations. The mean speed is the same for the two populations at any given *z* since we consider a uniform carpet without net drift, \(\langle {\boldsymbol{v}}\rangle =0\), as shown in Eq. (32). The particles switch directions at a mean rate of \(\frac{1}{2}{\tau }^{-1}\), set by the memory time of the active carpet. The total particle density is then given by *φ*(*z*, *t*) = *φ*_{up} + *φ*_{down}, and the diffusive flux of particles is *J*_{diff}(*z*, *t*) = *v*(*φ*_{up} − *φ*_{down}). Subsequently, using continuity of particle flow and conservation of particle number, the up and down populations evolve according to

In each equation, the first term describes spatial gradients in the moving populations, while the last two terms describe the switching between particles moving up and down, and vice versa. By adding and subtracting, the equations (57) can be rewritten in terms of the total density and diffusive flux, giving

Taking the time derivative of the first expression and combining with the second expression yields

Then, we assume that the high-frequency behaviour of particle movements can be neglected, so the second time derivative on the left-hand side vanishes,

Integrating this expression, we can solve for the diffusive flux *J*_{diff}. The constant of integration is set equal to zero to ensure that *J*_{diff} vanishes when the variance of the active fluctuations *v*^{2} are zero. If the particles are sedimenting with a constant drift velocity *v*_{d}, there is an additional flux *J*_{drift} = *v*_{d}*φ*, so the total flux is *J*_{z} = *J*_{diff} + *J*_{drift}. Then, combining all this information gives the first generalised Fick’s law in the vertical direction,

We find the other components of the flux by repeating the telegraph model analysis in the *x* and *y* directions, which gives

Since the variance tensor only has diagonal components, we can write the generalised flux in three dimensions as

which is written in the main text as Eq. (6). Note that the last term only has a vertical component for systems that obey translational invariance along the horizontal directions, when the variance tensor only depends on *z*. Finally, using the continuity equation we obtain the second generalised Fick’s law,

where repeated indices are summed over.

### Sedimentation towards an active carpet: simulation details

Here, we consider the motion of a sedimenting tracer particle above a carpet of moving parallel Stokeslets. For spherical particles, the sedimentation is described by a constant drift velocity \({\boldsymbol{v}}_{\text{d}}=-{v}_{\text{g}}\hat{\boldsymbol{z}}=\frac{{d}^{2}{{\Delta }}\rho }{18\mu }\boldsymbol{g}\), in terms of the gravitational acceleration \(\boldsymbol{g}\), the particle diameter *d*, its density difference with the medium Δ*ρ*, and the medium viscosity *μ*. The equations of motion of the moving forces are as described in ‘Methods: Simulating the diffusivity of particles near carpets made of moving actuators’, with the parameters *V* = 1, *D*_{r} = 10, *n* = 0.1, *L* = 500, *h* = 1 and *f*_{∥} = 10. For the tracer equation of motion we add the sedimentation, and we introduce a reflecting boundary condition at *z* = *h* to prevent the particles from crossing the active carpet. This system is integrated numerically for different sedimentation velocities, *v*_{g} ∈ [10^{−2}, 1]. We simulate over a long period of time, *t* ∈ [0, 10^{5}], to ensure that the sedimentation profile is well sampled. To clarify, we do not average this sedimentation profile over a statistical ensemble of independent carpet configurations. Therefore, since the only averaging is temporal, the results are informative about the dynamics of a given system. We then normalise this particle concentration profile,

where *N*_{p} = ∫*φ*d*z* is the number of tracers in the system, so ∫Φ(*z*)d*z* = 1. From these distributions we also evaluate the maximum \({z}_{\max }\) where \(\frac{{\mathrm{d}}\varphi }{{\mathrm{d}}z}=0\). These results are shown and compared with our analytical predictions in Fig. 3.

### Self-cleaning effect

These simulation results can be understood using the generalised Fick’s laws we discussed earlier. The first two terms on the RHS in Eq. (63) are identical to the case of constant diffusion (Eq. (55)), describing a flux of particles towards regions of low concentration. A third emerges, however, which describes diffusion towards regions of low fluid speed. Since the active fluctuations decay with distance (e.g. see Eq. (2)), this term leads to a flux directed away from the carpet.

To understand this better, we must quantify the contributions of the flux. Since the variance of all the active fluctuations in Table 1 feature an algebraic decay, we write

Similarly, for the memory time we write

because this algebraic form is common in natural systems. To name a few examples: *β* = 0 corresponds to a constant memory time. *β* = −1 corresponds to the time scale *τ* ~ *z*/*v*_{a} associated with an actuator moving underneath a tracer particle. *β* = −2 corresponds to the time scale *τ* ~ *z*^{2}/*D*_{a} associated with actuators diffusing underneath a tracer particle. Note that one could have inverted the exponent to keep *β* positive, but this does not change anything physically. Hence, we prefer to keep the expressions for *α* and *β* consistent.

Then the vertical diffusivity is \(D(z)=\tilde{D}/{z}^{\alpha +\beta }\), where the constant part is \(\tilde{D}={\tilde{{\mathcal{V}}}}_{zz}\tilde{\tau }\). For example, for slowly moving parallel Stokeslets we have \({\tilde{{\mathcal{V}}}}_{zz}=6\pi n{h}^{2}{f}_{\parallel }^{2}\) and \(\tilde{\tau }={D}_{\,\text{r}\,}^{-1}\) with *α* = 2 and *β* = 0 from Eq. (5). Inserting the expressions (70)–(71) into (63), one obtains

The first term is negative for sedimentation, and the second term still describes ordinary diffusion towards regions of low concentration. But the third term is always positive, repelling particles away from the carpet, which explains the self-cleaning effect.

To solve the sedimentation profile, we require that *J*_{z} = 0 at steady state, which yields

where *φ*_{0} is a normalisation factor. In the limit of a constant diffusivity (*α* = *β* = 0) we recover the Boltzmann distribution, \(\varphi (z)={\varphi }_{0}{\mathrm{e}}^{-{v}_{\text{g}}z/\tilde{D}}\). This is no longer true for active carpets with decaying fluctuations because of the *z*^{α/2} factor, where *α* = 2 for parallel Stokeslets. Therefore, the sedimentation profile features a maximum (Fig. 3c), which is located at

These results agree with our simulations (Fig. 3c, d), for different values of the sedimentation velocity.

### Sedimentation with active and thermal diffusion

Besides the fluctuating flows generated by the active carpet, the particles may also experience Brownian thermal fluctuations. This thermal diffusion *D*_{th} can be included explicitly in the generalised Fick’s law:

Then, the expression *J*_{z} = 0 can still be solved analytically to determine the steady-state sedimentation profile. For parallel Stokeslets, for example, with *α* = 2 and *β* = 0, we find the solution

It is important to note that this function has the same shape as the original solution (Eq. (73)). Indeed, particles are still repelled from the active surface, \({\mathrm{lim}\,}_{z\to 0}\varphi (z)=0\), and the function has a maximum at the same location as before, at \({z}_{\max }={(\tilde{D}/{v}_{\text{g}})}^{1/3}\) for all *D*_{th} ≥ 0. Thus, the self-cleaning effect is not affected by thermal diffusion. Of course, when the surface activity vanishes, \(\tilde{D}\to 0\), we recover the Boltzmann distribution.

### Diffusion from a source to an active carpet sink: simulation details

In this section, we consider the dynamics of particles that are spawned at a source and absorbed by an active sink. The particles are subject to active fluctuations due to a carpet of slowly moving parallel Stokeslets. The equations of motion of the *N*_{a} actuators are as described in ‘Methods: Simulating the diffusivity of particles near carpets made of moving actuators and Sedimentation towards an active carpet: simulation details’. For the tracer equation of motion we remove the sedimentation, we impose an absorbing boundary condition at *z*_{sink} = *h*, and a reflecting boundary condition at *z*_{source} = *H*. Whenever a particle is absorbed by the sink, we place it back at the source, at *x* = *y* = 0, and we redistribute all the actuators with new random positions and orientations to start a new fully independent trajectory. We run two separate types of simulations: First, the gap size is varied with different values of *H* ∈ [2, 20] with constant force *f*_{∥} = 10. Second, we vary *f*_{∥} = ∈ [1, 10] with constant source height *H* = 5. For each of these forces we measure the flux *J*_{z}, defined as the number of particles that diffuse from the source to the sink per unit time, *J*_{z} = −*N*_{p}/〈*t*_{mfp}〉, where 〈*t*_{mfp}〉 is the mean first-passage time and *N*_{p} is the number of particles. By symmetry, we expect the nutrient flux to scale quadratically with the force, because the diffusion equations are invariant under the transformation *f* → −*f*. Indeed, this is also observed in the simulations, as shown in Fig. 4.

### Diffusion from a source to an active carpet sink: theory details

To solve the system of non-equilibrium diffusion from a source to an active sink, we again consider the vertical flux given by Eq. (7). This time, the sedimentation velocity is equal to zero and we seek the steady-state solution (∂_{t}*φ* = 0) with fixed particle concentrations at the source and the sink. Hence, we must solve the continuity equation ∂_{z}*J*_{z} = 0 subject to the boundary conditions *φ*(*H*) = *φ*_{+} and *φ*(*h*) = 0. This gives the solution

for *α* ≥ 0 and *β* ≥ −1, or a slightly more complex function for other values. The corresponding solution for the flux, equivalent to the particle capture rate, is

When comparing this prediction with the simulated flux, care should be taken to account for the normalisation (Eq. (69)), because the number of particles *N*_{p} is coupled to concentration *φ*_{+} at the top boundary condition *φ*(*H*) = *φ*_{+}. This relationship can also be computed exactly,

which depends on *H*, so the power-law of *J*_{z}(*H*) should be rescaled based on whether *N*_{p} or *φ*_{+} is kept constant. In the limit *h* ≪ *H*, this simplifies to

Hence, using Eq. (78) for *h* ≪ *H* yields

For slowly moving parallel Stokeslets (*α* = 2 and *β* = 0), we then have *J*_{z} ∝ *z*^{−3} for a constant *φ*_{+} concentration, or *J*_{z} ∝ *z*^{−4} for a constant number of particles *N*_{p}. In Fig. 4c we show the latter.

### Diffusion from a source to an active carpet sink: comparison with thermal diffusion

We expect that the thermal diffusion will be more effective at transporting the particles if the distance between the source and the sink is large, because the thermal noise does not decay with *z*. To quantify this, we equate the active carpet flux (Eq. (81)) with the thermal flux,

The boundary concentration *φ*_{+} cancels out, so we find that the value of *H* for which the two are equal is

Inserting \(\tilde{D}=6\pi n{h}^{2}{f}_{\parallel }^{2}/{D}_{\text{r}}\) from Eq. (5) with *α* = 2 and *β* = 0 for parallel Stokeslets, and using the typical values *h* ~ 1 μm, *n* ~ 1 μm^{−2}, *D*_{r} ~ 1 s^{−1} and 8*π**μ**f*_{∥} ~ 1 pN, we find *H*^{*} ~ 350 and 11 μm, respectively, for micron-sized and molecular paxrticles of *D*_{th} ~ 0.5 and 500*μ*m^{2}/s.

The diffusive flux can also be computed in the presence of both thermal and active fluctuations. As before, we solve ∂_{z}*J*_{z} = 0 using the generalised Fick’s law that includes thermal diffusion (Eq. (75)) without gravity, with boundary conditions *φ*(0) = 0 and *φ*(*H*) = *φ*_{+}. For parallel Stokeslets this yields the concentration profile

and the corresponding flux

As expected, in the limit \(\tilde{D}\to 0\) we recover the thermal flux, *J*_{z} = −*D*_{th}*φ*_{+}/*H*. This corresponds to ~50 particles/second for molecular diffusion with *D*_{th} ~ 500 μm^{2}/s, and using *φ*_{+} ~1 particle/μm and *H* ~ 10 μm. Conversely, in the limit *D*_{th} → 0 we recover the original solution, \({J}_{z}=-2\tilde{D}{\varphi }_{+}/{H}^{3}\). This gives a ‘bare’ active flux of ~60 particles/second when inserting the same values as those below Eq. (84). Interestingly, these fluxes do not just add up because there is also a cross term. In fact, the total flux from Eq. (86) gives *J*_{z} ~ 128 particles/second. Therefore, the thermal diffusion can actually enhance the active diffusive flux, and vice versa, since they co-operate.

### Advective and diffusive transport

To investigate the relative importance of local advective and diffusive transport, we consider a carpet composed of perpendicular Stokeslets that fluctuate about a non-zero mean. Then, the Ornstein–Uhlenbeck process (Eq. (39)) becomes

where the mean force is \(\langle f\rangle =\bar{f}\) and its variance is Var(*f*) = *σ*^{2}*τ*/2 as before. The resulting flow is then described by an advective contribution, \({\boldsymbol{v}}_{\text{adv}}\) due to the mean force, \(\bar{f}\), and a diffusive contribution, \({\boldsymbol{v}}_{\text{diff}}\) due to its variance, Var(*f*). The mean of the diffusive contribution vanishes when averaging over the temporal noise but, at any one location, the advection does not.

Naturally, the advection is not significant in the limit of a small mean force, when \({\bar{f}}^{2}\ll \,\text{Var}\,(f)\). Even when the mean force is comparatively large, however, the active diffusion can still dominate far from the surface, depending on the structure of \(F({\boldsymbol{r}}_{a},{\boldsymbol{p}}_{a})\). This is explained in terms of the local heterogeneities becoming less important when *z* ≫ *r*_{nn}, where the typical nearest-neighbour distance between actuators is \({r}_{\text{nn}} \sim 1/\sqrt{n}\). To quantify this carefully, we consider a square lattice of perpendicular Stokeslets with lattice spacing *r*_{nn} = *λ*. That is, the forces are located at position (*i**λ*, *j**λ*, *h*) where *i* and *j* are integer numbers, so the number density *n* = 1/*λ*^{2}. The total advection generated by this active carpet is given by

where \({\boldsymbol{u}}_{\perp }\) is given by Eq. (15). This total flow is shown in Fig. 5 for different lattice spacings, where all Stokeslets have the same (negative) force \(\bar{f}\). In all cases, there is a down-welling region (downward flow) near the Stokeslets and, by incompressibility, up-welling regions between the Stokeslets. Perhaps counter-intuitively, at a given distance *z* from the surface, the sparse carpets (Fig. 5a, with large *λ*) drive stronger flows than the dense carpets (Fig. 5b, with small *λ*). This is highlighted in Fig. 5c, which plots the vertical flow velocity along the line *y* = 0 for different values of *λ*. These curves show the down-welling regions around *x* = 0, ±*λ*, and up-welling regions around *x* = ±*λ*/2, ±3*λ*/2, … , but their amplitude decreases strongly with decreasing *λ*, i.e. with increasing number density *n*. This is quantified further in Fig. 5d, which shows the vertical flow directly above a Stokeslet (*x* = *y* = 0). Using Eq. (22), we write the normalised total vertical flow as

where the dimensionless number \(\zeta =z/\lambda =z\sqrt{n}\) and the normalisation factor is \({u}_{\perp }^{z}(0,0,z)=12{h}^{2}\bar{f}/{z}^{3}\). Recall that \(\bar{f}\) has units m^{2}/s because forces are scaled with the fluid viscosity (see text under Eq. (15)). Then, in the limit *ζ* → 0 we recover the flow due to a single Stokeslet, Φ → 1, as expected. However, in the limit *ζ* → *∞* the flow tends to zero because the spatial gradients in the actuator density disappear. This decay is quite strong (Fig. 5d; black points), approximately like a Gaussian function, \({{\Phi }}(\zeta )\approx \exp (-{\zeta }^{2})\) (dashed blue line). Thus, the normalised advective transport decays rapidly with *ζ*, while the diffusive transport actually increases. Specifically, using \(\langle {v}_{\,\text{diff}\,,z}^{2}\rangle =15\pi n{h}^{4}\,\text{Var}\,(f)/{z}^{4}\), we can write the normalised diffusive transport as

which is shown in Fig. 5d as red lines. The relative importance of the diffusive and the advective transport is

Hence, the diffusion dominates over advection beyond a distance \({z}^{* }={\zeta }^{* }/\sqrt{n}\) from the carpet, where

in terms of the Lambert *W*_{0} function. This occurs at *ζ*^{*} ≈ 0.85 for \({\text{Var}}\,(f)={\bar{f}}^{2}\), which is fairly close to the active carpet. Even when their fluctuations are a thousand times weaker (Fig. 5d; red dotted line), for \({\bar{f}}^{2}=1{0}^{6}\,\text{Var}\,(f)\), the transition occurs at *ζ*^{*} ≈ 2.55, which is still not that far from the carpet. This is especially relevant for high actuator densities. Indeed, many organisms like *Vorticella* colonies can grow fairly dense, approaching close packing.

Another point to note is that the advective flows can average out in time: Consider a particle located in a down-welling region, slowly moving down towards an actuator. If the particle is also subject to active and/or passive fluctuations, it can diffuse horizontally into an up-welling region, so it can escape. To demonstrate this, we repeat our diffusion simulations (cf. Fig. 2a–d) for a carpet of perpendicular Stokeslets with very weak active fluctuations, Var(*f*) = 10^{−6}, compared to a strong mean force directed towards the carpet, \(\bar{f}=-1\). As shown in Fig. 5e, the MSD still transitions to diffusive motion when *t* > *τ*, and the ballistic advective motion disappears over time. Moreover, the resulting space-dependent diffusivity (Fig. 5f) still agrees with the theoretical prediction (Eq. (4)) for *z* ≳ 3 for all components of *D*_{ij}.

The nutrient flux that individual organisms receive therefore depends strongly on the distance of the nutrient source. If the source is located at *H* < *z*^{*}, then the flux to individual organisms can be large due to advection. If the source is located at *H* > *z*^{*}, then the flux is determined by diffusion. This spreads out the nutrients horizontally before they reach the surface, so on average all organisms receive the same global diffusive flux as discussed in ‘Methods: Diffusion from a source to an active carpet sink: simulation details’.

The theory can also be extended for situations where the relation \(\langle {v}_{\,\text{diff}\,}^{2}\rangle \gg {v}_{\,\text{adv}\,}^{2}\) does not hold. This may be important for scenarios in biology or synthetic carpets of intermediate actuator densities. As a first approximation, one could explicitly insert the advection flow as \({\boldsymbol{v}}_{\text{d}}={\boldsymbol{v}}_{\text{adv}}(\boldsymbol{r})\) into the three-dimensional generalised flux (Eq. (67)), assuming that the fluctuations are still uniformly distributed on the surface. This advection term could be written in terms of Stokeslets, or found with any other hydrodynamic technique such as the boundary-element method, a squirmer-like model, or computational fluid dynamics (CFD) simulations. The disadvantage of this formulation is that it is inherently system specific. The advantage is that any flow pattern of interest can be inserted (e.g. ciliary transport, filter feeding, bacteria on surfaces), so the resulting advection-diffusion equations can be solved accordingly.

### Quenched disorder

The connection between advective and diffusive transport is also related to quenched disorder, the notion that spatial heterogeneity can be frozen in place so a spatial average would not be equal to a local temporal average. In other words, a system features quenched disorder if it has random variables that are quenched (frozen) in time, so their dynamics cannot evolve as fast as the other time scales in the system. In general, active carpets can indeed feature quenched disorder. In that case, it would not be appropriate to model the tracer dynamics with the generalised Fick’s laws described here: This modelling approach is based on averaging with respect to a large ensemble of active carpet configurations, so it would not always be informative about the dynamics of a single specific system configuration.

However, the disorder need not necessarily be quenched for active carpets. The relevant random variables are the relative positions and orientations between the actuators and the tracer particles. Therefore, there is no quenched disorder when the actuators meander along the surface, like bacteria, as long as they move or turn rapidly. Similarly, even if the actuators themselves are fixed, the tracers may still diffuse in space under the right conditions, so the relative positions could still vary freely. They key question is what these conditions are.

The first requirement is that the advective transport (due to individual actuators locally) is much weaker than the diffusive transport (due all actuators together, possibly aided by thermal noise). If a tracer is caught in a local actuator current, then its dynamics are effectively quenched; however, if the particle can diffuse away and escape, the ballistic motion disappears and the advection flows tend to average out over time (Fig. 5e, f). Hence, as the particles explore space horizontally, their spreading over time becomes equivalent to spatial averaging, so the dynamics become annealed. As described in previous section, the diffusion dominates advection if the particles are located far away from the active surface (see Eq. (91)).

The second requirement is that the time scale of diffusive transport is slow compared to the time scale of the active fluctuations themselves. This ensures that the tracer motion is diffusive over time and not ballistic according to a specific carpet configuration. As described in ‘Methods: Derivation of the mean-squared displacement and space-dependent diffusivity’, this requirement is satisfied when Eq. (54) holds. Therefore, both requirements are fulfilled beyond a certain distance from the surface.

To verify the generalised Fick’s laws, we compared their predictions with detailed hydrodynamic simulations. These fully resolve all the actuator positions and orientations throughout time, so any quenched disorder is explicitly included. Importantly, all our main findings (enhanced diffusivities, sedimentation profiles, nutrient fluxes) are supported by this data. Indeed, we found that the simulations and the theory agree well with one another beyond a certain distance from the surface, when both requirements described above are fulfilled. Future theories could perhaps relax these conditions by taking the effects of quenched disorder into account. We expect this could be an exciting opportunity of further research in the field of non-equilibrium statistical mechanics and active matter systems.

### Extension to more complex geometries

In principle, our theory may be generalised for active carpets of more complex geometries by taking the following steps: First, one should find the hydrodynamic Green’s function (cf. the Blake tensor in Eq. (12)) that satisfies the Stokes equations and the boundary conditions of the geometry in question. Once this flow solution is known, one can start developing simulations to verify the following steps. Second, the mean flow \(\langle \boldsymbol{v}(\boldsymbol{r})\rangle\) should be determined by integrating this Green’s function over the carpet along with its force distribution, as in Eq. (32). This may tend to zero for homogeneously distributed carpets, depending on the surface shape and the distribution of actuator positions and orientations, but not necessarily. Third, one should determine the variance tensor \({{\mathcal{V}}}_{ij}(\boldsymbol{r})=\langle {v}_{i}{v}_{j}\rangle\) as in Eq. (33), which may in general be dependent on all three spatial coordinates. Fourth, the generalised flux may be extended by revisiting the telegraph model, as in ‘Methods: Generalised Fick’s laws’. These equations may then be solved numerically or analytically, but care should be taken that the conditions for the theory to be accurate are correctly translated to the new geometry.

## Data availability

All simulation data used for this paper are available from the corresponding authors upon request.

## Code availability

All simulation codes used for this paper are available from the corresponding authors upon request.

## References

- 1.
Chandrasekhar, S. Stochastic problems in physics and astronomy.

*Rev. Mod. Phys.***15**, 1–89 (1943). - 2.
Lauga, E. & Powers, T. R. The hydrodynamics of swimming microorganisms.

*Rep. Progr. Phys.***72**, 096601 (2009). - 3.
Ramaswamy, S. The mechanics and statistics of active matter.

*Annu. Rev. Condens. Matter Phys.***1**, 323–345 (2010). - 4.
Koch, D. L. & Subramanian, G. Collective hydrodynamics of swimming microorganisms: living fluids.

*Annu. Rev. Fluid Mech.***43**, 637–659 (2011). - 5.
Saintillan, D. & Shelley, M. J. Active suspensions and their nonlinear models.

*Compt. Rend. Phys.***14**, 497–517 (2013). - 6.
Marchetti, M. C. et al. Hydrodynamics of soft active matter.

*Rev. Mod. Phys.***85**, 1143–1189 (2013). - 7.
Elgeti, J., Winkler, R. G. & Gompper, G. Physics of microswimmers—single particle motion and collective behavior: a review.

*Rep. Progr. Phys.***78**, 056601 (2015). - 8.
Cates, M. & Tailleur, J. Motility-induced phase separation.

*Annu. Rev. Condens. Matter Phys.***6**, 219 (2015). - 9.
Bechinger, C. et al. Active particles in complex and crowded environments.

*Rev. Mod. Phys.***88**, 045006 (2016). - 10.
Zöttl, A. & Stark, H. Emergent behavior in active colloids.

*J. Phys. Condens. Matter***28**, 253001 (2016). - 11.
Needleman, D. & Dogic, Z. Active matter at the interface between materials science and cell biology.

*Nat. Rev. Mater.***2**, 17048 (2017). - 12.
Kanazawa, K., Sano, T. G., Cairoli, A. & Baule, A. Loopy Lévy flights enhance tracer diffusion in active suspensions.

*Nature***579**, 364–367 (2020). - 13.
Riedel, C. et al. The heat released during catalytic turnover enhances the diffusion of an enzyme.

*Nature***517**, 227 (2015). - 14.
Golestanian, R. Enhanced diffusion of enzymes that catalyze exothermic reactions.

*Phys. Rev. Lett.***115**, 108102 (2015). - 15.
Illien, P. et al. Exothermicity is not a necessary condition for enhanced diffusion of enzymes.

*Nano Lett.***17**, 4415–4420 (2017). - 16.
Wu, X.-L. & Libchaber, A. Particle diffusion in a quasi-two-dimensional bacterial bath.

*Phys. Rev. Lett.***84**, 3017 (2000). - 17.
Kim, M. J. & Breuer, K. S. Enhanced diffusion due to motile bacteria.

*Phys. Fluids***16**, L78–L81 (2004). - 18.
Thiffeault, J.-L. & Childress, S. Stirring by swimming bodies.

*Phys. Lett. A***374**, 3487–3490 (2010). - 19.
Pushkin, D. O., Shum, H. & Yeomans, J. M. Fluid transport by individual microswimmers.

*J. Fluid Mech.***726**, 5–25 (2013). - 20.
Jeanneret, R., Pushkin, D. O., Kantsler, V. & Polin, M. Entrainment dominates the interaction of microalgae with micron-sized objects.

*Nat. Commun.***7**, 12518 (2016). - 21.
Peng, Y. et al. Diffusion of ellipsoids in bacterial suspensions.

*Phys. Rev. Lett.***116**, 068303 (2016). - 22.
Vaccari, L., Molaei, M., Leheny, R. L. & Stebe, K. J. Cargo carrying bacteria at interfaces.

*Soft Matter***14**, 5643–5653 (2018). - 23.
Gilpin, W., Prakash, V. N. & Prakash, M. Rapid behavioral transitions produce chaotic mixing by a planktonic microswimmer. Preprint at https://arxiv.org/abs/1804.08773 (2018).

- 24.
Mathijssen, A. J. T. M., Jeanneret, R. & Polin, M. Universal entrainment mechanism controls contact times with motile cells.

*Phys. Rev. Fluids***3**, 033103 (2018). - 25.
Pedley, T., Hill, N. & Kessler, J. The growth of bioconvection patterns in a uniform suspension of gyrotactic micro-organisms.

*J. Fluid Mech.***195**, 223–237 (1988). - 26.
Hill, N. & Pedley, T. J. Bioconvection.

*Fluid Dyn. Res.***37**, 1–20 (2005). - 27.
Karimi, A. & Ardekani, A. Gyrotactic bioconvection at pycnoclines.

*J. Fluid Mech.***733**, 245–267 (2013). - 28.
Tuval, I. et al. Bacterial swimming and oxygen transport near contact lines.

*Proc. Natl Acad. Sci. USA***102**, 2277–2282 (2005). - 29.
Magar, V., Goto, T. & Pedley, T. J. Nutrient uptake by a self-propelled steady squirmer.

*Q. J. Mech. Appl. Math.***56**, 65–91 (2003). - 30.
Short, M. B. et al. Flows driven by flagella of multicellular organisms enhance long-range molecular transport.

*Proc. Natl Acad. Sci. USA***103**, 8315–8319 (2006). - 31.
Michelin, S. & Lauga, E. Optimal feeding is optimal swimming for all Péclet numbers.

*Phys. Fluids***23**, 101901 (2011). - 32.
Tam, D. & Hosoi, A. E. Optimal feeding and swimming gaits of biflagellated organisms.

*Proc. Natl Acad. Sci. USA***108**, 1001–1006 (2011). - 33.
Mathijssen, A. J. T. M., Culver, J., Bhamla, M. S. & Prakash, M. Collective intercellular communication through ultra-fast hydrodynamic trigger waves.

*Nature***571**, 560–565 (2019). - 34.
Shimmen, T. & Yokota, E. Cytoplasmic streaming in plants.

*Curr. Opin. Cell Biol.***16**, 68–72 (2004). - 35.
Juelicher, F., Kruse, K., Prost, J. & Joanny, J.-F. Active behavior of the cytoskeleton.

*Phys. Rep.***449**, 3–28 (2007). - 36.
Mizuno, D., Tardin, C., Schmidt, C. F. & MacKintosh, F. C. Nonequilibrium mechanics of active cytoskeletal networks.

*Science***315**, 370–373 (2007). - 37.
Wilhelm, C. Out-of-equilibrium microrheology inside living cells.

*Phys. Rev. Lett.***101**, 028101 (2008). - 38.
Fodor, É. et al. Activity-driven fluctuations in living cells.

*Europhys. Lett.***110**, 48005 (2015). - 39.
Goldstein, R. E. Fluid dynamics at the scale of the cell.

*J. Fluid Mech.***807**, 1–39 (2016). - 40.
Needleman, D. & Shelley, M. The stormy fluid dynamics of the living cell.

*Physics Today***72**, 32–38 (2019). - 41.
Bruot, N. & Cicuta, P. Realizing the physics of motile cilia synchronization with driven colloids.

*Annu. Rev. Condens. Matter Phys.***7**, 323–348 (2016). - 42.
Vilfan, A. & Jülicher, F. Hydrodynamic flow patterns and synchronization of beating cilia.

*Phys. Rev. Lett.***96**, 058102 (2006). - 43.
Brumley, D. R., Polin, M., Pedley, T. J. & Goldstein, R. E. Hydrodynamic synchronization and metachronal waves on the surface of the colonial alga

*Volvox carteri*.*Phys. Rev. Lett.***109**, 268102 (2012). - 44.
Elgeti, J. & Gompper, G. Emergence of metachronal waves in cilia arrays.

*Proc. Natl Acad. Sci. USA***110**, 4470–4475 (2013). - 45.
Nonaka, S. et al. De novo formation of left–right asymmetry by posterior tilt of nodal cilia.

*PLoS Biol.***3**, e268 (2005). - 46.
Lukens, S., Yang, X. & Fauci, L. Using Lagrangian coherent structures to analyze fluid mixing by cilia.

*Chaos***20**, 017511 (2010). - 47.
Ding, Y., Nawroth, J. C., McFall-Ngai, M. J. & Kanso, E. Mixing and transport by ciliary carpets: a numerical study.

*J. Fluid Mech.***743**, 124–140 (2014). - 48.
Ramirez-San Juan, G. R. et al. Multi-scale spatial heterogeneity enhances particle clearance in airway ciliary arrays.

*Nat. Phys.***16**, 958–964 (2020). - 49.
Berke, A. P., Turner, L., Berg, H. C. & Lauga, E. Hydrodynamic attraction of swimming microorganisms by surfaces.

*Phys. Rev. Lett.***101**, 038102 (2008). - 50.
Drescher, K., Dunkel, J., Cisneros, L. H., Ganguly, S. & Goldstein, R. E. Fluid dynamics and noise in bacterial cell–cell and cell–surface scattering.

*Proc. Natl Acad. Sci. USA***108**, 10940–10945 (2011). - 51.
Mathijssen, A. J. T. M., Pushkin, D. O. & Yeomans, J. M. Tracer trajectories and displacement due to a micro-swimmer near a surface.

*J. Fluid Mech.***773**, 498–519 (2015). - 52.
Mathijssen, A. J. T. M., Guzmán-Lastra, F., Kaiser, A. & Löwen, H. Nutrient transport driven by microbial active carpets.

*Phys. Rev. Lett.***121**, 248101 (2018). - 53.
Xu, H., Dauparas, J., Das, D., Lauga, E. & Wu, Y. Self-organization of swimmers drives long-range fluid transport in bacterial colonies.

*Nat. Commun.***10**, 1792 (2019). - 54.
Pepper, R. E., Roper, M., Ryu, S., Matsudaira, P. & Stone, H. A. Nearby boundaries create eddies near microscopic filter feeders.

*J. R. Soc. Interface***7**, 851–862 (2009). - 55.
Pepper, R. E. et al. A new angle on microscopic suspension feeders near boundaries.

*Biophys. J.***105**, 1796–1804 (2013). - 56.
Kirkegaard, J. B. & Goldstein, R. E. Filter-feeding, near-field flows, and the morphologies of colonial choanoflagellates.

*Phys. Rev. E***94**, 052401 (2016). - 57.
Roper, M., Dayel, M. J., Pepper, R. E. & Koehl, M. A. R. Cooperatively generated stresslet flows supply fresh fluid to multicellular choanoflagellate colonies.

*Phys. Rev. Lett.***110**, 228104 (2013). - 58.
Nielsen, L. T. et al. Hydrodynamics of microbial filter feeding.

*Proc. Natl Acad. Sci. USA***114**, 9373–9378 (2017). - 59.
Shapiro, O. H. et al. Vortical ciliary flows actively enhance mass transport in reef corals.

*Proc. Natl Acad. Sci. USA***111**, 13391–13396 (2014). - 60.
Shapiro, O. H., Kramarsky-Winter, E., Gavish, A. R., Stocker, R. & Vardi, A. A coral-on-a-chip microfluidic platform enabling live-imaging microscopy of reef-building corals.

*Nat. Commun.***7**, 1–10 (2016). - 61.
Durieux, D. M., Gemmell, B. & Du Clos, K. Benthic jellyfish dominate water mixing in mangrove ecosystems. Preprint at https://www.biorxiv.org/content/10.1101/784173v2.full (2019).

- 62.
Morad, M., Khalili, A., Roskosch, A. & Lewandowski, J. Quantification of pumping rate of

*Chironomus plumosus*larvae in natural burrows.*Aquatic Ecol.***44**, 143–153 (2010). - 63.
den Toonder, J. et al. Artificial cilia for active micro-fluidic mixing.

*Lab Chip***8**, 533–541 (2008). - 64.
Van Oosten, C. L., Bastiaansen, C. W. M. & Broer, D. J. Printed artificial cilia from liquid-crystal network actuators modularly driven by light.

*Nat. Mater.***8**, 677 (2009). - 65.
Bricard, A. et al. Emergent vortices in populations of colloidal rollers.

*Nat. Commun.***6**, 7470 (2015). - 66.
Maass, C. C., Krüger, C., Herminghaus, S. & Bahr, C. Swimming droplets.

*Annu. Rev. Condens. Matter Phys.***7**, 171–193 (2016). - 67.
Darnton, N., Turner, L., Breuer, K. & Berg, H. C. Moving fluid with bacterial carpets.

*Biophys. J.***86**, 1863–1870 (2004). - 68.
Jin, X. & Riedel-Kruse, I. H. Biofilm lithography enables high-resolution cell patterning via optogenetic adhesin expression.

*Proc. Natl Acad. Sci. USA***115**, 3698–3703 (2018). - 69.
Schaller, V., Weber, C., Semmrich, C., Frey, E. & Bausch, A. R. Polar patterns of driven filaments.

*Nature***467**, 73 (2010). - 70.
Gong, X., Mathijssen, A. J. T. M., Bryant, Z. & Prakash, M. Engineering reconfigurable flow patterns via surface-driven light-controlled active matter. https://arxiv.org/abs/2004.01368 (2020).

- 71.
Warner, M. & Terentjev, E. M.

*Liquid Crystal Elastomers*(Oxford University Press, 2007). - 72.
Stuart, M. A. C. et al. Emerging applications of stimuli-responsive polymer materials.

*Nat. Mater.***9**, 101 (2010). - 73.
Wang, E., Desai, M. S. & Lee, S.-W. Light-controlled graphene-elastin composite hydrogel actuators.

*Nano Lett.***13**, 2826–2830 (2013). - 74.
Blake, J. R. A note on the image system for a Stokeslet in a no-slip boundary.

*Math. Proc. Camb. Philos. Soc***70**, 303–310 (1971). - 75.
Chilvers, M. A. & O’Callaghan, C. Analysis of ciliary beat pattern and beat frequency using digital high speed imaging: comparison with the photomultiplier and photodiode methods.

*Thorax***55**, 314–317 (2000). - 76.
Zaid, I. M., Dunkel, J. & Yeomans, J. M. Lévy fluctuations and mixing in dilute suspensions of algae and bacteria.

*J. R. Soc. Interface***8**, 1314–1331 (2011). - 77.
Wang, B., Kuo, J., Bae, S. C. & Granick, S. When Brownian diffusion is not Gaussian.

*Nat. Mater.***11**, 481–485 (2012). - 78.
Schnitzer, M. J. Theory of continuum random walks and application to chemotaxis.

*Phys. Rev. E***48**, 2553–2568 (1993). - 79.
Bouchaud, J.-P. & Georges, A. Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications.

*Phys. Rep.***195**, 127–293 (1990). - 80.
Kiørboe, T. & Jackson, G. A. Marine snow, organic solute plumes, and optimal chemosensory behavior of bacteria.

*Limnol. Oceanogr.***46**, 1309–1318 (2001). - 81.
Stocker, R. Marine microbes see a sea of gradients.

*Science***338**, 628–633 (2012). - 82.
Rogers, S. A., Lisicki, M., Cichocki, B., Dhont, J. K. G. & Lang, P. R. Rotational diffusion of spherical colloids close to a wall.

*Phys. Rev. Lett.***109**, 098305 (2012). - 83.
Lisicki, M., Cichocki, B., Rogers, S. A., Dhont, J. K. & Lang, P. R. Translational and rotational near-wall diffusion of spherical colloids studied by evanescent wave scattering.

*Soft Matter***10**, 4312–4323 (2014). - 84.
Lisicki, M., Cichocki, B. & Wajnryb, E. Near-wall diffusion tensor of an axisymmetric colloidal particle.

*J. Chem. Phys.***145**, 034904 (2016). - 85.
Kiørboe, T., Andersen, A., Langlois, V. J., Jakobsen, H. H. & Bohr, T. Mechanisms and feasibility of prey capture in ambush-feeding zooplankton.

*Proc. Natl Acad. Sci. USA***106**, 12394–12399 (2009). - 86.
Miskin, M. Z. et al. Graphene-based bimorphs for micron-sized, autonomous origami machines.

*Proc. Natl Acad. Sci. USA***115**, 466–470 (2018). - 87.
Cui, J. et al. Nanomagnetic encoding of shape-morphing micromachines.

*Nature***575**, 164–168 (2019). - 88.
Ge, Q., Qi, H. J. & Dunn, M. L. Active materials by four-dimension printing.

*Appl. Phys. Lett.***103**, 131901 (2013). - 89.
Banani, S. F., Lee, H. O., Hyman, A. A. & Rosen, M. K. Biomolecular condensates: organizers of cellular biochemistry.

*Nat. Rev. Mol. Cell Biol.***18**, 285–298 (2017). - 90.
Yoshizawa, T., Nozawa, R.-S., Jia, T. Z., Saio, T. & Mori, E. Biological phase separation: cell biology meets biophysics.

*Biophys. Rev*.**12**, 1–21 (2020). - 91.
Deng, J., Molaei, M., Chisholm, N. G. & Stebe, K. J. Motile bacteria at oil-water interfaces:

*Pseudomonas aeruginosa*.*Langmuir***36**, 6888–6902 (2020). - 92.
Ahmadzadegan, A., Wang, S., Vlachos, P. P. & Ardekani, A. M. Hydrodynamic attraction of bacteria to gas and liquid interfaces.

*Phys. Rev. E***100**, 062605 (2019). - 93.
Kostka, J. E. et al. Hydrocarbon-degrading bacteria and the bacterial community response in Gulf of Mexico beach sands impacted by the Deepwater Horizon oil spill.

*Appl. Environ. Microbiol.***77**, 7962–7974 (2011). - 94.
Dubinsky, E. A. et al. Succession of hydrocarbon-degrading bacteria in the aftermath of the Deepwater Horizon oil spill in the Gulf of Mexico.

*Environ. Sci. Technol.***47**, 10860–10867 (2013).

## Acknowledgements

F.G.-L. acknowledges Millennium Nucleus Physics of Active Matter of ANID (Chile). H.L. acknowledges support from the Deutsche Forschungsgemeinschaft, DFG projects SPP 1726 and LO 418/23. A.J.T.M. acknowledges funding from the Human Frontier Science Program (Fellowship LT001670/2017) and the United States Department of Agriculture (USDA-NIFA AFRI grants 2020-67017-30776 and 2020-67015-32330). F.G.-L. and A.M. also acknowledge support from the American Physical Society (APS) for an International Research Travel Award (IRTAP).

## Author information

### Affiliations

### Contributions

F.G.-L. and A.J.T.M. contributed equally to this work and are joint first, last and corresponding authors. F.G.-L., H.L., and A.J.T.M. designed the research and wrote the manuscript. F.G.-L. and A.J.T.M. performed the simulations. A.J.T.M. developed the theory.

### Corresponding authors

## Ethics declarations

### Competing interests

The authors declare no competing interests.

## Additional information

**Peer review information** *Nature Communications* thanks Marco Polin and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Guzmán-Lastra, F., Löwen, H. & Mathijssen, A.J.T.M. Active carpets drive non-equilibrium diffusion and enhanced molecular fluxes.
*Nat Commun* **12, **1906 (2021). https://doi.org/10.1038/s41467-021-22029-y

Received:

Accepted:

Published:

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.