Analytical structure, dynamics, and coarse-graining of a kinetic model of an active fluid
Tong Gao, Meredith D. Betterton, An-Sheng Jhang, Michael J. Shelley
AAnalytical structure, dynamics, and coarse-graining of a kinetic model of anactive fluid
Tong Gao
Department of Mechanical EngineeringDepartment of Computational Mathematics,Science and Engineering,Michigan State University
Meredith D. Betterton
Department of Physics, University of Colorado at Boulder
An-Sheng Jhang
Courant Institute of Mathematical Sciences, New York University
Michael J. Shelley
Courant Institute of Mathematical Sciences, New York UniversityCenter for Computational Biology, Flatiron Institute, New York (Dated: March 6, 2017)We analyze one of the simplest active suspensions with complex dynamics: a suspensionof immotile “Extensor” particles that exert active extensile dipolar stresses on the fluid inwhich they are immersed. This is relevant to several experimental systems, such as recentlystudied tripartite rods that create extensile flows by consuming a chemical fuel. We firstdescribe 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 hydrody-namic instabilities, as well as the steric interactions of rod-like particles that lead to nematicalignment. This active nematic system yields complex flows and disclination defect dynam-ics very similar to phenomenological Landau-deGennes Q -tensor theories for active nematicfluids, as well as by more complex Doi-Onsager theories for polar microtubule/motor-proteinsystems. We apply the quasi-equilibrium Bingham closure, used to study suspensions of pas-sive microscopic rods, to develop a non-standard Q -tensor theory. We demonstrate throughsimulation that this “ BQ -tensor” theory gives an excellent analytical and statistical account-ing of the suspension’s complex dynamics, at a far reduced computational cost. Finally, weapply 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 weaknematic alignment, but for strong alignment find novel dynamics with activity-controlleddefect production and absorption at the boundaries of the domain. In bi-concave domains, aFredericks-like transition occurs as the width of the neck connecting the two disks is varied. I. INTRODUCTION
Active suspensions are non-equilibrium materials composed of suspended particles whose activ-ity, driven by consumption of a local fuel, lead to particle motions or local induced flows [1, 2].Examples include bacterial swarms [3, 4], collections of synthetic colloidal particles [5, 6], andmixtures of cytoskeletal filaments driven by molecular motors [7, 8]. One origin of the instabilitiesand complex dynamics in active suspension is the stresses created in the surrounding solvent bythe particles’ activity. The solvent is a coupling medium for multiscale dynamics: particle in-teractions through the solvent can manifest at the system scale as collective motion, which feeds a r X i v : . [ c ond - m a t . s o f t ] M a r back to alter their interactions. Theoretical and numerical studies have investigated various as-pects of active suspensions at different scales, from single particle dynamics to suspension rheology.While discrete-particle simulations that incorporate long-ranged hydrodynamic effects can captureobserved large-scale features of these systems [9–13], constructing accurate continuum models isessential to their characterization and analysis, and for making new predictions of macroscalebehaviors [2].Continuum descriptions of active suspensions of self-propelled particles, such as bacteria, havebeen based on a variety of approaches [14–20]. In all, coarse-grained variables, such as particleconcentration, concentration-dependent steric interactions, order parameters, and solvent stressesand 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 ofsuspensions of active rod-like particles [17, 21–25]. In this approach, a Smoluchowski equation de-scribes the evolution of a particle distribution function, coupled to a coarse-grained Stokes equationdriven by the extra stresses created by particle activity and interactions. The fluxes and stressesof the model are derived from modeling how such particles interact with mean-field flows andother coarse-grained variables. For suspensions of self-propelled rods which interact only hydrody-namically [17, 21] this approach identifies characteristic aspects of the dynamics: either uniformisotropic or globally aligned suspensions have hydrodynamic instabilities that depend upon theself-propulsion mechanism, particle concentration, and system size [20–22], leading to the emer-gence of complex, turbulent-like dynamics. When such particles also interact sterically throughconcentration-dependent alignment torques, this generates a polar “active-nematic” model with anisotropic-to-nematic ordering transition. The new nematically-aligned steady states can again beunstable to hydrodynamic instabilities, again leading to complex dynamics but modulated by aneffective liquid crystal elasticity [23]. Similar but more elaborate kinetic models have been derivedfor “bio-active” suspensions [7] composed of microtubules whose relative motions are driven byprocessive cross-linking molecular motors such as kinesin-1 [24, 25]. There, microtubule motionand motor activity create polarity-dependent fluxes and stresses that drive hydrodynamic insta-bility, and a complex dynamics characterized by nematic director fields with motile disclinationdefects that are continuously nucleated and annihilated.While having the advantage of being founded upon microscopic modeling, such kinetic theoriescan carry the price of complexity. For Doi-Onsager theories as above, the distribution functiondepends upon both particle position x and orientation vector p ( | p | = 1) and so, in three dimen-sions, 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 p variable, yielding more approximate models that describe the evolution of lower-order p -momentsof the distribution function, such as the concentration φ ( x , t ) (zeroth-moment), the polarity vector q ( x , t ) (vector of first-moments), and the unnormalized tensor order parameter D ( x , t ) (tensorof second-moments) [2, 26–28]. Tensorial moment-closure models can resemble, with importantdifferences, tensor models based upon phenomenological Landau-deGennes Q -tensor liquid crystaltheories [29–31].In this paper, we use a Doi-Onsager theory to study what is perhaps the simplest active sus-pension: a collection of non-motile, yet mobile, rod-like particles that exert active dipolar stresseson the solvent and interact with each other through hydrodynamic coupling and steric alignmenttorques. When the stresses are extensile, we call these particles “Extensors”. While this is a specialcase of the systems above, it is well-worth studying in its own right. Firstly, it presents all thebasic transitions and instabilities associated with motile suspensions but in a more revealing andsimplified form. Secondly, there are several types of active suspensions described at least qual-itatively by this theory, including particles that elongate through stretching or growth, creatingthe extensile flows associated with destabilizing dipolar stresses. Such particle dynamics occurs insome instances of bacterial cell division [32], as well as in intriguing liquid-crystal phase transitions[33, 34]. Such stretching is also thought to occur for bundles of microtubules (formed by a depletioninteraction) as they are “polarity-sorted” by kinesin-1 motor-proteins [7, 35]. A very different sys-tem is composed of tripartite Au-Pt-Au nanomotors [36, 37] for which catalytic surface reactionswith a surrounding aqueous hydrogen peroxide solution [38] create extensile surface flows. Theself-assembly of motile aggregates by small numbers of such nonmotile particles has been studiedexperimentally by Jewell et al. [37] and Wykes et al. [36], and numerically by Pandey et al. [39].Thirdly, we use this simple kinetic theory as a testing ground for constructing coarse-grainedmacro models for active suspensions using the Bingham moment closure [40], a closure schemeused successfully in the study of passive rod suspensions [41, 42]. The Bingham closure expressesfourth-moment tensors arising in the kinetic theory in terms of D ( x , t ) using the so-called Binghamdistribution, which is a quasi-equilibrium Ansatz . The Doi-Onsager system is then expressed solelyin terms of D . We call this BQ -tensor theory as D is an unnormalized form of the tensor orderparameter Q , itself the central evolved quantity in the Landau-deGennes theory of liquid crystaldynamics, which has been applied extensively to so-called active nematics [24, 25, 29–31]. Weshow that the BQ -tensor theory closely or exactly reproduces the instabilities of the kinetic theorynear 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 thekinetic theory model, both with and without steric interactions, of an active nematic suspensiondriven by extensile stresses.Finally, having established the fidelity of the BQ -tensor theory, we use it to examine the dynam-ics of active suspensions under confinement. This is inspired by several recent experiments showingthat confined collective suspensions of self-propelled particles can organize into auto-circulatingstates [43, 44], and develop emergent ordering and density shocks [6]. The dynamics of activesuspensions has been studied previously using both discrete particle simulations and continuummodels in straight channels [27, 45], circular disks [26, 46], annuli, and racetracks [28]. For sus-pensions in a circular chamber we find perfect agreement with the previous results for suspensionsthat interact only hydrodynamically [26], and focus instead on the active nematic case where stericinteractions 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 atthe boundaries. In biconcave geometries, which are essentially two disks smoothly connected by abridge, we find evidence for a Friedricks transition in which the elasticity of the ordered fluid pre-vents 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 rodsthat produce extensile surface flows. Except for the sign of the active dipolar stress, this theoryis nearly identical to the classical Doi-Onsager models for passive liquid crystal polymers wherethe suspended molecules have fore-aft symmetry [42, 47]. We discuss its analytical and stabilityproperties, as well as introduce simplified models, particularly coarse graining through the Binghamclosure. In Sect. III, we study numerically the dynamics of Extensor suspensions. In periodic opendomains, 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, flowpatterns, and statistical structure. We then simulate the dynamics of the BQ -tensor model inboth circular and biconcave domains using a finite element method. Following a discussion inSect. IV, we present some additional results in appendices. (a)(b) (c) Au Au Pt b l s u p u p p x FIG. 1: (a) A schematic of an Extensor particle that produces extensile surface flows which yield extensiledipolar stresses. Related examples: (b) bacterial filamentation; (c) Au-Pt-Au nanomotor; (c) bundle of MTsmixed with motile crosslinking motors.
II. DOI-ONSAGER THEORY FOR A SIMPLE ACTIVE SUSPENSIONA. A micro-mechanical model
Consider a suspension of N rigid rod-like active particles in a volume V , each of length l anddiameter b . The rods are considered to be slender so that the aspect ratio r = b/l <<
1. Fordefiniteness we assume a simple form of activity: each rod produces a symmetric surface flow v ( s ) = sgn ( s ) u p with − l/ ≤ s ≤ l/ p its unitorientation vector, and u the signed surface flow speed; see the schematic in Fig. 1. Because ofsurface flow symmetry such active rods are not motile. If u > u <
0. These flows are associated with dipolarextra stresses. Directly following the derivation for Pusher suspensions [48], assuming that activeparticles are small relative to the large-scale suspension flows they create in ensemble , and usingslender 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 activesurface flow is σ a pp where σ a = − ηl u / η = 8 πµ/ | ln( er ) | with µ the fluid viscosity. Note that σ a < σ a > σ a < B. Doi-Onsager kinetic theory
As an expression of the conservation of particle number, one can derive a Smoluchowski equationfor the (dimensional) distribution function ψ ( x , p , t ) for number density of particles at position x with orientation p : ∂ψ∂t + ∇ · ( ˙ x ψ ) + ∇ p · ( ˙ p ψ ) = 0 . (1)The distribution function is chosen so that (cid:82) V dV (cid:82) S dS p ψ = N . Here and for the remainderof this paper, derivative operators without subscripts (e.g., ∇ ) denote spatial derivatives. Theorientational gradient operator on the unit sphere is ∇ p = ( I − pp ) · ∂/∂ p . Moments of ψ withrespect to p arise naturally in the theory: the particle concentration φ , (unnormalized) polarity field q , second-moment tensor D (unnormalized tensor-order parameter), and fourth-moment tensor S are given by φ ( x , t ) = (cid:90) S dS p ψ ( x , p , t ) , (2) q ( x , t ) = (cid:90) S dS p p ψ ( x , p , t ) (3) D ( x , t ) = (cid:90) S dS p pp ψ ( x , p , t ) , (4) S ( x , t ) = (cid:90) S dS p pppp ψ ( x , p , t ) . (5)The trace-free (normalized) tensor order parameter, the so-called Q -tensor, is defined by Q ( x , t ) = φ ( x , t ) − D − I /d , with d = 2 , Q has a maximal nonnegativeeigenvalue λ max satisfying 0 ≤ λ max ≤ ( d − /d . Assuming that λ max is isolated, then we callits associated unit eigenvector the director m , and 0 ≤ s = λ max d/ ( d − ≤ x and ˙ p arise from thedynamics of the particle’s center-of-mass position x and orientation p , calculated using slender-body theory for the Stokes equations [49], and Maier-Saupe theory for steric alignment [50] ofrod-like particles. This yields˙ x = u − D T ∇ ln ψ, (6)˙ p = ( I − pp ) · ( ∇ u + 2 ζ D ) · p − D R ∇ p ln ψ, (7)where u is the background fluid velocity field. Equation (6) states that the rod’s positions movewith the background fluid flow and undergo translational diffusion (with coefficient D T ), whileEq. (7) extends Jeffrey’s equation for the rotation of rods by the fluid velocity gradient to includea local alignment torque arising from a rod-rod interaction potential with strength ζ . The lastterm of Eq. (7) models rotational diffusion with coefficient D R .The equations are closed by relating the incompressible background velocity u and pressureΠ to the extra-stress tensor σ produced by the particles’ presence [51], through a forced Stokesequation ∇ Π − µ ∆ u = ∇ · σ , (8) ∇ · u = 0 . (9)The stress tensor σ has three contributions: σ = σ a + σ c + σ s . The dipolar active stress tensor σ a = σ a D arises from the active surface flows generated by particles. For Extensor particles, σ a is negative, and positive for “Contractors” (see Fig. 1). The so-called constraint stress tensor σ c = σ c S : E , where E = 1 / ∇ u + ∇ u T ) is the symmetric rate-of-strain tensor and σ c = ηl / σ s = − σ s ( D · D − S : D ) with σ s = πηl ζ / ν = nbl where n = N/V is the meannumber density [23, 47]. To nondimensionalize, we choose a characteristic length scale l c = b/ν ,velocity scale | u | , and stress scale µ | u | /l c . The distribution function is normalized such that(1 /V ) (cid:82) Ω dV (cid:82) S dS p Ψ( x , p , t ) = 1. Equation (1) retains its form, where the translational androtational fluxes are now ˙ x = u − d T ∇ ln Ψ , (10)˙ p = ( I − pp ) · ( ∇ u + 2 ζ D ) · p − d R ∇ p ln Ψ . (11)The Stokes equation becomes ∇ Π − ∆ u = ∇ · σ , (12) ∇ · u = 0 , (13)where the dimensionless stress tensor is given by σ = α D + β S : E − ζβ ( D · D − S : D ) . (14)The dimensionless control parameters are α = σ a µ | u | l , ζ = ζ | u | l , β = πrν r ) , d T = νb | u | D T , d R = bν | u | D R . (15)Given the normalization of Ψ, the average density is ¯ φ = 1, and the state of uniform isotropyis given by Ψ ≡ Ψ = 1 / π in 3D, and 1 / π in 2D. In the case of spatially uniform density, D = Q + I /d . Note that α inherits the sign of σ a , and that the nondimensional translationaland rotational diffusion coefficients are proportional and inversely proportional, respectively, tothe effective volume fraction ν . We will sometimes refer to α as the activity and treat it as afree parameter though, in fact, it is a geometric parameter associated with particle shape and therelative placement of active and passive stresses upon it [17, 22, 24]. For the specific case of thetripartite rods considered by Wykes et al. [36], it is estimated that α ≈ −
1. Alternatively, if u isinterpreted as the velocity scale for a base level of activity, then increases in α (either positively ornegatively) can be interpreted directly as multiplicative increases in surface velocity over this baselevel.Rotational thermal fluctuations of rod-like particles also give rise to an extra stress of form K D with K > α of either sign, is a much larger contribution.A very similar model describes suspensions of stretching filaments where the distribution func-tion now depends upon particle length l ( t ). Assuming that particle length growth, ˙ l = f ( l ), isindependent of particle position and orientation, one can derive a closed evolution equation for themarginal distribution by integrating out l . The main differences are that u is replaced by f ( l ),and l dependencies are replaced by distributional averages. As expected, ˙ l >
1. Concentration Fluctuations.
Unlike suspensions of Pusher swimmers, which also produceextensile flows and where fluctuations can grow through nonlinear effects [21, 48], concentrationfluctuations in the nonmotile case decay to zero. Denote the fluid domain as Ω, and assume either (i)
Ω is a cubic domain [ − L/ , L/ over which Ψ is periodic, or (ii) there is neither mass norparticle flux across the boundary ∂ Ω. Integrating Eq. (1) over the unit p -sphere, and using Eq. (6)and velocity incompressibility, then gives ∂φ∂t + u · ∇ φ = d T ∆ φ. (16)Integrating Eq. (16) over Ω, integrating by parts and applying boundary conditions (periodic orzero flux) yields ddt (cid:90) Ω dV φ = − d T (cid:90) Ω dV |∇ φ | . (17)Hence, concentration fluctuations from equilibrium will decay to zero, regardless of particle type.
2. An Energy Identity.
Like the motile case [21], the kinetic theory has an “energy” identitygoverning the evolution of a total energy composed of the configurational entropy and the steric-interaction energy [52]. Let E = S + κ P = (cid:90) Ω dV (cid:90) S dS p ΨΨ ln ΨΨ − κ (cid:90) Ω dV ζ ( D : D ) , (18)where κ = dβ/α . The conformational entropy S is non-negative, and zero if and only if Ψ ≡ Ψ (the state of uniform isotropy). The contribution P arises from the Maier-Saupe approximation tothe classical Onsager potential [53]. We can show that E satisfies the identity˙ E = − dα (cid:18) (cid:90) Ω dV E : E + β (cid:90) Ω dV E : S : E (cid:19) − dζ βα (cid:90) Ω dV D : ( D · D − S : D )+ 2 dζ (cid:18) βd R α (cid:19) (cid:90) Ω dV D : D + dζβd T α (cid:90) Ω dV |∇ D | − (cid:90) Ω dV (cid:90) S dS p (cid:16) d T |∇ ln Ψ | + d R |∇ p ln Ψ | (cid:17) , (19)where every integral density is non-negative (in particular, D : ( D · D − S : D ) ≥ ζ ), and the last term is a negative-definite dissipationterm arising from thermal fluctuations.Apparently, the above identity reveals that the nature of the active stress, i.e., contractile orextensile, is a central determinant of system behavior. For a contractile active stress ( α > ζ = 0), the energy E is driven to zero, corresponding toa uniform isotropic state. The same result holds for the motile case [21] and is consistent withthe simulations of semi-dilute active particle suspensions [11]. When steric interactions are moredominant, ordering effects can drive active contractile systems into complex dynamics [23]. When α <
3. The stability of uniform isotropic suspensions.
The linear theory of Extensor suspensionshas a simple and evocative structure that is obscured in the motile case. We give below the resultsof a stability calculation using the stream function Φ (i.e., u = ∇ × Φ ), and the details are givenin Appendix A.For an isotropic steady state we have Ψ = 1 / π , D = I / u = , and S ijkl = ( δ ik δ jl + δ il δ jk + δ ij δ kl ). Perturbing this steady state as D = I / ε D (cid:48) , S = S + ε S (cid:48) , u = ε u (cid:48) ,where ε <<
1, we find the 3D Smoluchowski equation (1) can be linearized at O ( ε ) as [17, 22, 23]: ∂ Ψ (cid:48) ∂t − π pp : E (cid:48) = 3 ζ π pp : D (cid:48) + d T ∆Ψ (cid:48) + d R ∇ p Ψ (cid:48) . (20)The momentum-balance equation in (12,13), together with Eq. (14), can be linearized as ∇ Π − (cid:18) β (cid:19) ∆ u (cid:48) = (cid:18) α − βζ (cid:19) ∇ · D (cid:48) , (21) ∇ · u (cid:48) = 0 . (22)Rewriting Eq. (21) in terms of Φ allows significant simplification. After some calculation (seederivations in Appendix A), and defining g = ∇ Φ , the linearized equations reduce to ∂ g ∂t = ( − C ( β ) α + C ( β ) ζ − d R ) g + d T ∆ g , (23)where C ( z ) = (3 / / (1 + z/ C ( z ) = (4 /
5) + (6 z/ / (1 + z/ σ = − C ( β ) α + C ( β ) ζ − d R − d T k . (24)The revealing aspect of this linear analysis is that the effect of activity, steric interactions, androtational diffusion all compete at the level of scale-independent dissipation and growth. Stericinteractions associated with ordering drive the system away from isotropy, which is abetted byactivity if α < k → + σ ( k = 0 + ) = − (cid:18) β + 15 (cid:19) α + 2 (cid:18) β + 6 β + 15 (cid:19) ζ − d R . (25)Hence, for periodic boundary conditions imposed on sufficiently large domains, the fastest growingscale is directly associated with the domain’s first periodic mode (which agrees well with oursimulation results in Sect. III.)As discussed below, another accessible steady state when steric interactions are strong is thenematic state. We will turn to a numerical treatment to study its stability except in the case ofsharply aligned suspensions. C. Reduced Models
In this section we will discuss two useful reductions of the kinetic theory. In the first, weinvestigate a special case of suspensions where the active particles are strictly aligned at eachpoint in space. This is an exact reduction of the kinetic theory, and is useful for constructingphenomenological models illustrating the instabilities of aligned suspensions. Since dynamicalsimulation of the kinetic theory is demanding given the number of independent variable (e.g., threein space, two on the unit sphere, plus time in 3D), we investigate the Bingham closure whichhas been successfully applied to passive rod suspensions. This is a coarse graining that expressesthe fourth-order S tensor in terms of the second-order D tensor, leading to a so-called Q -tensordynamics which we term BQ -tensor theory. We find that it gives an excellent accounting of thecomplex dynamics we observe in the kinetic theory, and at a far lower computational cost byremoving the orientational degrees of freedom.
1. Sharply aligned suspensions.
One analytically useful approximation is to assume negligible translational and rotational diffu-sion, and that the suspension is sharply aligned [21]. That is, we write Ψ( x , p , t ) = ¯ φδ [ p − q ( x , t )]where δ is the Dirac delta function on the unit sphere, and further assume that the concentrationis uniform, i.e., φ ≡ ¯ φ = 1. Integrating the Smoluchowski equation (1) against p , and making useof ∇ · u = 0 and the identity (cid:82) S dS p p ∇ p · ( ˙ p Ψ) = − (cid:82) S dS p ˙ p Ψ [17] yields ∂ q ∂t = − u · ∇ q + ( I − qq ) · ∇ u · q . (26)The extra stress in Eq. (14) simplifies to σ = α qq + β qqqq : E , and u is solved by ∇ Π − ∆ u = ∇ · ( α qq + β qqqq : E ) , (27) ∇ · u = 0 . (28)Equations (26)-(28) are apolar, and also invariant under the transformation q ↔ − q .Here we have a few comments: Equations (26)-(28) provide a useful comparison with the Bingham closure discussed below. Onthis point, we can derive an equation for the evolution of the dipolar tensor D = qq , using that ∂∂t ( qq ) = − ( u · ∇ ) ( qq ) + ∇ u · ( qq ) + ( qq ) · ∇ u T − qq : E ) qq . (29)This evolution equation is D ∇ + 2 ( D : E ) D = , (30)where D ∇ = ∂ D ∂t + u · ∇ D − (cid:0) ∇ u · D + D · ∇ u T (cid:1) is the upper-convected time derivative. Therefore,it becomes a closed evolution equation when taken together with the momentum balance equationand incompressible constraint: ∇ Π − ∆ u = ∇ · ( α D + β D ( D : E )) , (31) ∇ · u = 0 . (32) The sharply aligned dynamics are interesting to consider in the Lagrangian frame of the back-ground flow. As a reminder, this is defined by the flow map X ( α , t ) from the initial position α (i.e., the Lagrangian variable) to the current position X through ∂ X ∂t ( α , t ) = u ( X ( α , t ) , t ) (33)where X ( α ,
0) = α . The associated deformation tensor is then defined as F = ∂ X /∂ α ( F ij = ∂X i /∂α j ). When defining the field of unit vectors r ( α , t ) = F · q | F · q | , and using the identity D | F · q | Dt = | F · q | ( ∇ u : rr ), it is then straightforward to show that r satisfies Jeffery’s equation: ∂ r ∂t = ( I − rr ) · ∇ u · r . (34)Noting that F ( α , ≡ I and setting q ( α ,
0) = q ( α ), then the uniqueness of solutions to initialvalue problems gives immediately that, in the Lagrangian frame, the solution to Eq. (26) is q ( α , t ) = F · q | F · q | . (35)Equation (35) is simply a normalized version of the Result of Cauchy that ω = F · ω with ω thevorticity field evolved by the 3D incompressible Euler equations [54]. Naturally, this also yields aLagrangian expression for the tensor D = qq : D = ( F · q ) ( F · q ) | F · q | = F · D · F T tr ( F · D · F T ) . (36)To be consistent with its origins, D must be initialized by the dyadic product of a unit vector withitself. These expressions can be used to derive an integro-differential version in the Lagrangian0frame akin to the vorticity-based Lagrangian formulation of the incompressible Euler equations[54]. The sharply aligned dynamics, described either by Eqs. (26)-(28) or Eqs. (30)-(32), is partic-ularly useful for analyzing the stability of aligned suspensions. A globally-aligned steady state isgiven by q ≡ ˆz (or D = ˆzˆz ) and u ≡ . Following [21] by linearizing about this state, and seekingplane-wave solutions with wave-vector k , yields two growth rates: σ , = − αH , (Φ) , (37)where H (Φ) = cos Φ and H (Φ) = cos Φ (cid:0) Φ − (cid:1) , with Φ ∈ [0 , π ) the angle between thenormalized wave-vector ˆk and ˆz (i.e., cos Φ = ˆk · ˆz ). For α <
0, there is exponential growth at allplane-wave angles, independently of k , with the exception of Φ = π/ 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 ( I − qq ) · ∆ q on the right-hand-side of Eq. (26) [56]. This simply yields the modified growthrates ˜ σ , = − αH , (Φ) − d T k . (38)
2. Coarse-graining via the Bingham closure
The kinetic theory is high dimensional, having three spatial dimensions and two angles on theunit sphere (in 3D), and thus expensive to simulate. Alternatively, we can consider approximationsbased on moment-closure schemes to find approximate evolution equations for low-order p -averagedquantities. Many such schemes have been proposed [42, 47] but we consider here the Binghamclosure which is based on a quasi-equilibrium Ansatz and has been previously applied and comparedto kinetic theories of passive rod suspensions in simple rheological flows [42] (essentially the sameas ours but with α > p of theSmoluchowski equation yields the concentration equation (16), while integrating the tensor pp against the Smoluchowski equation yields a dynamical equation for D : D ∇ + 2 E : S = 4 ζ ( D · D − S : D ) + d T ∆ D − dd R (cid:18) D − φd I (cid:19) . (39)The left-hand-side of Eq. (39) is a tensor transport operator with its last term being generated byJeffery’s equation [57]. This evolution equation for D will be closed if the fourth-moment tensor S is expressed in terms of lower order moments, in particular if S = S [ D ]. Here we use the Binghamclosure [40, 41] which constructs an approximate local orientation distribution function Ψ B giventhe second-moment tensor D . The Bingham distribution takes the axisymmetric form:Ψ B [ T ] = exp ( T : pp ) Z [ T ] , (40)where T is a traceless symmetric tensor and Z is a normalization constant. For given D , the tensor T is determined by solving the relation D = (cid:90) S dS p pp Ψ B [ T ] (41)1at each point in space. Computationally, one can take advantage of the fact that the eigenvectorsof T and D are co-aligned and so both can be rotated into diagonal forms in a common frame.Consequently, only the eigenvalues of T need to be determined [41]. Given T , and hence Ψ B , thefourth-moment tensor S is then approximated by S B [ D ] = (cid:82) S dS p Ψ B pppp . The final equationstake the following forms: D ∇ + 2 E : S B [ D ] = 4 ζ ( D · D − S B [ D ] : D ) + d T ∆ D − dd R (cid:18) D − φd I (cid:19) , (42)together with ∇ Π − ∆ u = ∇ · σ B [ D ] , (43) ∇ · u = 0 . (44)And the extra stress is expressed as: σ B [ D ] = α D + β S B [ D ] : E − ζβ ( D · D − S B [ D ] : D ) . (45)We refer to Eqs. (42)-(45) as the BQ -tensor model. We can also show that both Eq. (39) andEq. (42) conserve tr( D ).For the above BQ − tensor model, we have several comments: Previous work [41, 42] has shown that in potential flows and in the weak flow limit (i.e., theequilibrium state), the Bingham distribution function yields the exact solution of the Smoluchowskiequation (with the rotational flux only). The Bingham distribution also arises naturally as describ-ing spatially uniform, nematically-ordered steady states. Such steady states can be obtained bysetting ˙ p = in Eq. (11), leading to a balance between the rotational diffusion and the Maier-Saupealignment torque (see, for example, Ezhilan et al. [23] and Gao et al. [25]): ∇ p ln Ψ = ξ ( I − pp ) · D · p , (46)where ξ = 2 ζ/d R . Consider the unit sphere with angle coordinates ( θ, Φ) and polar axis along ˆ z ,so that p = sin Φ cos θ ˆ x + sin Φ sin θ ˆ y + cos Φˆ z with θ ∈ [0 , π ) and Φ ∈ [0 , π ). In seeking a solutionΨ(Φ), we find D has the form D [Ψ] = A [Ψ] ˆzˆz + B [Ψ] I (47)where A [Ψ] = π (cid:82) π d ΦΨ(3 cos Φ −
1) sin Φ (see [23]). Then Eq. (46) can be integrated toΨ(Φ) = exp ( δ ( ξ ) cos 2Φ)2 π (cid:82) π exp ( δ ( ξ ) cos 2Φ) sin Φ d Φ = exp( T : pp T ) Z [ T ] , (48)where T = δ ( ξ )diag ( − / , − / , / δ ( ξ ) is the solution of the equation: δ = ξ (cid:82) π d Φ sin Φ(3 cos Φ −
1) exp ( δ cos 2Φ) (cid:82) π d Φ sin Φ exp ( δ cos 2Φ) , (49)which captures a bifurcation, with increasing ξ , from an isotropic ( δ ( ξ ) = 0, Ψ ≡ / π , and A [Ψ] = 0) to a nematically aligned state ( δ ( ξ ) (cid:54) = 0 and A (cid:54) = 0).Thus the Bingham distribution has a microscopic origin with respect to the kinetic theory, andcaptures the the isotropic-to-nematic transition. However, to our knowledge there is no asymptoticanalysis, rigorous or formal, that establishes error estimates between the Bingham closure and the2original kinetic theory. We suspect that such a connection could be established in the particularlimit of strong alignment torques at large ξ . 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, thekinetic theory itself is derived under several modeling and separation-of-scale assumptions). Whilestructurally similar to the active and passive nematic Q -tensor theories that follow the phenomeno-logical Landau-deGennes (LdG) approach for liquid-crystalline fluids (see, for example, [31]), thereare significant differences. In the LdG approach, the ordering and relaxation dynamics, and elasticstresses, follow from the LdG free-energy. The BQ -tensor theory has different nonlinearities gov-erning the steric ordering, in both the D dynamics and the associated stress, especially due to thepresence of S B . While in the LdG theory the spatial diffusion coefficient of Q is a Frank elasticcoefficient, in the BQ -tensor theory for D it is the spatial diffusion coefficient, d T , of the activeparticles. In the LdG theory, the transport operator for Q is corotational, while in the BQ -tensortheory 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 fromthe constraint of particle rigidity. Note that while the original kinetic theory has an energy law (i.e., Eqs. (18,19)), the BQ -tensortheory apparently does not. However, Li et al. [58] have recently developed more complicatedBingham closure models for passive rod suspensions that preserve the energetic structure. In theirtheory, the configurational entropy is approximated using the Bingham distribution, leading toa theory where the T tensor appears directly in the dynamics, rather than only indirectly in thecalculation of S B . This approach leads to a closer concordance of the BQ -tensor and LdG Q -tensortheories, though we have not pursued that here in the context of active suspensions. In the absence of spatial diffusion, Eq. (42) evolves D using only local information. Hence, likethe sharply aligned case, this equation could be reformulated naturally to describe the dynamicsin the Lagrangian frame. One simple and direct comparison of the kinetic and BQ theories is their linear behavior nearsteady state. We examine this in order, first for the uniform isotropic steady state, and second fora nematically aligned steady state. i. Stability near uniform isotropy.
We will show that the kinetic and BQ theories give identicalresults. This is easily seen by considering the evolution of the perturbation D (cid:48) . Multiplying thelinearized Smoluchowski equation (20) by pp , and integrating over the unit p -sphere gives theevolution equation for D (cid:48) : ∂ D (cid:48) ∂t − E (cid:48) = 4 ζ D (cid:48) + d T ∆ D (cid:48) − d R D (cid:48) . (50)The above equation is closed since E (cid:48) is a linear functional of D (cid:48) through the linearized Stokesequations (21). Inverting Eq. (21) in Fourier space for u (cid:48) (with u (cid:48) ( x ) = ˜ u ( k ) exp( i k · x )), forming˜ E (cid:48) [ ˜ D (cid:48) ], and substituting into Eq. (20) yields˜ D (cid:48) t = − γ (cid:16) ( I − ˆ k ˆ k ) · ˜ D (cid:48) · ˆ k ˆ k + ˆ k ˆ k · ˜ D (cid:48) · ( I − ˆ k ˆ k ) (cid:17) − ω ˜ D (cid:48) . (51)where γ = 15 (cid:32) α − ζβ β (cid:33) and ω k = − (cid:18) k d T − ζ d R (cid:19) , (52)and ˆ k = k /k .3Linearizing the BQ -theory Eqs. (42)-(44) gives an identical result, which is mainly due to thefact that the term D · D − S : D appears commonly in both the D transport and momentum-balanceequations (here we suppress the B subscript in the Bingham case). Also linearization of S : D produces S : D (cid:48) + S (cid:48) : D , and in both models we have S : D (cid:48) = 2 D (cid:48) /
15 and S (cid:48) : D = D (cid:48) / R as a rotation matrix such that ˆ k = R · ˆ z , and defineˆ D (cid:48) = R T · ˜ D (cid:48) · R . Then Eq. (51) can be rewritten asˆ D (cid:48) t = − γ (cid:16) ( I − ˆ z ˆ z ) · ˆ D (cid:48) · ˆ z ˆ z + ˆ z ˆ z · ˆ D (cid:48) · ( I − ˆ z ˆ z ) (cid:17) − ω k ˆ D (cid:48) . (53)In the rotated frame, the off-diagonal shear components ˆ D (cid:48) and ˆ D (cid:48) evolve asˆ D (cid:48) t = − ( γ + ω k ) ˆ D (cid:48) , (54)which yields the growth rate, σ = − ( γ + ω k ), identical to that in Eq. (24). The remainingcomponents of ˆ D (cid:48) instead satisfy ˆ D (cid:48) t = − ω k ˆ D (cid:48) . (55)In short, only ˆ D (cid:48) and ˆ D (cid:48) reflect the contributions of fluid flows, which can be stabilizing ordestabilizing. The remainder are decoupled from flow and show decay or growth only due tothermal fluctuations and steric interactions, respectively. ii. Stability near a nematically aligned homogeneous state.
As discussed above, with the inclusionof steric interactions via the Maier-Saupe potential there are additional homogeneous solutions thatarise as a balance between the steric alignment torque and the rotational diffusion. These solutionsare parametrized by ξ = 2 ζ/d R , and only the isotropic solutions exist when ξ is below a certaincritical value ξ c . In 3D there is a transcritical bifurcation at ξ = ξ c ≈ .
46 beyond which theisotropic 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 alignedas ξ → ∞ [23]. In 2D the bifurcation is instead a supercritical pitchfork bifurcation occurring at ξ c = 8 [25]. As discussed in Comment ( ) the form of distribution function for the nematicallyaligned 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 theoriesfor Pusher and Extensor suspensions [23], and for microtubule/motor-protein suspensions [24, 25].These show that activity can drive a bending instability where the wave-vector of maximal growth isaligned with the base orientation of the nematic state, as already suggested by the sharp-alignmentanalysis (see Eq. (38)). Here we examine the bending instability for Extensor suspensions in 2D forboth the kinetic and BQ -tensor theories. Solving Eq. (49) for δ ( ξ ) gives the Bingham distributionΨ , the base-state tensors T , D and S , and the base-state normalization factor Z . Perturbationsof these base quantities are introduced as ( ε (cid:28) T = T + ε T (cid:48) , Z = Z + εZ (cid:48) , D = D + ε D (cid:48) , S = S + ε S (cid:48) . (56)The perturbations can be further expressed in terms of T (cid:48) as Z (cid:48) = Z D : T (cid:48) , D (cid:48) = ( S − D D ) : T (cid:48) , S (cid:48) = ( R − S D ) : T (cid:48) , (57)with R = Z (cid:82) pppppp exp (cid:0) T : pp T (cid:1) d p the sixth-order p -moment tensor. These relationsallow us to express the high-order moments in terms of D (cid:48) through the intermediate variable T (cid:48) .The linearized equations are then transformed to Fourier space where we solve for the plane-wavegrowth rates σ numerically.4 k σ ζ /d R = 5 ζ /d R = 10 ζ /d R = 20 ζ /d R = 40sharply aligned FIG. 2: growth rates as a function of wave-number k as alignment torques are increased (i.e., increasing ζ/d R ). The black dotted curve is the growth rate predicted from the sharply-aligned case. Here α = − . β = 0 . d T = d R = 0 . Figure 2 shows the growth rates σ , for both the kinetic and BQ -tensor theories, of a plane-waveperturbation to the homogeneous nematic state. Here the plane-wave vector is in the direction ofnematic alignment, which has the greatest growth rates. Similar to the isotropic case, the fastestgrowth rate is obtained as k → + . These features are in agreement with the sharply alignedanalysis, and for motile and microtubule suspensions. We show the effect of increasing alignmenttorque while also comparing the results of the kinetic theory with those of the BQ -tensor theory.As ζ is increased for fixed α , we see an increase in the maximal growth rate and the band ofunstable modes. We also see an increasingly good correspondence between the full kinetic theoryand 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 featuresof the growth rate curves, it is apparently not the asymptotic limit of the BQ -tensor model as ζ tends to infinity. III. NUMERICAL STUDIES AND COMPARISONS
In this section we study numerically the dynamics of Extensor suspensions in both periodic andclosed geometries. We use simulations in 2D periodic geometries to compare the kinetic modelwith its coarse-grained BQ -tensor model at different levels of activity (i.e., varying α ). Herewe 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 ζ = β = 0). Model I makes a connection with our earlier studies of active suspensions [17, 21],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 thecomplex dynamics, particularly for Model II. We then use the BQ -tensor theory to study thedynamics of active suspensions under confinement. Flows in a disk have been studied previouslyusing other closures of Model I to show bifurcations, with increasing activity, from isotropic toauto-circulating flows, and thence to more complex dynamics. We recapitulate this using a BQ -5tensor version of Model I, and find excellent agreement (see Appendix B for details). We thenfocus on the full model (i.e., Model II) to study confined suspensions when steric interactions arestrong. These steric interactions make a nematically ordered suspension the natural state againstwhich activity competes, and we find both new bifurcations and dynamical phenomena. A. Numerical comparison of the kinetic and BQ -tensor theories In 2D, we evolve Ψ( x , θ ) where θ ∈ [0 , π ) is the particle orientation in the xy -plane (i.e., p = (cos θ, sin θ )). For the kinetic theory we use Eqs. (1,10,11) to evolve the distribution functionΨ, with the requisite background velocity field obtained by solving Eq. (8). The numerical methodis pseudo-spectral with 256-512 Fourier modes in each spatial direction, and 64 modes in theorientation θ . Equation (8) is inverted via Fourier transform as˜ u ( k ) = ik (cid:16) I − ˆ k ˆ k (cid:17) · ˜ σ · k . (58)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[24, 59]. The nondimensional diffusion coefficients, d T and d R , are chosen between 0.01-0.05. Weevolve the BQ − tensor system in a similar fashion, again using 256-512 Fourier modes in eachspatial direction.As a first comparison, Fig. 3 shows the late-time evolution of Model I for both the kinetic theoryand its BQ -tensor closure (again, setting β = ζ = 0). These simulations are performed from initialdata near uniform isotropy. In particular, the initial condition for Ψ has the formΨ( x, θ,
0) = 12 π (cid:18) i (cid:15) i cos( k i · x + ϕ i ) P i ( θ ) (cid:19) , (59)where the amplitudes (cid:15) i ∈ [ − .
01; 0 .
01] and phases ϕ i ∈ [0 , π ) are sampled from uniform distribu-tions, with P i a low-order trigonometric polynomial which satisfies Ψ ( θ + π ) = Ψ ( θ ). We includeonly the 10 longest-wavelength spatial modes in Eq. (59). The initial data for the BQ system issimply found by forming D from Ψ( x , θ, m ( x , t )overlaying the scalar order parameter s ( x , t ). The system evolves quickly away from the initiallyisotropic state to a temporally fluctuating state with high local order. Panels (b) and (d) show thefluid velocity field u = ( u, v ) overlaying the planar (scalar) vorticity ( ω = v x − u y ). The temporaland spatial dynamics for both models are complex and seemingly chaotic, and qualitatively similarto those seen in models of bacterial Pusher suspensions [15, 21] and of biofilament/molecular-motorassemblies [24, 25, 29, 30].While both the kinetic and BQ -tensor models start with the same data for D , the systemdynamics being chaotic means that the two descriptions diverge in their details. However, bothdescriptions do evolve towards a similar statistical structure, as shown by the time evolution in themean nematic order S ( t ) = (cid:82) Ω s ( x , t ) dA/ area(Ω), from near zero to a temporally fluctuating statewith mean near 0.6 (note the BQ theory slightly overestimating the order). The linear theory inSection A predicts that the kinetic and BQ -tensor models have the identical linear growth ratesfrom a uniform isotropic base state. Exponential growth to saturation is seen in S ( t ) following ashort equilibration period where high wave-number modes rapidly decay. Fitting S ( t ) during thegrowth period yields a growth rate σ ≈ .
16 for both models. This is close to the growth rate of6 a b c d FIG. 3: Collective motion of Extensor suspensions at late times ( t = 100) when using Model I, starting froman initial condition of a small perturbation about the uniform state by kinetic theory (a,b) and BQ − tensortheory (c,d). Panels (a,c) show the nematic director field, superimposed on the colormap of the scalar orderparameter s ; right panels (b,d) show the background fluid velocity vector field superposed upon the colormapof the associated vorticity. The simulation parameters are chosen as L = 15, α = − . d T = d R = 0 . the first fundamental mode in the simulation box, σ ≈ . BQ -tensor theories produce the same magnitude of flow is seen in the verysimilar 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 functionCor[ u , u ]( R ) = (cid:10) u ( x ) · u ( x + R ) / | u ( x ) | (cid:11) where R = | R | . We see a close match in this measurebetween the two models, as well as a slight negative minimum (indicating oppositely directed flowat R ≈
13, which is near the box size). This is consistent with the first fundamental spatial modein 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)), theseregions are associated with likewise oriented extensional straining flows (which are necessarily ofsmall vorticity). As found in simulations of active nematics [24, 25], both models show ± / m field, inhabiting regions of small nematic order, and co-located withfluid jets between oppositely signed vortical regions.7 a crack b c d FIG. 4: Active nematic states of Extensor suspensions at late times ( t = 100) when using Model II, startingfrom near-uniform state by kinetic theory (a,b) and BQ − tensor theory (c,d). Panels (a,c) show the nematicdirector field, superimposed on the colormap of the scalar order parameter s ; right panels (b,d) show thebackground fluid velocity vector field superposed upon the colormap of the associated vorticity. In (a), typicaldisclination defects are marked by open circles (+1 /
2) and squares ( − / L = 15, α = − . ζ = 1 . β = 0 . d T = d R = 0 . 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 containsmotile disclination defect pairs of order ± /
2. The induced active flows in panels (b) and (d) looksimilar to those of Model I, but with stronger vorticity. The defects shown here exist in regionsof small nematic order (dark blue), and are born as opposing pairs in elongated “incipient crack”regions [24, 25]. These crack structures locally decrease nematic order, and increase curvature ofdirector field lines. Characteristically, we find that the +1 / − / / BQ -tensor models, evolving from the same D data, show statistically similar long-time dynamics,8 Time M ea n n e m a ti c o r d e r Model IIModel I R C o r [ u , u ] FIG. 5: The mean (i.e., spatially averaged) nematic order S as a function of time. Inset: Cor [ u , u ]. Solidlines and open circles represent the results of the kinetic theory and the BQ − tensor model, respectively. and nearly identical growth from isotropy. As expected, we observe a more rapid growth in ModelII 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 fityielding σ ≈ .
2, close to the linear growth rate of the first fundamental mode in the box σ ≈ . BQ -tensor closure (andimproved over Model I). B. Active Flows in Confined Geometries
Having established an excellent correspondence between the kinetic and BQ -tensor theories, wenow 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 [46], Quinckerollers [6], and MT/motor-protein assemblies [7] to study the effects of confinement. Previous work[26, 27] has studied bacterial suspensions using models closely related to the BQ -tensor model ofModel I. To establish correspondence with these works, in Appendix B we compare the dynamicsof Model I in a circular chamber to previous analytical and numerical results by Goldstein andWoodhouse [26]. Overall, we find excellent agreement, showing a bifurcation with decreasing α < 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 Ω with boundary ∂ Ω; for example, see the disk in Fig. 6. Herewe assume the no-slip condition u | ∂ Ω = for the velocity field when solving the Stokes equation.For the D advection-diffusion equation we return to the kinetic description in Eqs. (1,10,11),9and find that conservation of total particle number is enforced by the boundary condition ( u Ψ − d T ∇ Ψ) | ∂ Ω · n = 0, which reduces to d T ∇ Ψ | ∂ Ω · n = 0. Taking the second moment with respect to p then gives the zero-flux boundary condition ( n · ∇ ) D = on ∂ Ω. Note that this condition doesnot enforce any particular alignment direction at the boundary.
Numerical Method.
We solve this finite-domain problem using a finite element method withan unstructured triangular grid. The governing equations are discretized using the standard 2DGalerkin formulation: (cid:90) Ω (cid:18) Re ∂ u ∂t · ˜u − p ∇ · ˜u + σ [ D ] : ∇ ˜u (cid:19) dA = 0 (60) (cid:90) Ω ( ∇ · u )˜ p dA = 0 (61) (cid:90) Ω (cid:18) D ∇ + 2 ∇ u : S − ζ ( D · D − D : S ) − d T ∇ D + 4 d R (cid:18) D − I (cid:19)(cid:19) : ˜D d Ω = 0 (62)where ˜ u , ˜ p , and ˜ D are test functions. Here we are solving an unsteady Stokes equation with apartial derivative term Re (cid:0) ∂ u ∂ t (cid:1) with an estimated Reynolds number of Re = 10 − . In the finiteelement solver, we choose different interpolation functions for different unknown variables in orderto satisfy the LBB condition [60, 61]. The fluid velocity is approximated by piecewise quadraticfunctions which are continuous over Ω (P2). The pressure and stress are approximated by piecewiselinear functions (P1). For this type of mixed finite elements, it is known that when Re is small,the numerical solution of the flow field can achieve third order accuracy in space [62]. For thisstudy, the temporal discretization is a second-order difference scheme. The governing equationsare reduced to a nonlinear system of algebraic equations which is solved by a modified Newton-Raphson algorithm [60, 61]. Specifically, in each iteration of the algorithm, the discretized linearsystem is solved by the GMRES method [63]. An incomplete LU preconditioner [64] is used toaccelerate convergence. In the following, we mainly vary parameters α and R but fix ζ = 0 . β = 0 . d T = d R = 0 .
1. Active flows in a circular chamber
To start, in Fig. 6 we consider an Extensor suspension with weak activity, α = − .
3. Froma nearly isotropic state, the orientation field rapidly orders (note the increase in the scalar orderparameter as compared to the cases in Appendix B where ζ = β = 0 are used in Model I) and two+1 / D .Increasing activity (i.e., decreasing α ) leads to persistent dynamics, as shown in the accompany-ing supplementary videos. In Fig. 7, α has been decreased to -0.7; see video S1 in the supplementaryinformation. Initially it follows the course of the previous case: the system moves quickly awayfrom isotropy with the appearance of two +1 / ± /
2; panels (a,e)). The − / / (a) (b) (c) (d) -0.1 0.1 (e) t = 30 -0.08 0.06 (f) t = 50 -0.04 0.03 (g) t = 60 (h) t = 200 FIG. 6: Evolution of an Extensor suspension from near isotropy, for small but negative α ( α = − . R = 2 .
0. The upper row shows the nematic field m overlaying the scalar orderparameter s (color field), while the lower row shows the velocity vector field u overlaying its scalar vorticity(color field). (a) (b) (c) (d) -0.2 0.4 (e) t = 100 -0.15 0.2 (f) t = 110 -0.15 0.15 (g) t = 120 -0.25 0.25 (h) t = 130 FIG. 7: Evolution of an Extensor suspension from near isotropy with α = − .
7, confined to a disk ofradius R = 2 .
0. The upper row shows the nematic field m overlaying the scalar order parameter s (colorfield), while the lower row shows the velocity vector field u overlaying its scalar vorticity (color field); seesupplementary video S1. (a) (b) (c) -0.4 0.8 (d) t = 170 -0.4 0.8 (e) t = 172 -0.4 0.8 (f) t = 176 FIG. 8: Evolution of an Extensor suspension from near isotropy with α = − .
0, confined to a disk ofradius R = 2 .
0. The upper row shows the nematic field m overlaying the scalar order parameter s (colorfield), while the lower row shows the velocity vector field u overlaying its scalar vorticity (color field); seesupplementary video S2. formation of two low nematic order cracks (panels (d,h)). Accompanying the nematic structurevariation, we observe the persistence of background rotational flows which strengthen and weakenthrough the various stages of this periodic dynamics (lower row of panels).The dynamics becomes quasi-periodic or chaotic at yet higher activity. For α = −
2, Fig. 8shows a short period in time that illustrates the complex system dynamics (see also supplementaryvideo S2). In panels(a,d) we see that as α further decreases, there is still a basic dynamics of tworotating +1 / − / / α = −
5, which shows the yet morecomplex interplays between defects, cracks, and the defect pairs thereby produced. The associatedvortical flows are no longer dominated by a single-signed vortex but is much more dipolar, with afluctuating predominance of one signed vortex over the other. While complicated both spatiallyand temporally, the dynamics nonetheless maintains a quasi-periodicity.2
Time M e a n o r d er α = -1.0 α = -2.5 α = -3.0 α = -4.0 Time M e a n o r d er α = -0.3 α = -0.7 α = -2.0 α = -5.0 - α σ analyticalnumerical - α σ numerical FIG. 9: Evolution of the mean nematic order S ( t ) at different levels of activity with R = 2, correspondingto the cases shown in Figs. 6-8, and supplementary videos S1-S3. The initial growth of S is found to beexponential in time, and the inset shows a fit to the growth rates σ (open squares) as a function of − α . Thedashed linear fit suggests that the growth rate scales approximately with | α | . Figure 9 shows the evolution of the mean nematic order parameter S ( t ) for these various simula-tions. In all cases the rapid growth of S ( t ) away from zero reflects the transition from isotropy to astate of higher nematic order. For the simulation with the lowest activity, α = − . S = 0 . < t <
60 (see Fig. 6), leads to relaxation to a homogeneously orderedstate with S slightly less than 0.9. The emergence of persistent and complex defect dynamics inthe other cases is reflected in the saturated behaviors of S ( t ) following growth from isotropy. Inall cases the initial growth from zero can be well-fitted by an exponential in time. The calculatedgrowth rates are plotted in the inset, and show a linear increase with increasing activity α .We have also examined the dynamics in circular chambers of different radii and found generallysimilar behaviors. That being said, in a smaller confinement we have also discovered a fascinatingperiodic coherent structure. Figure 10 shows the dynamics of the system for α = − cf. Fig. 8)but with confinement radius R = 1 .
25, as it enters into a periodic dynamics. The defect dynamicsstarts from the apparent splitting of a +1-order defect into two +1 / (a) (b) (c) -0.2 0.4 (d) t = 100 -0.2 0.4 (e) t = 180 -0.2 0.4 (f) t = 330 FIG. 10: Evolution of an Extensor suspension from near isotropy with α = − .
0, when confined to a smallerdisk of radius R = 1 .
25. The upper row shows the nematic field m overlaying the scalar order parameter s (color field), while the lower row shows the velocity vector field u overlaying its scalar vorticity (color field);see supplementary video S4.
2. Active flows in a biconcave chamber
The coarse-grained BQ -tensor model and its finite-element discretization greatly facilitates theexploration of active nematic flows in more complex geometries. In Fig. 11, we study the dynamicsof an Extensor suspension confined to a biconcave chamber where two circular chambers, each orradius R = 2, are connected smoothly by a bridge of width d . Here we fix α = −
2. Supplementaryvideos S5-S8 show the dynamics as d is successively increased ( d =0.2, 0.8, 1.2, and 2). We findthat when the neck is thin, the dynamics in the two chambers are very similar to the case shownin Fig. 7 for α = −
2, and appear to be evolving independently of each other; see supplementaryvideo S5. This lack of interchamber communication is consistent with the velocity stagnation zonein the bridge (panel (b)). While a thin neck would be sufficient to prevent interplay between thetwo chambers, this effect is likely reinforced by the alignment of the nematic field lines across theneck 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, apparentlyblocking any substantial breaching flows. As d is further increased, beyond d crit ≈ .
0, the dynamicsin the two chambers begins to couple. Figure 7(a) and supplementary video S7 for d = 1 . (a) -0.3 0.4 (b) d FIG. 11: For two circular chambers connected by a narrow neck of width d = 0 .
2, and filled with an Extensorsuspension with α = −
2, (a) the nematic vector field m overlaying the scalar order parameter s (color field),and (b) the fluid velocity field u overlaying its scalar vorticity ω . See supplementary video S5. from one side to the other. For the largest neck width, d = 2 .
0, Figure 7(b) and supplementaryvideo S8 show that flows now move freely between the two chambers. This transition is reminiscentof a Fredricks transition where an electric or magnetic field needs to exceed a critical intensity todeform a uniformly aligned nematic liquid crystal material spanning between two plates. Frederickstransitions have been studied in other contexts in the field of active matter [65, 66].
IV. DISCUSSION
We explored a multi-scale model, and its simplifications, for suspensions of rod-like particlesthat exert active dipolar extensile stresses on the immersing solvent. This model is relevant tosuspensions or bundles of microtubules that undergo polarity sorting driven by crosslinking motorproteins [67–70], and directly models suspensions of chemically-active particles whose consumptionof a chemical fuel creates extensile flows along the particle [36, 37, 39]. Aside from analysisof the model, a central goal here was to investigate the Bingham closure [40–42] as applied toactive suspensions. We found that the resultant BQ -tensor theory gave an excellent accountingof the chaotic self-driven flows of active suspensions (both with and without strong alignmentinteractions), and allowed us to investigate with relative ease the behavior of an active-nematicfluid under confinement. We found several fascinating behaviors, such as defects being absorbedat (and perhaps birthed by) boundaries, complex rotational flows, and an apparent Frederickstransition associated with nematic elasticity.5 (a) (b) -0.3 0.4 FIG. 12: Two circular chambers connected by a neck of width d , and filled with an Extensor suspensionwith α = −
2. (a) For d = 1 .
2, the nematic vector field m overlaying the scalar order parameter s (colorfield). The arrows show the direction of shearing of the nematic field lines. The inset shows the fluid velocityfield u overlaying its scalar vorticity ω in the neck region. The arrows also correspond to the location anddirection of fluid jets between the two chambers. See supplementary video S7. (b) The nematic vector fieldand scalar order parameter for d = 2 .
0. See supplementary video S8.
We remark that the BQ -tensor theory has no unfixed constants with respect to the full kinetictheory, and is far cheaper to simulate (essentially by a factor inversely proportional to the numberof 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 suchas those based on the Landau-deGennes approach (e.g. [30, 31, 58]). We are currently working onimproving such closure approaches, from a variational perspective, and extending them to polarkinetic theories such as for motile suspensions [17, 20] and microtubule/motor-protein suspensions[24, 25]. A. Length scale determination
We close with a short digression into the determination of length scales of instability in activesuspensions, as a suspension of Extensors is a particularly simple case to discuss. As previouslydiscussed, it seems generic that bulk models of active suspensions lack an intrinsic length scaledetermined by a most unstable mode [14, 17, 24]. This is seen most directly in Eq. (24) for thegrowth rates of plane-wave perturbations from isotropy, and in our stability analyses of alignedstates. It is also consistent with our simulations in periodic domains where observed growth rates6 k σ ζ /d R = 5 ζ /d R = 10 ζ /d R = 20 ζ /d R = 40sharply aligned FIG. 13: Linear stability analysis for a thin layer of homogenous nematically-aligned Extensor suspensionconfined between two fluids. growth rates are plotted as a function of wave-number k at various ζ/d R usingboth the BQ -tensor theory (dashed lines) and the kinetic theory (solid lines). The black dotted curve isthe growth rate predicted from the sharply-aligned case. Computation parameters are chosen as α = − . β = 0 . d T = d R = 0 . from isotropy are well-described by that for the first fundamental mode with k = 2 π/L of theperiodic domain. Nonetheless, some studies have sought to identify an intrinsic length scale ofinstability [71]. 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. [24] inmodeling the experiments of Sanchez et al. [7], as well as Leoni and Liverpool [72]), or sits atop asubstrate [73].To be specific, consider a thin active layer evolving in the xy -plane at z = 0, across which layeractivity 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 Stokesflows can be solved analytically under the conditions of zero velocity as | z | → ∞ and continuity ofvelocity at z = 0. This determines the interfacial velocity field u ( x, y ) as a function of the activelayer extra stress σ [24, 25, 74]: ˜ u ( k ) = i k (cid:16) I − ˆ k ˆ k (cid:17) · ˜ σ · k (63)where ˆk = k / | k | is the 2D unit wavevector. While not strictly necessary, we have assumed that theinterfacial flow is incompressible (as seems consistent with the experiments of Sanchez et al. [7],private communication with Z. Dogic). Thus, activity at the interface drives flows in the two outerfluids, which in turn provides a drag on the active surface. Simulating this situation amounts toreplacing 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/ γ and ω k (Eq. (52)) that govern the growth/decayrate σ k = − ( γ + ω k ) (Eq. (54)). For the 2D (or 3D) bulk fluid case, γ has the form γ = ˜ α/ (1 +˜ β ); while for the immersed interface case considered here, it is modified to γ k = ( ˜ αk/ / (1 +( ˜ βk/ k from zero to ˜ α/ ˜ β . For ˜ α <
0, competition of activity withtranslational diffusion in ω k yields a unique critical wavenumber k c > σ k is maximized.This same effect manifests itself in the instability of a homogenous nematically-aligned state, andyields an intrinsic scale that is observed in nonlinear simulations (see Gao et al. [24] for the more7complicated case of microtubule/motor-protein suspensions). For a simple Extensor suspension,Fig. 13 shows the numerically calculated growth rates σ k (red curve) of a plane-wave perturbation tothe homogeneous nematic state, with wave-vector taken in the direction of maximal growth, θ = 0.This being the direction of maximal growth is in agreement with the sharply-aligned analysis ofSect. II C 1, and has been found in other studies of active suspensions [17, 21, 23–25]. The computedcurve also shows an intrinsic wave-number of maximal growth, as was found analytically for theisotropic state. As ζ is increased for fixed α , we see an increase in maximal growth rate anddecrease in its associated length scale. In that limit we find an increasing correspondence betweenthe full kinetic theory and the BQ -theory.The parabolic-like structure of the growth rate curve is in qualitative agreement with our phe-nomenological sharply-aligned model when modified to the immersed interface case. This modifi-cation alters the plane-wave growth rates ˜ σ , in Eqs. (38)(where β = 0) to the formˆ σ , = − αH , (Φ) k/ − d T k . (64)For α <
0, this expression yields at each Φ (with H maximal at Φ = 0) a maximal growth rate σ max = α H (Φ) / d T , occurring at k c = − αH (Φ) / d T . To emphasize an important point: thequadratic increase in maximal growth rate with activity follows from k c increasing linearly withactivity (i.e., decreasing negative α ). For the bulk case, where the system size sets the mostunstable scale, maximal growth rates should scale linearly with activity (as in Eqs. (24), (38), andthe 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 grantK25GM110486 (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, itdoes have an especially simple and evocative structure that is obscured in the motile case. For anisotropic steady state we have Ψ = 1 / π , D = I / u = , and S ijkl = ( δ ik δ jl + δ il δ jk + δ ij δ kl ).By perturbing this steady state as D = I / ε D (cid:48) , S = S + ε S (cid:48) , u = ε u (cid:48) , the 3D linearizedSmoluchowski equation (1) becomes [17, 22, 23] ∂ Ψ (cid:48) ∂t − π pp : E (cid:48) = 3 ζ π pp : D (cid:48) + d T ∆Ψ (cid:48) + d R ∇ p Ψ (cid:48) . (A1)This result uses that, for any p -independent tensor A , ∇ p · (( I − pp ) · A · p ) = tr( A ) − pp : A ,and that tr( D (cid:48) ) = 0 for a constant density flow. Multiplying Eq. (A1) by pp , and integrating overthe unit p -sphere gives the evolution equation for D (cid:48) : ∂ D (cid:48) ∂t − E (cid:48) = 4 ζ D (cid:48) + d T ∆ D (cid:48) − d R D (cid:48) . (A2)This result uses that S : A = 2 A /
15 for any symmetric trace-free tensor A . An identical result isobtained through linearization of Eq. (39). The momentum balance equation and incompressiblecondition, Eq. (12)-(14), linearize to ∇ Π − (cid:18) β (cid:19) ∆ u (cid:48) = (cid:18) α − βζ (cid:19) ∇ · D (cid:48) , (A3) ∇ · u (cid:48) = 0 . (A4)8Rewriting Eq. (A3) in terms of the stream function allows significant simplification. By assumingthat the fluid domain is simply connected, we can define the vector stream function Φ such that u = ∇× Φ and ∆ Φ = −∇×∇× Φ = − ω . Taking the curl of Eq. (A3) then yields ∇ Φ = ˜ α ∇ × ∇ · D (cid:48) , (A5)where ˜ α = ( α − βζ ) / (1 + β/ D (cid:48) we can express the right-hand-sideof this equation, component-wise, in terms of Ψ (cid:48) :( ∇ × ∇ · D (cid:48) ) i = (cid:90) S dS p (cid:15) ijk ∂∂x j ∂∂x l ( p k p l ) Ψ (cid:48) = (cid:90) S dS p L i Ψ (cid:48) . (A6)This defines the vector operator L , given component-wise by L i = ε ijk p k p l ∂ /∂x j ∂x l . Equation(A3) then becomes ∇ Φ = ˜ α (cid:90) S dS p L Ψ (cid:48) . (A7)Now we turn to the evolution equation (A1). The contracted term pp : ∇ u can be rewritten as pp : ∇ u = − L · Φ (this uses that pp : E = pp : ∇ u , and that ε ijk = − ε kji , followed by an indexrelabeling). Then, Eq. (A1) can be rewritten as ∂ Ψ (cid:48) ∂t = − π L · Φ + 3 ζ π pp : D (cid:48) + d T ∆Ψ (cid:48) + d R ∇ p Ψ (cid:48) . (A8)Now we take a time derivative of Eq. (A7) (i.e. ∇ Φ t = ˜ α (cid:82) S dS p L Ψ (cid:48) t ), and use Eqs. (A8) and(A7) to find ∂ ∇ Φ ∂t = − α π (cid:90) S dS p L ( L · Φ ) + 3 ζ ˜ α π (cid:90) S dS p L ( pp : D (cid:48) ) + d T ∇ Φ + ˜ αd R (cid:90) S dS p L ∇ p Ψ (cid:48) . (A9)The terms on the right-hand-side of Eq. (A9) can be calculated as (cid:90) S dS p L ( L · Φ ) = −∇× (cid:18) ∇ · (cid:90) S dS p ( pp ( ∇ u : pp )) (cid:19) = 4 π ∇ Φ , (A10) (cid:90) S dS p L (cid:0) pp T : D (cid:48) (cid:1) = 8 π (cid:90) S dS p L Ψ (cid:48) , (A11) (cid:90) S dS p L ∇ p Ψ (cid:48) = − (cid:90) S dS p L Ψ (cid:48) . (A12)Finally, by defining g = ∇ Φ and using the definition for ˜ α the linearized dynamics reduces to ∂ g ∂t = ( − C ( β ) α + C ( β ) ζ − d R ) g + d T ∆ g , (A13)where C ( z ) = (3 / / (1 + z/ C ( z ) = (4 /
5) + (6 z/ / (1 + z/ σ = C ( β ) α + C ( β ) ζ − d R − d T k . (A14)9 (a) α = -0.6 -0.04 0.08 (b) α = -0.7 r/R v / V m a x (c) analyticalnumerical FIG. 14: Panels (a) and (b) show the fluid velocity vector field u and its scalar vorticity ω (color field)for two Extensor suspensions confined to a circular chamber of radius R = 2 using Model I. Panel (a) uses α = − .
6, while (b) uses α = − .
7, which spans the critical value α ≈ − .
63 calculated by Woodhouseand Goldstein [26]. Panel (a) shows the result of relaxation to a quiescent flow, while panel (b) shows theemergence of an auto-circulating flow. Panel (c) shows a comparison of the analytical calculation for theazimuthal velocity in (b), with the numerically determined value taken along the white line in (b).
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 simulationnear uniform isotropy in the weakly aligned limit with ζ = β = 0, which corresponds to neglectingsteric interactions between filaments. This limit corresponds to the model previously studied byWoodhouse and Goldstein [26], who discovered a symmetry-breaking bifurcation at which thestationary state transitions to a stable circulating state. For higher activity, another bifurcationproduces an oscillatory state with a pair of mobile defects.We expected that our Bingham closure would reproduce the Woodhouse and Goldstein resultsthat used the Hinch and Leal closure [75] when neglecting steric interactions. Indeed, we recover anidentical system of linearized equations for perturbation about an isotropic state. Numerically, wecapture the bifurcation to a circulatory state (Fig. 14), and the growth rates agree quantitativelywith linear stability analysis (Fig. 16). The shape of the flow profile matches that determinedanalytically [26] (Fig. 14).For increasing activity (i.e., decreasing α ), in Fig. 15 we reproduce the oscillatory two-defectstate found by Woodhouse and Goldstein [26]. For sufficiently high activity, this oscillatory stateis unstable to the production of cracks and defects in the nematic field, and multiple vortices. [1] S. Ramaswamy. The mechanics and statistics of active matter. Ann. Rev. Cond. Matt. Phys. , 1:323–345,2010.[2] D. Saintillan and M. Shelley. Active suspensions and their nonlinear models.
C. R. Phys. , 14:497–517,2013.[3] Christopher Dombrowski, Luis Cisneros, Sunita Chatkaew, Raymond E. Goldstein, and John O.Kessler. Self-Concentration and Large-Scale coherence in bacterial dynamics.
Physical Review Let-ters , 93(9):098103, August 2004.[4] H. P. Zhang, Avraham Be’er, E.-L. Florin, and Harry L. Swinney. Collective motion and densityfluctuations in bacterial colonies.
Proceedings of the National Academy of Sciences , 107(31):13626– (a) (b) (c) (d) -0.1 0.2 (e) α = -1.0 -0.4 0.9 (f) α = -2.5 -0.4 1.2 (g) α = -3.0 -2 1 (h) α = -4.0 FIG. 15: Snapshots of the characteristic nematic and flow structures of Extensor suspensions of increasingactivity(i.e., decreasing α ), when confined to a circular chamber of radius R = 2. The upper row shows thenematic field m overlaying the scalar order parameter s (color field), while the lower row shows the velocityvector field u overlaying its scalar vorticity (color field).13630, August 2010.[5] Jeremie Palacci, Stefano Sacanna, Asher Preska Steinberg, David J Pine, and Paul M Chaikin. Livingcrystals of light-activated colloidal surfers. Science , 339(6122):936–940, 2013.[6] A. Bricard, J. Caussin, N. Desreumaux, O. Dauchot, and D. Bartolo. Emergence of macroscopicdirected motion in populations of motile colloids.
Nature , 503:95–98, 2013.[7] T. Sanchez, D. Chen, S. DeCamp, M. Heymann, and Z. Dogic. Spontaneous motion in hierarchicallyassembled active matter.
Nature , 491:431–434, 2012.[8] M. Shelley. The dynamics of microtubule/motor-protein assemblies in biology and physics.
Annu. Rev.Fluid Mech. , accepted.[9] Juan P. Hernandez-Ortiz, Christopher G. Stoltz, and Michael D. Graham. Transport and collec-tive dynamics in suspensions of confined swimming particles.
Physical Review Letters , 95(20):204501,November 2005.[10] David Saintillan and Michael J. Shelley. Orientational order and instabilities in suspensions of self-locomoting rods.
Physical Review Letters , 99, August 2007.[11] D. Saintillan and M. J. Shelley. Emergence of coherent structures and large-scale flows in motilesuspensions.
Journal of The Royal Society Interface , 9(68):571–585, August 2011.[12] K. Yeo, E. Lushi, and P. Vlahovska. Collective dynamics in a binary mixture of hydrodynamicallycoupled microrotors.
Physical Review Letters , 114:188301, 2015.[13] B. Delmotte, E. Keaveny, F. Plourabou´o, and E. Climent. Large-scale simulation of steady and time-dependent active suspensions with the force-coupling method.
J. Comp. Phys. , 114:312, 1994.[14] R. A. Simha and S. Ramaswamy. Statistical hydrodynamics of ordered suspensions of self-propelledparticles: waves, giant number fluctuations and instabilities.
Physica A: Statistical Mechanics and itsApplications , 306:262–269, 2002.[15] Igor S. Aranson, Andrey Sokolov, John O. Kessler, and Raymond E. Goldstein. Model for dynamicalcoherence in thin films of self-propelled microorganisms.
Physical Review E , 75(4):040901, April 2007.[16] Patrick T. Underhill, Juan P. Hernandez-Ortiz, and Michael D. Graham. Diffusion and spatial corre-lations in suspensions of swimming particles.
Physical Review Letters , 100(24):248101, 2008.[17] D. Saintillan and M. J. Shelley. Instabilities and pattern formation in active particle suspensions: Time M e a n o r d er α = -1.0 α = -2.5 α = -3.0 α = -4.0 Time M e a n o r d er α = -0.3 α = -0.7 α = -2.0 α = -5.0 - α σ analyticalnumerical - α σ numerical FIG. 16: Corresponding to the simulations shown in Fig. 15, the evolution of the mean nematic order S ( t )at various values of α . Again, the initial growth of S is found to be exponential in time, and the inset showsa fit to the growth rate σ (open squares) as a function of − α . The solid line is the growth rate predicted bythe analytical axisymmetric mode: σ = − α − d R − (cid:0) λ R (cid:1) d T .kinetic theory and continuum simulations. Physical Review Letters , 100(17):178103, 2008.[18] Charles W. Wolgemuth. Collective swimming and the dynamics of bacterial turbulence.
BiophysicalJournal , 95(4):1564–1574, August 2008.[19] Aparna Baskaran and M. Cristina Marchetti. Statistical mechanics and hydrodynamics of bacterialsuspensions.
Proceedings of the National Academy of Sciences , 106(37):15567–15572, September 2009.[20] Ganesh Subramanian and Donald L. Koch. Critical bacterial concentration for the onset of collectiveswimming.
Journal of Fluid Mechanics , 632:359–400, 2009.[21] David Saintillan and Michael J. Shelley. Instabilities, pattern formation, and mixing in active suspen-sions.
Physics of Fluids , 20:123304–123304–16, 2008.[22] Christel Hohenegger and Michael J. Shelley. Stability of active suspensions.
Physical Review E ,81(4):046311, April 2010.[23] B. Ezhilan, M. Shelley, and D. Saintillan. Instabilities and nonlinear dynamics of concentrated activesuspensions.
Phys. Fluids , 25:070607, 2013.[24] T. Gao, R. Blackwell, M. A. Glaser, M. D. Betterton, and M. J. Shelley. Multiscale polar theory ofmicrotubule and motor-protein assemblies.
Phys. Rev. Lett. , 114:048101, 2015.[25] T. Gao, R. Blackwell, M. Glaser, M. Betterton, and M. Shelley. Multiscale modeling and simulation ofmicrotubule/motor protein assemblies.
Phys. Rev. E , 92:062709, 2015.[26] F. Woodhouse and R. Goldstein. Spontaneous circulation of confined active suspensions.
PhysicalReview Letters , 109, 2012.[27] B. Ezhilan and D. Saintillan. Transport of a dilute active suspension in pressure-driven channel flow.
Journal of Fluid Mechanics , 777:482–522, 2015.[28] M. Theillard, R. Alonso-Matilla, and D. Saintillan. Geometric control of active collective motion.
SoftMatter , to appear.[29] S. Thampi, R. Golestanian, and J. Yeomans. Velocity correlations in an active nematic.
Phys. Rev.Lett. , 111:118101, 2013.[30] L. Giomi, M. Bowick, X. Ma, and M. Marchetti. Defect annihilation and proliferation in active nematics.
Phys. Rev. Lett. , 110:228101, 2013.[31] L. Giomi. Geometry and topology of turbulence in active nematics.
Phys. Rev. X , 5:031003, 2015.[32] D. W. Adams and J. Errington. Bacterial cell division: assembly, maintenance and disassembly of thez ring.
Nat. Re.v Micro. , 7:642–653, 2009. [33] A. Buka, P. Palffy-Muhoray, and Z. Racz. Viscous fingering in liquid crystals. Phys. Rev. A , 36:3984– 3989, 1987.[34] M. Shelley and T. Ueda. The Stokesian hydrodynamics of flexing, stretching filaments.
Phys. D , 146:221– 245, 2000.[35] Felix C. Keber, Etienne Loiseau, Tim Sanchez, Stephen J. DeCamp, Luca Giomi, Mark J. Bowick,M. Cristina Marchetti, Zvonimir Dogic, and Andreas R. Bausch. Topology and dynamics of activenematic vesicles.
Science , 345:1135–1139, 2014.[36] M. Davies Wykes, J. Palacci, T. Adachi, L. Ristroph, X. Zhong, M. Ward, J. Zhang, and M. Shelley.Dynamic self-assembly of microscale rotors and swimmers.
Soft Matter , 12:4584–4589, 2016.[37] E. Jewell, W. Wang, and T. Malloukl. Catalytically driven assembly of trisegmented metallic nanorodsand polystyrene tracer particles.
Soft Matter , 12:2501–2504, 2016.[38] Walter F. Paxton, Kevin C. Kistler, Christine C. Olmeda, Ayusman Sen, Sarah K. St. Angelo, YanyanCao, Thomas E. Mallouk, Paul E. Lammert, and Vincent H. Crespi. Catalytic nanomotors: Au-tonomous movement of striped nanorods.
Journal of the American Chemical Society , 126(41):13424–13431, October 2004.[39] A. Pandey, P. B. Sunil Kumar, and R. Adhikari. Flow-induced nonequilibrium self-assembly in suspen-sions of stiff, apolar, active filaments.
Soft Matter , 12:9068–9076, 2016.[40] C. Bingham. An antipodally symmetric distribution on the sphere.
Ann. Stat. , 2:1201?225, 1974.[41] C. Chaubal and L. G. Leal. A closure approximation for liquid-crystalline polymer models based onparametric density estimation.
J. Rheol. , 42:177–201, 1998.[42] J. Feng, C. Chaubal, and L. G. Leal. Closure approximations for the doi theory: which to use insimulating complex flows of lcps.
J. Rheol. , 42:1095–1119, 1998.[43] H. Wioland, F. Woodhouse, J. Dunkel, J. Kessler, and R. Goldstein. Confinement stabilizes a bacterialsuspension into a spiral vortex.
Phys. Rev. Lett. , 110:268102, 2013.[44] H. Wioland, F. Woodhouse, J. Dunkel, and R. Goldstein. Ferromagnetic and antiferromagnetic orderin bacterial vortex lattices.
Nature Physics , 4, 2016.[45] A.C.H. Tsang and E. Kanso. Density shock waves in con ned microswimmers.
Physical Review Letters ,116:048101, 2016.[46] E. Lushi, H. Wioland, and R. Goldstein. Fluid flows generated by swimming bacteria drive self-organization in confined fluid suspensions.
Proceedings of the National Academy of Sciences , 111:9733–9738, 2014.[47] M. Doi and S. F. Edwards.
The theory of polymer dynamics , volume 73. Oxford University Press, USA,1986.[48] C. Hohenegger and M. Shelley.
Dynamics of complex bio-fluids . Oxford University Press, 2011.[49] J. Keller and S. Rubinow. Slender-body theory for slow viscous flow.
Journal of Fluid Mechanics ,75:705–714, 1976.[50] W. Maier and A. Saupe. Eine einfache molekulare theorie des nematischen kristallinfl¨ussigen zustandes.
Zeit. Nat. Teil A , 13:564, 1958.[51] G. K. Batchelor. The stress system in a suspension of force-free particles.
Journal of Fluid Mechanics ,41(3):545–570, 1970.[52] J. Han, Y. Luo, W. Wang, P. Zhang, and Z. Zhang. From microscopic theory to macroscopic theory:a systematic study on modeling for liquid crystals.
Arch. Rational Mech. Anal , 215:741–809, 2015.[53] L. Onsager. The effects of shapes on the interaction of colloidal particles.
Ann. N.Y. Acad. Sci. , 51:627,1949.[54] S. Childress.
An introduction to theoretical fluid mechanics , volume 19. Courant Lecture Notes, 2009.[55] L. Giomi, L. Mahadevan, B. Chakraborty, and M. Hagan. Excitable patterns in active nematics.
Physical Review Letters , 106, 2011.[56] F. Lin and C. Liu. Nonparabolic dissipative systems modeling the flow of liquid crystals.
Communica-tions on Pure and Applied Mathematics , 48:501–537, 1995.[57] G. Jeffery. The motion of ellipsoidal particles immersed in a viscous fluid.
Proc. Roy. Soc. Lond. Ser.A , 102:161–179, 1922.[58] S. Li, W. Wang, and P. Zhang. Local well-posedness and small deborah limit of a molecular-basedq-tensor system.
DCDS-B , 20:2611–2655, 2015.[59] T. Hou, J. Lowengrub, and M. J. Shelley. Removing the stiffness from interfacial flows with surfacetension.
J. Comp. Phys. , 114:312, 1994. [60] T. Gao and H. H. Hu. Deformation of elastic particles in viscous shear flow. J. Comput. Phys. ,228:2132–2151, 2009.[61] H. H. Hu, M. Y. Zhu, and N. Patankar. Direct numerical simulations of fluid solid systems using thearbitrary lagrangian eulerian technique.
J. Comput. Phys. , 169:427–462, 2001.[62] F. Thomasset.
Implementation of Finite Element Methods for Navier-Stokes Equations . SpringerVerlag, New York, 1981.[63] Y. Saad and Y. Schultz. Gmres: A generalized minimal residual algorithm for solving nonsymmetriclinear systems.
SIAM J. Sci. Stat. Comput. , 7:856–869, 1986.[64] E. Chow and Y. Saad. Experimental study of ilu preconditioners for indefinite matrices.
J. Comput.Appl. Math. , 86:387–414, 1997.[65] Luca Giomi, M. Cristina Marchetti, and Tanniemola B. Liverpool. Complex spontaneous flows andconcentration banding in active polar films.
Physical Review Letters , 101:198101, 2008.[66] R. Voituriez, J. Joanny, and J. Prost. Spontaneous flow transition in active polar gels.
EurophysicsLetters , 70:404, 2005.[67] Thomas Surrey, Fran¸cois N´ed´elec, Stanislas Leibler, and Eric Karsenti. Physical properties determiningSelf-Organization of motors and microtubules.
Science , 292(5519):1167–1171, May 2001.[68] Poul M. Bendix, Gijsje H. Koenderink, Damien Cuvelier, Zvonimir Dogic, Bernard N. Koeleman,William M. Brieher, Christine M. Field, L. Mahadevan, and David A. Weitz. A quantitative analysis ofcontractility in active cytoskeletal protein networks.
Biophysical Journal , 94(8):3126–3136, April 2008.[69] Christian Hentrich and Thomas Surrey. Microtubule organization by the antagonistic mitotic motorskinesin-5 and kinesin-14.
The Journal of Cell Biology , 189(3):465–480, May 2010.[70] Daniel Gordon, Anne Bernheim-Groswasser, Chen Keasar, and Oded Farago. Hierarchical self-organization of cytoskeletal active networks.
Physical Biology , 9(2):026005, April 2012.[71] S. Thampi, R. Golestanian, and J. Yeomans. Instabilities and topological defects in active nematics.
Europhys. Lett. , 105:18001, 2014.[72] M. Leoni and T. B. Liverpool. Swimmers in thin films: from swarming to hydrodynamic instabilities.
Phys. Rev. Lett. , 105:238102, 2010.[73] S. Thampi, R. Golestanian, and J. Yeomans. Active nematic materials with substrate friction.
Phys.Rev. E , 90:062307, 2014.[74] H. Masoud and M. Shelley. Collective surfing of chemically active particles.
Phys. Rev. Lett. , 112:128304,2014.[75] E. Hinch and L. G. Leal. Constitutive equations in suspension mechanics. part 2. approximate forms fora suspension of rigid particles affected by brownian rotations.