# Asymptotic transport and dispersion of active particles in periodic porous media

###### Abstract

We elucidate the long-time transport of active Brownian particles flowing through a porous lattice using generalized Taylor dispersion theory and Langevin simulations. The asymptotic spreading of a dilute cloud of microswimmers is shown to obey an obstacle-free advection-diffusion equation, which we use to unravel the effects of motility, lattice geometry and fluid flow. Our model suggests that the interplay of particle self-propulsion, entropic trapping by obstacles and shear-induced dispersion can be exploited for controlling the transport of active swimmers through microstructured environments.

###### pacs:

…The transport of active self-propelled particles through complex microstructured media has important implications in microbial ecology as well as human health. Examples include the spreading of contaminants in soils and groundwater aquifers, bacterial filtering, biodegradation and bioremediation processes, and the transport of motile cells inside the body. Engineering areas such as medical diagnostics and biochemical analysis also often hinge on the manipulation and control of active particle motions through complex geometries C. Bechinger et al. (2016).

While the effective transport of passive tracers in porous media flows has been analyzed minutely in the past, the case of active particles remains largely unexplored. In free space, a microswimmer with constant speed and translational and rotational diffusivities and performs a spatial random walk with long-time dispersivity in dimensions Berg (1993). In a porous matrix and without flow, this active dispersion should intuitively be reduced, not only by geometric obstruction but also by the tendency of microswimmers to accumulate at boundaries Ezhilan and Saintillan (2015); Ezhilan et al. (2015); Yan and Brady (2015). Trapping near obstacles is indeed seen in experiments Takagi et al. (2014) and is thought to be enhanced by hydrodynamic interactions Spagnolie et al. (2015). When an external flow is applied, the coupling of particle orientations with the fluid shear further complicates transport. First, diffusion and swimming across streamlines alter long-time spreading by the classic phenomenon of shear-induced Taylor dispersion Taylor (1953); Aris (1956). Secondly, the rotation of particles in the strong shear near boundaries is expected to result in their alignment against the flow E. Altshuler et al. (2013); Creppy et al. (2018), possibly leading to a net reduction in mean transport by an effect similar to upstream swimming in pressure-driven flows Kaya and Koser (2012); Ezhilan and Saintillan (2015). The situation is yet more complex in the presence of asymmetric obstacles, where the gliding of microswimmers along curved boundaries can drive a net migration even in quiescent conditions Galajda et al. (2007); Yariv and Schnitzer (2014); M. S. Davies Wykes et al. (2017); Tong and Shelley (2017).

Some of these trends have been confirmed in experiments A. T. Brown et al. (2016); Sosa-Hernández et al. (2017); Morin et al. (2017); Creppy et al. (2018) as well as various computational models Chepizhko and Peruani (2013); Potiguar et al. (2014); Chamolly et al. (2017); Khalilian and Fazli (2016), yet a general theoretical framework for predicting active transport and dispersion in even simple porous matrices remains lacking. In this Letter, we present a macro-continuum transport model based on generalized Taylor dispersion theory Brenner (1980) that predicts the asymptotic transport velocity and dispersivity of active Brownian particles in periodic porous media. Starting from the Fokker-Planck equation for the configuration of a single particle in a periodic lattice, we derive theoretical expressions for its long-time mean velocity and dispersivity dyadic in an applied flow, where both quantities simply involve the solution of a boundary-value problem on the pore scale. We also propose a macroscale Eulerian interpretation for these coefficients, and use it to unravel the roles of particle motility, fluid flow, and lattice geometry on asymptotic transport properties.

Problem definition.—We analyze the long-time asymptotic transport of a dilute collection of self-propelled active particles dispersed in a viscous solvent and moving through the interstices of a doubly-periodic 2D porous material. The porous medium is modeled as an infinite square lattice comprised of rigid obstacles of characteristic dimension and generated through the discrete translation of a geometric unit cell of linear dimension (Fig. 1). Each cell in the lattice is labeled by two integers , which identify its horizontal and vertical positions with respect to origin , chosen as the centroid of an arbitrary cell. The array is characterized by its porosity , or ratio of the interstitial fluid area of a cell over its total area.

In a dilute system, we can neglect interparticle interactions and need only consider the transport of a single swimmer. Its instantaneous configuration is described by its position with respect to point and by its swimming direction in the plane of motion, where . The global coordinate can be decomposed as , where is the position of the unit cell where the particle is located and is a local coordinate with respect to the center of that cell. The particle is released at at position with orientation . We introduce the following gradient operators: global , local and orientational .

Particle motion results from self-propulsion, from advection and rotation by the fluid flow, and from translational and rotational diffusion. We model its dynamics by two coupled Langevin equations

(1) | ||||

(2) |

where and are Gaussian random vectors with zero mean and variance and , respectively. The velocity field is a 2D Stokes flow driven by a macroscopic pressure gradient applied across the array, with uniform speed and incoming angle upstream of the array; we solve for it numerically by the boundary integral method (see Supplemental Material ^{1}^{1}1Supplemental material can be found at XXXX.). The fluid vorticity also causes rotation of the particle in Eq. (2). Our aim is to describe the long-time statistics of , which we proceed to model within a continuum framework.

Fokker-Planck description.—We introduce the conditional probability density of finding the particle at a position with orientation at time , given that it was initially released at . It satisfies the Fokker-Planck equation Chandrasekhar (1943); Doi and Edwards (1986)

(3) |

where the initial condition is incorporated by the source term on the right-hand side. The spatial and orientational probability fluxes are easily inferred from the Langevin equations (1)–(2) as Note1 ()

(4) | ||||

(5) |

For a given cell (), Eq. (3) can be recast solely in terms of local variables as

(6) |

where . The probability density function is normalized as

(7) |

which implies convergence of the previous sum and the subsequent convergence of the probability moments. A no-flux condition is enforced on the obstacle boundaries, as well as continuity at the edges of the unit cells Note1 (). We henceforth use dimensionless equations where we scale variables using length scale , time scale , and velocity scale , and we also normalize by . In addition to , non-dimensionalization of the governing equations yields three dimensionless groups: , , and Ezhilan and Saintillan (2015).

Macrotransport model.—Following generalized Taylor dispersion theory Brenner (1980), we seek the mean asymptotic transport properties in terms of the th-order polyadic local moments of defined as

(8) |

The mean asymptotic transport velocity can be obtained from the steady-state zeroth-order local moment , which describes the long-time probability of finding the particle at local position with orientation irrespective of which unit cell it is traversing. It satisfies the simplified steady Fokker-Planck equation , with fluxes

(9) | ||||

(10) |

subject to normalization condition , to the no-flux condition on the surface of the pillar, and to jump conditions and on the unit cell edges and . We solve for numerically using a finite-volume method. is then obtained by simple quadrature as Note1 ()

(11) |

where is the outward unit normal along the cell edges.

Determining the long-time dispersivity tensor requires consideration of the first-order moment , whose asymptotic form can be derived as

(12) |

The fluctuation field encapsulates information about dispersion and obeys the steady problem

(13) |

subject to the no-flux condition at the obstacle walls, and to jump conditions and on and . We numerically solve for , and its knowledge provides as a second quadrature

(14) |

where denotes the spatial and orientational average of with weighting factor Note1 ().

On long length and time scales, a dilute cloud of active particles is thus expected to be transported with mean velocity and to spread with dispersivity as provided by Eqs. (11) and (14). This suggests seeking a coarse-grained Eulerian interpretation of and at the so-called Darcy scale, where the fluid and solid obstacles that comprise the porous material become indistinguishable and distances shorter than the characteristic size of the unit cell are irrelevant. To this end, we introduce the discrete conditional probability density Brenner (1980) for a particle to be located in cell () at time . On large length scales, we can assimilate to a continuous variable and define a corresponding macro-scale gradient operator . The Darcy-scale probability density , which is the continuous counterpart of , is then expected to formally satisfy an obstacle-free advection-diffusion equation Note1 ()

(15) |

As shown below, the analytical solution of Eq. (15) is in exquisite agreement with simulation results at long times, thereby solidifying the interpretation of and .

Results and discussion.—Results from Brownian Dynamics (BD) simulations of the spreading of a concentrated point source of active swimmers based on Eqs. (1)–(2) are compared to the analytical solution of Eq. (15) in Fig. 2. On short time scales, the cloud develops a complex shape that is continually distorted by the pillars and exhibits wakes, boundary layers, and other features typical of advective-diffusive transport around obstacles. At later times, the cloud shape becomes increasingly smooth and is very well captured by the theoretical model, which predicts spreading into a translating anisotropic Gaussian. This comparison asserts the strength of the Darcy-scale interpretation for describing asymptotic transport, but also highlights key differences. In particular, we observe tailed distributions at pre-asymptotic times in Fig. 2(b), a phenomenon also recently reported in experiments using E. coli Creppy et al. (2018). The distributions are found to be negatively skewed in the flow direction in strong flows of weakly swimming particles (high and low ), whereas they become positively skewed in the opposite case Note1 (). At longer times, the tails lose their asymmetry and the distributions converge to the asymptotic diffusive solution as shown in Fig. 2(,).

Asymptotic transport characteristics can be attributed to pore-scale features of the fluid flow and of the steady-state zeroth-order moment as illustrated in Fig. 2(). In the absence of flow, the interplay of self-propulsion and translational diffusion near obstacle boundaries creates a net swimmer polarization against the obstacles Ezhilan and Saintillan (2015); Ezhilan et al. (2015); Yan and Brady (2015), which vanishes at the cell edges by symmetry. This polarization is accompanied by a net increase in density , which is azimuthally symmetric near a circular obstacle and enhanced in regions of high curvature for non-circular pillars Yan and Brady (2018). When a flow is applied, the fore-aft symmetry is broken and particles now preferentially accumulate near stagnation points, with a maximum density reached in the rear of circular obstacles. This accumulation, which is observed in experiments G. Miño et al. (2018) and is akin to that occurring past funnel constrictions E. Altshuler et al. (2013), results from the combination of advection by the flow and of shear-rotation and swimming in the fast-flowing high-shear regions above and below the pillars. Advection and rotation also induce a density minimum ahead of the obstacle, which coincides with a stagnation point in the polarization field. As explained below, these features play a central role in determining the effective spreading behavior.

The dependence of dispersion coefficients on external flow strength and swimming activity is shown in Fig. 3(,) for a lattice comprised of circular pillars. In weak flows, and have similar magnitudes, with faster-swimming particles (high ) spreading more efficiently by active dispersion. Quite remarkably, increasing flow strength first causes a decrease in spreading. This decrease is unique to self-propelled particles, and indeed disappears at low values of . It is attributed to a reduction in active dispersion by the flow, which prevents particles from freely sampling all swimming directions. At higher values of , starts increasing again and displays an asymptotic scaling of , with an exponent of at intermediate flow strengths and of in very strong flows, which is the classical scaling for the porous dispersion of passive tracers Edwards et al. (1991). Interestingly, activity is seen to hinder in strong flows in Fig. 3(), as swimming enhances cross-stream migration, which reduces axial spreading within the framework of shear-induced dispersion. Unlike , the transverse dispersivity in Fig. 3() monotonically decays with as stronger flows cause faster rotations of the swimmers in the local shear and thus limit their ability to travel in the -direction.

Dispersion can be further controlled by tuning the incoming flow angle in Fig. 3(,). For arbitrary , the dispersion tensor is no longer diagonal and maximum dispersion occurs along the eigendirection associated with its maximum eigenvalue . The magnitude of is maximum for , which allows for easy gliding of the swimmers between rows of pillars with minimal collisions. Other flow directions result in more severe obstruction by the pillars and thus display much lower values of , which reaches its minimum for at high flow rates. While in strong flows the direction of maximum dispersion is close to that of the imposed flow, such is not the case in weak flows where and have opposite signs, indicating significant spreading transverse to the flow.

The effect of lattice porosity is analyzed in Fig. 4(,) and further underscores the subtle interplay of active vs shear-induced dispersion. In the absence of flow and at low porosity, dispersion is weak as the pillars act as entropic barriers that hinder active transport Laachi et al. (2007). Increasing porosity thus causes a monotonic increase in and as geometric obstruction plays a lesser role. In strong flows (), porosity instead causes a monotonic decrease in , as the obstacles act as regions of shear production, thus enhancing shear-induced dispersion. In weaker flows (), shows a non-monotonic behavior: as porosity is increased, axial dispersion first decreases due to a drop in local shear rate, before increasing again at low obstacle area fractions due to active swimming. Curiously, the limit of differs from the case of swimmers in free-space, as infinitesimal obstacles still induce fluid shear in an imposed flow while also modifying the swimmer distribution in their vicinity.

Pillar aspect ratio and fore-aft asymmetry both also affect transport in non-trivial ways. As shown in Fig. 4(), elliptical obstacles aligned with a weak flow cause higher streamwise dispersion than if aligned perpendicular to it, as the former allow swimmers to glide past while the latter act as entropic barriers that obstruct transport Laachi et al. (2007). These trends reverse at high flow rates, as perpendicular obstacles induce stronger velocity gradients thereby enhancing shear-induced dispersion. Fore-aft asymmetric obstacles also lead to complex trends in Fig. 4(). In the absence of flow, decreases as the shape parameter deviates from (circle) by the aforementioned entropic obstruction mechanism. When a flow is applied, fluid shear causes particles to swim along the walls and against the flow, resulting in an enhancement of particle accumulation around obstacles for and in a reduction for . Consequently, the effects of shear-alignment (in weak flows) and shear-induced dispersion (in strong flows) on transport are magnified for obstacle shapes that promote accumulation (), as seen in Fig. 4(). This asymmetry becomes negligible at very high values of , as activity becomes subdominant and the magnitude of the shear gradients that drive dispersion is independent of the sign of by reversibility of Stokes flow.

Concluding remarks.—We have analyzed the long-time transport properties of active particles in a porous lattice under both quiescent and flow conditions. The predictions of our continuum model, which agree perfectly with Langevin simulations at long times, highlight the complex interplay of particle motility, alignment under shear, cross-stream migration, lattice porosity and pillar geometry on asymptotic dispersion. Our theory provides a simple framework for analyzing microorganismal transport in natural structured environments and for designing engineered porous media in applications involving microswimmers. The fundamental premise of non-interacting active Brownian particles glosses over many details that may become relevant in specific systems. Future work should address the roles of hydrodynamic interactions with obstacles Spagnolie et al. (2015), swimmer-specific scattering dynamics Kantsler et al. (2013), rheological effects due to activity Saintillan (2018), and the potential emergence of spontaneous flows and active turbulence in denser systems Theillard et al. (2017).

###### Acknowledgements.

D.S. gratefully acknowledges funding from NSF Grant DMS-1463965.## References

- C. Bechinger et al. (2016) C. Bechinger et al., Rev. Mod. Phys. 88, 045006 (2016).
- Berg (1993) H. C. Berg, Random Walks in Biology (Princeton University Press, 1993).
- Ezhilan and Saintillan (2015) B. Ezhilan and D. Saintillan, J. Fluid Mech. 777, 482 (2015).
- Ezhilan et al. (2015) B. Ezhilan, R. Alonso-Matilla, and D. Saintillan, J. Fluid Mech. 781, R4 (2015).
- Yan and Brady (2015) W. Yan and J. F. Brady, J. Fluid Mech. 785, R1 (2015).
- Takagi et al. (2014) D. Takagi, J. Palacci, A. B. Braunschweig, M. J. Shelley, and J. Zhang, Soft Matter 10, 1784 (2014).
- Spagnolie et al. (2015) S. E. Spagnolie, G. R. Moreno-Flores, D. Bartolo, and E. Lauga, Soft Matter 11, 3396 (2015).
- Taylor (1953) G. I. Taylor, Proc. R. Soc. Lond. A 219, 186 (1953).
- Aris (1956) R. Aris, Proc. R. Soc. Lond. A 235, 67 (1956).
- E. Altshuler et al. (2013) E. Altshuler et al., Soft Matter 9, 1864 (2013).
- Creppy et al. (2018) A. Creppy, E. Clément, C. Douarche, M. V. D’Angelo, and H. Auradou, arXiv:1802.01879 (2018).
- Kaya and Koser (2012) T. Kaya and H. Koser, Biophys. J. 102, 1514 (2012).
- Galajda et al. (2007) P. Galajda, J. Keymer, P. Chaikin, and R. Austin, J. Bacteriol. 189, 8704 (2007).
- Yariv and Schnitzer (2014) E. Yariv and O. Schnitzer, Phys. Rev. E 90, 032115 (2014).
- M. S. Davies Wykes et al. (2017) M. S. Davies Wykes et al., Soft Matter 13, 4681 (2017).
- Tong and Shelley (2017) J. Tong and M. J. Shelley, to appear in SIAM J. Appl. Math. (2017).
- A. T. Brown et al. (2016) A. T. Brown et al., Soft Matter 12, 131 (2016).
- Sosa-Hernández et al. (2017) J. E. Sosa-Hernández, M. Santillán, and J. Santana-Solano, Phys. Rev. E 95, 032404 (2017).
- Morin et al. (2017) A. Morin, D. L. Cardozo, V. Chikkadi, and D. Bartolo, Phys. Rev. E 96, 042611 (2017).
- Chepizhko and Peruani (2013) O. Chepizhko and F. Peruani, Phys. Rev. Lett. 111, 160604 (2013).
- Potiguar et al. (2014) F. Q. Potiguar, G. Farias, and W. Ferreira, Phys. Rev. E 90, 012307 (2014).
- Chamolly et al. (2017) A. Chamolly, T. Ishikawa, and E. Lauga, New J. Phys. 19, 115001 (2017).
- Khalilian and Fazli (2016) H. Khalilian and H. Fazli, J. Chem. Phys. 145, 164909 (2016).
- Brenner (1980) H. Brenner, Phil. Trans. R. Soc. Lond. A 297, 81 (1980).
- (25) Supplemental material can be found at XXXX.
- Avron et al. (2004) J. Avron, O. Gat, and O. Kenneth, Phys. Rev. Lett. 93, 186001 (2004).
- Chandrasekhar (1943) S. Chandrasekhar, Rev. Mod. Phys. 15, 1 (1943).
- Doi and Edwards (1986) M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Oxford University Press, 1986).
- Yan and Brady (2018) W. Yan and J. F. Brady, Soft Matter 14, 279 (2018).
- G. Miño et al. (2018) G. Miño et al., Adv. Microbiol. 8, 451 (2018).
- Edwards et al. (1991) D. Edwards, M. Shapiro, H. Brenner, and M. Shapira, Trans. Porous Media 6, 337 (1991).
- Laachi et al. (2007) N. Laachi, M. Kenward, E. Yariv, and K. D. Dorfman, Europhys. Lett. 80, 50009 (2007).
- Kantsler et al. (2013) V. Kantsler, J. Dunkel, M. Polin, and R. E. Goldstein, Proc. Natl. Acad. Sci. USA 110, 1187 (2013).
- Saintillan (2018) D. Saintillan, Annu. Rev. Fluid Mech. 50, 563 (2018).
- Theillard et al. (2017) M. Theillard, R. Alonso-Matilla, and D. Saintillan, Soft Matter 13, 363 (2017).