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 $$ {\mathit{r}}_a=\left(x_a,y_a,h\right)$$ and orientation unit vector $$ {\mathit{p}}_a$$. Each actuator a drives an individual flow, $$ \mathit{u}\left(\mathit{r},{\mathit{r}}_a,{\mathit{p}}_a\right)$$, evaluated at position $$ \mathit{r}$$. In general, this flow can be written in terms of the Blake tensor $$ \mathcal{B}\left(\mathit{r},{\mathit{r}}_a\right)$$⁷⁴ and a multipole expansion thereof (see ‘Methods: Individual flow fields’).

Fig. 1

Flows driven by different types of active carpets.

Diagrams of a few examples: a Molecular motors walking along a surface, described by parallel Stokeslets. b Sessile suspension feeders that draw in feeding currents, described by perpendicular Stokeslets. c E. coli bacteria swimming along a surface, described by Stokes dipoles. d Actuators that rotate about the z-axis, such as nodal cilia, described by Stokes rotlets. e–h Flow fields generated by these different actuator types, shown in the xy plane at z₀ = 2h. i–l Probability distribution functions, PDF(v_z), of the total flow velocity due to N_a = 10⁵ actuators, at different heights z₀. Insets show rescaled histograms that highlight skewness and kurtosis, especially for small z₀ values. m–p The variance of these distributions, $$ ⟨v_i^2⟩$$, corresponding to the strength of active fluctuations in different directions i, as a function of distance z₀ from the active carpet. Insets show the skewness for n and o, but this is zero for m and p. Symbols show simulation results and the lines are the theoretical predictions of Eq. (1) and Table 1.

For example, molecular motors walking or cells crawling along a substrate can entrain the surrounding fluid (Fig. 1a)^{34–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 $$ {\mathit{u}}_{\parallel }=f_{\parallel }\mathcal{B}\cdot {\mathit{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, $$ {\mathit{u}}_{\perp }=f_{\perp }\mathcal{B}\cdot {\mathit{e}}_z$$. Similarly, the flow generated by an E. coli bacterium swimming parallel to a wall (Fig. 1c) is described well⁵⁰ by a ‘Stokes dipole’ flow, given by $$ {\mathit{u}}_{\mathrm{D}}=\kappa \left({\mathit{p}}_a\cdot {\mathit{\nabla }}_a\right)\left(\mathcal{B}\cdot {\mathit{p}}_a\right)$$, with dipole moment κ. Finally, a torque-generating actuator can be described by a Stokes rotlet (Fig. 1d), given by $$ {\mathit{u}}_{\mathrm{R}}=\varrho \left(\left({\mathit{e}}_x\cdot {\mathit{\nabla }}_a\right)\left(\mathcal{B}\cdot {\mathit{e}}_y\right)-\left({\mathit{e}}_y\cdot {\mathit{\nabla }}_a\right)\left(\mathcal{B}\cdot {\mathit{e}}_x\right)\right)/2$$, with strength ϱ. This could represent rotors on a surface, such as bacteria with tethered flagella, or a carpet of nodal cilia⁴⁵ that move around in circles in the xy 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 $$ \mathit{u}$$ are established, the total flow $$ \mathit{v}={\sum }_a\mathit{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\left({\mathit{r}}_a,{\mathit{p}}_a\right)$$ of finding an actuator at position $$ {\mathit{r}}_a$$ and orientation $$ {\mathit{p}}_a$$, the average total flow evaluated at position $$ \mathit{r}$$ is $$ ⟨\mathit{v}\left(\mathit{r}\right)⟩=\int \mathit{u}F\mathrm{d}{\mathit{r}}_a\mathrm{d}{\mathit{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⁵². 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 $$ ⟨\mathit{v}⟩=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 xy plane. We then evaluate the total flow $$ \mathit{v}$$ evaluated at particle position $$ {\mathit{r}}_0=\left(0,0,z_0\right)$$. By repeating this for a large ensemble of N_e independent carpet configurations, we evaluate the distribution of the total flow, $$ \mathrm{PDF}\left(\mathit{v}\right)$$, analogous to the Holtsmark distribution¹. 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₀. 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,

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 $$ ⟨v_x^2⟩=⟨v_y^2⟩=\frac{1}{2}⟨v_z^2⟩$$ for parallel Stokeslets.

Table 1

Properties of fluctuations driven by an active carpet.

	γ	Parallel Stokeslets	Perpendicular Stokeslets	Parallel dipoles	Rotlets
$$ \mathit{u}\left(\mathit{r};{\mathit{r}}_a,{\mathit{p}}_a\right)$$		$$ {\mathit{u}}_{\parallel }$$, Eq. (14)	$$ {\mathit{u}}_{\perp }$$, Eq. (15)	$$ {\mathit{u}}_{\mathrm{D}}$$, Eq. (16)	$$ {\mathit{u}}_{\mathrm{R}}$$, Eq. (17)
Mean	1	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)	(0, 0, 0)
Variance	2	$$ \left(1,1,2\right)\frac{3\pi n{h}^2{f}_{\parallel }^2}{z_0^2}$$	$$ \left(1,1,5\right)\frac{3\pi n{h}^4{f}_{\perp }^2}{z_0^4}$$	$$ \left(7,7,30\right)\frac{3\pi n{h}^2{\kappa }^2}{4z_0^4}$$	$$ \left(1,1,0\right)\frac{3\pi n{h}^2{\varrho }^2}{2z_0^4}$$
Skewness	3	(0, 0, 0)	$$ \left(0,0,\frac{483840}{4199}\right)\frac{\pi n{h}^6{f}_{\perp }^3}{z_0^7}$$	$$ \left(0,0,-\frac{552960}{4199}\right)\frac{\pi n{h}^3{\kappa }^3}{z_0^7}$$	(0, 0, 0)
Kurtosis	4	$$ \left(3,3,20\right)\frac{54\pi n{h}^4{f}_{\parallel }^4}{35z_0^6}$$	$$ \left(267,267,12,880\right)\frac{81\pi n{h}^8{f}_{\perp }^4}{1001z_0^{10}}$$	$$ \left(2797,2797,107,680\right)\frac{27\pi n{h}^4{\kappa }^4}{2288z_0^{10}}$$	$$ \left(1,1,0\right)\frac{27\pi n{h}^4{\varrho }^4}{14z_0^{10}}$$

Moreover, when repeating this algebra for Stokeslets oriented perpendicular to the surface, we find that $$ ⟨v_z^2⟩$$ is five times stronger than $$ ⟨v_x^2⟩$$ (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 xy plane, so their resulting variance only has horizontal components. This fact could be exploited to tune the relative magnitude of $$ ⟨v_x^2⟩$$ and $$ ⟨v_z^2⟩$$. 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, $$ \mathrm{PDF}\left(\mathit{v}\right)$$, 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. Otherwise, Lévy flights must be considered^12,16,76,77. Additionally, these higher-order moments also decay with z₀, 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 \left(t\right)$$, 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 $$ \mathit{v}\left(t\right)$$ 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’).

Fig. 2

Diffusion driven by an active carpet.

For a–d the surface is covered with sessile perpendicular Stokeslets with independent random forces. Also see Supplementary Movie 1. For e–h the surface is covered by parallel Stokeslets moving along a surface with a constant velocity V and rotational diffusion D_r. For both cases we simulate the motion of tracer particles subject to these active fluctuations, ensemble-averaged over N_e = 100 independent tracer trajectories, as a function of their initial distance from the surface, z₀. a Diagram and typical time course of a random force f(t) described by an Ornstein–Uhlenbeck process. b Velocity correlation function (VCF) of the tracer velocity over time. Blue lines show different values of z₀, each ensemble-averaged, and the red dashed line is the prediction of Eq. (48). The inset shows that the resulting correlation time τ is independent of z₀. c Mean-squared displacement (MSD) for different values of z₀, each ensemble-averaged, transitioning from ballistic to diffusive motion. The red dashed lines show the theoretical approximation of Eq. (3) for z₀ = 2, 20. d Anisotropic diffusivity in the directions parallel (blue, green) and perpendicular (purple) to the surface. The solid lines show the prediction of Eq. (4). e Diagram and typical time course of the orientation angle ϕ_a(t) described by rotational diffusion. f–h are equivalent to b–d for moving actuators.

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²〉 = 2D_xxt, and similarly in the other directions. Hence, we determine the components of the diffusion tensor as a function of z₀ (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, $$ ⟨\mathrm{\Delta }{\mathit{r}}^2⟩\ll {\mathit{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,

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)², τ ~1 s and 8πμf_⊥ ~500 pN⁵⁵, the active diffusion can be dominant up to distances of $$ z_{\mathrm{th}}=\sqrt[4]{15\pi n{h}^4⟨f^2⟩\tau /D_{\mathrm{th}}}~1$$ mm for small nutrient molecules of D_th ~ 500 μm²/s. That is much larger than the cell size. Moreover, for micron-sized prey of D_th ~0.5 μm²/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₀/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₀ (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

Again, the active diffusion is significant compared to thermal diffusion. Inserting some typical values into Eq. (5), say h ~ 1 μm, n ~ 1 μm⁻², D_r ~ 1 s⁻¹ and 8πμf ~ 1 pN, we have $$ z_{\mathrm{th}}=\sqrt{6\pi n{h}^2{f}_{\parallel }^2/D_{\mathrm{r}}{D}_{\mathrm{th}}}~7.7\mathrm{\mu }\mathrm{m}$$ for molecules of D_th ~ 500 μm²/s. Similarly, for micron-sized particles with D_th ~ 0.5 μm²/s we have z_th ~ 240 μm, so again much larger than h. Since nh² 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 $$ \phi \left(\mathit{r},t\right)$$ and an external drift flow $$ {\mathit{v}}_{\mathrm{d}}$$ to the total flux, $$ \mathit{J}={\mathit{v}}_{\mathrm{d}}\phi -\tilde{D}\mathit{\nabla }\phi $$. Fick’s second law describes how the particle concentration evolves in time, $$ {\partial }_t\phi =-\mathit{\nabla }\cdot \left({\mathit{v}}_{\mathrm{d}}\phi \right)+\tilde{D}{\nabla }^2\phi $$. This follows directly from the first law and the continuity equation, $$ {\partial }_t\phi =-\mathit{\nabla }\cdot \mathit{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⁷⁸. Indeed, as we saw earlier (Eqs. (4) and (5)), the diffusion tensor $$ D_{ij}\left(z\right)={\mathcal{V}}_{ij}\left(z\right)\tau \left(z\right)$$ 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⁷⁸. The particle flux is then described by the first generalised Fick’s law,

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.

Fig. 3

Sedimentation towards an active carpet.

Particles sink with velocity v_g due to gravity, but are pushed up by the active fluctuations generated by the carpet (dashed purple line), which in this case is composed of moving parallel Stokeslets. a Diagram showing the steady-state sedimentation profile. b Typical trajectories z(t) for a range of sedimentation velocities v_g. c Steady-state distributions of particle density for different v_g. Histograms show simulations and dashed lines are the prediction of Eq. (8). d Position of the maximum particle density as a function of v_g. Markers show simulations and the dashed line is Eq. (9).

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}\left(z\right)=\tilde{\mathcal{V}}/z^{\alpha }$$ as in Table 1, and similarly for the memory time, $$ \tau \left(z\right)=\tilde{\tau }/z^{\beta }$$, the vertical component of the generalised flux (Eq. (6)) becomes

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.

Fig. 4

Diffusion from a source to an active sink.

Particles are spawned at a source (orange line) and subsequently diffuse to a sink (purple line) due to active fluctuations generated by the carpet, which is composed in this case of moving parallel Stokeslets. a Typical trajectory z(t) of a particle. Every time it hits the sink, we place it back at the source (events marked with *). b Steady-state distributions of the particle density for different source positions. Histograms show simulations and dashed lines are Eq. (10). c Diffusive flux as a function of distance. d Diffusive flux as a function of active fluctuation strength. For both c and d the markers show simulations and dashed lines show the prediction of Eq. (82).

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 ∂_zJ_z = 0 subject to the boundary conditions φ(H) = φ₊ and φ(h) = 0. This gives the solution

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 -{\phi }_{+}\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⁻¹. That is more than the thermal flux for micron-sized and molecular particles, ∣J_z∣ ~ 0.05 and 50 s⁻¹, respectively. Interestingly, the thermal and active fluctuations can co-operate to give an even larger flux of ~128 s⁻¹ 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\left({\mathit{r}}_a,{\mathit{p}}_a\right)$$ is constant so the average flow vanishes, $$ ⟨\mathit{v}⟩=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 $$ \overline{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, $$ \mathit{v}={\mathit{v}}_{\mathrm{adv}}+{\mathit{v}}_{\mathrm{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.

Fig. 5

Comparison of advective and diffusive transport.

Here the active carpet is a square lattice of perpendicular Stokeslets. a, b Flows (Eq. (88)) produced by a carpet of λ = 10 (sparse) and λ = 2.5 (dense), respectively, in the plane y = 0. Colours shows the flow magnitude and black arrows are streamlines. c The total advection flow, v_adv,z (x, 0, 2.5), normalised by the flow of a single actuator, u_⊥,z (0, 0, 2.5), for different lattice spacings: λ = 10 (blue), λ = 2.5 (green), λ = 1 (red). The flows vanish as λ decreases. d Comparison of advective and active diffusive transport. Black points show the normalised advection flow, Φ(ζ) (see Eq. (89)). The dashed lines are the decaying functions e^−ζ (dashed green) and $$ e^{-{\zeta }^2}$$ (dashed blue), so the advection vanishes for dense carpets and large distances from the surface. However, the normalised diffusive transport (Eq. (90)) increases with ζ (red lines). e Mean square displacement (MSD) of a particle near a carpet with very weak active fluctuations compared to a strong mean force, Var(f) = 10⁻⁶ and $$ \overline{f}=-1$$. The other parameters used are density n = 1, h = 1, τ = 0.1. The red dashed lines show the prediction of Eq. (3). f Corresponding space-dependent diffusivity. The solid lines show the prediction of Eq. (4). Despite the strong advection currents near the carpet, the theory still holds beyond a certain distance from the surface.

Individual flow fields

The flow generated by an individual actuator is written in terms of the Blake tensor⁷⁴, given by

The flow due to a point force f_∥ oriented parallel to the surface is given by

Note, we have scaled the forces as f = f_N/(8πμ), where μ ~ 10⁻³ Pa s is the viscosity of water, so f_N are in Newtons and f have units [m²/s]. The Stokes dipole flow with prefactor κ in [m³/s] is found by taking the derivative along the $$ {\mathit{p}}_a$$ direction,

Finally, the Stokes rotlet describes an actuator that generates a torque about the z-axis. Its flow with prefactor ϱ in [m³/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,

To be explicit, for the parallel Stokeslets evaluated at $$ \mathit{r}=\left(0,0,z_0\right)$$, and using $$ {\mathit{r}}_{\mathit{a}}=\left({\rho }_a\cos {\theta }_a,{\rho }_a\sin {\theta }_a,h\right)$$ and $$ {\mathit{p}}_{\mathit{a}}=\left(\cos {\varphi }_a,\sin {\varphi }_a,0\right)$$, this gives the Cartesian components

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

For Stokes dipoles oriented parallel to the surface, with orientation $$ {\mathit{p}}_{\mathit{a}}=\left(\cos {\varphi }_a,\sin {\varphi }_a,0\right)$$, we have

Characterising fluctuations: simulation details

To characterise the strength of the active fluctuations, we simulate the distribution $$ \mathrm{PDF}\left(\mathit{v}\right)$$ of the total flow evaluated at particle position $$ {\mathit{r}}_0=\left(0,0,z_0\right)$$. The carpet consists of N_a = 4nL² 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 $$ \mathit{v}={\sum }_a\mathit{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 $$ \mathrm{PDF}\left(\mathit{v}\right)$$ and its moments as a function of distance from the carpet, z₀. In non-dimensionalised simulation units, we use the parameters h = 1, n = 0.1, N_e = 10⁴, and L = 500 so N_a = 10⁵. 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, $$ ⟨\mathit{v}⟩=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

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 = 4nL = 10⁵. 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₀). 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

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

In our simulations, we use an average force magnitude of 〈f²〉 = (3/4)² and τ = 1/10 so that σ² = 2〈f²〉/τ, 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

By averaging over the ensemble of N_e trajectories we compute the VCFs, $$ C\left(\mid t^{\prime }-t^{″}\mid \right)=⟨v_i\left(\mathit{r},t^{\prime }\right)v_j\left(\mathit{r},t^{″}\right)⟩$$, shown in Fig. 2b for different values of z₀ ∈ [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 2D_ijt = 〈Δ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}\mathit{r}}{\mathrm{d}t}\approx \mathit{v}\left({\mathit{r}}_0,t\right)$$, plus higher-order terms that can be expanded as a power series in 1/z₀⁵¹. Consequently, the MSD is given by

The relationship between this ensemble-averaged flow correlation and the Ornstein–Uhlenbeck force correlations is 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² ≫ 2D_zzt = 30πnh⁴〈f²〉τt/z⁴ using Eq. (53), we can rearrange this for a temporal condition, t/τ ≪ z⁶/(30πnh⁴〈f²〉τ²). 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

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 $$ \mathit{V}=V{\mathit{p}}_a$$ along their director, $$ {\mathit{p}}_a$$, with a constant speed V. They generate a flow with force $$ \mathit{f}=f_{\parallel }{\mathit{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}{\varphi }_a}{\mathrm{d}t}=\sqrt{2D_{\mathrm{r}}}\eta \left(t\right)$$ where η is white Gaussian noise with 〈η(t)〉 = 0 and 〈η(t₁)η(t₂)〉 = δ(t₁ − t₂). 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₀ ∈ [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₀ ∈ [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φ = −∂_zJ_z, this gives Fick’s second law,

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⁷⁸, 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\left(z\right)=\sqrt{⟨v_z^2⟩}=\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²(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’⁷⁸ 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, $$ ⟨\mathit{v}⟩=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² 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

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 $$ {\mathit{v}}_{\mathrm{d}}=-v_{\mathrm{g}}\widehat{\mathit{z}}=\frac{d^2\mathrm{\Delta }\rho }{18\mu }\mathit{g}$$, in terms of the gravitational acceleration $$ \mathit{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⁻², 1]. We simulate over a long period of time, t ∈ [0, 10⁵], 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,

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

Then the vertical diffusivity is $$ D\left(z\right)=\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_{\mathrm{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

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, $$ {\lim }_{z\to 0}\phi \left(z\right)=0$$, and the function has a maximum at the same location as before, at $$ z_{\max }={\left(\tilde{D}/v_{\mathrm{g}}\right)}^{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 ∂_zJ_z = 0 subject to the boundary conditions φ(H) = φ₊ and φ(h) = 0. This gives the solution

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,

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

For slowly moving parallel Stokeslets (α = 2 and β = 0), we then have J_z ∝ z⁻³ for a constant φ₊ concentration, or J_z ∝ z⁻⁴ 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_{\mathrm{r}}$$ from Eq. (5) with α = 2 and β = 0 for parallel Stokeslets, and using the typical values h ~ 1 μm, n ~ 1 μm⁻², D_r ~ 1 s⁻¹ 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²/s.

The diffusive flux can also be computed in the presence of both thermal and active fluctuations. As before, we solve ∂_zJ_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

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²/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}{\phi }_{+}/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

Naturally, the advection is not significant in the limit of a small mean force, when $$ {\overline{f}}^2\ll \mathrm{Var}\left(f\right)$$. 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\left({\mathit{r}}_a,{\mathit{p}}_a\right)$$. 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_{\mathrm{nn}}~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/λ². The total advection generated by this active carpet is given by

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

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⁻⁶, compared to a strong mean force directed towards the carpet, $$ \overline{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 $$ ⟨v_{\mathrm{diff}}^2⟩\gg v_{\mathrm{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 $$ {\mathit{v}}_{\mathrm{d}}={\mathit{v}}_{\mathrm{adv}}\left(\mathit{r}\right)$$ 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 $$ ⟨\mathit{v}\left(\mathit{r}\right)⟩$$ 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}\left(\mathit{r}\right)=⟨v_i{v}_j⟩$$ 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.