^{1}, Michael J. Shelley

^{2}and David Saintillan

^{1,a)}

### Abstract

Suspensions of active particles, such as motile microorganisms and artificial microswimmers, are known to undergo a transition to complex large-scale dynamics at high enough concentrations. While a number of models have demonstrated that hydrodynamic interactions can in some cases explain these dynamics, collective motion in experiments is typically observed at such high volume fractions that steric interactions between nearby swimmers are significant and cannot be neglected. This raises the question of the respective roles of steric vs hydrodynamic interactions in these dense systems, which we address in this paper using a continuum theory and numerical simulations. The model we propose is based on our previous kinetic theory for dilute suspensions, in which a conservation equation for the distribution function of particle configurations is coupled to the Stokes equations for the fluid motion [D. Saintillan and M. J. Shelley, “Instabilities, pattern formation, and mixing in active suspensions,” Phys. Fluids20, 123304 (Year: 2008)]10.1063/1.3041776. At high volume fractions, steric interactions are captured by extending classic models for concentrated suspensions of rodlike polymers, in which contacts between nearby particles cause them to align locally. In the absence of hydrodynamic interactions, this local alignment results in a transition from an isotropic base state to a nematic base state when volume fraction is increased. Using a linear stability analysis, we first investigate the hydrodynamic stability of both states. Our analysis shows that suspensions of pushers, or rear-actuated swimmers, typically become unstable in the isotropic state before the transition occurs; suspensions of pullers, or head-actuated swimmers, can also become unstable, though the emergence of unsteady flows in this case occurs at a higher concentration, above the nematic transition. These results are also confirmed using fully nonlinear numerical simulations in a periodic cubic domain, where pusher and puller suspensions are indeed both found to exhibit instabilities at sufficiently high volume fractions; these instabilities lead to unsteady chaotic states characterized by large-scale correlated motions and strong density fluctuations. While the dynamics in suspensions of pushers are similar to those previously reported in the dilute regime, the instability of pullers is novel and typically characterized by slower dynamics and weaker hydrodynamic velocities and active input power than in pusher suspensions at the same volume fraction.

The authors thank E. Lushi for illuminating conversations on sterically induced stresses, and acknowledge funding from National Science Foundation (NSF) Grant Nos. DMS-0930930 and DMS-0930931 and from the (U.S.) Department of Energy (DOE) Grant No. DE-FG02-88ER25053. Computational resources were provided by Teragrid Grant No. TG-CTS100007.

I. INTRODUCTION

II. KINETIC MODEL

A. Conservation equation

B. Steric interactions

C. Diffusion coefficients

D. Mean-field fluid velocity

E. Non-dimensionalization

III. STABILITY ANALYSES

A. Isotropic and nematic base states

B. Linearized equations and eigenvalue problem

C. Stability of the isotropic base state

1. Eigenvalue problem

2. Spherical harmonic expansion

D. Stability of the nematic base states

E. Summary of stability analyses and stability diagrams

IV. NUMERICAL SIMULATIONS

A. Simulation method and parameter selection

B. Results and discussion

V. CONCLUDING REMARKS

### Key Topics

- Suspensions
- 119.0
- Hydrodynamics
- 37.0
- Diffusion
- 29.0
- Tensor methods
- 27.0
- Eigenvalues
- 21.0

## Figures

(a) Solutions of the equation g(δ) = 0, where g is defined in Eq. (31) , as functions of ξ = 2U 0ν/d. Full lines show the branches that, at a given value of ξ, minimize the steric interaction energy. (b) Steric interaction energy E(ξ) along each of the three branches found in (a). (c) Orientation distributions for ξ = 20 corresponding to the three solution branches.

(a) Solutions of the equation g(δ) = 0, where g is defined in Eq. (31) , as functions of ξ = 2U 0ν/d. Full lines show the branches that, at a given value of ξ, minimize the steric interaction energy. (b) Steric interaction energy E(ξ) along each of the three branches found in (a). (c) Orientation distributions for ξ = 20 corresponding to the three solution branches.

Marginal stability curve in the long-wave limit (k → 0) in terms of effective volume fraction ν and dimensionless dipole strength αν/30d, for various values of 2U 0/d [see Eq. (68) ]. In this plot, β = 1.75, corresponding to a particle aspect ratio of r ≈ 10.

Marginal stability curve in the long-wave limit (k → 0) in terms of effective volume fraction ν and dimensionless dipole strength αν/30d, for various values of 2U 0/d [see Eq. (68) ]. In this plot, β = 1.75, corresponding to a particle aspect ratio of r ≈ 10.

Numerical solutions of the dispersion relations Eqs. (60) and (61) for the complex growth rates λ0, λ1, and λ2 of azimuthal modes m = 0, 1, 2 as functions of the wavenumber k. (a)–(c) show the real parts, while (d)–(f) show the imaginary parts. In (b) and (e), we set β = 0 to isolate the leading-order effect of the steric torque, and α/4U 0 = −0.2 (pushers).

Numerical solutions of the dispersion relations Eqs. (60) and (61) for the complex growth rates λ0, λ1, and λ2 of azimuthal modes m = 0, 1, 2 as functions of the wavenumber k. (a)–(c) show the real parts, while (d)–(f) show the imaginary parts. In (b) and (e), we set β = 0 to isolate the leading-order effect of the steric torque, and α/4U 0 = −0.2 (pushers).

Unstable range of wavenumbers as a function of ξ = 2U 0ν/d for azimuthal modes m = 0 (a), 1 (b), and 2 (c). For mode 1, we set β = 0 to isolate the leading-order effect of the steric torque.

Unstable range of wavenumbers as a function of ξ = 2U 0ν/d for azimuthal modes m = 0 (a), 1 (b), and 2 (c). For mode 1, we set β = 0 to isolate the leading-order effect of the steric torque.

Dependence of the maximum reduced growth rate Re(λ) governing the stability of the nematic base states on the dimensionless active stress magnitude α, for ξ = 20: (a) branch 2 and (b) branch 3. Results for two different wave orientations Θ are shown. Insets show close-ups on the region near the origin.

Dependence of the maximum reduced growth rate Re(λ) governing the stability of the nematic base states on the dimensionless active stress magnitude α, for ξ = 20: (a) branch 2 and (b) branch 3. Results for two different wave orientations Θ are shown. Insets show close-ups on the region near the origin.

Dependence of the maximum reduced growth rate Re(λ) on the parameter ξ, along both nematic branches, for two different wave directions Θ: (a) branch 2, Θ = 0; (b) branch 2, Θ = π/2; (c) branch 3, Θ = 0; and (d) branch 3, Θ = π/2.

Dependence of the maximum reduced growth rate Re(λ) on the parameter ξ, along both nematic branches, for two different wave directions Θ: (a) branch 2, Θ = 0; (b) branch 2, Θ = π/2; (c) branch 3, Θ = 0; and (d) branch 3, Θ = π/2.

Dependence of the maximum reduced growth rate Re(λ) on: (a) wave direction Θ (in the limit of k → 0), and (b) wavenumber k (for a wave with orientation Θ = 0). Both plots were obtained on branch 2, with ξ = 20.

Dependence of the maximum reduced growth rate Re(λ) on: (a) wave direction Θ (in the limit of k → 0), and (b) wavenumber k (for a wave with orientation Θ = 0). Both plots were obtained on branch 2, with ξ = 20.

Stability diagrams for: (a) movers (α = 0), (b) pushers (α < 0), and (c) pullers (α > 0). A branch is labeled unstable if there exists a positive growth rate Re(λ) > 0. In the case of movers, branch 2 is only weakly unstable, as the growth rates on that branch are very low (two orders of magnitude lower than on other unstable branches).

Stability diagrams for: (a) movers (α = 0), (b) pushers (α < 0), and (c) pullers (α > 0). A branch is labeled unstable if there exists a positive growth rate Re(λ) > 0. In the case of movers, branch 2 is only weakly unstable, as the growth rates on that branch are very low (two orders of magnitude lower than on other unstable branches).

Simulation results for pushers and pullers. Panels (a)–(d) show the concentration fields c(x, t) at an arbitrary time after the initial transient for the following cases: (a) pushers at ν = 0.05, (b) pushers at ν = 0.07, (c) pushers at ν = 0.2, and (d) pullers at ν = 0.2. Panels (e)–(h) show the corresponding nematic parameter fields N(x, t) defined in Eq. (70) , for the same cases. Panels (i)–(l) show the nematic orientation fields in a two-dimensional slice, obtained by taking the eigenvector of with largest eigenvalue (enhanced online). [URL: http://dx.doi.org/10.1063/1.4812822.1]doi: 10.1063/1.4812822.1.

Simulation results for pushers and pullers. Panels (a)–(d) show the concentration fields c(x, t) at an arbitrary time after the initial transient for the following cases: (a) pushers at ν = 0.05, (b) pushers at ν = 0.07, (c) pushers at ν = 0.2, and (d) pullers at ν = 0.2. Panels (e)–(h) show the corresponding nematic parameter fields N(x, t) defined in Eq. (70) , for the same cases. Panels (i)–(l) show the nematic orientation fields in a two-dimensional slice, obtained by taking the eigenvector of with largest eigenvalue (enhanced online). [URL: http://dx.doi.org/10.1063/1.4812822.1]doi: 10.1063/1.4812822.1.

Spatially averaged nematic parameter ⟨N(x, t)⟩, defined in Eq. (70) , as a function of time in suspensions of pushers (α = −1), pullers (α = +1), as well as shakers (V 0 = 0, α = ±1) at various volume fractions.

Spatially averaged nematic parameter ⟨N(x, t)⟩, defined in Eq. (70) , as a function of time in suspensions of pushers (α = −1), pullers (α = +1), as well as shakers (V 0 = 0, α = ±1) at various volume fractions.

Time evolution of the spatial averages of (a) |n(x, t)|, (b) |u(x, t)|, and (c) |n(x, t) + u(x, t)|, in suspensions of pushers at ν = 0.05, 0.07, and 0.2, and pullers at ν = 0.2.

Time evolution of the spatial averages of (a) |n(x, t)|, (b) |u(x, t)|, and (c) |n(x, t) + u(x, t)|, in suspensions of pushers at ν = 0.05, 0.07, and 0.2, and pullers at ν = 0.2.

Time evolution of the correlation defined in Eqs. (72)–(74) : (a) C 1(t), (b) C 2(t), and (c) C 3(t), in suspensions of pushers at ν = 0.05, 0.07, and 0.2, and pullers at ν = 0.2.

Time evolution of the correlation defined in Eqs. (72)–(74) : (a) C 1(t), (b) C 2(t), and (c) C 3(t), in suspensions of pushers at ν = 0.05, 0.07, and 0.2, and pullers at ν = 0.2.

Mean active power P(t) defined in Eq. (75) as a function of time, in suspensions of pushers at ν = 0.05, 0.07, and 0.2, and pullers at ν = 0.2. The plot also shows results for shakers (V 0 = 0, α = ±1) at ν = 0.2.

Mean active power P(t) defined in Eq. (75) as a function of time, in suspensions of pushers at ν = 0.05, 0.07, and 0.2, and pullers at ν = 0.2. The plot also shows results for shakers (V 0 = 0, α = ±1) at ν = 0.2.

Simulation results for shakers (V 0 = 0) at an effective volume fraction of ν = 0.2 (ξ = 26.92). Panels (a)–(c) are for pushers (α = −1) and (d)–(f) for pullers (α = +1). (a) and (d) show the nematic parameter N(x, t) defined in Eq. (70) ; (b) and (c) show the nematic orientation fields in a two-dimensional slice, obtained by taking the eigenvector of with largest eigenvalue; and (c) and (f) show the hydrodynamic velocity fields in a two-dimensional slice.

Simulation results for shakers (V 0 = 0) at an effective volume fraction of ν = 0.2 (ξ = 26.92). Panels (a)–(c) are for pushers (α = −1) and (d)–(f) for pullers (α = +1). (a) and (d) show the nematic parameter N(x, t) defined in Eq. (70) ; (b) and (c) show the nematic orientation fields in a two-dimensional slice, obtained by taking the eigenvector of with largest eigenvalue; and (c) and (f) show the hydrodynamic velocity fields in a two-dimensional slice.

Standard deviation of the concentration field c(x, t) as a function of time, in suspensions of pushers (α = −1), pullers (α = +1), and shakers (V 0 = 0) of both types, at ν = 0.2.

Standard deviation of the concentration field c(x, t) as a function of time, in suspensions of pushers (α = −1), pullers (α = +1), and shakers (V 0 = 0) of both types, at ν = 0.2.

## Multimedia

Article metrics loading...

Full text loading...

Commenting has been disabled for this content