Analytical structure, dynamics, and coarse-graining of a kinetic model of an active fluid

TL;DR

This paper develops a simplified yet accurate kinetic and Q-tensor model for active suspensions of extensile particles, capturing complex flow and defect dynamics with reduced computational effort.

Contribution

It introduces the BQ-tensor theory, a novel coarse-grained model that efficiently describes active nematic suspensions and their boundary-driven defect dynamics.

Findings

The BQ-tensor model accurately reproduces complex suspension dynamics.

Novel defect production and absorption observed in strong nematic alignment regimes.

Boundary effects induce a Fredericks-like transition in bi-concave domains.

Abstract

We analyze one of the simplest active suspensions with complex dynamics: a suspension of immotile "Extensor" particles that exert active extensile dipolar stresses on the fluid in which they are immersed. This is relevant to several experimental systems, such as recently studied tripartite rods that create extensile flows by consuming a chemical fuel. We first describe the system through a Doi-Onsager kinetic theory based on microscopic modeling. This theory captures the active stresses produced by the particles that can drive hydrodynamic instabilities, as well as the steric interactions of rod-like particles that lead to nematic alignment. This active nematic system yields complex flows and disclination defect dynamics very similar to phenomenological Landau-deGennes Q-tensor theories for active nematic fluids, as well as by more complex Doi-Onsager theories for polar microtubule/motor-protein systems. We apply the quasi-equilibrium Bingham closure, used to study suspensions of passive microscopic rods, to develop a non-standard Q-tensor theory. We demonstrate through simulation that this "BQ-tensor" theory gives an excellent analytical and statistical accounting of the suspension's complex dynamics, at a far reduced computational cost. Finally, we apply the BQ-tensor model to study the dynamics of Extensor suspensions in circular and bi-concave domains. In circular domains, we reproduce previous results for systems with weak nematic alignment, but for strong alignment find novel dynamics with activity-controlled defect production and absorption at the boundaries of the domain. In bi-concave domains, a Fredericks-like transition occurs as the width of the neck connecting the two disks is varied.

Figures

Equations

$$\frac{\partial \psi }{\partial t}+\nabla \cdot (\dot{\mathbf{x}}\psi )+{\nabla }_p\cdot (\dot{\mathbf{p}}\psi )=0.$$

$$\varphi (\mathbf{x},t)$$

$$\mathbf{q}(\mathbf{x},t)$$

$$\mathbf{D}(\mathbf{x},t)$$

$$\mathbf{S}(\mathbf{x},t)$$

$$\dot{\mathbf{x}}$$

$$\dot{\mathbf{p}}$$

$$\nabla \Pi -\mu \Delta \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$\dot{\mathbf{x}}$$

$$\dot{\mathbf{p}}$$

$$\nabla \Pi -\Delta \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$\mathbf{\sigma }=\alpha \mathbf{D}+\beta \mathbf{S}:\mathbf{E}-2\zeta \beta \left(\mathbf{D}\cdot \mathbf{D}-\mathbf{S}:\mathbf{D}\right).$$

$$\alpha =\frac{{\sigma }_a}{\mu |u_0|l^2},\zeta =\frac{{\zeta }_0}{|u_0|l^2},\beta =\frac{\pi r\nu }{6\ln (2r)},d_T=\frac{\nu }{b|u_0|}{D}_T,d_R=\frac{b}{\nu |u_0|}{D}_R.$$

$$\frac{\partial \varphi }{\partial t}+\mathbf{u}\cdot \nabla \varphi =d_T\Delta \varphi .$$

$$\frac{d}{dt}{\int }_{\Omega }dV\frac{1}{2}{\varphi }^2=-d_T{\int }_{\Omega }dV|\nabla \varphi {|}^2.$$

$$\mathcal{E}=\mathcal{S}+\kappa \mathcal{P}={\int }_{\Omega }dV{\int }_{S}d{S}_p\frac{\Psi }{{\Psi }_0}\ln \frac{\Psi }{{\Psi }_0}-\kappa {\int }_{\Omega }dV\zeta (\mathbf{D}:\mathbf{D}),$$

$$\dot{\mathcal{E}}=$$

$$+2d\zeta \left(1+\frac{2\beta d_R}{\alpha }\right){\int }_{\Omega }dV\mathbf{D}:\mathbf{D}+\frac{d\zeta \beta d_T}{\alpha }{\int }_{\Omega }dV{\left|\nabla \mathbf{D}\right|}^2$$

$$-{\int }_{\Omega }dV{\int }_{S}d{S}_p\left(d_T{\left|\nabla \ln \Psi \right|}^2+d_R{\left|{\nabla }_p\ln \Psi \right|}^2\right),$$

$$\frac{\partial {\Psi }^{\prime }}{\partial t}-\frac{3}{4\pi }\mathbf{p}\mathbf{p}:{\mathbf{E}}^{\prime }=\frac{3\zeta }{2\pi }\mathbf{p}\mathbf{p}:{\mathbf{D}}^{\prime }+d_T\Delta {\Psi }^{\prime }+d_R{\nabla }_p^2{\Psi }^{\prime }.$$

$$\nabla \Pi -\left(1+\frac{\beta }{15}\right)\Delta {\mathbf{u}}^{\prime }$$

$$\nabla \cdot {\mathbf{u}}^{\prime }$$

$$\frac{\partial \mathbf{g}}{\partial t}=\left(-C_1\left(\beta \right)\alpha +C_2\left(\beta \right)\zeta -6d_R\right)\mathbf{g}+d_T\Delta \mathbf{g},$$

$$\sigma =-C_1(\beta )\alpha +C_2(\beta )\zeta -6d_R-d_T{k}^2.$$

$$\sigma (k=0^{+})=-\left(\frac{3}{\beta +15}\right)\alpha +2\left(\frac{\beta +6}{\beta +15}\right)\zeta -6d_R.$$

$$\frac{\partial \mathbf{q}}{\partial t}=-\mathbf{u}\cdot \nabla \mathbf{q}+(\mathbf{I}-\mathbf{q}\mathbf{q})\cdot \nabla \mathbf{u}\cdot \mathbf{q}.$$

$$\nabla \Pi -\Delta \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$\frac{\partial }{\partial t}\left(\mathbf{q}\mathbf{q}\right)=-\left(\mathbf{u}\cdot \nabla \right)\left(\mathbf{q}\mathbf{q}\right)+\nabla \mathbf{u}\cdot \left(\mathbf{q}\mathbf{q}\right)+\left(\mathbf{q}\mathbf{q}\right)\cdot \nabla {\mathbf{u}}^T-2\left(\mathbf{q}\mathbf{q}:\mathbf{E}\right)\mathbf{q}\mathbf{q}.$$

$${\mathbf{D}}^{\nabla }+2\left(\mathbf{D}:\mathbf{E}\right)\mathbf{D}=0,$$

$$\nabla \Pi -\Delta \mathbf{u}$$

$$\nabla \cdot \mathbf{u}$$

$$\frac{\partial \mathbf{X}}{\partial t}(\mathbf{\alpha },t)=\mathbf{u}(\mathbf{X}(\mathbf{\alpha },t),t)$$

$$\frac{\partial \mathbf{r}}{\partial t}=\left(\mathbf{I}-\mathbf{r}\mathbf{r}\right)\cdot \nabla \mathbf{u}\cdot \mathbf{r}.$$

Full text

I Introduction

Active suspensions are non-equilibrium materials composed of suspended particles whose activity, driven by consumption of a local fuel, lead to particle motions or local induced flows ramaswamy10 ; SS2013 . Examples include bacterial swarms dombrowski04 ; zhang10 , collections of synthetic colloidal particles PalacciEtAl2013 ; bricard13 , and mixtures of cytoskeletal filaments driven by molecular motors SanchezEtAl2012 ; Shelley2016 . One origin of the instabilities and complex dynamics in active suspension is the stresses created in the surrounding solvent by the particles’ activity. The solvent is a coupling medium for multiscale dynamics: particle interactions through the solvent can manifest at the system scale as collective motion, which feeds back to alter their interactions. Theoretical and numerical studies have investigated various aspects of active suspensions at different scales, from single particle dynamics to suspension rheology. While discrete-particle simulations that incorporate long-ranged hydrodynamic effects can capture observed large-scale features of these systems hernandez-ortiz05 ; SS2007 ; SS2011 ; Yeo2015 ; Delmotte2015 , constructing accurate continuum models is essential to their characterization and analysis, and for making new predictions of macroscale behaviors SS2013 .

Continuum descriptions of active suspensions of self-propelled particles, such as bacteria, have been based on a variety of approaches simha02 ; aranson07 ; underhill08 ; SS2008a ; wolgemuth08 ; baskaran09 ; subramanian09 . In all, coarse-grained variables, such as particle concentration, concentration-dependent steric interactions, order parameters, and solvent stresses and velocity, are used to capture salient features of the macroscopic dynamics. In our own work, we and our collaborators have developed Doi-Onsager kinetic models to describe the evolution of suspensions of active rod-like particles SS2008a ; SS2008b ; HS2010 ; ESS2013 ; GBGBS2015a ; GBGBS2015b . In this approach, a Smoluchowski equation describes the evolution of a particle distribution function, coupled to a coarse-grained Stokes equation driven by the extra stresses created by particle activity and interactions. The fluxes and stresses of the model are derived from modeling how such particles interact with mean-field flows and other coarse-grained variables. For suspensions of self-propelled rods which interact only hydrodynamically SS2008a ; SS2008b this approach identifies characteristic aspects of the dynamics: either uniform isotropic or globally aligned suspensions have hydrodynamic instabilities that depend upon the self-propulsion mechanism, particle concentration, and system size SS2008b ; subramanian09 ; HS2010 , leading to the emergence of complex, turbulent-like dynamics. When such particles also interact sterically through concentration-dependent alignment torques, this generates a polar “active-nematic” model with an isotropic-to-nematic ordering transition. The new nematically-aligned steady states can again be unstable to hydrodynamic instabilities, again leading to complex dynamics but modulated by an effective liquid crystal elasticity ESS2013 . Similar but more elaborate kinetic models have been derived for “bio-active” suspensions SanchezEtAl2012 composed of microtubules whose relative motions are driven by processive cross-linking molecular motors such as kinesin-1 GBGBS2015a ; GBGBS2015b . There, microtubule motion and motor activity create polarity-dependent fluxes and stresses that drive hydrodynamic instability, and a complex dynamics characterized by nematic director fields with motile disclination defects that are continuously nucleated and annihilated.

While having the advantage of being founded upon microscopic modeling, such kinetic theories can carry the price of complexity. For Doi-Onsager theories as above, the distribution function depends upon both particle position ${\bf x}$ and orientation vector ${\bf p}$ ($|{\bf p}|=1$) and so, in three dimensions, has five independent variables plus time $t$. This makes their simulation challenging. Hence, simplifications have been sought, often following from moment-closure schemes that eliminate the ${\bf p}$ variable, yielding more approximate models that describe the evolution of lower-order ${\bf p}$-moments of the distribution function, such as the concentration $\phi({\bf x},t)$ (zeroth-moment), the polarity vector ${\bf q}({\bf x},t)$ (vector of first-moments), and the unnormalized tensor order parameter ${\bf D}({\bf x},t)$ (tensor of second-moments) woodhouse2012 ; SS2013 ; ezhilin2015 ; Theillard2017 . Tensorial moment-closure models can resemble, with important differences, tensor models based upon phenomenological Landau-deGennes $Q$-tensor liquid crystal theories TGY2013 ; giomi13 ; Giomi2015 .

In this paper, we use a Doi-Onsager theory to study what is perhaps the simplest active suspension: a collection of non-motile, yet mobile, rod-like particles that exert active dipolar stresses on the solvent and interact with each other through hydrodynamic coupling and steric alignment torques. When the stresses are extensile, we call these particles “Extensors”. While this is a special case of the systems above, it is well-worth studying in its own right. Firstly, it presents all the basic transitions and instabilities associated with motile suspensions but in a more revealing and simplified form. Secondly, there are several types of active suspensions described at least qualitatively by this theory, including particles that elongate through stretching or growth, creating the extensile flows associated with destabilizing dipolar stresses. Such particle dynamics occurs in some instances of bacterial cell division Adams09 , as well as in intriguing liquid-crystal phase transitions BPR1987 ; SU2000 . Such stretching is also thought to occur for bundles of microtubules (formed by a depletion interaction) as they are “polarity-sorted” by kinesin-1 motor-proteins SanchezEtAl2012 ; keber14 . A very different system is composed of tripartite Au-Pt-Au nanomotors Wykes16 ; Jewell2016 for which catalytic surface reactions with a surrounding aqueous hydrogen peroxide solution paxton04 create extensile surface flows. The self-assembly of motile aggregates by small numbers of such nonmotile particles has been studied experimentally by Jewell et al. Jewell2016 and Wykes et al. Wykes16 , and numerically by Pandey et al. Pandey2016 .

Thirdly, we use this simple kinetic theory as a testing ground for constructing coarse-grained macro models for active suspensions using the Bingham moment closure Bingham74 , a closure scheme used successfully in the study of passive rod suspensions Chaubal98 ; Feng98 . The Bingham closure expresses fourth-moment tensors arising in the kinetic theory in terms of ${\bf D}({\bf x},t)$ using the so-called Bingham distribution, which is a quasi-equilibrium Ansatz. The Doi-Onsager system is then expressed solely in terms of ${\bf D}$. We call this $BQ$-tensor theory as ${\bf D}$ is an unnormalized form of the tensor order parameter ${\bf Q}$, itself the central evolved quantity in the Landau-deGennes theory of liquid crystal dynamics, which has been applied extensively to so-called active nematics TGY2013 ; giomi13 ; Giomi2015 ; GBGBS2015a ; GBGBS2015b . We show that the $BQ$-tensor theory closely or exactly reproduces the instabilities of the kinetic theory near homogeneous isotropic or nematic steady states. We then show through simulations that the $BQ$-tensor theory gives an excellent statistical accounting of the complex dynamics resolved by the kinetic theory model, both with and without steric interactions, of an active nematic suspension driven by extensile stresses.

Finally, having established the fidelity of the $BQ$-tensor theory, we use it to examine the dynamics of active suspensions under confinement. This is inspired by several recent experiments showing that confined collective suspensions of self-propelled particles can organize into auto-circulating states Wioland2013 ; wioland2016 , and develop emergent ordering and density shocks bricard13 . The dynamics of active suspensions has been studied previously using both discrete particle simulations and continuum models in straight channels ezhilin2015 ; tsang2016 , circular disks woodhouse2012 ; Lushi2014 , annuli, and racetracks Theillard2017 . For suspensions in a circular chamber we find perfect agreement with the previous results for suspensions that interact only hydrodynamically woodhouse2012 , and focus instead on the active nematic case where steric interactions are strong. We observe a plethora of dynamical states as the activity is increased. In all, defects are produced in the bulk, propagate, and annihilate, both with each other and at the boundaries. In biconcave geometries, which are essentially two disks smoothly connected by a bridge, we find evidence for a Friedricks transition in which the elasticity of the ordered fluid prevents flows between the two sides. When the neck is wide enough, activity overcomes the elasticity, and material (and defects) move between the two sides.

The paper is organized as follows. In Sect. II we develop the kinetic theory for microscopic rods that produce extensile surface flows. Except for the sign of the active dipolar stress, this theory is nearly identical to the classical Doi-Onsager models for passive liquid crystal polymers where the suspended molecules have fore-aft symmetry DE1986 ; Feng98 . We discuss its analytical and stability properties, as well as introduce simplified models, particularly coarse graining through the Bingham closure. In Sect. III, we study numerically the dynamics of Extensor suspensions. In periodic open domains, we implement a pseudo-spectral method to solve both the kinetic theory and the $BQ$-tensor theory, and compare the two descriptions in terms of emergent coherent structures, flow patterns, and statistical structure. We then simulate the dynamics of the $BQ$-tensor model in both circular and biconcave domains using a finite element method. Following a discussion in Sect. IV, we present some additional results in appendices.

II Doi-Onsager theory for a simple active suspension

II.1 A micro-mechanical model

Consider a suspension of $N$ rigid rod-like active particles in a volume $V$, each of length $l$ and diameter $b$. The rods are considered to be slender so that the aspect ratio $r=b/l<<1$. For definiteness we assume a simple form of activity: each rod produces a symmetric surface flow ${\bf v}(s)=sgn(s)u_{0}{\bf p}$ with $-l/2\leq s\leq l/2$ the signed arclength along the rod center-line, ${\bf p}$ its unit orientation vector, and $u_{0}$ the signed surface flow speed; see the schematic in Fig. 1. Because of surface flow symmetry such active rods are not motile. If $u_{0}>0$ there is an extensional straining flow along the rod, and a compressive one if $u_{0}<0$. These flows are associated with dipolar extra stresses. Directly following the derivation for Pusher suspensions HS2011 , assuming that active particles are small relative to the large-scale suspension flows they create in ensemble, and using slender body theory for the Stokes equations, one can calculate the rod center-of-mass velocity, its rate of rotation, and the force/length exerted by the rod upon the fluid which, consequently, gives the rod’s contributions to the added stress. The single particle contribution from the active surface flow is $\sigma_{a}{\bf p}{\bf p}$ where $\sigma_{a}=-\eta l^{2}u_{0}/8$ has units of force times length (i.e., a “stresslet”), and $\eta=8\pi\mu/|\ln(er^{2})|$ with $\mu$ the fluid viscosity. Note that $\sigma_{a}<0$ for extensional surface flows, and $\sigma_{a}>0$ for compressive ones. Overall, when $\sigma_{a}<0$ the system can exhibit activity-driven hydrodynamic instabilities simha02 ; SS2008a .

II.2 Doi-Onsager kinetic theory

As an expression of the conservation of particle number, one can derive a Smoluchowski equation for the (dimensional) distribution function $\psi({\bf x},{\bf p},t)$ for number density of particles at position ${\bf x}$ with orientation ${\bf p}$:

[TABLE]

The distribution function is chosen so that $\int_{V}dV\int_{S}dS_{p}\psi=N$. Here and for the remainder of this paper, derivative operators without subscripts (e.g., $\nabla$) denote spatial derivatives. The orientational gradient operator on the unit sphere is $\nabla_{p}=({\bf I}-{\bf p}{\bf p})\cdot\partial/\partial{\bf p}$. Moments of $\psi$ with respect to ${\bf p}$ arise naturally in the theory: the particle concentration $\phi$, (unnormalized) polarity field ${\bf q}$, second-moment tensor ${\bf D}$ (unnormalized tensor-order parameter), and fourth-moment tensor ${\bf S}$ are given by

[TABLE]

The trace-free (normalized) tensor order parameter, the so-called $Q$-tensor, is defined by ${\bf Q}({\bf x},t)=\phi({\bf x},t)^{-1}{\bf D}-{\bf I}/d$, with $d=2,3$ the spatial dimension. The tensor ${\bf Q}$ has a maximal nonnegative eigenvalue $\lambda_{max}$ satisfying $0\leq\lambda_{max}\leq(d-1)/d$. Assuming that $\lambda_{max}$ is isolated, then we call its associated unit eigenvector the director ${\bf m}$, and $0\leq s=\lambda_{max}d/(d-1)\leq 1$ the scalar order parameter.

As described in previous work SS2008a ; SS2008b ; ESS2013 ; GBGBS2015a ; GBGBS2015b , the conformational fluxes $\dot{\bf x}$ and $\dot{\bf p}$ arise from the dynamics of the particle’s center-of-mass position ${\bf x}$ and orientation ${\bf p}$, calculated using slender-body theory for the Stokes equations Keller1976 , and Maier-Saupe theory for steric alignment maier58 of rod-like particles. This yields

[TABLE]

where ${\bf u}$ is the background fluid velocity field. Equation (6) states that the rod’s positions move with the background fluid flow and undergo translational diffusion (with coefficient $D_{T}$), while Eq. (7) extends Jeffrey’s equation for the rotation of rods by the fluid velocity gradient to include a local alignment torque arising from a rod-rod interaction potential with strength $\zeta_{0}$. The last term of Eq. (7) models rotational diffusion with coefficient $D_{R}$.

The equations are closed by relating the incompressible background velocity ${\bf u}$ and pressure $\Pi$ to the extra-stress tensor ${\bm{\sigma}}$ produced by the particles’ presence batchelor70 , through a forced Stokes equation

[TABLE]

The stress tensor ${\bm{\sigma}}$ has three contributions: ${\bm{\sigma}}={\bm{\sigma}}_{a}+{\bm{\sigma}}_{c}+{\bm{\sigma}}_{s}$. The dipolar active stress tensor ${\bm{\sigma}}_{a}=\sigma_{a}{\bf D}$ arises from the active surface flows generated by particles. For Extensor particles, $\sigma_{a}$ is negative, and positive for “Contractors” (see Fig. 1). The so-called constraint stress tensor ${\bm{\sigma}}_{c}=\sigma_{c}{\bf S}:{\bf E}$, where ${\bf E}=1/2(\nabla{\bf u}+\nabla{\bf u}^{T})$ is the symmetric rate-of-strain tensor and $\sigma_{c}=\eta l^{3}/24$ arises from particle rigidity DE1986 ; ESS2013 . The stress tensor generated by steric interactions is ${\bm{\sigma}}_{s}=-\sigma_{s}({\bf D}\cdot{\bf D}-{\bf S}:{\bf D})$ with $\sigma_{s}=\pi\eta l^{3}\zeta_{0}/3$ DE1986 ; ESS2013 .

Following ESS2013 , we introduce an effective volume fraction $\nu=nbl^{2}$ where $n=N/V$ is the mean number density DE1986 ; ESS2013 . To nondimensionalize, we choose a characteristic length scale $l_{c}=b/\nu$, velocity scale $|u_{0}|$, and stress scale $\mu|u_{0}|/l_{c}$. The distribution function is normalized such that $(1/V)\int_{\Omega}dV\int_{S}dS_{p}~{}\Psi({\bf x},{\bf p},t)=1$. Equation  (1) retains its form, where the translational and rotational fluxes are now

[TABLE]

The Stokes equation becomes

[TABLE]

where the dimensionless stress tensor is given by

[TABLE]

The dimensionless control parameters are

[TABLE]

Given the normalization of $\Psi$, the average density is $\bar{\phi}=1$, and the state of uniform isotropy is given by $\Psi\equiv\Psi_{0}=1/4\pi$ in 3D, and $1/2\pi$ in 2D. In the case of spatially uniform density, ${\bf D}={\bf Q}+{\bf I}/d$. Note that $\alpha$ inherits the sign of $\sigma_{a}$, and that the nondimensional translational and rotational diffusion coefficients are proportional and inversely proportional, respectively, to the effective volume fraction $\nu$. We will sometimes refer to $\alpha$ as the activity and treat it as a free parameter though, in fact, it is a geometric parameter associated with particle shape and the relative placement of active and passive stresses upon it SS2008a ; HS2010 ; GBGBS2015a . For the specific case of the tripartite rods considered by Wykes et al. Wykes16 , it is estimated that $\alpha\approx-1$. Alternatively, if $u_{0}$ is interpreted as the velocity scale for a base level of activity, then increases in $\alpha$ (either positively or negatively) can be interpreted directly as multiplicative increases in surface velocity over this base level.

Rotational thermal fluctuations of rod-like particles also give rise to an extra stress of form $K{\bf D}$ with $K>0$ DE1986 . Our assumption and expectation is that the deterministic active stress, with $\alpha$ of either sign, is a much larger contribution.

A very similar model describes suspensions of stretching filaments where the distribution function now depends upon particle length $l(t)$. Assuming that particle length growth, $\dot{l}=f(l)$, is independent of particle position and orientation, one can derive a closed evolution equation for the marginal distribution by integrating out $l$. The main differences are that $u_{0}$ is replaced by $f(l)$, and $l$ dependencies are replaced by distributional averages. As expected, $\dot{l}>0$ yields extensile stresses.

Here we discuss some analytical properties of the kinetic model.

1. Concentration Fluctuations. Unlike suspensions of Pusher swimmers, which also produce extensile flows and where fluctuations can grow through nonlinear effects SS2008b ; HS2011 , concentration fluctuations in the nonmotile case decay to zero. Denote the fluid domain as $\Omega$, and assume either (i) $\Omega$ is a cubic domain $[-L/2,L/2]^{3}$ over which $\Psi$ is periodic, or (ii) there is neither mass nor particle flux across the boundary $\partial\Omega$. Integrating Eq. (1) over the unit ${\bf p}$-sphere, and using Eq. (6) and velocity incompressibility, then gives

[TABLE]

Integrating Eq. (16) over $\Omega$, integrating by parts and applying boundary conditions (periodic or zero flux) yields

[TABLE]

Hence, concentration fluctuations from equilibrium will decay to zero, regardless of particle type.

2. An Energy Identity. Like the motile case SS2008b , the kinetic theory has an “energy” identity governing the evolution of a total energy composed of the configurational entropy and the steric-interaction energy Han15 . Let

[TABLE]

where $\kappa=d\beta/\alpha$. The conformational entropy ${\cal S}$ is non-negative, and zero if and only if $\Psi\equiv\Psi_{0}$ (the state of uniform isotropy). The contribution ${\cal P}$ arises from the Maier-Saupe approximation to the classical Onsager potential Onsager1949 . We can show that ${\cal E}$ satisfies the identity

[TABLE]

where every integral density is non-negative (in particular, ${\bf D}:({\bf D}\cdot{\bf D}-{\bf S}:{\bf D})\geq 0$). The first term in braces is a generalized rate of viscous dissipation, the next three terms are generated by steric interactions (and hence are multipled by $\zeta$), and the last term is a negative-definite dissipation term arising from thermal fluctuations.

Apparently, the above identity reveals that the nature of the active stress, i.e., contractile or extensile, is a central determinant of system behavior. For a contractile active stress ($\alpha>0$) in the absence of steric interactions ($\zeta=0$), the energy $\mathcal{E}$ is driven to zero, corresponding to a uniform isotropic state. The same result holds for the motile case SS2008b and is consistent with the simulations of semi-dilute active particle suspensions SS2011 . When steric interactions are more dominant, ordering effects can drive active contractile systems into complex dynamics ESS2013 . When $\alpha<0$ (i.e., Extensor particles), the particle activity increases the energy in the system through the first term, while its effect upon the steric terms is somewhat complex.

3. The stability of uniform isotropic suspensions. The linear theory of Extensor suspensions has a simple and evocative structure that is obscured in the motile case. We give below the results of a stability calculation using the stream function ${\bm{\Phi}}$ (i.e., ${\bf u}=\nabla\times{\bm{\Phi}}$ ), and the details are given in Appendix A.

For an isotropic steady state we have $\Psi^{0}=1/4\pi$, ${\bf D}^{0}={\bf I}/3$, ${\bf u}^{0}={\bf 0}$, and $S^{0}_{ijkl}=\frac{1}{15}\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}+\delta_{ij}\delta_{kl}\right)$. Perturbing this steady state as ${\bf D}={\bf I}/3+\varepsilon{\bf{D^{\prime}}},\;{\bf{S}}={{\bf{S}}^{0}}+\varepsilon{\bf{S^{\prime}}},\;{\bf{u}}=\varepsilon{\bf{u^{\prime}}}$, where $\varepsilon<<1$, we find the 3D Smoluchowski equation (1) can be linearized at $O\left(\varepsilon\right)$ as SS2008a ; HS2010 ; ESS2013 :

[TABLE]

The momentum-balance equation in (12,13), together with Eq. (14), can be linearized as

[TABLE]

Rewriting Eq. (21) in terms of ${\bm{\Phi}}$ allows significant simplification. After some calculation (see derivations in Appendix A), and defining ${\bf g}=\nabla^{4}{\bm{\Phi}}$, the linearized equations reduce to

[TABLE]

where $C_{1}(z)=(3/5)/(1+z/15))$ and $C_{2}(z)=(4/5)+(6z/25)/(1+z/15)$. Thus, growth or decay rates of the plane-wave solutions are given by

[TABLE]

The revealing aspect of this linear analysis is that the effect of activity, steric interactions, and rotational diffusion all compete at the level of scale-independent dissipation and growth. Steric interactions associated with ordering drive the system away from isotropy, which is abetted by activity if $\alpha<0$ (i.e. for Extensor particles). Further, as is generic for active suspension theories in the absence of an external scale SS2008a ; SS2013 ; GBGBS2015a ; GBGBS2015b , we find the fastest growing mode occurs as $k\rightarrow 0^{+}$

[TABLE]

Hence, for periodic boundary conditions imposed on sufficiently large domains, the fastest growing scale is directly associated with the domain’s first periodic mode (which agrees well with our simulation results in Sect. III.)

As discussed below, another accessible steady state when steric interactions are strong is the nematic state. We will turn to a numerical treatment to study its stability except in the case of sharply aligned suspensions.

II.3 Reduced Models

In this section we will discuss two useful reductions of the kinetic theory. In the first, we investigate a special case of suspensions where the active particles are strictly aligned at each point in space. This is an exact reduction of the kinetic theory, and is useful for constructing phenomenological models illustrating the instabilities of aligned suspensions. Since dynamical simulation of the kinetic theory is demanding given the number of independent variable (e.g., three in space, two on the unit sphere, plus time in 3D), we investigate the Bingham closure which has been successfully applied to passive rod suspensions. This is a coarse graining that expresses the fourth-order ${\bf S}$ tensor in terms of the second-order ${\bf D}$ tensor, leading to a so-called $Q$-tensor dynamics which we term $BQ$-tensor theory. We find that it gives an excellent accounting of the complex dynamics we observe in the kinetic theory, and at a far lower computational cost by removing the orientational degrees of freedom.

II.3.1 Sharply aligned suspensions.

One analytically useful approximation is to assume negligible translational and rotational diffusion, and that the suspension is sharply aligned SS2008b . That is, we write $\Psi({\bf x},{\bf p},t)={\bar{\phi}}\delta[{\bf p}-{\bf q}({\bf x},t)]$ where $\delta$ is the Dirac delta function on the unit sphere, and further assume that the concentration is uniform, i.e., $\phi\equiv{\bar{\phi}}=1$. Integrating the Smoluchowski equation (1) against ${\bf p}$, and making use of $\nabla\cdot{\bf u}=0$ and the identity $\int_{S}dS_{p}\ {\bf p}\nabla_{p}\cdot(\dot{\bf p}\Psi)=-\int_{S}dS_{p}\ \dot{\bf p}\Psi$ SS2008a yields

[TABLE]

The extra stress in Eq. (14) simplifies to ${\bm{\sigma}}=\alpha{\bf q}{\bf q}+\beta{\bf q}{\bf q}{\bf q}{\bf q}:{\bf E}$, and ${\bf u}$ is solved by

[TABLE]

Equations (26)-(28) are apolar, and also invariant under the transformation ${\bf q}\leftrightarrow-{\bf q}$.

Here we have a few comments:

1. Equations (26)-(28) provide a useful comparison with the Bingham closure discussed below. On this point, we can derive an equation for the evolution of the dipolar tensor ${\bf D}={\bf q}{\bf q}$, using that

[TABLE]

This evolution equation is

[TABLE]

where ${\bf D}^{\nabla}=\frac{\partial{\bf D}}{\partial t}+{\bf u}\cdot\nabla{\bf D}-\left(\nabla{\bf u}\cdot{\bf D}+{\bf D}\cdot\nabla{\bf u}^{T}\right)$ is the upper-convected time derivative. Therefore, it becomes a closed evolution equation when taken together with the momentum balance equation and incompressible constraint:

[TABLE]

2. The sharply aligned dynamics are interesting to consider in the Lagrangian frame of the background flow. As a reminder, this is defined by the flow map ${\bf X}({\bm{\alpha}},t)$ from the initial position ${\bm{\alpha}}$ (i.e., the Lagrangian variable) to the current position ${\bf X}$ through

[TABLE]

where ${\bf X}({\bm{\alpha}},0)={\bm{\alpha}}$. The associated deformation tensor is then defined as ${\bf F}=\partial{\bf X}/\partial{\bm{\alpha}}$ ($F_{ij}=\partial X_{i}/\partial\alpha_{j}$). When defining the field of unit vectors ${\bf r}({\bm{\alpha}},t)=\frac{{\bf F}\cdot{\bf q}_{0}}{\left|{\bf F}\cdot{\bf q}_{0}\right|}$, and using the identity $\frac{D\left|{\bf F}\cdot{\bf q}_{0}\right|}{Dt}=\left|{\bf F}\cdot{\bf q}_{0}\right|\left(\nabla{\bf u}:{\bf r}{\bf r}\right)$, it is then straightforward to show that ${\bf r}$ satisfies Jeffery’s equation:

[TABLE]

Noting that ${\bf F}({\bm{\alpha}},0)\equiv{\bf I}$ and setting ${\bf q}({\bm{\alpha}},0)={\bf q}_{0}({\bm{\alpha}})$, then the uniqueness of solutions to initial value problems gives immediately that, in the Lagrangian frame, the solution to Eq. (26) is

[TABLE]

Equation (35) is simply a normalized version of the Result of Cauchy that ${\bm{\omega}}={\bf F}\cdot{\bm{\omega}}_{0}$ with ${\bm{\omega}}$ the vorticity field evolved by the 3D incompressible Euler equations Childress2009 . Naturally, this also yields a Lagrangian expression for the tensor ${\bf D}={\bf q}{\bf q}$:

[TABLE]

To be consistent with its origins, ${{\bf D}}_{0}$ must be initialized by the dyadic product of a unit vector with itself. These expressions can be used to derive an integro-differential version in the Lagrangian frame akin to the vorticity-based Lagrangian formulation of the incompressible Euler equations Childress2009 .

3. The sharply aligned dynamics, described either by Eqs. (26)-(28) or Eqs. (30)-(32), is particularly useful for analyzing the stability of aligned suspensions. A globally-aligned steady state is given by ${\bf q}\equiv{\bf{\hat{z}}}$ (or ${\bf D}={\bf{\hat{z}}}{\bf{\hat{z}}}$) and ${\bf u}\equiv{\bf 0}$. Following SS2008b by linearizing about this state, and seeking plane-wave solutions with wave-vector ${\bf k}$, yields two growth rates:

[TABLE]

where $H_{1}(\Phi)=\cos^{2}\Phi$ and $H_{2}(\Phi)=\cos^{2}\Phi\left(2\cos^{2}\Phi-1\right)$, with $\Phi\in[0,\pi)$ the angle between the normalized wave-vector ${\bf{\hat{k}}}$ and ${\bf{\hat{z}}}$ (i.e., $\cos\Phi={\bf{\hat{k}}}\cdot{\bf{\hat{z}}}$). For $\alpha<0$, there is exponential growth at all plane-wave angles, independently of $k$, with the exception of $\Phi=\pi/2$ which is neutrally stable. On both branches, the direction of maximal growth is $\Phi=0$ for which the wave-vector is aligned with the suspension. This is in agreement with the more complex expression for motile suspensions, and also seems to be a generic result for active matter theories SS2008b ; giomi2011 ; SS2013 ; GBGBS2015a . Maximal growth rates for any plane-wave vector ${\bf k}$ are again, like the isotropic case, obtained as $k$ tends to zero.

Spatial diffusion can be incorporated in the model phenomenologically by including the term $d_{T}({\bf I}-{\bf q}{\bf q})\cdot\Delta{\bf q}$ on the right-hand-side of Eq. (26) lin1995 . This simply yields the modified growth rates

[TABLE]

II.3.2 Coarse-graining via the Bingham closure

The kinetic theory is high dimensional, having three spatial dimensions and two angles on the unit sphere (in 3D), and thus expensive to simulate. Alternatively, we can consider approximations based on moment-closure schemes to find approximate evolution equations for low-order ${\bf p}$-averaged quantities. Many such schemes have been proposed DE1986 ; Feng98 but we consider here the Bingham closure which is based on a quasi-equilibrium Ansatz and has been previously applied and compared to kinetic theories of passive rod suspensions in simple rheological flows Feng98 (essentially the same as ours but with $\alpha>0$ and externally forced). The zeroth moment with respect to ${\bf p}$ of the Smoluchowski equation yields the concentration equation (16), while integrating the tensor ${\bf p}{\bf p}$ against the Smoluchowski equation yields a dynamical equation for ${\bf D}$:

[TABLE]

The left-hand-side of Eq. (39) is a tensor transport operator with its last term being generated by Jeffery’s equation Jeffery22 . This evolution equation for ${\bf D}$ will be closed if the fourth-moment tensor ${\bf S}$ is expressed in terms of lower order moments, in particular if ${\bf S}={\bf S}[{\bf D}]$. Here we use the Bingham closure Bingham74 ; Chaubal98 which constructs an approximate local orientation distribution function $\Psi_{B}$ given the second-moment tensor ${\bf D}$. The Bingham distribution takes the axisymmetric form:

[TABLE]

where ${\bf T}$ is a traceless symmetric tensor and $Z$ is a normalization constant. For given ${\bf D}$, the tensor ${\bf T}$ is determined by solving the relation

[TABLE]

at each point in space. Computationally, one can take advantage of the fact that the eigenvectors of ${\bf T}$ and ${\bf D}$ are co-aligned and so both can be rotated into diagonal forms in a common frame. Consequently, only the eigenvalues of ${\bf T}$ need to be determined Chaubal98 . Given ${\bf T}$, and hence $\Psi_{B}$, the fourth-moment tensor ${\bf S}$ is then approximated by ${\bf S}_{B}[{\bf D}]=\int_{S}dS_{p}{\Psi_{B}}{\bf p}{\bf p}{\bf p}{\bf p}$. The final equations take the following forms:

[TABLE]

together with

[TABLE]

And the extra stress is expressed as:

[TABLE]

We refer to Eqs. (42)-(45) as the $BQ$-tensor model. We can also show that both Eq. (39) and Eq. (42) conserve $\rm{tr}({\bf D})$.

For the above $BQ-$tensor model, we have several comments:

1. Previous work Chaubal98 ; Feng98 has shown that in potential flows and in the weak flow limit (i.e., the equilibrium state), the Bingham distribution function yields the exact solution of the Smoluchowski equation (with the rotational flux only). The Bingham distribution also arises naturally as describing spatially uniform, nematically-ordered steady states. Such steady states can be obtained by setting $\dot{\bf p}={\bf 0}$ in Eq. (11), leading to a balance between the rotational diffusion and the Maier-Saupe alignment torque (see, for example, Ezhilan et al. ESS2013 and Gao et al. GBGBS2015b ):

[TABLE]

where $\xi=2\zeta/d_{R}$. Consider the unit sphere with angle coordinates $(\theta,\Phi)$ and polar axis along $\hat{\bf z}$, so that ${\bf p}=\sin\Phi\cos\theta\hat{{\bf x}}+\sin\Phi\sin\theta\hat{\bf{y}}+\cos\Phi\hat{\bf{z}}$ with $\theta\in[0,2\pi)$ and $\Phi\in[0,\pi)$. In seeking a solution $\Psi(\Phi)$, we find ${\bf D}$ has the form

[TABLE]

where $A[\Psi]=\pi\int_{0}^{\pi}d\Phi\Psi(3\cos^{2}\Phi-1)\sin\Phi$ (see ESS2013 ). Then Eq. (46) can be integrated to

[TABLE]

where ${\bf{T}}=\delta(\xi){\rm{diag}}\left(-2/3,-2/3,4/3\right)$. In the above, $\delta(\xi)$ is the solution of the equation:

[TABLE]

which captures a bifurcation, with increasing $\xi$, from an isotropic ($\delta(\xi)=0$, $\Psi\equiv 1/4\pi$, and $A[\Psi]=0$) to a nematically aligned state ($\delta(\xi)\neq 0$ and $A\neq 0$).

Thus the Bingham distribution has a microscopic origin with respect to the kinetic theory, and captures the the isotropic-to-nematic transition. However, to our knowledge there is no asymptotic analysis, rigorous or formal, that establishes error estimates between the Bingham closure and the original kinetic theory. We suspect that such a connection could be established in the particular limit of strong alignment torques at large $\xi$.

2. The $BQ$-tensor system is a direct closure of the kinetic theory for an active or passive liquid-crystal polymer suspension, and each of its parameters has a clear origin (though admittedly, the kinetic theory itself is derived under several modeling and separation-of-scale assumptions). While structurally similar to the active and passive nematic $Q$-tensor theories that follow the phenomenological Landau-deGennes (LdG) approach for liquid-crystalline fluids (see, for example, Giomi2015 ), there are significant differences. In the LdG approach, the ordering and relaxation dynamics, and elastic stresses, follow from the LdG free-energy. The $BQ$-tensor theory has different nonlinearities governing the steric ordering, in both the ${\bf D}$ dynamics and the associated stress, especially due to the presence of ${\bf S}_{B}$. While in the LdG theory the spatial diffusion coefficient of ${\bf Q}$ is a Frank elastic coefficient, in the $BQ$-tensor theory for ${\bf D}$ it is the spatial diffusion coefficient, $d_{T}$, of the active particles. In the LdG theory, the transport operator for ${\bf Q}$ is corotational, while in the $BQ$-tensor theory it is, in part, an upper-convected derivative arising from the microscopic modeling. Finally, in the $BQ$-tensor theory there are additional terms, in both dynamics and stresses, that arise from the constraint of particle rigidity.

3. Note that while the original kinetic theory has an energy law (i.e., Eqs. (18,19)), the $BQ$-tensor theory apparently does not. However, Li et al. LWZ2015 have recently developed more complicated Bingham closure models for passive rod suspensions that preserve the energetic structure. In their theory, the configurational entropy is approximated using the Bingham distribution, leading to a theory where the ${\bf T}$ tensor appears directly in the dynamics, rather than only indirectly in the calculation of ${\bf S}_{B}$. This approach leads to a closer concordance of the $BQ$-tensor and LdG $Q$-tensor theories, though we have not pursued that here in the context of active suspensions.

4. In the absence of spatial diffusion, Eq. (42) evolves ${\bf D}$ using only local information. Hence, like the sharply aligned case, this equation could be reformulated naturally to describe the dynamics in the Lagrangian frame.

5. One simple and direct comparison of the kinetic and $BQ$ theories is their linear behavior near steady state. We examine this in order, first for the uniform isotropic steady state, and second for a nematically aligned steady state.

i. Stability near uniform isotropy. We will show that the kinetic and $BQ$ theories give identical results. This is easily seen by considering the evolution of the perturbation ${\bf D}^{\prime}$. Multiplying the linearized Smoluchowski equation (20) by ${\bf p}{\bf p}$, and integrating over the unit ${\bf p}$-sphere gives the evolution equation for ${\bf D}^{\prime}$:

[TABLE]

The above equation is closed since ${\bf E}^{\prime}$ is a linear functional of ${\bf D}^{\prime}$ through the linearized Stokes equations (21). Inverting Eq. (21) in Fourier space for ${\bf u}^{\prime}$ (with ${\bf u}^{\prime}({\bf x})=\tilde{\bf u}({\bf k})\exp(i{\bf k}\cdot{\bf x})$), forming $\tilde{\bf E}^{\prime}[\tilde{\bf D}^{\prime}]$, and substituting into Eq. (20) yields

[TABLE]

where

[TABLE]

and $\hat{\bf k}={\bf k}/k$.

Linearizing the $BQ$-theory Eqs. (42)-(44) gives an identical result, which is mainly due to the fact that the term ${\bf D}\cdot{\bf D}-{\bf S}:{\bf D}$ appears commonly in both the ${\bf D}$ transport and momentum-balance equations (here we suppress the $B$ subscript in the Bingham case). Also linearization of ${\bf S}:{\bf D}$ produces ${\bf S}^{0}:{\bf D}^{\prime}+{\bf S}^{\prime}:{\bf D}^{0}$, and in both models we have ${\bf S}^{0}:{\bf D}^{\prime}=2{\bf D}^{\prime}/15$ and ${\bf S}^{\prime}:{\bf D}^{0}={\bf D}^{\prime}/3$.

To finish the analysis, we employ ${\bf R}$ as a rotation matrix such that $\hat{\bf k}={\bf R}\cdot\hat{\bf z}$, and define $\hat{\bf D}^{\prime}={\bf R}^{T}\cdot\tilde{\bf D}^{\prime}\cdot{\bf R}$. Then Eq. (51) can be rewritten as

[TABLE]

In the rotated frame, the off-diagonal shear components $\hat{D}_{13}^{\prime}$ and $\hat{D}_{23}^{\prime}$ evolve as

[TABLE]

which yields the growth rate, $\sigma=-(\gamma+\omega_{k})$, identical to that in Eq. (24). The remaining components of $\hat{\bf D}^{\prime}$ instead satisfy

[TABLE]

In short, only $\hat{D}_{13}^{\prime}$ and $\hat{D}_{23}^{\prime}$ reflect the contributions of fluid flows, which can be stabilizing or destabilizing. The remainder are decoupled from flow and show decay or growth only due to thermal fluctuations and steric interactions, respectively.

ii. Stability near a nematically aligned homogeneous state. As discussed above, with the inclusion of steric interactions via the Maier-Saupe potential there are additional homogeneous solutions that arise as a balance between the steric alignment torque and the rotational diffusion. These solutions are parametrized by $\xi=2\zeta/d_{R}$, and only the isotropic solutions exist when $\xi$ is below a certain critical value $\xi_{c}$. In 3D there is a transcritical bifurcation at $\xi=\xi_{c}\approx 13.46$ beyond which the isotropic state becomes energetically unstable (in terms of the Maier-Saupe interaction energy). The system will be driven towards a nematic state where particles become increasingly aligned as $\xi\rightarrow\infty$ ESS2013 . In 2D the bifurcation is instead a supercritical pitchfork bifurcation occurring at $\xi_{c}=8$ GBGBS2015b . As discussed in Comment (1.) the form of distribution function for the nematically aligned states is of Bingham type, and hence yields steady states for the $BQ$-tensor theory.

The plane-wave stability of nematically aligned states has been examined in kinetic theories for Pusher and Extensor suspensions ESS2013 , and for microtubule/motor-protein suspensions GBGBS2015a ; GBGBS2015b . These show that activity can drive a bending instability where the wave-vector of maximal growth is aligned with the base orientation of the nematic state, as already suggested by the sharp-alignment analysis (see Eq. (38)). Here we examine the bending instability for Extensor suspensions in 2D for both the kinetic and $BQ$-tensor theories. Solving Eq. (49) for $\delta(\xi)$ gives the Bingham distribution $\Psi_{0}$, the base-state tensors ${\bf T}_{0}$, ${\bf D}_{0}$ and ${\bf S}_{0}$, and the base-state normalization factor $Z_{0}$. Perturbations of these base quantities are introduced as ($\varepsilon\ll 1$)

[TABLE]

The perturbations can be further expressed in terms of ${\bf T}^{\prime}$ as

[TABLE]

with ${{\bf{R}}^{0}}=\frac{1}{Z^{0}}\int{{\bf{pppppp}}\exp\left({{{\bf{T}}^{0}}:{\bf{pp}}^{T}}\right)}d{\bf{p}}$ the sixth-order ${\bf p}$-moment tensor. These relations allow us to express the high-order moments in terms of ${\bf D}^{\prime}$ through the intermediate variable ${\bf T}^{\prime}$. The linearized equations are then transformed to Fourier space where we solve for the plane-wave growth rates $\sigma$ numerically.

Figure 2 shows the growth rates $\sigma$, for both the kinetic and $BQ$-tensor theories, of a plane-wave perturbation to the homogeneous nematic state. Here the plane-wave vector is in the direction of nematic alignment, which has the greatest growth rates. Similar to the isotropic case, the fastest growth rate is obtained as $k\rightarrow 0^{+}$. These features are in agreement with the sharply aligned analysis, and for motile and microtubule suspensions. We show the effect of increasing alignment torque while also comparing the results of the kinetic theory with those of the $BQ$-tensor theory. As $\zeta$ is increased for fixed $\alpha$, we see an increase in the maximal growth rate and the band of unstable modes. We also see an increasingly good correspondence between the full kinetic theory and the $BQ$-theory. In addition, when examining the prediction from the sharply-aligned analysis (black-dotted curve), we find while it gives an excellent accounting for the overall scale and features of the growth rate curves, it is apparently not the asymptotic limit of the $BQ$-tensor model as $\zeta$ tends to infinity.

III Numerical studies and comparisons

In this section we study numerically the dynamics of Extensor suspensions in both periodic and closed geometries. We use simulations in 2D periodic geometries to compare the kinetic model with its coarse-grained $BQ$-tensor model at different levels of activity (i.e., varying $\alpha$). Here we consider two versions of the model (both kinetic and its $BQ$-tensor closure): the full model, termed Model II, and Model I where we omit both steric interactions and constraint stresses (i.e., set $\zeta=\beta=0$). Model I makes a connection with our earlier studies of active suspensions SS2008a ; SS2008b , and emphasizes the fundamental differences arising from steric interactions and nematic ordering. In either case we find that the $BQ$-tensor theory gives an excellent statistical accounting of the complex dynamics, particularly for Model II. We then use the $BQ$-tensor theory to study the dynamics of active suspensions under confinement. Flows in a disk have been studied previously using other closures of Model I to show bifurcations, with increasing activity, from isotropic to auto-circulating flows, and thence to more complex dynamics. We recapitulate this using a $BQ$-tensor version of Model I, and find excellent agreement (see Appendix B for details). We then focus on the full model (i.e., Model II) to study confined suspensions when steric interactions are strong. These steric interactions make a nematically ordered suspension the natural state against which activity competes, and we find both new bifurcations and dynamical phenomena.

III.1 Numerical comparison of the kinetic and $BQ$-tensor theories

In 2D, we evolve $\Psi({\bf x},\theta)$ where $\theta\in[0,2\pi)$ is the particle orientation in the $xy$-plane (i.e., ${\bf p}=(\cos\theta,\sin\theta)$). For the kinetic theory we use Eqs. (1,10,11) to evolve the distribution function $\Psi$, with the requisite background velocity field obtained by solving Eq. (8). The numerical method is pseudo-spectral with 256-512 Fourier modes in each spatial direction, and 64 modes in the orientation $\theta$. Equation (8) is inverted via Fourier transform as

[TABLE]

We use pseudo-spectral collocation to evaluate the nonlinear advection terms. For time-stepping, we use a second-order Adams-Bashforth scheme of the Fourier-transformed distribution function, combined with an integrating factor method for handling the diffusion terms accurately and stably HLS1994 ; GBGBS2015a . The nondimensional diffusion coefficients, $d_{T}$ and $d_{R}$, are chosen between 0.01-0.05. We evolve the $BQ-$tensor system in a similar fashion, again using 256-512 Fourier modes in each spatial direction.

As a first comparison, Fig. 3 shows the late-time evolution of Model I for both the kinetic theory and its $BQ$-tensor closure (again, setting $\beta=\zeta=0$). These simulations are performed from initial data near uniform isotropy. In particular, the initial condition for $\Psi$ has the form

[TABLE]

where the amplitudes $\epsilon_{i}\in[-0.01;0.01]$ and phases $\varphi_{i}\in[0,2\pi)$ are sampled from uniform distributions, with $P_{i}$ a low-order trigonometric polynomial which satisfies $\Psi\left(\theta+\pi\right)=\Psi\left(\theta\right)$. We include only the 10 longest-wavelength spatial modes in Eq. (59). The initial data for the $BQ$ system is simply found by forming ${\bf D}$ from $\Psi({\bf x},\theta,0)$.

In agreement with analytical prediction, we find that the initial concentration fluctuations decay, and the concentration becomes uniform. Figure 3(a) and (c) show the nematic director field ${\bf m}({\bf x},t)$ overlaying the scalar order parameter $s({\bf x},t)$. The system evolves quickly away from the initially isotropic state to a temporally fluctuating state with high local order. Panels (b) and (d) show the fluid velocity field ${\bf u}=(u,v)$ overlaying the planar (scalar) vorticity ($\omega=v_{x}-u_{y}$). The temporal and spatial dynamics for both models are complex and seemingly chaotic, and qualitatively similar to those seen in models of bacterial Pusher suspensions aranson07 ; SS2008b and of biofilament/molecular-motor assemblies TGY2013 ; giomi13 ; GBGBS2015a ; GBGBS2015b .

While both the kinetic and $BQ$-tensor models start with the same data for ${\bf D}$, the system dynamics being chaotic means that the two descriptions diverge in their details. However, both descriptions do evolve towards a similar statistical structure, as shown by the time evolution in the mean nematic order $S(t)=\int_{\Omega}s({\bf x},t)dA/{\rm area}(\Omega)$, from near zero to a temporally fluctuating state with mean near 0.6 (note the $BQ$ theory slightly overestimating the order). The linear theory in Section A predicts that the kinetic and $BQ$-tensor models have the identical linear growth rates from a uniform isotropic base state. Exponential growth to saturation is seen in $S(t)$ following a short equilibration period where high wave-number modes rapidly decay. Fitting $S(t)$ during the growth period yields a growth rate $\sigma\approx 0.16$ for both models. This is close to the growth rate of the first fundamental mode in the simulation box, $\sigma\approx 0.167$, predicted by the linear theory. This figure also suggests that the growth rate of the first mode sets the time scale to saturation.

That the kinetic and $BQ$-tensor theories produce the same magnitude of flow is seen in the very similar scales of their vorticity dynamics. To quantify the flows’ length scales, the inset to Fig. 5 plots the normalized space- and time-averaged velocity-velocity spatial autocorrelation function ${\rm{Cor}}[{\bf u},{\bf u}](R)=\left\langle{\bf u}({\bf x})\cdot{\bf u}({\bf x}+{\bf R})/|{\bf u}({\bf x})|^{2}\right\rangle$ where $R=|{\bf R}|$. We see a close match in this measure between the two models, as well as a slight negative minimum (indicating oppositely directed flow at $R\approx 13$, which is near the box size). This is consistent with the first fundamental spatial mode in the box growing the most rapidly.

In comparing Figs. 3(a) and (c), we find in both models large regions of high nematic order (light to dark red) and director alignment. In the corresponding velocity fields ((b) and (d)), these regions are associated with likewise oriented extensional straining flows (which are necessarily of small vorticity). As found in simulations of active nematics GBGBS2015a ; GBGBS2015b , both models show $\pm 1/2$-order disclination defects in the ${\bf m}$ field, inhabiting regions of small nematic order, and co-located with fluid jets between oppositely signed vortical regions.

In Fig. 4, we show a late-time comparison for Model II using the kinetic theory and its $BQ$-tensor closure. Panels (a) and (c) show clearly that, in both models, the nematic field contains motile disclination defect pairs of order $\pm 1/2$. The induced active flows in panels (b) and (d) look similar to those of Model I, but with stronger vorticity. The defects shown here exist in regions of small nematic order (dark blue), and are born as opposing pairs in elongated “incipient crack” regions GBGBS2015a ; GBGBS2015b . These crack structures locally decrease nematic order, and increase curvature of director field lines. Characteristically, we find that the $+1/2$ defects propagate approximately along their central axis and have a much higher velocity than those of $-1/2$ order. The relatively higher flow velocity in the neighborhood of a $+1/2$ defect appears as a well-localized jet, in the direction of defect motion, between two oppositely signed vortices. These dynamics are similar to those observed previously in studies of polar kinetic models of microtuble/motor-protein suspensions GBGBS2015a ; GBGBS2015b , and other active nematic models TGY2013 ; giomi13 ; Giomi2015 .

Again, as shown by the evolution of the mean nematic order parameter in Fig. 5, the kinetic and $BQ$-tensor models, evolving from the same ${\bf D}$ data, show statistically similar long-time dynamics, and nearly identical growth from isotropy. As expected, we observe a more rapid growth in Model II away from isotropy due to the Maier-Saupe alignment torques; see Eq. (24). The kinetic and $BQ$-tensor models again show nearly the same exponential growth to saturation of $S$, with a fit yielding $\sigma\approx 1.2$, close to the linear growth rate of the first fundamental mode in the box $\sigma\approx 1.08$ (which gives a time scale for growth to saturation). Finally, the velocity-velocity autocorrelation function shows excellent agreement between the kinetic model and its $BQ$-tensor closure (and improved over Model I).

III.2 Active Flows in Confined Geometries

Having established an excellent correspondence between the kinetic and $BQ$-tensor theories, we now study the dynamical behavior of an active suspension confined to an enclosure, using the $BQ$-tensor version of Model II only. Experiments have used motile suspensions of bacteria Lushi2014 , Quincke rollers bricard13 , and MT/motor-protein assemblies SanchezEtAl2012 to study the effects of confinement. Previous work woodhouse2012 ; ezhilin2015 has studied bacterial suspensions using models closely related to the $BQ$-tensor model of Model I. To establish correspondence with these works, in Appendix B we compare the dynamics of Model I in a circular chamber to previous analytical and numerical results by Goldstein and Woodhouse woodhouse2012 . Overall, we find excellent agreement, showing a bifurcation with decreasing $\alpha<0$ from isotropy to an axisymmetric swirling state, further to a symmetry breaking associated with +1/2-order disclination singularities, and thence to complex aperiodic dynamics. Here, within the context of Model II, we focus on the interaction of active extensile stresses and steric ordering with the geometry of the confining chamber.

Boundary Conditions. To solve the confined equations, we must supplement Eqs. (42)-(45) with appropriate boundary conditions. We choose the simplest that conserve total particle number. Consider a bounded 2D fluid domain $\Omega$ with boundary $\partial\Omega$; for example, see the disk in Fig. 6. Here we assume the no-slip condition ${\bf u}|_{\partial\Omega}={\bf 0}$ for the velocity field when solving the Stokes equation. For the ${\bf D}$ advection-diffusion equation we return to the kinetic description in Eqs. (1,10,11), and find that conservation of total particle number is enforced by the boundary condition $({\bf u}\Psi-d_{T}\nabla\Psi)|_{\partial\Omega}\cdot{\bf n}=0$, which reduces to $d_{T}\nabla\Psi|_{\partial\Omega}\cdot{\bf n}=0$. Taking the second moment with respect to ${\bf p}$ then gives the zero-flux boundary condition $({\bf n}\cdot\nabla){\bf D}={\bf 0}$ on $\partial\Omega$. Note that this condition does not enforce any particular alignment direction at the boundary.

Numerical Method. We solve this finite-domain problem using a finite element method with an unstructured triangular grid. The governing equations are discretized using the standard 2D Galerkin formulation:

[TABLE]

where $\tilde{\textbf{u}}$, $\tilde{p}$, and $\tilde{\bf D}$ are test functions. Here we are solving an unsteady Stokes equation with a partial derivative term $\rm{Re}\left(\frac{\partial\bf u}{\partial t}\right)$ with an estimated Reynolds number of $\rm{Re}=10^{-3}$. In the finite element solver, we choose different interpolation functions for different unknown variables in order to satisfy the LBB condition Gao09 ; Hu01 . The fluid velocity is approximated by piecewise quadratic functions which are continuous over $\Omega$ (P2). The pressure and stress are approximated by piecewise linear functions (P1). For this type of mixed finite elements, it is known that when $\rm{Re}$ is small, the numerical solution of the flow field can achieve third order accuracy in space Thomasset81 . For this study, the temporal discretization is a second-order difference scheme. The governing equations are reduced to a nonlinear system of algebraic equations which is solved by a modified Newton-Raphson algorithm Hu01 ; Gao09 . Specifically, in each iteration of the algorithm, the discretized linear system is solved by the GMRES method Saad86 . An incomplete LU preconditioner Chow97 is used to accelerate convergence. In the following, we mainly vary parameters $\alpha$ and $R$ but fix $\zeta=0.5$, $\beta=0.874$, and $d_{T}=d_{R}=0.05$.

III.2.1 Active flows in a circular chamber

To start, in Fig. 6 we consider an Extensor suspension with weak activity, $\alpha=-0.3$. From a nearly isotropic state, the orientation field rapidly orders (note the increase in the scalar order parameter as compared to the cases in Appendix B where $\zeta=\beta=0$ are used in Model I) and two $+1/2$-order disclination defects appear (panels(a,e)), moving about the center (as seen in Fig. 15b), while a vortical field appears concomitantly. However, these two defects spiral outwards to the disk boundary where they are absorbed (panels (b,f)). No new defects are thereafter created and the vortical flow weakens (panels (c,g)) as the nematic orientation field relaxes to a homogeneously aligned state (panels (d,h)), which is allowed by the zero-flux boundary condition for ${\bf D}$.

Increasing activity (i.e., decreasing $\alpha$) leads to persistent dynamics, as shown in the accompanying supplementary videos. In Fig. 7, $\alpha$ has been decreased to -0.7; see video S1 in the supplementary information. Initially it follows the course of the previous case: the system moves quickly away from isotropy with the appearance of two $+1/2$-order disclination defects and a vortical flow. The defects spiral outwards, are absorbed by the walls, and the system relaxes into a state of nearly homogeneous nematic alignment. However, this state is unstable to hydrodynamic instabilities driven by activity GBGBS2015a . As seen in video S1, two cracks of low nematic order appear, with each birthing a pair of disclination defects of opposing sign (order $\pm 1/2$; panels (a,e)). The $-1/2$-order defects are now absorbed by the wall (panels (b,f)), followed shortly thereafter by the $+1/2$-order defects (panels (c,g). The process begins anew with the bending of nematic field lines and the formation of two low nematic order cracks (panels (d,h)). Accompanying the nematic structure variation, we observe the persistence of background rotational flows which strengthen and weaken through the various stages of this periodic dynamics (lower row of panels).

The dynamics becomes quasi-periodic or chaotic at yet higher activity. For $\alpha=-2$, Fig. 8 shows a short period in time that illustrates the complex system dynamics (see also supplementary video S2). In panels(a,d) we see that as $\alpha$ further decreases, there is still a basic dynamics of two rotating $+1/2$-order defects. The rotational flow is more finely scaled as there is now a secondary vortex of the opposite sign to the primary (red) counter-clockwise vortex. The dipolar form of the vortical distribution is associated with a strong fluid jet, marked as the white arrow in panel (a). The primary vortex is more compact, leading to a somewhat smaller circulatory flow. Shortly thereafter (panels (b,e)), a low-order crack appears near one of the defects and births a defect pair. The newly formed $-1/2$-order defect merges with and annihilates the older $+1/2$-order defect nearby (or sometimes both defects are absorbed into the wall), and the process begins anew; see panels (c,f).

Finally, we refer the reader to the supplementary video S3 for $\alpha=-5$, which shows the yet more complex interplays between defects, cracks, and the defect pairs thereby produced. The associated vortical flows are no longer dominated by a single-signed vortex but is much more dipolar, with a fluctuating predominance of one signed vortex over the other. While complicated both spatially and temporally, the dynamics nonetheless maintains a quasi-periodicity.

Figure 9 shows the evolution of the mean nematic order parameter $S(t)$ for these various simulations. In all cases the rapid growth of $S(t)$ away from zero reflects the transition from isotropy to a state of higher nematic order. For the simulation with the lowest activity, $\alpha=-0.3$ (black curve), this initial rise yields a transitory plateau slightly above $S=0.7$ during which the dynamics is characterized by formation and outward spiraling of two +1/2-order defects. Their absorption by the chamber walls, when $40<t<60$ (see Fig. 6), leads to relaxation to a homogeneously ordered state with $S$ slightly less than 0.9. The emergence of persistent and complex defect dynamics in the other cases is reflected in the saturated behaviors of $S(t)$ following growth from isotropy. In all cases the initial growth from zero can be well-fitted by an exponential in time. The calculated growth rates are plotted in the inset, and show a linear increase with increasing activity $\alpha$.

We have also examined the dynamics in circular chambers of different radii and found generally similar behaviors. That being said, in a smaller confinement we have also discovered a fascinating periodic coherent structure. Figure 10 shows the dynamics of the system for $\alpha=-2$ (cf. Fig. 8) but with confinement radius $R=1.25$, as it enters into a periodic dynamics. The defect dynamics starts from the apparent splitting of a $+1$-order defect into two $+1/2$ defects (panel (a); see supplementary video S4) co-rotating with the associated nearly axisymmetric vortical field (panel (d)). The two defects then destabilize from a common rotation centered about the origin (panels (b,e)), and begin to, apparently, move into and out of the wall (i.e., the wall both absorbs and produces defects). The apparent production of defects at confinement boundaries has been observed in experiments involving microtubule/motor-protein assemblies (private communication with Z. Dogic). Eventually, the system relaxes to having one defect seemingly trapped outside of the physical flow domain, leaving only a vestigial low-order crack which moves clockwise around the boundary in the direction opposite to the counter-clockwise circulation of the flow (panels (c,f)).

III.2.2 Active flows in a biconcave chamber

The coarse-grained $BQ$-tensor model and its finite-element discretization greatly facilitates the exploration of active nematic flows in more complex geometries. In Fig. 11, we study the dynamics of an Extensor suspension confined to a biconcave chamber where two circular chambers, each or radius $R=2$, are connected smoothly by a bridge of width $d$. Here we fix $\alpha=-2$. Supplementary videos S5-S8 show the dynamics as $d$ is successively increased ($d=$0.2, 0.8, 1.2, and 2). We find that when the neck is thin, the dynamics in the two chambers are very similar to the case shown in Fig. 7 for $\alpha=-2$, and appear to be evolving independently of each other; see supplementary video S5. This lack of interchamber communication is consistent with the velocity stagnation zone in the bridge (panel (b)). While a thin neck would be sufficient to prevent interplay between the two chambers, this effect is likely reinforced by the alignment of the nematic field lines across the neck which lends elastic rigidity to the material contained therein.

The essentially independent dynamics of the two chambers persists to quite wide necks, e.g. $d=0.8$. As seen in supplementary video S6, the nematic field lines still span the neck, apparently blocking any substantial breaching flows. As $d$ is further increased, beyond $d_{crit}\approx 1.0$, the dynamics in the two chambers begins to couple. Figure 7(a) and supplementary video S7 for $d=1.2$ show significant bending and shearing of the nematic field lines, with the panel(a) inset showing that these deformations are associated with a vortex bounded above and below by two fluid jets that penetrate from one side to the other. For the largest neck width, $d=2.0$, Figure 7(b) and supplementary video S8 show that flows now move freely between the two chambers. This transition is reminiscent of a Fredricks transition where an electric or magnetic field needs to exceed a critical intensity to deform a uniformly aligned nematic liquid crystal material spanning between two plates. Fredericks transitions have been studied in other contexts in the field of active matter giomi08 ; Voituriez2005 .

IV Discussion

We explored a multi-scale model, and its simplifications, for suspensions of rod-like particles that exert active dipolar extensile stresses on the immersing solvent. This model is relevant to suspensions or bundles of microtubules that undergo polarity sorting driven by crosslinking motor proteins surrey01 ; bendix08 ; hentrich10 ; gordon12 , and directly models suspensions of chemically-active particles whose consumption of a chemical fuel creates extensile flows along the particle Wykes16 ; Jewell2016 ; Pandey2016 . Aside from analysis of the model, a central goal here was to investigate the Bingham closure Bingham74 ; Chaubal98 ; Feng98 as applied to active suspensions. We found that the resultant $BQ$-tensor theory gave an excellent accounting of the chaotic self-driven flows of active suspensions (both with and without strong alignment interactions), and allowed us to investigate with relative ease the behavior of an active-nematic fluid under confinement. We found several fascinating behaviors, such as defects being absorbed at (and perhaps birthed by) boundaries, complex rotational flows, and an apparent Fredericks transition associated with nematic elasticity.

We remark that the $BQ$-tensor theory has no unfixed constants with respect to the full kinetic theory, and is far cheaper to simulate (essentially by a factor inversely proportional to the number of points of the sphere (3D) or circle (2D) used to resolve the distribution function in orientation). This also establishes a connection between Doi-Onsager kinetic theories and $Q$-tensor theories such as those based on the Landau-deGennes approach (e.g. giomi13 ; Giomi2015 ; LWZ2015 ). We are currently working on improving such closure approaches, from a variational perspective, and extending them to polar kinetic theories such as for motile suspensions SS2008a ; subramanian09 and microtubule/motor-protein suspensions GBGBS2015a ; GBGBS2015b .

IV.1 Length scale determination

We close with a short digression into the determination of length scales of instability in active suspensions, as a suspension of Extensors is a particularly simple case to discuss. As previously discussed, it seems generic that bulk models of active suspensions lack an intrinsic length scale determined by a most unstable mode simha02 ; SS2008a ; GBGBS2015a . This is seen most directly in Eq. (24) for the growth rates of plane-wave perturbations from isotropy, and in our stability analyses of aligned states. It is also consistent with our simulations in periodic domains where observed growth rates from isotropy are well-described by that for the first fundamental mode with $k=2\pi/L$ of the periodic domain. Nonetheless, some studies have sought to identify an intrinsic length scale of instability thampi14b . Such a scale is present when the system is subject to an external drag or damping, as when an active material is confined to a thin layer between two fluids (e.g., Gao et al. GBGBS2015a in modeling the experiments of Sanchez et al. SanchezEtAl2012 , as well as Leoni and Liverpool leoni10 ), or sits atop a substrate Thampi2014c .

To be specific, consider a thin active layer evolving in the $xy$-plane at $z=0$, across which layer activity gives a tangential stress jump between two outer 3D Stokesian fluids. By taking a 2D $xy$ Fourier transform of the 3D Stokes equations, the two half-space problems for the 3D Stokes flows can be solved analytically under the conditions of zero velocity as $|z|\rightarrow\infty$ and continuity of velocity at $z=0$. This determines the interfacial velocity field ${\bf u}(x,y)$ as a function of the active layer extra stress ${\bm{\sigma}}$ MS2014 ; GBGBS2015a ; GBGBS2015b :

[TABLE]

where ${\bf{\hat{k}}}={{\bf{k}}}/{{\left|{\bf{k}}\right|}}$ is the 2D unit wavevector. While not strictly necessary, we have assumed that the interfacial flow is incompressible (as seems consistent with the experiments of Sanchez et al. SanchezEtAl2012 , private communication with Z. Dogic). Thus, activity at the interface drives flows in the two outer fluids, which in turn provides a drag on the active surface. Simulating this situation amounts to replacing the inversion of the 2D Stokes operator in Eq. (58) with Eq. (63).

The inversion formula in Eq. (63) differs by a factor of $k/2$ from that in Eq. (58) where the Stokes equation is forced by a bulk stress. The propagation of that difference through the stability analysis is found by examining the coefficients $\gamma$ and $\omega_{k}$ (Eq. (52)) that govern the growth/decay rate $\sigma_{k}=-(\gamma+\omega_{k})$ (Eq. (54)). For the 2D (or 3D) bulk fluid case, $\gamma$ has the form $\gamma=\tilde{\alpha}/(1+\tilde{\beta})$; while for the immersed interface case considered here, it is modified to $\gamma_{k}=(\tilde{\alpha}k/2)/(1+(\tilde{\beta}k/2))$ which rises monotonically in $k$ from zero to $\tilde{\alpha}/\tilde{\beta}$. For ${\tilde{\alpha}}<0$, competition of activity with translational diffusion in $\omega_{k}$ yields a unique critical wavenumber $k_{c}>0$ at which $\sigma_{k}$ is maximized.

This same effect manifests itself in the instability of a homogenous nematically-aligned state, and yields an intrinsic scale that is observed in nonlinear simulations (see Gao et al. GBGBS2015a for the more complicated case of microtubule/motor-protein suspensions). For a simple Extensor suspension, Fig. 13 shows the numerically calculated growth rates $\sigma_{k}$ (red curve) of a plane-wave perturbation to the homogeneous nematic state, with wave-vector taken in the direction of maximal growth, $\theta=0$. This being the direction of maximal growth is in agreement with the sharply-aligned analysis of Sect. II.3.1, and has been found in other studies of active suspensions SS2008a ; SS2008b ; ESS2013 ; GBGBS2015a ; GBGBS2015b . The computed curve also shows an intrinsic wave-number of maximal growth, as was found analytically for the isotropic state. As $\zeta$ is increased for fixed $\alpha$, we see an increase in maximal growth rate and decrease in its associated length scale. In that limit we find an increasing correspondence between the full kinetic theory and the $BQ$-theory.

The parabolic-like structure of the growth rate curve is in qualitative agreement with our phenomenological sharply-aligned model when modified to the immersed interface case. This modification alters the plane-wave growth rates $\tilde{\sigma}_{1,2}$ in Eqs. (38)(where $\beta=0$) to the form

[TABLE]

For $\alpha<0$, this expression yields at each $\Phi$ (with $H$ maximal at $\Phi=0$) a maximal growth rate $\sigma_{max}=\alpha^{2}H^{2}(\Phi)/16d_{T}$, occurring at $k_{c}=-\alpha H(\Phi)/4d_{T}$. To emphasize an important point: the quadratic increase in maximal growth rate with activity follows from $k_{c}$ increasing linearly with activity (i.e., decreasing negative $\alpha$). For the bulk case, where the system size sets the most unstable scale, maximal growth rates should scale linearly with activity (as in Eqs. (24), (38), and the inset to Fig. 9).

Acknowledgements. This work was funded by NSF grants DMR-0820341 (NYU MRSEC: TG, MS), DMS-1463962 and DMS-1620331 (MS), DMS 1619960 (TG, MDB, MS), DMR-0847685 (MDB), DMR-1551095 (MDB), DOE grant DE-FG02-88ER25053 (TG, MS), NIH grant K25GM110486 (MDB), R01 GM104976 (MDB, MS).

Appendix A Stability of 3D uniform isotropic suspension

While the linear theory of the Extensor suspension is a special case of that for a motile one, it does have an especially simple and evocative structure that is obscured in the motile case. For an isotropic steady state we have $\Psi^{0}=1/4\pi$, ${\bf D}^{0}={\bf I}/3$, ${\bf u}^{0}={\bf 0}$, and $S^{0}_{ijkl}=\frac{1}{15}\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}+\delta_{ij}\delta_{kl}\right)$. By perturbing this steady state as ${\bf D}={\bf I}/3+\varepsilon{\bf{D^{\prime}}},\;{\bf{S}}={{\bf{S}}^{0}}+\varepsilon{\bf{S^{\prime}}},\;{\bf{u}}=\varepsilon{\bf{u^{\prime}}}$, the 3D linearized Smoluchowski equation (1) becomes SS2008a ; HS2010 ; ESS2013

[TABLE]

This result uses that, for any ${\bf p}$-independent tensor ${\bf A}$, $\nabla_{p}\cdot\left(({\bf I}-{\bf p}{\bf p})\cdot{\bf A}\cdot{\bf p}\right)=\mbox{tr}({\bf A})-3{\bf p}{\bf p}:{\bf A}$, and that $\mbox{tr}({\bf D}^{\prime})=0$ for a constant density flow. Multiplying Eq. (65) by ${\bf p}{\bf p}$, and integrating over the unit ${\bf p}$-sphere gives the evolution equation for ${\bf D}^{\prime}$:

[TABLE]

This result uses that ${\bf S}_{0}:{\bf A}=2{\bf A}/15$ for any symmetric trace-free tensor ${\bf A}$. An identical result is obtained through linearization of Eq. (39). The momentum balance equation and incompressible condition, Eq. (12)-(14), linearize to

[TABLE]

Rewriting Eq. (67) in terms of the stream function allows significant simplification. By assuming that the fluid domain is simply connected, we can define the vector stream function ${\bm{\Phi}}$ such that ${\bf u}={\nabla\times}{\bm{\Phi}}$ and $\Delta{\bm{\Phi}}=-{\nabla\times}{\nabla\times}{\bm{\Phi}}=-{\bm{\omega}}$. Taking the curl of Eq. (67) then yields

[TABLE]

where $\tilde{\alpha}=(\alpha-\frac{2}{5}\beta\zeta)/(1+\beta/15)$. Through the definition of ${\bf D}^{\prime}$ we can express the right-hand-side of this equation, component-wise, in terms of $\Psi^{\prime}$:

[TABLE]

This defines the vector operator ${\bf L}$, given component-wise by $L_{i}=\varepsilon_{ijk}p_{k}p_{l}{\partial}^{2}/{\partial x_{j}}{\partial x_{l}}$. Equation (67) then becomes

[TABLE]

Now we turn to the evolution equation (65). The contracted term ${\bf p}{\bf p}:\nabla{\bf u}$ can be rewritten as ${\bf p}{\bf p}:\nabla{\bf u}=-{\bf L}\cdot{\bm{\Phi}}$ (this uses that ${\bf p}{\bf p}:{\bf E}={\bf p}{\bf p}:\nabla{\bf u}$, and that $\varepsilon_{ijk}=-\varepsilon_{kji}$, followed by an index relabeling). Then, Eq. (65) can be rewritten as

[TABLE]

Now we take a time derivative of Eq. (71) (i.e. $\nabla^{4}{\bm{\Phi}}_{t}=\tilde{\alpha}\int_{S}dS_{p}\ {\bf L}\Psi^{\prime}_{t}$), and use Eqs. (72) and (71) to find

[TABLE]

The terms on the right-hand-side of Eq. (73) can be calculated as

[TABLE]

Finally, by defining ${\bf g}=\nabla^{4}{\bm{\Phi}}$ and using the definition for ${\tilde{\alpha}}$ the linearized dynamics reduces to

[TABLE]

where $C_{1}(z)=(3/5)/(1+z/15))$ and $C_{2}(z)=(4/5)+(6z/25)/(1+z/15)$. Thus, growth or decay rates of plane-wave solutions are given by

[TABLE]

Appendix B The dynamics of Model I in a circular domain

We confine an Extensor suspension to a circular chamber of radius $R$, starting the simulation near uniform isotropy in the weakly aligned limit with $\zeta=\beta=0$, which corresponds to neglecting steric interactions between filaments. This limit corresponds to the model previously studied by Woodhouse and Goldstein woodhouse2012 , who discovered a symmetry-breaking bifurcation at which the stationary state transitions to a stable circulating state. For higher activity, another bifurcation produces an oscillatory state with a pair of mobile defects.

We expected that our Bingham closure would reproduce the Woodhouse and Goldstein results that used the Hinch and Leal closure hinch1976 when neglecting steric interactions. Indeed, we recover an identical system of linearized equations for perturbation about an isotropic state. Numerically, we capture the bifurcation to a circulatory state (Fig. 14), and the growth rates agree quantitatively with linear stability analysis (Fig. 16). The shape of the flow profile matches that determined analytically woodhouse2012 (Fig. 14).

For increasing activity (i.e., decreasing $\alpha$), in Fig. 15 we reproduce the oscillatory two-defect state found by Woodhouse and Goldstein woodhouse2012 . For sufficiently high activity, this oscillatory state is unstable to the production of cracks and defects in the nematic field, and multiple vortices.