Active carpets drive non-equilibrium diffusion and enhanced molecular fluxes (2024)

Individual flow fields

The flow generated by an individual actuator is written in terms of the Blake tensor⁷⁴, 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. ⁵¹, 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⁻³ 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 \({\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³/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 = 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 \(\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₀. 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, \(\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 Table1.

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 Table1 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

where τ is a relaxation time, σ is a constant that relates to the force strength, and η is Gaussian white noise with 〈η(t)〉 = 0 and 〈η(t₁)η(t₂)〉 = δ(t₁ − t₂). 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²〉 = (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

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₀.

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₀ ∈ [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}}\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₀⁵¹. 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²〉 = 〈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 Table1.

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

and similarly for other actuator types. Inserting the values used for simulations in Fig.2a–d, being n = 0.1, h = 1, 〈f²〉 = (3/4)², τ = 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₁)η(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,

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⁷⁸, 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²(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, \(\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² 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⁻², 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,

where N_p = ∫φdz is the number of tracers in the system, so ∫Φ(z)dz = 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 Table1 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²/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 φ₀ 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 ∂_zJ_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⁻³ 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}_{\text{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

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²/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 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/λ². 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²/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₀ 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⁻⁶, 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 hom*ogeneously 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.