Active matter invasion of a viscous fluid: unstable sheets and a no-flow theorem

TL;DR

This paper studies how active particles invade viscous fluids, revealing instabilities in particle sheets and proving a no-flow theorem that describes isotropic distributions causing no fluid movement.

Contribution

It introduces a new instability analysis for active particle invasion and proves a surprising no-flow theorem for isotropic distributions in viscous fluids.

Findings

Aligned pusher sheets are always unstable.

Puller sheets can be stable or unstable depending on motility.

Isotropic distributions cause no fluid flow at any time.

Abstract

We investigate the dynamics of a dilute suspension of hydrodynamically interacting motile or immotile stress-generating swimmers or particles as they invade a surrounding viscous fluid. Colonies of aligned pusher particles are shown to elongate in the direction of particle orientation and undergo a cascade of transverse concentration instabilities, governed at small times by an equation which also describes the Saffman-Taylor instability in a Hele-Shaw cell, or Rayleigh-Taylor instability in two-dimensional flow through a porous medium. Thin sheets of aligned pusher particles are always unstable, while sheets of aligned puller particles can either be stable (immotile particles), or unstable (motile particles) with a growth rate which is non-monotonic in the force dipole strength. We also prove a surprising "no-flow theorem": a distribution initially isotropic in orientation loses isotropy immediately but in such a way that results in no fluid flow everywhere and for all time.

Figures

Equations

$${\Psi }_t+\nabla \cdot \left(\dot{\mathbf{x}}\Psi \right)+{\nabla }_{\mathbf{p}}\cdot \left(\dot{\mathbf{p}}\Psi \right)=0,$$

$$\dot{\mathbf{x}}=V_0\mathbf{p}+\mathbf{u}(\mathbf{x})-d_t\nabla \left(\ln \Psi \right),$$

$$\dot{\mathbf{p}}=(\mathbf{I}-\mathbf{pp})\cdot (\mathbf{p}\cdot \nabla \mathbf{u})-d_r{\nabla }_{\mathbf{p}}\left(\ln \Psi \right),$$

$$-\nabla q+\mu {\nabla }^2\mathbf{u}+\nabla \cdot {\mathbf{\Sigma }}_a=0,\nabla \cdot \mathbf{u}=0,$$

$$c_t+\nabla \cdot (c\mathbf{u})-d_t{\nabla }^{2}c=-V_0\nabla \cdot ⟨\mathbf{p}⟩,$$

$$h_t+V_0{h}_x=v+V_0\varphi ,{\varphi }_t+V_0{\varphi }_x=v_x,$$

$$v=-\frac{\alpha \kappa }{4L}\mathcal{H}[h_x],$$

$$\mathcal{H}[f](x)=\frac{1}{\pi }\text{vbox}-{\int }_{-\infty }^{\infty }\frac{f(y)}{x-y}dy.$$

$${\lambda }_{\pm }=\frac{\pi }{4L^4}\left(-\alpha \kappa \left|k\right|-8ikL{V}_0\pm {\gamma }_k\right),$$

$$h_t=\frac{-\alpha \kappa }{4L}\mathcal{H}[h_x].$$

$${\Psi }_t^{\ast }+V_0\mathbf{p}\cdot \nabla {\Psi }^{\ast }-d_t{\nabla }_{\mathbf{x}}^2{\Psi }^{\ast }-d_r{\nabla }_{\mathbf{p}}^2{\Psi }^{\ast }=0,$$

$${\hat{\mathbf{u}}}_{\mathbf{k}}[{\Psi }^{\ast }]=(8\pi |\mathbf{k}{|}^2{)}^{-1}(\mathbf{I}-\hat{\mathbf{k}}\hat{\mathbf{k}})\cdot {\hat{\mathbf{\Sigma }}}_a\cdot \mathbf{k},$$

$${\hat{\mathbf{\Sigma }}}_a\cdot \mathbf{k}={\int }_D\mathbf{p}(\mathbf{p}\cdot \mathbf{k}){\hat{\Psi }}_{\mathbf{k}}^{\ast }(\mathbf{p},t)d\mathbf{p},$$

Full text

The last decade has seen an explosion of interest in the collective dynamics of active particles immersed in fluids, from swimming microorganisms to magnetically driven and phoretic colloidal particles Pedley and Kessler (1992); Dombrowski et al. (2004); Cisneros et al. (2007); Underhill et al. (2008); Baskaran and Cristina (2009); Ramaswamy (2010); Koch and Subramanian (2011); Vicsek and Zafeiris (2012); Marchetti et al. (2013); Saintillan and Shelley (2015); Elgeti et al. (2015); Zöttl and Stark (2016); Yeo et al. (2016) to kinesin-driven microtubule assemblies Helmke et al. (2013); Sanchez et al. (2012); Keber et al. (2014); Prost et al. (2015); Foster et al. (2015); Gao et al. (2015); Shelley (2016); Maryshev et al. (2018). A first-principles model of active suspensions is a mean-field kinetic theory which tracks the distribution of particle positions and orientations and which may include hydrodynamic interactions Saintillan and Shelley (2007, 2008); Subramanian and Koch (2009); Saintillan and Shelley (2015); Koch and Subramanian (2011) and short-range physics Mehandia and Nott (2008); Subramanian and Koch (2009). Constituent particles are classified as either “pushers” or “pullers” depending on the sign of the generated stresslet flow, which in turn depends on the geometry of the body and the mechanism of stress-generation Saintillan and Shelley (2008); Lauga and Powers (2009); Spagnolie and Lauga (2012); Ghose and Adhikari (2014); Lauga and Michelin (2016); Nasouri and Elfring (2018). Other models range from Landau-de Gennes “Q tensor” theories to moment-closure theories Ramaswamy et al. (2003); Woodhouse and Goldstein (2012); Forest et al. (2013); Giomi et al. (2013); Thampi et al. (2013); Brotto et al. (2013). Generic features in these systems include long-range coherence, topological defects, and instability Simha and Ramaswamy (2002); Saintillan and Shelley (2008); Hohenegger and Shelley (2010); Bertin et al. (2013); Ezhilan et al. (2013); Brotto et al. (2013); Shi et al. (2014).

Much is known about active suspensions which cover the entire available physical domain. Far less is known about the invasion of a surrounding particle-free environment, though this is of considerable importance in the dynamic self assembly of swarms Shapiro et al. (1997); Copeland and Weibel (2009); Darnton et al. (2010), and in the formation of biofilms, mycelia, and fruiting bodies Claessen et al. (2014). Novel means of bringing bacteria into a confined region using external flows have allowed for a closer look at rapid expansion, including acoustic-trapping Takatori et al. (2016); Gutiérrez-Ramos et al. (2018), UV-light exposure Patteson et al. (2018) and vortical flows Sokolov and Aranson (2016); Sokolov et al. (2018). The effects of confinement by soft boundaries with surface tension has seen theoretical treatment Gao and Li (2017); Gao et al. (2017), and unstable bands of active particles have been studied in a dry system Ngo et al. (2014).

In this Letter we investigate the dynamics of colonies (a coherent collective) of either motile or immotile active particles as they invade a surrounding viscous fluid. Colonies of aligned pushers are shown to elongate in the direction of particle orientation and then undergo a cascade of transverse concentration instabilities. The initial instability in two-dimensions is shown to be governed at small times by an equation which also describes the Saffman-Taylor instability in a Hele-Shaw cell (flow through a small gap between two nearby plates), or the Rayleigh-Taylor instability in two-dimensional flow through a porous medium, respectively. Linear stability analysis offers approximations that match the results of full numerical simulations. We close with a proof and demonstration of a counter-intuitive “no-flow theorem,” that an isotropically oriented distribution with any initial concentration profile results in no fluid flow everywhere and for all time.

Mathematical model: Following Refs. Saintillan and Shelley (2008); Subramanian and Koch (2009), we describe a dilute suspension of $N$ self-propelled rod-like particles in a viscous fluid by the particle distribution function, $\Psi(\mathbf{x},\mathbf{p},t)$, where $\mathbf{x}$ is the particle position in a periodic spatial domain $D$ while $\mathbf{p}$ is the particle orientation vector on the unit ball $S$ ($|\mathbf{p}|=1$). The number of particles is conserved, $N=\int_{D}\int_{S}\Psi(\mathbf{x},\mathbf{p},t)d\mathbf{p}d\mathbf{x}$, resulting in a Smoluchowski equation,

[TABLE]

where $\nabla=\nabla_{\mathbf{x}}$ and $\nabla_{\mathbf{p}}=(\mathbf{I}-\mathbf{p}\mathbf{p})\cdot\partial/\partial\mathbf{p}$. Neglecting collisions Ezhilan et al. (2013), the fluxes $\dot{\mathbf{x}}$ and $\dot{\mathbf{p}}$ are given by

[TABLE]

with $V_{0}$ the swimming speed, $d_{t}$ ($d_{r}$) the translational (rotational) diffusivity, $\mathbf{u}(\mathbf{x},t)$ the fluid velocity, and $\mathbf{pp}$ a dyadic product.

The environment is assumed to be a viscous Newtonian fluid, driven by stresses generated by the suspended particles. The flow field $\mathbf{u}$ satisfies the Stokes equations, consisting of momentum and mass conservation,

[TABLE]

with $q$ the pressure, $\mu$ the dynamic viscosity, and $\bm{\Sigma}_{a}=\sigma\langle\mathbf{pp}\rangle$ the active stress (proportional to the second orientational moment of $\Psi$, see below). The coefficient $\sigma$ is the force dipole (or stresslet) strength, with $\sigma<0$ for pushers and $\sigma>0$ for pullers Saintillan and Shelley (2015), which has been computed for ellipsoidal “Janus” particles Spagnolie and Lauga (2012); Lauga and Michelin (2016) and for more general particle types Pak and Lauga (2015); Reinken et al. (2018); Nasouri and Elfring (2018), and has been obtained experimentally for a few types of swimming cells Drescher et al. (2010, 2011). Orientational moments will be denoted by $\langle h(\mathbf{p})\rangle=\int_{S}h(\mathbf{p})\Psi\,d\mathbf{p}$. For example, integrating Eq. (1) gives an evolution equation for the particle concentration, $c=\langle 1\rangle$, of the form

[TABLE]

where $\langle\mathbf{p}\rangle$ is the polarity Saintillan and Shelley (2015).

With $\ell$ the active particle length, we scale velocities by the swimming speed, $V_{0}$, and lengths by the mean free path $\ell_{c}=(V/V_{p})\ell$, where $V$ is the total fluid volume and $V_{p}=N\ell^{3}$ is an effective volume of particles, hence $\ell_{c}=V(N\ell^{2})^{-1}$. Time is scaled by $\ell_{c}/V_{0}$, force densities are scaled by $\mu V_{0}/\ell_{c}^{2}$, and $\Psi$ is normalized by the particle number density, $N/V=(\ell^{2}\ell_{c})^{-1}=(\ell/\ell_{c})\ell^{-3}$. The dimensionless dipole strength is defined as $\alpha=\sigma/[\mu V_{0}\ell_{c}^{2}]$. With all variables now taken to be dimensionless, particle conservation is written as $\tilde{\kappa}^{-1}\int_{D}\int_{S}\Psi\,d\mathbf{p}d\mathbf{x}=1$, where $\tilde{\kappa}^{-1}=\ell_{c}^{3}/V$ is proportional to the particle volume fraction. In the case of immotile particles a different velocity scale is required [] (2018).

The far-field velocity due to an individual swimmer at the origin, oriented in the direction $\mathbf{p}$, is $\mathbf{u}=\alpha(8\pi)^{-1}\mathbf{p}\mathbf{p}:\nabla\mathbf{G}(\mathbf{x})$ where $G_{ij}=\delta_{ij}/|\mathbf{x}|+x_{i}x_{j}/|\mathbf{x}|^{3}$ is the Stokeslet singularity Pozrikidis (1992). The active force density is then given by $\mathbf{f}_{a}=\alpha(\ell_{c}/\ell)^{2}\nabla\cdot\langle\mathbf{pp}\rangle$. Following Ref. Saintillan and Shelley (2008) for the sake of comparison, we set $\ell_{c}/\ell=1$. The swimming speed, $V_{0}$, now taken as dimensionless, is unity for motile systems and zero otherwise.

We will consider the case of confinement to motion in two-dimensions in a periodic domain $(x,y)\in[0,L)\times[0,H)$, and invariance in the $\hat{\mathbf{z}}$ direction, writing $\mathbf{p}=(\cos\theta,\sin\theta,0)$ and $\Psi(\mathbf{x},\mathbf{p},t)=\Psi(x,y,\theta,t)$. It is expedient to then define $\kappa=HL$ so that $\kappa^{-1}\int_{0}^{L}\int_{0}^{H}\int_{0}^{2\pi}\Psi\,d\theta\,dy\,dx=1$. Numerical solution of Eqs. (1)-(4) using $L=H=50$ and a pseudospectral method with $256^{3}$ gridpoints and dealiasing (3/2 rule) Fornberg (1998); an integrating factor method along with a second-order accurate Adams-Bashforth scheme is used for time-stepping.

Dynamics of thin active sheets: To motivate the study to come we first consider the invasion of a concentrated cylindrical colony, Gaussian in cross-section, of motile particles initially aligned in the $\mathbf{\hat{x}}$ direction into an empty viscous fluid, shown in Fig. 1a. The associated global flow field is exactly that of a single regularized force dipole, resulting in the case of pullers in a stable concentration elongation in a direction orthogonal to the original swimmer orientation [] (2018). For pushers the colony-induced velocity field changes sign and elongation is parallel to the swimming direction, but if slightly perturbed a transverse concentration instability ensues. Fore-aft symmetry is broken due to particle motility; the colony splays at the leading tip on the right, while undergoing a periodic folding at the rear reminiscent of the buckling of planar viscous jets Cruickshank and Munson (1981) and extruded beams Gosselin et al. (2014) into viscous fluids.

To better understand this dramatic evolution we turn to the behavior of an infinite sheet of particles which are initially in alignment. Figure 1b shows the evolution of a distribution of immotile (top) and motile (bottom) pusher particles, initially confined to a thin band and with a small transverse concentration perturbation. Early stages show rapid growth of the wave amplitude. At the same time individual particles are rotated towards the principal direction of the velocity gradient, so that they remain nearly tangent to the concentration band which results in a secondary instability and self-folding. The same structures are observed again and again on smaller length scales, though particle motility breaks left-right symmetry and significantly alters the structure of subsequent folding events. At longer times the system is finally drawn to an unsteady roiling state, with uniform concentration for immotile particles ($c$ satisfies a pure advection-diffusion equation in Eq. (5) in this case) or concentration bands described by Saintillan & Shelley Saintillan and Shelley (2008) for motile particles.

To analyze the instability, let $h(x,t)$ and $\phi(x,t)$ represent the vertical displacement and polarity of the line distribution, respectively, with $\mathbf{n}=\langle\mathbf{p}\rangle/c=(\cos\phi,\sin\phi,0)$ the normalized polarity. We study the dynamics of this line distribution through its far-field self-influence. For small $h$ and $\phi$, and solving (4) for the velocity field [] (2018) we find

[TABLE]

where $v$ is the vertical component of velocity evaluated on the flat surface $h(x,t)=0$, and $\mathcal{H}[\cdot]$ is the Hilbert transform,

[TABLE]

The Hilbert transform is diagonalized by changing to a Fourier basis, with $\mathcal{H}[e^{iqx}]=-i\,\mbox{sign}(q)e^{iqx}$. The Ansatz $(h(x,t),\phi(x,t))=\sum_{k}(\hat{h}_{k}(t),\hat{\phi}_{k}(t))\exp(2\pi ikx/L)$ therefore results in a quadratic eigenvalue problem, and the eigenvalues

[TABLE]

where $\gamma_{k}=\sqrt{\alpha\kappa(\alpha\kappa k^{2}-16iLV_{0}k|k|)}$. A comparison to numerics is shown in Fig. 2 for motile pushers with three negative dipole strengths along with the theory for a wider range of $\alpha$ and $k$. The analytical predictions are accurate for the entire range studied, with discrepancies owing to the vertical periodic boundary condition and the non-vanishing thickness of the line distribution in the simulations.

In the immotile case, $V_{0}=0$ (or in the limit as $|\alpha|/V_{0}\to\infty$), sheets of pushers are all unstable and sheets of pullers are all stable, with respective growth and decay rates both given by $-\pi\alpha\kappa|k|/(2L^{2})$. This behavior owes to the velocity field created by the active stress, illustrated in Fig. 2 (see Supplementary Movie M2), which either amplifies or damps the initial concentration perturbation. The linearized dynamics are now governed solely by the equation

[TABLE]

This expression establishes an unexpected connection to well-studied phenomena in entirely different settings: interfacial instabilities in gravity or pressure-driven Hele-Shaw problems, or two-dimensional flows in porous media, without surface tension, whose flow is governed by Darcy’s Law Tryggvason and Aref (1983); Hou et al. (1994) (also known as the Muskat problem Siegel et al. (2004); Córdoba et al. (2011)). There, as in the present setting, the classical Saffman-Taylor or Rayleigh-Taylor instabilities are modified to exponential growth rate dependence which is linear in $|k|$ Saffman and Taylor (1958). Such interfacial instabilities are associated with the formation of singularities in free-surface flows, for example the finite-time “Moore singularity” development on a vortex sheet in an inviscid fluid with no surface tension described by the Kelvin-Helmholtz instability, a higher-order system that shares linear growth rate dependence on $|k|$ Moore (1979); Meiron et al. (1982); Krasny (1986); Shelley (1992); Cowley et al. (1999). We thus observe an identical initial growth behavior, but nonlinear terms for large amplitude waves result in a unique folding event in $t\approx 2.5$ in Fig. 1b and very different long-time behavior.

Meanwhile, in the motile case, $V_{0}=1$, sheets of pushers remain unstable for any dipole strength. Sheets of pullers, however, excite a positive-real-part eigenvalue in (9). Unlike in the case of pushers the maximal eigenvalue is not monotonic in the force dipole strength (Fig. 3). Expanding about small $\alpha$, the largest eigenvalue is found when $\alpha=2L/\kappa$, at which point $\mbox{Re}[\lambda_{+}]=\pi k(2L)^{-1}$ (since $\mbox{Re}[\lambda_{+}]\sim\pi k\alpha\kappa(4L^{2})^{-1}$). For either motile or immotile pullers, the velocity field (oppositely signed to that illustrated in Fig. 2) rapidly damps the initial displacement, $h(x,t)$, and now rotates particles towards the direction perpendicular to the concentration band. Motility, however, allows the displacement and orientation fields to synchronize, leading ultimately to rapid growth of the concentration band amplitude and a large departure away from the initial profile, as shown in Fig. 3.

For the motile suspensions above with non-zero dipole strength $\alpha$ there can be competing effects; in particular, if $\max(\mbox{Re}[\lambda_{\pm}])<0$ all solutions to the linear system eventually arrive at the stable base state, but if the system departs from the linearized region of phase space fast enough such solutions may not be realized in the fully nonlinear dynamics. This potential for departure is seen most clearly if the particles are not stress-generating: with $\alpha=0$ the wave amplitude grows linearly in time since any particles with nonzero initial orientation angle drift off without resistance along characteristic curves.

Isotropic suspensions remain velocity free: a “no-flow theorem”: Assuming uniqueness of solutions to Eqs. (1)-(4), active suspensions of motile or immotile particles modeled by Eqs. (1)-(4) which are initially isotropic in orientation, $\Psi(\mathbf{x},\mathbf{p},t=0)=\Psi_{0}(\mathbf{x})$, result in no fluid flow, $\mathbf{u}(\mathbf{x},t)=\mathbf{0}$, everywhere and for all time

Sketch of the proof: The proof assumes uniqueness of solutions for Eqs. (1)-(4), which was shown for two-dimensions by Chen & Liu Chen and Liu (2013). Consider first the solution $\Psi^{*}$ to the Smoluchowski equation without velocity terms,

[TABLE]

with an initial condition which is isotropic in orientation. The velocity field generated by this solution, $\mathbf{u}[\Psi^{*}]$, is given in Fourier space by

[TABLE]

where $\mathbf{\hat{k}}=\mathbf{k}/|\mathbf{k}|$. Writing $\mathbf{k}$ in a spherical (3D) or polar (2D) coordinate system about $\mathbf{p}$ we find $\hat{\bm{\Sigma}}_{a}\cdot\mathbf{k}=\lambda_{\mathbf{k}}(t)\mathbf{k}$ for some scalar function $\lambda_{\mathbf{k}}(t)$. Hence $\hat{\mathbf{u}}_{\mathbf{k}}[\Psi^{*}]=\mathbf{0}$ and then $\mathbf{u}[\Psi^{*}]=\mathbf{0}$. Since $\mathbf{u}[\Psi^{*}]=\mathbf{0}$, $\Psi^{*}$ also solves Eqs. (1)-(4) with velocity terms included. By the uniqueness assumption we finally have that $\mathbf{u=0}$ everywhere and for all time for any initially isotropic distribution. A more detailed proof is included in the Supplementary Material.

The result is surprising since the system immediately loses orientational isotropy (see Fig. 4), which would suggest the quick onset of a non-trivial flow field but this is not observed. Physically, the active force $\mathbf{f}_{a}=\nabla\cdot\bm{\Sigma}_{a}$ is non-trivial for any $t>0$ but it is curl-free, so by the Helmholtz decomposition theorem $\mathbf{f}_{a}=\nabla\lambda(t)$ for some scalar field $\lambda(t)$, which thus only modifies the pressure. As time progresses the force distribution evolves with the local active particle alignment, illustrated for two initially uniform colonies in Fig. 4, but the expanding colonies simply pass through each other as linear waves. This behavior can be inferred even when including two-particle correlations Stenhammar et al. (2017).

Moreover, any distributions which result in $\mathbf{u}=\mathbf{0}$ for all time may be superimposed without generating a velocity field, for instance a random isotropic distribution may be perturbed by another distribution which has the property that $\nabla\cdot\langle\mathbf{pp}\rangle=\nabla\chi(t)$ for any scalar $\chi(t)$, and still $\mathbf{u=0}$ for all time. Physics which introduce nonlinearities in Eq. (1), such as near-field steric repulsion, are expected to nullify the theorem.

The stability of the theorem to an initial localized alignment is not simply determined, as the initially isotropic state is not a stable base state. However, in light of the stability of the isotropic state of uniform concentration to large wavenumber perturbations Saintillan and Shelley (2008) we expect an initial damping back towards isotropy. But on an extremely long timescale in a sufficiently large domain, the low wavenumber residue of the initial disruption is expected to lead to eventual growth along with a non-trivial flow. We have verified this prediction in at least one setting by numerical simulation, placing an aligned colony as in Fig. 1a into a random concentration field which is orientationally isotropic. Persistent nematic alignment, for instance due to a boundary, may result in a more immediate transition to a global mean flow.

Conclusion: We have investigated colonies of active particles in the dilute regime as they invade a quiescent fluid. Colony-scale elongation depends on the sign of the active stress and can result in a self-buckling and self-folding cascade. Exponential growth at small times, linear in $|k|$, is mathematically equivalent to the Saffman-Taylor instability in a Hele-Shaw cell, or Rayleigh-Taylor instability in two-dimensional flow through a porous medium. The stability of sheets of pullers depends on particle motility with a growth rate which is non-monotonic in the dipole strength. Strikingly, a suspension modeled by pure far-field hydrodynamic interactions which is initially isotropic in orientation, even though isotropy is not preserved, results in no mean-field fluid flow everywhere and for all time.

Acknowledgements: This project was initiated at the Woods Hole Oceanographic Institute as part of the Geophysical Fluid Dynamics summer program. Financial support is acknowledged by MJS (National Science Foundation Grants DMR-0820341 [NYU MRSEC], DMS-1463962, and DMS-1620331) and SES (DMR 767-1121288 [UW MRSEC] and DMS-1661900).