Boundary design regulates the diffusion of active matter in heterogeneous environments

Abstract

Physical boundaries play a key role in governing the overall transport properties of nearby self-propelled particles. In this work, we develop dispersion theories and conduct Brownian dynamics simulations to predict the coupling between surface accumulation and effective diffusivity of active particles in boundary-rich media. We focus on three models that are well-understood for passive systems: particle transport in (i) an array of fixed volume-excluding obstacles; (ii) a pore with spatially heterogeneous width; and (iii) a tortuous path with kinks and corners. While the impact of these entropic barriers on passive particle transport is well established, we find that these classical models of porous media flows break down due to the unique interplay between activity and the microstructure of the internal geometry. We study the activity-induced slowdown of effective diffusivity by formulating a Smoluchowski description of long-time self diffusivity which contains contributions from the density and fluctuation fields of the active particles. Particle-based and finite element simulations corroborate this perspective and reveal important nonequilibrium considerations of active transport.

---

1 Introduction

Motile microorganisms self-propel and enhance their transport through complex, boundary-rich environments. The diffusion of living systems is crucial to public health; chronic infections can occur when pathogenic bacteria invade wound sites,^1–3 whereas bioremediation requires bacterial penetration into soil sediments.^4,5 In agriculture, rhizosphere bacteria move through soils and protect plant roots.⁶ Unfortunately, many biologically or industrially relevant porous materials are opaque and the details of the microstructure are challenging to resolve. Only recently have experiments been able to directly observe and track individual cell motion through a 3D porous domain;⁷ therefore, the development of theoretical models to predict boundary effects on self-propelled constituents is critical to our understanding of self-propelled (active) transport in complex physical environments.

Porous media guides species transport through a combination of enthalpic and entropic interactions.^8,9 In this work, we will focus our attention on entropic, excluded-volume interactions between diffusing particles and hard external boundaries. Introducing even a dilute concentration of immobile obstacles with purely excluded-volume interactions can provide a substantial slowdown to diffusive flux.^10–12 While self-propelled bacteria move via nonequilibrium forcing, recent research has found success in describing swimmers as Brownian particles under a higher effective temperature and swim energy k_sT_s. This is then translated to a bulk diffusivity using the drag coefficient ζ by means of a modification to the Stokes–Einstein–Sutherland relation: D₀ = (k_BT + k_sT_s)/ζ.

This effective temperature approach is an approximation of the true nonequilibrium behavior, but it has shown success in predicting swim pressure, bulk self-diffusivity, and the energy transferred to a tracer in an active bath.^13–16 For an active Brownian particle (ABP), one can define the swim energy in 2D:¹⁷k_sT_s = ζU₀²τ_R/2. This corresponds to a long time self diffusivity of D₀ = D_T + U₀²τ_R/2, where D_T = k_BT/ζ is the thermal diffusivity, U₀ is the self-propulsive velocity, τ_R is the reorientation time, and D_R = 1/τ_R is the rotational diffusion coefficient.¹⁸ The effective temperature ansatz has been successful in describing the diffusivity of free suspensions; however, it may break down in the presence of excluded volume interactions or external boundaries. Boundaries are especially important for active matter systems due to their activity-induced enrichment along surfaces. Unlike passive (equilibrium) Brownian particles, active particles accumulate at boundaries due to their persistent motion, even in the absence of particle-boundary attraction.^14,19–23

Many researchers have examined how confinement influences the behavior of self-propelled colloidal particles.^7,12,24–29 In arrays of fixed obstacles, specific swimmer models and geometries allow one to observe subdiffusive^30,31 and superdiffusive³² transport. Increased concavity and surface area enhances the partitioning of ABPs at boundaries and perturbs the microstructure further from Boltzmann statistics.^33,34 We hypothesize that the transport properties of active matter can be tuned by the careful control of porous media microstructure, just as the steady-state properties can. We have recently determined that variations in boundary curvature reduces the mobility of nearby swimmers,³⁴ but it is unclear how the design of porous materials reduces the macroscopic transport of swimmers from their well-known bulk diffusivity. Classical theories of diffusive transport in porous media rely on averaging over local degrees of freedom to find the effective diffusivity, D_E. These theories are accurate for a dilute suspension of passive Brownian particles that exhibit no accumulation on surfaces; however, their applicability to self-propelled swimmers is unclear.

In this work, we demonstrate that existing theories of porous media transport fail due to the accumulation of active swimmers on boundaries. Boundaries create significant, long-ranged disturbances to the distribution of active particles; we show that these disturbances correspond to a reduced scaled diffusivity in boundary-rich media.

2 Model

We design three representative examples of distinct porous media to interrogate the accuracy of the corresponding equilibrium-based model with the transport of active Brownian particles (Fig. 1).

Fig. 1 Diagrams of the three sample systems used to probe diffusivity models. (A) An infinite lattice of hard spheres, used to study generalized Taylor dispersion theory. (B) A narrow pore of varying constriction, used to probe Fick–Jacobs approximations. (C) A winding path, used to test the applicability of tortuosity relations.

First, in Fig. 1A, we consider a lattice array of circular 2D inclusions. For a dilute concentration of inclusions, the swimmers collide with the hard boundaries and then escape without difficulty, leading to a minor slowdown that depends on the obstacle density and size. Second, in Fig. 1B, we consider the transport across a narrow constriction. This constriction can be considered as the limit of dense packing for a lattice array of obstacles. Finally, in Fig. 1C, we consider the diffusion along a tortuous path that represents a simple model of a porous network. Each model system is designed to test how the geometry-dependent accumulation of active particles cause deviations in classical theories developed for passive Brownian systems. We focus on the transport of dilute ABPs in rigid, 2D porous media, such as swimmers in confined microfluidic devices, but our qualitative results and framework apply in 3D as well.

An ABP particle i follows the overdamped Langevin equation with a swim force of constant magnitude and white noise Brownian forces and torques.^35–38 To simulate a dilute system, the ABPs interact with obstacles via a purely repulsive potential, but do not interact with each other (i.e., the particles are “ideal”).

Here, F_i,rep is the repulsive force on particle i from the boundary excluded volume, ζ is the drag coefficient, and F_i,swim is the swim force. The swim force is set as F_i,swim = U₀ζq_i, where q_i = [cos θ_i, sin θ_i] is the unit vector describing particle i's orientation in 2D. Finally, the thermal and rotational diffusion coefficients are D_T and D_R = 1/τ_R, with (η_i, ξ_i) as random variables obeying the zero mean and variance consistent with the fluctuation-dissipation theorem. The obstacle-particle repulsive force F_i,rep is implemented using the Weeks–Chandler–Andersen (WCA)³⁹ potential on overlapping spheres defining the obstacle surface (see ESI† for details). We find the effective diffusivity tensor (D_E) of an ABP in the porous medium from the slope of the mean squared displacement (MSD) at long times.As a complementary description of the Langevin equation, the probability distribution P(R, q, t) of a single ABP at position R and orientation q in porous media can be studied using a dimensionless Smoluchowski equation:with the nondimensional variables t̃ = t/τ_R, R̃ = R/L, and . We impose a no-flux boundary condition on all obstacle surfaces. Time is scaled by the reorientation time τ_R, and distance is scaled by a geometric feature of length L. In principle, F_rep(R) should also appear, but we are using excluded volume interactions which can be implemented as boundary conditions. This form of nondimensionalization introduces two important length scales. The run length [small script l] = U₀τ_R is the average distance an ABP self-advects before reorienting and the microscopic length is the root mean squared distance an ABP travels by thermal motion before reorienting. These represent the strength of active and thermal forces, respectively.

Taking orientational moments of the full Smoluchowski equation turns eqn (2) into an infinite series of coupled conservation equations.⁴⁰ The density field depends on the polar order, , which depends on the nematic order, , and so on. We solve eqn (3) numerically and take orientational moments of the full solution to calculate the effective translational diffusivity using methods discussed in the following section. We will examine these equations using the list of relevant length scales shown in Table 1.

Table 1 Table of relevant length scales for active matter in boundary-rich materials

Length	Formula	Description
[small script l]	U 0 τ R	Run length of a free ABP
δ		Thermal diffusive distance moved in one reorientation time
λ ^−1		Screening length of ABP boundary layer
σ	N/A	Diameter of the ABP σ ≪ R, L, H, b
R, L, H, b	N/A	Material length scales

---

3 Case studies

3.1 Diffusion around rigid inclusions

3.1.1 Background. Porous materials can be modelled as a continuous fluid interspersed with discrete rigid inclusions. As a model system, we study a repeating square lattice of disks of radius R and excluded area fraction ϕ = ρπR², where ρ is the disk number density.

For passive Brownian particles in a high porosity lattice (ϕ ≪ 1), as in Fig. 2A, the effective diffusivity decreases linearly with the area density of obstacles D_E/D_T = (1 − ϕ) + [scr O, script letter O](ϕ²),^41,42 and has been solved to arbitrary order using series solutions.¹⁰ If the active particles are equivalent to passive Brownian particles at a higher effective temperature, their reduced diffusivity should track the series solution after replacing D_T with D₀ = D_T + U₀²τ_R/2.

Fig. 2 Diffusivity of swimmers in a lattice array depends strongly on the obstacle density. (A) At high porosity, diffusivity decreases linearly with obstacle area density ϕ. The radius of the disk is defined by the region excluded to the ABP R_c = R_disk + σ/2. In all our results, we keep σ constant and simplify notation by R = R_c. (B) At dense packing, ABPs are trapped in small pockets with narrow escape domains of size ε. The lattice is made of disks centered in a unit cell of length L.

In addition to Brownian dynamics simulations, we numerically solve eqn (3) using generalized Taylor dispersion theory.^{10,12,43–47} We obtain the effective velocity vector u_E and diffusivity tensor D_E of the active particles by splitting the full distribution into local and global contributions. The ABPs position vector R is defined as the sum of the global unit cell vector (X) and the local position inside the unit cell (r). For a particle that enters an L × L unit cell containing a central obstacle, the normalized probability density of finding a particle at position r and orientation q at time t is denoted as g₀(r, q, t). Particle density fluctuations due to the presence of obstacles give rise to an effective diffusivity that is distinct from the bulk diffusivity D₀. The strength and direction of these density fluctuations is measured by the fluctuation field d(r, q, t). Equations are kept in dimensional form for clarity. (See ESI† for full derivation):

where ∇_r indicates the gradient operator over the local position variable r.

We solve these equations at steady state subject to no-flux boundary conditions along the surface of the inclusion, periodic boundary conditions at the unit cell edges, and normalization constraints 〈g₀〉_q,r = 1,〈d〉_q,r = 0 using the finite element method. Where . We use the steady state solutions to compute the average effective velocity u_E and diffusivity D_E of the system:

u_E = U₀〈qg₀〉_q,r − D_T〈∇_rg₀〉_q,r (5a)

D_E = D_TI − U₀〈qd〉_q,r + D_T〈∇_rd〉_q,r. (5b)

3.1.2 Results. In Fig. 3A, we summarize the results of our Brownian dynamics simulations and dispersion theory calculations. In the limit of high porosity (ϕ → 0), the microstructure of the ABPs is weakly perturbed from that of a uniform suspension. The self-diffusivity decreases linearly with obstacle volume fraction ϕ as described by effective temperature theories: D_E = D₀(1 − ϕ). However, as demonstrated in Fig. 3A, the initial slope of the scaled diffusivity varies for highly active particles ([small script l] ≫ δ or [small script l] ≫ R). The effective temperature expression only holds for obstacles spaced far enough apart that the ABPs interact with at most one obstacle over many reorientation times. Therefore, as the activity increases, the effective temperature is restricted to a smaller and smaller ϕ.

Fig. 3 The effective diffusivity of active swimmers depends non-monotonically on the run length ([small script l]) and the microscopic length (δ). (A) Scaled diffusivity for an infinite square lattice as a function of obstacle volume density. The 4 sets represent example leaflets generated by varying particle activity. Points are Brownian dynamics simulations data, dashed lines are numerical solutions to dispersion theory. Dispersion theory calculations are not included for [small script l]/R = 0.25, δ²/R² = 2.5 × 10⁻⁴ because the large relative activity [small script l] ≫δ generates a sharp boundary layer along the obstacle surface that is numerically challenging to resolve. (B) Activity has a nonmonotonic effect on scaled diffusivity due to boundary layer coupling. Scaled diffusivity deviation as a function of run length and for various microscopic length at constant area fraction ϕ = 0.35. Error bars in (A and B) are the standard deviation from 5 independent simulations, and when not visible are smaller than the marker size.

As the obstacle density increases, the diffusivities diverge from those predicted by the effective temperature theory. At intermediate obstacle densities, ϕ ≈ 0.1–0.7, the relative microscopic length (δ) and run length ([small script l]) become important predictors of scaled diffusivity. While each curve monotonically decreases with higher obstacle density, at high activities there can be large deviations from the effective temperature theory.

Finally, as the system approaches close packing, ϕ → π/4, the diffusivities rapidly decrease to zero due to particle transport being limited by a narrow escape through a small pore.

To further examine the effect of run length and microscopic length on active transport at intermediate obstacle densities, ϕ ≈ 0.1–0.7, we define the normalized deviation in measured diffusivity from the model's expected effective temperature diffusivity D_Model as the quantity ΔD/D₀:

This quantity measures the deviation in diffusivity the real ABPs experience compared to that predicted for a Brownian particle in the same system at a higher effective temperature. Therefore, eqn (6) isolates the true coupling between activity-induced boundary accumulation and the obstacle-driven slowdown.

In Fig. 3B, we show the diffusivity deviation for a fixed obstacle density, ϕ = 0.35, while varying both the run length and microscopic length of the ABPs. For weak activity, [small script l]/R ≪ 1, the scaled diffusivity matches the expected results from the passive limit, leading to minimal departure from D_Model. As activity increases beyond [small script l]/R ≈ 1, the effective temperature theories no longer hold. Reductions in transport coefficients are reflected through the increased accumulation within the boundary layer. Strongly active particles are “trapped” within the boundary layer until they reorient or slide off. In contrast, when activity is weak, translational Brownian fluctuations dominate and carry the ABP away from the boundary layer before reorientation occurs. Inactive or weakly active particles “forget” the obstacle surface much faster than their reorientation time or sliding time. As [small script l]/R increases, more ABPs accumulate on the surface. While increasing (δ/R)² may decrease the magnitude of surface accumulation, it also increases the boundary layer thickness, resulting in a net slowdown of scaled ABP diffusivity over the range of parameters studied.

Interestingly, for very small values of the microscopic length [(δ/R)² = 0.01 in Fig. 3B], we see an enhancement in scaled diffusivity beyond the passive theory. We suspect this is due to the unique geometry of the square lattice. When Brownian fluctuations are small in a square lattice, the particles exhibit sustained runs in the interstitial space, effectively lowering the number of obstacles an ABP runs into in a manner similar to the enhancement described by Pattanayak et al.⁴⁸

The boundary-induced accumulation of ABPs occurs over a boundary layer¹⁹ of size λ⁻¹:

In Fig. 4, we show contour plots to connect the ABP microstructure and the diffusivity reduction discussed previously. Fig. 4A and D are the passive case, and exhibit a uniform density field in space. Moving across the top row of Fig. 4A–C, both the boundary layer width and strength are growing, leading to higher relative number of active particles that are influenced by the obstacle before moving to the bulk. The second row of Fig. 4D–F shows the corresponding local diffusivity deviation determined via Taylor dispersion theory that is then averaged over the free space to obtain the normalized effective diffusivity deviation ΔD/D₀. As activity increases across Fig. 4D–F, the local diffusivity decreases in regions near the obstacle surface due to ABPs accumulating within the boundary layer. When an active swimmer is moving unobstructed through the bulk fluid, it spreads with its swim diffusivity D₀. The local diffusivity is quadrupolar near the obstacle, showing reductions at the surface along the axis of motion (in this case the left and right for D_xx), and enhancements above and below due to sliding. The increased surface accumulation in 2D along a weakly curved obstacle with curvature 1/R is described by:³³At a constant density, larger obstacles and activities increase the surface accumulation, which perturbs the behavior of the ABPs away from the effective temperature theory described by k_sT_s.

Fig. 4 Surface accumulation leads to slowdowns in scaled diffusivity. (A–C) Local number density of particles for different activity parameters. As the activity increases in strength, the particles experience larger accumulation near the obstacle surface as shown through the contours. (D–F) Pointwise diffusivity deviation along the x direction calculated numerically via dispersion theory (see ESI† for more details). The pointwise diffusivity deviation is ΔD_xx/D₀ = (D_xx(x,y) − D_Model)/D₀. The pointwise diffusivity has a quadrupolar form due to particles getting trapped in the leading and trailing wakes of the obstacles. (G) Cartoon demonstrating that ABPs within the boundary layer move with a reduced diffusivity compared to the ABPs in the bulk. Subplots (A and D) correspond to [small script l]/R = 0 and δ²/R² = N/A. Subplots (B and E) correspond to [small script l]/R = 1 and δ²/R² = 1. Subplots (C and F) correspond to [small script l]/R = 10 and δ²/R² = 10. All data collected for ϕ = 0.38.

Fig. 4G is a schematic demonstrating the relationship between surface accumulation and transport reduction. For ABPs inside the boundary layer, their self-propulsive force is directed inward against the wall. These ABPs are trapped against the surface until they reorient, and are therefore only diffusing by thermal motion until they escape. This indicates that there are three domains relevant to ABP transport through a packed bed: the occupied area of obstacles ϕ, the boundary layer around obstacles ϕ_BL (which must be weighted by the fraction of ABPs partitioned inside the boundary layer), and then the bulk space ϕ_Free.

Within each of these domains, the ABP moves with a different effective diffusivity. Outside the boundary layer, an ABP in ϕ_Free does not feel any boundary effects and moves with bulk diffusivity D₀. Inside the boundary layer ϕ_BL, an ABP is slowed by collisions with the obstacles, and moves with a local diffusivity approximately equal to D_T. Finally, for impenetrable obstacles no transport can occur and therefore the diffusivity is zero.

3.2 Diffusion through a narrow pore

3.2.1 Background. Taylor dispersion theory is a powerful tool to predict transport properties for dilute systems; however, analytical treatment becomes intractable as the packing density approaches the percolation threshold. Additionally, due to the close packing, the active boundary layers become numerically challenging to resolve. For all of these reasons, we use a simpler Fick–Jacobs^49,50 model of the diffusion reduction due to narrow constriction sites.

The Fick–Jacobs equation describes diffusion in a channel with a slowly varying cross section of width w(y), . Utilizing the methods described by Rubí and Reguera,^51,52 we consider the flux of Brownian particles through a thin pore as a 1D problem, as shown in Fig. 5, with local variations in transverse width as an “entropic force”. For a Brownian particle, the Fick–Jacobs equation is:

See ESI† for the full derivation. The Fick–Jacobs model connects the geometry of the pore with thermodynamic potentials; the constrictions and extensions of the boundaries create entropic traps of the form V_eff = −k_BT ln(w) that reduce the axial flux. Rubí and Reguera^51,52 show that the applicability of the 1D Fick–Jacobs equation can be extended via the introduction of a heuristically defined position dependent diffusion coefficientwith the scaling exponent 1/3 used to empirically match a series solution developed by Zwanzig.⁵⁰

As shown earlier, if the self-propelled particles behave as Brownian particles at a higher effective temperature, one should be able to replace their thermal diffusivity with their bulk diffusivity, D₀ = D_T + U₀τ_R/2. This expression can be used to find the one-dimensional axial effective diffusivity by using the Lifson–Jackson formula⁵³ for motion in a periodic potential described by our entropic barrier

3.2.2 Results. The accuracy of the Fick–Jacobs (FJ) model is predicated on the assumption of rapid local equilibrium in the transverse direction compared to the axial direction. For passive systems, this assumption holds when variations in channel width are slow compared to diffusive motion (|w′(y)| ≪ 1). For an active system, particles move ballistically for their run length [small script l] = U₀τ_R; therefore, we expect the FJ model to be the most accurate when the pore width varies slowly compared to the run length.

Fig. 5 Schematic of the periodic narrow pore used to describe the Fick–Jacobs approximation.

In Fig. 6A, we show the scaled diffusivity of several example active particle systems as a function of the nondimensional constriction. As the pore constricts [(H − b)/H → 1], transport through the small opening is blocked, rapidly reducing the diffusivity down to zero.

Fig. 6 Fick–Jacobs (FJ) theory overestimates the ABP scaled diffusivity. (A) ABP diffusivity decreases below the Brownian reference due to high accumulation along the concave regions of the boundary. (H − b)/H is the degree of constriction due to increased amplitude of pore perturbations (sinusoidal amplitude). The case of H = b corresponds to two flat plates, whereas at b = 0 the pore walls are touching and passage is impossible. The black dashed line is the FJ theory estimate. The data points are from simulations of ABPs. (B) Diffusivity deviation of the ABPs from the Brownian FJ theory. For all simulations, we keep the max wall–wall distance constant at H/Δy = 1 with minimum distance controlled by changing the amplitude of the sinusoid. Lines are to guide the eye. Error bars are the standard deviation from 3 independent simulations, and when not visible are smaller than the marker size.

In Fig. 6B, we plot the scaled diffusivity deviation from the FJ theory as a function of run length and microscopic length. Over the range of parameters studied, FJ theories consistently overestimate the effective diffusivity past [small script l]/Δy ≈ 0.2, with the exact deviation point depending on the pore amplitude and wavelength. Taylor dispersion theory is not used in Fig. 6 due to the high number of finite element meshing points needed to converge the solution at high activities. For example Taylor dispersion theory calculations of weakly active systems, see the ESI† for 3 contour plots of ABP concentration and local diffusivity deviation in a sinusoidal pore.

Our results corroborate the previous findings by Sandoval and Dagdug,⁵⁴ who found that a similar overestimation of diffusivity occurred for a zig–zag and semicircular cavity as the swim speed increased. Their work focused on the impact of varying swim speed U₀. However, our findings reveal this deviation is a result of geometric factors with [small script l] and δ providing the true measure of activity.

Axial diffusivity is reduced primarily through the high accumulation of active particles along the concave regions of the pore. As demonstrated in eqn (8), surface accumulation depends on the radius of curvature. For concave surfaces, the curvature κ is negative, causing surface accumulation to increase rapidly with more negative values of κ: n^surf ≈ 1 + [small script l]²/2δ² − [small script l]²λκ + [scr O, script letter O]([small script l]²λκ)². As the amplitude increases, the radius of curvature shrinks, leading to high accumulation at the concave valleys, away from the small convex pore, thus lowering scaled diffusivity further.

3.3 Diffusion through a tortuous path

3.3.1 Background. Porous materials often contain spatial and/or temporal heterogeneities in their microstructure that preclude detailed characterization. To differentiate between multiple materials with the same porosity, catalyst beds and soil sediments are often described with an empirically measured dimensionless tortuosity, [scr T, script letter T].

Assuming all motion is diffusive, the tortuosity can be determined using a tracer particle of known diffusivity. The relative reduction in measured diffusivity compared to bulk provides an estimate of the reduction factor. Following our definition provided in Fig. 7, the ratio of measured diffusivity in the pore to bulk diffusivity is the ratio in paths traveled in the same amount of time; therefore, one can define a nondimensional tortuosity:

In general, the tortuosity [scr T, script letter T] is system dependent and measured experimentally or via image analysis; however, many models and correlations exist for simple structures.⁵⁵ In order to test the validity of the tortuosity effective temperature approach for active systems, we have designed a winding path that allows us to define the tortuosity directly from the geometry. Our system is made of 5 individual segments of length L_i and width H. Covering a vertical displacement Δy = 3L_i.

Fig. 7 The tortuosity is defined as [scr T, script letter T] = ΔL/Δy, a ratio of the total path length divided by the vertical displacement. For the tortuous paths considered in this work, the center-line path is ΔL = 5L_i and Δy = 3L_i, which gives a tortuosity of [scr T, script letter T] = 5/3.

3.3.2 Results. In Fig. 8A, we show the scaled effective diffusivity measured from active Brownian particle (ABP) simulations. The hatched region between two horizontal lines represents the predictions from the tortuosity theory. The lower bound estimates ΔL as the path through the maze center-line. The upper bound estimates ΔL as the shortest possible path through the maze. While the shortest path changes depending on the geometry, in Fig. 8A we derive the diffusivity from the shortest path of the L_i/H = 4 geometry. Using the L_i/H = 8 geometry would lower the upper bound by 20%; therefore, this choice gives the broadest scaled diffusivity range.

Fig. 8 Tortuosity scalings fail as run length increases beyond the channel width. (A) Scaled effective diffusivity of active particles confined to the tortuous path from Brownian dynamics simulations. The scaled diffusivity decreases as the run length increases above the pore width. The different symbols correspond to paths of different aspect ratios. The hatched region represents the expected values described by tortuosity theory. The top line estimates ΔL as the shortest possible path through the maze, and the bottom line estimates ΔL as the path along the pore center-line. (B) ABPs accumulate significantly at concave corners. The local number density was found via a histogram of particle positions in the Brownian simulation of bin size 0.04σ². That histogram was then divided by the max value. When [small script l]/H > 1, strong corner accumulation decreases the diffusivity far below the tortuosity predictions. This image is from simulation data at L_i/H = 4, [small script l]/H = 20, and δ²/H² = 4. Error bars in (A) are the standard deviation from 3 independent simulations, and when not visible are smaller than the marker size. Lines are to guide the eye.

At low activity, the scaled diffusivity matches the prediction from the center-line path D_E/D₀ ≈ 1/[scr T, script letter T]² = 0.36. This indicates that the diffusion of ABPs roughly follows that of the center-line path, without extreme deviations due to backtracking or strong accumulation.

Our geometry informed tortuosity begins to fail when the run length [small script l] = U₀τ_R approaches the same size as either the segment length L_i or pore width H. Motion is still diffusive over the unit cell, but the transport through each individual segment becomes increasingly ballistic. ABPs move along a segment until getting trapped in a corner. Eventually, the ABP reorients enough to hop across the next segment.

At high activities ([small script l]/H) ≫ 1, the ABPs move ballistically along a segment until colliding with the corner and then spend a significant amount of time trapped before turning. This ballistic motion leads to the concentration profile shown in Fig. 8B. Particles deplete from the convex regions, and accumulate in the highly concave corners. The ballistic motion also increases the effective path length, creating a long-tailed reduction in scaled diffusivity, as shown in Fig. 8A. Further exploration of the parameter space (such as L_i/H ≫ 1 and full sweeps of [small script l]/H at constant (δ/H)²) can be found in the ESI.† Those additional simulations justify our scaling by H instead of L_i, as well as demonstrate that the solutions are relatively insensitive to (δ/H)².

Unlike the narrow pore and the fixed obstacle array, the tortuous path never exhibits a complete blockage of flow, meaning that for networks of the type described in this section, a nonzero diffusivity is possible, even for an arbitrarily large number of segments. Changing the channel boundary patterning from flat walls is expected to further decrease diffusivity. Small pockets of local curvature create an effect similar to that of the sinusoidal pore, providing an additional barrier to whatever maze the ABPs must diffuse through. While this work focuses on a design with only one possible path, the addition of concave “dead ends” should partition the ABPs away from the optimal path, resulting in a larger relative reduction in diffusivity.

4 Discussion

In this work, we test the validity of classical geometric scaling laws for diffusion in heterogeneous environments using analytical theory and Brownian dynamics simulations. We find that even mildly active particles have a scaled diffusivity that deviates from simple effective temperature arguments. The nonequilibrium nature of self-propelled particles needs to be carefully considered in transport calculations, just as it must be for phase behavior and surface accumulation.

Through this work, we have examined the applicability (or lack thereof) of traditional porous media estimates for effective diffusivity in active matter. While Taylor dispersion theory is rigorous and accurate for all activities, numerical instabilities and analytic intractability create difficulty in practical use. Particle-based simulations allow for a micromechanical explanation of the decrease in scaled diffusivity due to accumulation. Both the Fick–Jacobs and tortuosity relationships are accurate in the weakly active limit (provided that one includes the swim diffusivity when calculating the bulk diffusivity), but the strong accumulation at boundary layers and ballistic transport on length scales commensurate with the geometry reduce scaled diffusivity beyond predicted amounts. Interestingly, for a packed bed/lattice, it is possible for ABPs to achieve scaled diffusivity above the passive predictions, and future work aims to understand this intriguing behavior.

Future work on this topic may look at modifications to the effective temperature theory when incorporating it into boundary rich environments. As a first correction, we can utilize the decreased local diffusivity in the boundary layer as D_E ≈ D₀ϕ_Free + αD_Tϕ_BL. With an unknown coefficient α related to the increased partitioning ABPs experience due to boundary curvature.

Previous research by Khatami et al. found that different models of active matter with the same bulk diffusivity have different transport rates when traversing through a maze.⁵⁶ They report that run and tumble models have a smaller mean first passage time than active Brownian particles. Additionally, Kurzthaler et al.⁵⁷ found that the introduction of a reorientation mechanism may lead to greater absolute diffusivity when the reversal run length is commensurate with the maximal chord length of a 3D porous medium. These results indicate that porous media may be used as a novel sorting method for mixtures of active particles utilizing different self-propulsion mechanisms.

This work focuses on the introduction of rigid, immobile obstacles and walls, similar to those found in soil sediments and etched microfluidic devices. However, soft porous materials present a rich opportunity for future study. Boundary fluctuations are important for passive transport in mucus, hydrogels, and other polymeric networks.^8,58 For example, active particles can push through pores smaller than their diameter, if the pore or the particle can deform under thermal or swim forces.

Polymer networks are well-known to reduce the passive (non-motile) particle transport. Steric hindrance, nonspecific interactions (hydrophobicity, electrostatics), and specific interactions (ligand-receptor binding), all play a part to trap foreign debris.^8,9 While this work focused only on the effect of the obstacle excluded volume, other forms of interactions can lead to controllable transport based on the swimmer and obstacle chemistry.

External torques also influence transport in a nontrivial matter. Rod-shaped bacteria can align with the boundaries and with other bacteria to create nematic ordering/flocking.⁵⁹ Chemoattractants, repellents, light sources, and other external fields can influence the favored direction of bacteria in the absence of hard-wall interactions. For example, a random array of point-source repellents can lead to subdiffusive transport in closed spirals.^30,31

Hydrodynamic interactions can provide qualitative or quantitative changes to active transport.^60,61 Hydrodynamic coupling to walls leads to strong interactions between swimmers and surfaces, which would further alter the boundary-accumulation effect on transport.⁶² Pusher type active matter (like bacteria and sperm) are hydrodynamically driven to align parallel to the walls, while puller types align perpendicular to the walls.⁶³ In addition, the surface-modified flow field creates a constant torque leading to motion via counterclockwise (for E. Coli.) spirals.⁶⁴

Finally, in 3D there are additional degrees of freedom a swimmer can use to avoid obstacles. In 2D, close packing of disks precludes transport; however, a close packed bed of monodispersed spheres in 3D is permeable, resulting in a reduced but nonzero diffusivity. Future work in this space would benefit from mean field approximations for the network, in addition to a detailed analysis of the variance of pore sizes.

Author contributions

K. J. M. and S. C. T. conceived of the study; all authors designed research; K. J. M. performed simulations and analytical calculations; S. C. T. supervised the study; and all authors wrote the paper.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-21-1-0287. K. J. M. is supported by the National Science Foundation Graduate Research Fellowship under Grant No. 1650114. S. C. T. is supported by the Packard Fellowship in Science and Engineering. Use was made of computational facilities purchased with funds from the National Science Foundation (OAC-1925717) and administered by the Center for Scientific Computing (CSC). The CSC is supported by the California NanoSystems Institute and the Materials Research Science and Engineering Center (MRSEC; NSF DMR 1720256) at UC Santa Barbara.

References

1. D. Ribet and P. Cossart, Microbes Infect., 2015, 17, 173–183.
2. S. X. Gu and S. R. Lentz, J. Clin. Invest., 2018, 128, 3243–3245.
3. C. Y. Kao, W. H. Lin, C. C. Tseng, A. B. Wu, M. C. Wang and J. J. Wu, Eur. J. Clin. Microbiol. Infect. Dis., 2014, 33, 2157–2162.
4. J. Gannon, U. Mingelgrin, M. Alexander and R. Wagenet, Soil Biol. Biochem., 1991, 23, 1155–1160.
5. J. S. T. Adadevoh, S. Triolo, C. A. Ramsburg and R. M. Ford, Environ. Sci. Technol., 2016, 50, 181–187.
6. O. O. Babalola, Biotechnol. Lett., 2010, 32, 1559–1570.
7. T. Bhattacharjee and S. S. Datta, Nat. Commun., 2019, 10, 2075.
8. L.-H. Cai, S. Panyukov and M. Rubinstein, Macromolecules, 2015, 48, 847–862.
9. J. Witten and K. Ribbeck, Nanoscale, 2017, 9, 8080–8095.
10. M. Mangeat, T. Guérin and D. S. Dean, J. Chem. Phys., 2020, 152, 234109.
11. M. Saxton, Biophys. J., 1994, 66, 394–401.
12. R. Alonso-Matilla, B. Chakrabarti and D. Saintillan, Phys. Rev. Fluids, 2019, 4, 043101.
13. D. Loi, S. Mossa and L. F. Cugliandolo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 051111.
14. S. C. Takatori and J. F. Brady, Curr. Opin. Colloid Interface Sci., 2016, 21, 24–33.
15. G. Szamel, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 012111.
16. E. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco and F. van Wijland, Phys. Rev. Lett., 2016, 117, 038103.
17. S. C. Takatori and J. F. Brady, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 032117.
18. H. C. Berg, Random Walks in Biology, Princeton University Press, 1984.
19. W. Yan and J. F. Brady, J. Fluid Mech., 2015, 785, R1.
20. F. Smallenburg and H. Löwen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 032304.
21. S. C. Takatori, W. Yan and J. F. Brady, Phys. Rev. Lett., 2014, 113, 028103.
22. G. Volpe, S. Gigan and G. Volpe, Am. J. Phys., 2014, 82, 659–664.
23. C. M. Kjeldbjerg and J. F. Brady, Soft Matter, 2022, 18, 2757–2766.
24. L. J. Perez, T. Bhattacharjee, S. S. Datta, R. Parashar and N. L. Sund, Phys. Rev. E, 2021, 103, 012611.
25. N. A. Licata, B. Mohari, C. Fuqua and S. Setayeshgar, Biophys. J., 2016, 110, 247–257.
26. V. K. Truong, D. E. Mainwaring, P. Murugaraj, D. H. K. Nguyen and E. P. Ivanova, J. Mater. Chem. B, 2015, 3, 8704–8710.
27. M. Brun-Cosme-Bruny, E. Bertin, B. Coasne, P. Peyla and S. Rafaï, J. Chem. Phys., 2019, 150, 104901.
28. M. Kumar, J. S. Guasto and A. M. Ardekani, J. Rheol., 2022, 66, 375–397.
29. A. Martínez-Calvo, C. Trenado-Yuste and S. S. Datta, Active transport in complex environments, 2021, arxiv:2108.07011.
30. O. Chepizhko and F. Peruani, Phys. Rev. Lett., 2013, 111, 160604.
31. A. Morin, D. Lopes Cardozo, V. Chikkadi and D. Bartolo, Phys. Rev. E, 2017, 96, 042611.
32. M. Zeitz, K. Wolff and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2017, 40, 23.
33. W. Yan and J. F. Brady, Soft Matter, 2018, 14, 279–290.
34. K. J. Modica, Y. Xi and S. C. Takatori, Frontiers Physics, 2022, 10, 869175.
35. M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gompper and H. Rieger, Nat. Rev. Phys., 2020, 2, 181–199.
36. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
37. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143–1189.
38. A. P. Solon, M. E. Cates and J. Tailleur, Eur. Phys. J.: Spec. Top., 2015, 224, 1231–1262.
39. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
40. D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2008, 100, 178103.
41. J. C. Maxwell, A Treatise on Electricity and Magnetism, Clarendon Press, 1881.
42. G. K. Batchelor, J. Fluid Mech., 1976, 74, 1–29.
43. S. C. Takatori and J. F. Brady, Soft Matter, 2014, 10, 9433–9445.
44. R. N. Zia and J. F. Brady, J. Fluid Mech., 2010, 658, 188–210.
45. S. C. Takatori and J. F. Brady, Phys. Rev. Lett., 2017, 118, 018003.
46. Z. Peng and J. F. Brady, Phys. Rev. Fluids, 2020, 5, 073102.
47. E. W. Burkholder and J. F. Brady, Phys. Rev. E, 2017, 95, 052605.
48. S. Pattanayak, R. Das, M. Kumar and S. Mishra, Eur. Phys. J. E: Soft Matter Biol. Phys., 2019, 42, 62.
49. M. H. Jacobs, Diffusion Processes, Springer Berlin Heidelberg, Berlin, Heidelberg, 1967.
50. R. Zwanzig, J. Phys. Chem. C, 1992, 96, 3926–3930.
51. D. Reguera and J. M. Rubí, Phy. Rev. E, 2001, 64, 061106.
52. D. Reguera, G. Schmid, P. S. Burada, J. M. Rubí, P. Reimann and P. Hänggi, Phys. Rev. Lett., 2006, 96, 130603.
53. S. Lifson and J. L. Jackson, J. Chem. Phys., 1962, 36, 2410–2414.
54. M. Sandoval and L. Dagdug, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 062711.
55. L. Shen and Z. Chen, Chem. Eng. Sci., 2007, 62, 3748–3755.
56. M. Khatami, K. Wolff, O. Pohl, M. R. Ejtehadi and H. Stark, Sci. Rep., 2016, 6, 37670.
57. C. Kurzthaler, S. Mandal, T. Bhattacharjee, H. Löwen, S. S. Datta and H. A. Stone, Nat. Commun., 2021, 12, 7088.
58. I. V. Bodrenko, S. Salis, S. Acosta-Gutierrez and M. Ceccarelli, J. Chem. Phys., 2019, 150, 211102.
59. M. Bär, R. Großmann, S. Heidenreich and F. Peruani, Annu. Rev. Condens. Matter Phys., 2020, 11, 441–466.
60. E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601.
61. K. Drescher, J. Dunkel, L. H. Cisneros, S. Ganguly and R. E. Goldstein, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 10940–10945.
62. S. Rode, J. Elgeti and G. Gompper, New J. Phys., 2019, 21, 013016.
63. A. P. Berke, L. Turner, H. C. Berg and E. Lauga, Phys. Rev. Lett., 2008, 101, 038102.
64. R. G. Winkler and G. Gompper, Handbook of Materials Modeling, Springer International Publishing, Cham, 2018, pp. 1–21.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2sm01421a

---