Inertia drives a flocking phase transition in viscous active fluids

TL;DR

This paper investigates how inertia influences phase transitions in viscous active fluids, revealing a critical dimensionless parameter that governs the transition from turbulence to flocking, with implications for mesoscale swimmer suspensions.

Contribution

It introduces a dimensionless parameter R that predicts a continuous flocking transition driven by inertia in viscous active suspensions.

Findings

Inertia induces a flocking phase transition in viscous active fluids.

The transition is continuous with a growing correlation length near the threshold.

Numerical simulations confirm the transition and steady states in 2D and 3D.

Abstract

How fast must an oriented collection of extensile swimmers swim to escape the instability of viscous active suspensions? We show that the answer lies in the dimensionless combination $R=\rho v_0^2/2\sigma_a$, where $\rho$ is the suspension mass density, $v_0$ the swim speed and $\sigma_a$ the active stress. Linear stability analysis shows that for small $R$ disturbances grow at a rate linear in their wavenumber $q$, and that the dominant instability mode involves twist. The resulting steady state in our numerical studies is isotropic hedgehog-defect turbulence. Past a first threshold $R$ of order unity we find a slower growth rate, of $O(q^2)$; the numerically observed steady state is {\it phase-turbulent}: noisy but {\it aligned} on average. We present numerical evidence in three and two dimensions that this inertia driven flocking transition is continuous, with a correlation length that grows on approaching the transition. For much larger $R$ we find an aligned state linearly stable to perturbations at all $q$. Our predictions should be testable in suspensions of mesoscale swimmers [D Klotsa, Soft Matter \textbf{15}, 8946 (2019)].

Tables

Table 1. Table 1: Spatial dimension D 𝐷 D of the domain and parameters ℒ , N , v 0 , K ℒ 𝑁 subscript 𝑣 0 𝐾 \mathcal{L},~{}N,~{}v_{0},~{}K , and R 𝑅 R used in our direct numerical simulations. The suspension density ρ = 1 𝜌 1 \rho=1 , λ = 0.1 𝜆 0.1 \lambda=0.1 , μ = 0.1 𝜇 0.1 \mu=0.1 and the rotational mobility Γ = 1 Γ 1 \Gamma=1 are kept fixed for all the runs. Note that as R 𝑅 R approaches R 2 subscript 𝑅 2 R_{2} the range of linearly unstable modes shrinks and is restricted to small wave-numbers, i.e., large length-scales. To resolve these unstable modes as well as the small-scale fluctuations that arise because of the nonlinear couplings, for R = 8 𝑅 8 R=8 we use a square domain with each side of length 160 ​ π 160 𝜋 160\pi and discretize it with 4096 2 superscript 4096 2 4096^{2} collocation points.

	D	$ℒ$	$N$	$v_0(\times {10}^{-2})$	$K(\times {10}^{-3})$	$R\equiv \rho v_0^2/2{\sigma }_a$
SPP1	3	$2\pi$	$128$	$3.16$	$1$	$0.02,0.0625$
SPP2	3	$10\pi$	$160$	$0.7,13.4$	$2$	$0.01,4$
SPP3	3	$10\pi$	$320$	$3.16$	$1$	$0.1-2$
SPP4	2	$20\pi$	$1024$	$3.16$	$1$	$0.05-2.0$
SPP5	2	$32\pi$	$1024$	$3.16$	$1$	$0.15,0.20$
SPP6	2	$40\pi$	$2048$	$3.16$	$1$	$0.25,0.30,0.35$
SPP7	2	$64\pi$	$3072$	$3.16$	$1$	$0.4-0.6$
SPP8	$2$	$80\pi$	$4096$	$3.16$	$1$	$0.7$
SPP9	$2$	$80\pi$	$2048$	$3.16$	$1$	$1.25$
SPP10	$2$	$160\pi$	$4096$	$3.16$	$1$	$8.0$
SPP11	$2$	$128\pi$	$8192$	$3.16$	$1$	$0.01$

Equations

$$\rho ({\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u})$$

$${\partial }_t\mathbf{p}+(\mathbf{u}+v_0\mathbf{p})\cdot \nabla \mathbf{p}$$

$${\partial }_{t}c+\nabla \cdot \left[(\mathbf{u}+v_1\mathbf{p})c\right]$$

$${\mathbf{\Sigma }}^a\equiv -{\sigma }_a\mathbf{p}\mathbf{p}-{\gamma }_a(\nabla \mathbf{p}+\nabla {\mathbf{p}}^T)$$

$${\mathbf{\Sigma }}^r=\frac{1-\lambda }{2}\mathbf{p}\mathbf{h}-\frac{1+\lambda }{2}\mathbf{h}\mathbf{p}-\ell (\nabla \mathbf{h}+\nabla {\mathbf{h}}^T)$$

$$F=\int d^{3}r\left[\frac{1}{4}(\mathbf{p}\cdot \mathbf{p}-1{)}^2+\frac{K}{2}(\nabla \mathbf{p}{)}^2-E\mathbf{p}\cdot \nabla c\right]$$

$${\text{T}}_{\mathbf{q}}\equiv \text{I}-\hat{\mathbf{q}}\hat{\mathbf{q}}$$

$$\left(\rho {\partial }_t+\mu q^2\right)\delta {\mathbf{u}}_{\perp \mathbf{q}}$$

$${\partial }_t\delta {\mathbf{p}}_{\perp \mathbf{q}}$$

$$\omega ={\omega }_{\pm }^s=\frac{1}{2}{v}_{0}q\cos \varphi -i\frac{{\mu }_{+}}{2\rho }{q}^2\pm {\left(\frac{{\sigma }_a}{2\rho }\right)}^{1/2}{\left[A(\varphi )q^2+iB(\varphi )q^3+G(\varphi )q^4\right]}^{1/2}$$

$${\omega }_{\pm }^t=\frac{1}{2}{v}_{0}q\cos \varphi -i\frac{{\mu }_{+}}{2\rho }{q}^2\pm {\left(\frac{{\sigma }_a}{2\rho }\right)}^{1/2}{\left[A(0){\cos }^2\varphi q^2+iB(0)\cos \varphi q^3+\bar{G}(\varphi )q^4\right]}^{1/2}$$

$$R\equiv \rho v_0^2/2{\sigma }_a,$$

$$\omega ={\omega }_{\pm }^s=\frac{q}{2}\left\{v_0\cos \varphi \pm {\left[\frac{2{\sigma }_a}{\rho }A(\varphi )\right]}^{1/2}\right\}-\frac{i}{2}\frac{\mu }{\rho }{q}^2\left\{1+\beta \mp (1-\beta ){\left[\frac{R{\cos }^2\varphi }{A(\varphi )}\right]}^{1/2}\right\}$$

$$\omega ={\omega }_{\pm }^t=\frac{q}{2}\cos \varphi \left\{v_0\pm {\left[\frac{2{\sigma }_a}{\rho }A(0)\right]}^{1/2}\right\}-\frac{i}{2}\frac{\mu }{\rho }{q}^2\left\{1+\beta \mp (1-\beta ){\left[\frac{R}{A(0)}\right]}^{1/2}\right\}$$

$$R_2=\frac{{\mu }_{+}^2}{{\mu }_{+}^2-{\mu }_{-}^2}{R}_1=\frac{1+\lambda }{4\beta }(1+\beta {)}^2.$$

$${\ell }_v\equiv \mu /v_0\rho \text{mbox}and{\ell }_{\sigma }\equiv \mu /\sqrt{\rho {\sigma }_a}=R^{1/2}{\ell }_v$$

$${\ell }_K\equiv \sqrt{K/{\sigma }_a}$$

$${\omega }^s$$

$${\omega }^t$$

$$E_u(q)$$

$$E_p(q)$$

$$\left(\rho {\partial }_t+\mu q^2\right)\delta {\mathbf{u}}_{\perp \mathbf{q}}$$

$${\partial }_t\delta {\mathbf{p}}_{\perp \mathbf{q}}$$

$$({\partial }_t+i{v}_1{q}_x)\delta c_{\mathbf{q}}$$

$$\rho ({\partial }_t\mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{u})$$

$${\partial }_t\mathbf{p}+(\mathbf{u}+v_0\mathbf{p})\cdot \nabla \mathbf{p}$$

$${\omega }_{\pm }^t=\frac{1}{2}{v}_{0}q\cos \varphi -i\frac{{\mu }_{+}}{2\rho }{q}^2\pm {\left(\frac{{\sigma }_a}{2\rho }\right)}^{1/2}{\left[A(0){\cos }^2\varphi q^2+iB(0)\left(1-\frac{{\gamma }_a(1+\lambda )+2\ell {\sigma }_a}{v_0{\mu }_{-}}\right)\cos \varphi q^3+\left(\bar{G}(\varphi )+\frac{2{\gamma }_a\ell }{{\sigma }_a}\right)q^4\right]}^{1/2}$$

$${\omega }_{\pm }^t=\frac{1}{2}{v}_{0}q\cos \varphi -i\frac{{\mu }_{+}}{2\rho }{q}^2\pm {\left(\frac{{\sigma }_a}{2\rho }\right)}^{1/2}{\left[A(0){\cos }^2\varphi q^2+iB(0)\cos \varphi q^3+\bar{G}(\varphi )q^4\right]}^{1/2}.$$

$$D_1={\gamma }_a/\mu v_0,D_2=\ell {\sigma }_a/\mu v_0$$

Full text

I Introduction

The theory of active matter [1, 2, 3, 4, 5, 6, 7, 8, 9] – systems whose constituents convert a sustained supply of fuel into movement – is the framework of choice for understanding the collective behaviour of motile particles. Like condensed matter in general, active systems display many types of order and operate in a variety of dynamical regimes. Our interest in this paper is in groups of motile organisms in a bulk fluid medium, spontaneously organized into a flock in which their tail-to-head vectors on average point in a common direction. In now-standard terminology [2], we consider polar, wet active matter, described by a vector order parameter characterizing the degree and direction of common orientation and movement.

In the world of Stokesian hydrodynamics, where inertia is absent and viscosity holds sway, an ordered flock in bulk fluid is impossible [10]. Indeed this limit has shaped the defining image of active suspensions as inescapably unstable, dissolving via spontaneous flow [10, 11] and defect proliferation [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27] into a kind of turbulence [28, 29, 30, 19, 31, 32, 33, 34, 35, 36, 37]. We remind the reader of the instability mechanism [10, 1, 2]: an aligned state of active particles is a state of uniform uniaxial stress; perturbing this state – through bend or splay respectively for “pusher” or “puller” particles – creates spatially varying stresses; force balance requires that these are accompanied by flow; this flow rotates the alignment further in the direction of the perturbation. Note that this description refers neither to the directed motion of the particles nor to acceleration. Stresses and flows are uniaxial but apolar, that is, fore-aft symmetric, and flow responds instantaneously to active stress in the Stokesian approximation. Improved descriptions including polar order alone [10, 15, 16, 38], or inertia alone [10], do not mitigate the instability.

Stable flocks in bulk fluid are of course widely observed in the form of fish schools [39, 40], which are very polar and very far from Stokesian. We do not venture into the regime of schooling at high Reynolds number, governed by purely inviscid hydrodynamic interactions [41, 42], but consider weak inertial effects, which are know to alter significantly the viscous hydrodynamic interaction between slow swimmers [43, 44]. A recent Perspective [45] makes a persuasive case for the study of active fluids with small but non-negligible inertia, the regime we explore here. A result from [10] is relevant in this context: a linearized treatment retaining only acceleration and active stresses finds a parameter domain in which flocks in fluid are neutrally stable to first order in wavenumber, with a wavelike dynamic response. Interestingly, such waves of bend excitations have an analog in models without momentum conservation, i.e., “dry” flocking models, when rotational inertia is taken into account [46, 47, 48, 49]; the coupled dynamics of classical spin angular momentum and orientation in such models is formally similar to that of hydrodynamic vorticity and orientation in [10] and the present work. Staying close to viscous hydrodynamics and the force-dipole picture of swimmers, we ask: can inertia and polar order together defeat the Stokesian instability of flocks?

I.1 Summary of results

In this article, focusing on extensile or “pusher” [50] suspensions, we show that the introduction of inertia qualitatively alters our understanding of the viscous hydrodynamics of polar active matter, that is, flocking in fluids. Here are our main results. We show that speed matters: the dimensionless combination $R\equiv\rho v_{0}^{2}/2\sigma_{a}$, where $\rho$ is the suspension mass density, $\sigma_{a}$ the scale of active stress, and $v_{0}$ the self-advection speed, governs the stability of active suspensions. Flocks in fluid are stable for large $R$, and their inexorable Stokesian instability [10] is the $R=0$ limit of a far richer picture, Fig. 1. For small $R$ perturbations about the aligned state grow at a rate $\propto q$ for wavenumber $q\to 0$ while for moderate $R$ the linear instability persists but with a growth rate $\propto q^{2}$. Crucially, direct numerical simulations of the hydrodynamic equations reveal that the two regimes correspond to qualitatively distinct statistical steady states separated by a nonequilibrium phase transition. The small-$R$ regime is isotropic hedgehog-defect turbulence while that at moderate $R$ is a phase-turbulent [51, 52, 53, 54] but ordered flock. Our numerical results suggest a continuous order-parameter onset and a growing correlation length upon approaching the transition.

This paper is organized as follows. In § II we present the equations of hydrodynamics for polar active suspensions and investigate the linear stability of the uniaxially ordered state. § III describes our numerical studies, the flocking transition from defect to phase turbulence, and the properties of the turbulent states. We close in § IV with a summary, suggestions for experiment, and open questions.

II Governing equations and stability analysis

II.1 Hydrodynamics of active suspensions

We begin by constructing, from general principles, the hydrodynamic equations of motion for a flock in fluid. We do not employ the language of forces and fluxes or display the dependence of “active” coefficients on a maintained chemical driving force [2, 5, 55]. We adopt the general, symmetry-based approach of  [10] but our treatment is self-contained and does not presuppose familiarity with that work. The reader will see the results of [10] emerge as a limiting case in section II.2.2. We emphasize that our equations constitute a general effective description on length-scales much larger than a swimmer (as we shall call our self-propelled particles hereafter). They contain parameters such as viscosity and elastic constants; these are phenomenological coefficients in our coarse-grained description of this internally driven system, and are named based on the form of the terms they govern. Their values cannot in principle be estimated from a near-equilibrium hydrodynamic theory of the suspension. For example, we would imagine the viscosity in our equations receives “eddy” contributions from flows on scales of a few swimmers, and we expect that the elastic constants encoding the aligning tendency are at least partly behavioural rather than mechanical. Provided interactions are local in space and time, these features do not limit the validity of our approach, which depends only on conservation laws and symmetries.

For a steady state such as a flock, which spontaneously breaks a continuous invariance, the slow or hydrodynamic variables [56] are the local densities of conserved quantities and the broken-symmetry or Nambu-Goldstone [57, 58] fields. At a continuous transition to such an ordered state, the amplitude of the order parameter is an additional slow variable. In the absence of reproduction, death, and external forces, the conserved quantities are the total number of swimmers, the total amount of fluid, and the total momentum of swimmers plus fluid. Energy conservation does not play a role as each swimmers is endowed with a built-in power source. The slow variables corresponding to these conservation laws are then the densities $\rho$ and $\bm{g}=\rho\bm{u}$ of mass and momentum of swimmers plus fluid (defining the suspension velocity field $\bm{u}$), and the number density $c$ of swimmers. The broken-symmetry modes and the magnitude of order are jointly contained in the polar order parameter field $\bm{p}$, which is the local average of the orientation unit vectors of the particles [10]. It is interesting to note that an equilibrium liquid crystal with macroscopic vectorial, i.e., polar, order has only very recently been discovered [59].

As the order parameter is a space vector, our description is invariant under the joint inversion of $\bm{p}$ and the spatial coordinate $\bm{r}$, but not $\bm{p}$ alone. The absence of $\bm{p}\to-\bm{p}$ symmetry is central to our narrative. We therefore include at this stage, at leading order in gradients, all terms that break this symmetry in the equations of motion [60, 61] (see also [62, 63, 2, 15, 16, 64]), although we will shortly pass to a more economical description. The equations read

[TABLE]

Equation (3) expresses number conservation; the active particles self-propel with velocity $v_{1}\bm{p}$ in the frame of the suspension, and hence $\bm{u}+v_{1}\bm{p}$ in the laboratory frame. In (2), the polar order parameter $\bm{p}$ is carried by the hydrodynamic velocity $\bm{u}$ and by its self-advection $v_{0}$ [not related by any symmetry [65, 49] to $v_{1}$ in (3)], ${\bm{S}}$ and ${\bm{\Omega}}$ are the symmetric and antisymmetric parts of the velocity gradient tensor $\nabla{\bm{u}}$ which couple orientation to flow as in ordinary nematic liquid crystals [66], $\Gamma$ is the kinetic coefficient governing relaxation in the local molecular field $\bm{h}$ to be discussed below, and $\ell$ is the polar flow coupling at leading order in a gradient expansion. The $v_{0}$, $v_{1}$ and $\ell$ terms are polar: their presence implies that the equations are not invariant under $\bm{p}\to-\bm{p}$. In (1), the hydrodynamic pressure $P$ enforces incompressibility $\nabla\cdot{\bm{u}}=0$.

[TABLE]

is the intrinsic stress associated with swimming activity, which we display up to first subleading order in gradients. In (4) $\sigma_{a}$ is the force-dipole density [10, 67, 2]. In microscopic terms, the forces exerted by a swimmer and the ambient fluid on each other add to zero, so the associated force density has zero monopole moment. The minimal model for a swimmer is thus a point dipolar force density. A collection of such swimmers, each with dipole strength $W$, local concentration $c$, and mean local alignment given by the polar order parameter $\bm{p}$ can readily be seen [10, 67, 2] to have force density $-W\nabla\cdot(c\bm{p}\bm{p})$. $W>0$ and $W<0$ correspond respectively to extensile swimmers, that push fluid back with their tails and move forward, and contractile swimmers, that advance by pulling fluid toward themselves from the front. Thus, in (4), $\sigma_{a}=Wc$. The polar contribution to the active stress, given at leading order in a gradient expansion by the $\gamma_{a}$ term in (4), arises [62, 67] if the force dipole on each particle is displaced with respect to the center of drag of the particle, as it must be to achieve locomotion. From the foregoing it is plausible that $\gamma_{a}$ should be proportional to $\sigma_{a}$, with the proportionality factor being a length that measures the fore-aft asymmetry of the active particles. In principle all parameters in our equations should be functions of the local concentration $c$. Phenomenologically, $\sigma_{a},\,\gamma_{a}$ are tied to the presence of active particles and should therefore be proportional to $c$ for $c\to 0$. In the simple microscopic picture discussed above, if $W$ is treated as an intrinsic single-particle property and therefore independent of $c$, the proportionality is exact. The contribution

[TABLE]

is the reversible thermodynamic stress for an equilibrium polar liquid crystal. The expression (5) extends the form found in  [68, 2] to include the leading-order polar $\ell$ term, which is the Onsager counterpart of the polar flow-coupling term $\ell\nabla^{2}\bm{u}$ in (2), discussed in [60] and [61]. $\bm{h}=-\delta F/\delta{\bm{p}}$ is the molecular field conjugate to ${\bm{p}}$, derived from a free-energy functional

[TABLE]

favouring a ${\bm{p}}$-field of uniform magnitude [69] which we have rescaled to unity. A single Frank constant [70, 71, 72] $K$ penalizes gradients in ${\bm{p}}$, and $E$ promotes alignment of ${\bm{p}}$ up or down gradients of $c$, according to its sign. $\mu$ is the shear viscosity of the suspension, and $\Gamma$, the collective rotational mobility for the relaxation of the polar order parameter field, is expected to be of order $1/\mu$. $\lambda$ is the nematic flow-alignment parameter [73, 74] and $\ell$, with units of length, governs the lowest-order polar flow-coupling term [60, 61]. If our equations were derived from a microscopic model of particles in a fluid, we expect that both $\ell$ and the length $\gamma_{a}/\sigma_{a}$ would be related to a fore-aft asymmetry in the dimensions of the active particles.

II.1.1 Essential and incidental polar contributions

Equations (1)-(3) are endowed with a surfeit of parameters originating in the polar character of our system – the speeds $v_{0}$ and $v_{1}$ at which the orientation advects itself and the concentration respectively, the polar active stress coefficient $\gamma_{a}$ and the passive polar flow-coupling length-scale $\ell$. In the work [10] that initiated the study of the hydrodynamics of active liquid crystals, the polar character of the order parameter of a flock played an important role, combining with inertia to yield a propagative mode-structure. As this was a leading-order feature in a gradient expansion, neither $\gamma_{a}$ nor $\ell$, which enter at next-to-leading order relative to $\sigma_{a}$ in (4) and $1\pm\lambda$ in (5) and (2) respectively, were considered, and polar effects thus entered [10] only through $v_{0}$ and $v_{1}$. The instability of active liquid crystals in the Stokesian domain – the other major finding of [10] – commanded much greater interest in the field thanks to its connection to experimental realizations in cellular and microbial settings. An analysis that ignores polarity altogether and works only with the axis of orientation offers a satisfactory conceptual understanding of that instability [2, 30], though interesting complexities arise [15] in a Stokesian setting through the polar parameters $\gamma_{a}$, $v_{0}$ and $v_{1}$. A final remark in this context is that polarity will assert its presence in any formulation in terms of a vector order parameter $\bm{p}$, even if the equations of motion are invariant under $\bm{p}\to-\bm{p}$, through the nature of topological defects [20, 75].

For the purposes of the present work what matters is that the self-advection speed $v_{0}$ plays a distinct – and crucial – role. The other polar parameters contribute in an incidental manner. First, we are concerned here only with the extensile case $\sigma_{a}>0$, for which the instability mode is bend, which decouples from concentration in the linear theory. In what follows we therefore ignore the concentration field, and hence $v_{1}$ drops out of our analysis. Next, as we show in the Appendix, it is only through $v_{0}$ – in the form $\rho v_{0}^{2}$ and its competition with $\sigma_{a}$ – that the stabilizing effects of inertia enter our treatment. $\gamma_{a}$ and $\ell$ leave unaltered both the coefficient of the $O(q)$ contribution to the mode frequency and the parameter value at which the instability growth rate changes from $O(q)$ to $O(q^{2})$. They simply shift the coefficients of the $O(q^{2})$ piece of the mode frequency by amounts of relative order unity. We therefore work with an economical description in which $\gamma_{a}$ and $\ell$ are zero and polar effects enter only through $v_{0}$ and, of course, the nature of the allowed topological defects. Crucially, $v_{0}$ and $\sigma_{a}$ are independent quantities in our coarse-grained treatment, a point we will return to later in the paper.

II.2 Linear Stability analysis

Defining the ordering direction to be $\hat{\bm{x}}$ and directions in the $yz$ plane as $\perp$, we have investigated the stability of a uniform ordered flock ($c=c_{0},{\bm{u}}={\bm{0}}$, and ${\bm{p}}=\hat{\bm{x}}$, which is a stationary solution of Eqs.(1) - (3) to small perturbations ($\delta{\bm{u}}_{\perp},\delta{\bm{p}}_{\perp},\delta c$), where the presence of only the $\perp$ components is a result of incompressibility and the “fast” nature of $p_{x}$. We present here the results for the case where the concentration field $c$ is removed from the analysis. This is sufficient for our purposes, because $c$ does not participate significantly in the linear instabilities of relevance, as we now argue. Taking the curl with respect to $\nabla_{\perp}$ eliminates $c$ from the $\perp$ component of Eq. (2). A similar curl removes it from the $\perp$ component of Eq. (1) as well. Thus concentration does not participate in the linear dynamics of the twist-bend mode [10]. The 3-dimensional instability of extensile active fluids is known, numerically, to be twist-dominated [76], an observation for which our linear stability analysis below provides the natural explanation. A description without a concentration field should thus be a reasonable guide to instabilities and active turbulence in our system. It is important to note that the neglect of the concentration field in our treatment does not amount to an incompressibility constraint on the polar order parameter field. A formal connection between the complete equations and those without a concentration field can be achieved by introducing birth and death of particles so that $c$ becomes “fast” [77] and can be eliminated in favour of the slow variables $\bm{p}_{\perp}$ and $\bm{u}_{\perp}$, with at most a finite shift in parameter values in the equations for the slow variables. Changes in the linear stability analysis upon inclusion of the concentration are quantitative, not qualitative, and can be found in the Appendix. Defining the projector

[TABLE]

transverse to $\bm{q}$ and linearizing Eqs. (1) and (2) about the ordered state we find

[TABLE]

As in [10], the divergence and curl of Eqs. (8) and (9) describe respectively the dynamics of splay and twist, with an admixture of bend in each case for $q_{x}\neq 0$. Defining $\phi$ to be the angle between the wavevector $\bm{q}$ and the alignment ($\hat{\bm{x}}$) direction, the resulting dispersion relations for the frequency $\omega$, valid for all $\bm{q}$, for modes of the form $e^{i(\bm{q}\cdot{\bm{r}}-\omega t)}$, are

[TABLE]

for the splay-bend modes and

[TABLE]

for the twist-bend modes. In Eqs. (10) and (11) we have defined $A(\phi)\equiv R\cos^{2}\phi-\cos 2\phi(1+\lambda\cos 2\phi)$, $B(\phi)\equiv(v_{0}\mu_{-}/\sigma_{a})\cos\phi$, $G(\phi)=-(\mu_{-}^{2}/2\rho\sigma_{a})+(K/2\sigma_{a})(1+\lambda\cos 2\phi)^{2}$, and $\bar{G}(\phi)\equiv-(\mu_{-}^{2}/2\rho\sigma_{a})+(K/2\sigma_{a})(1+\lambda)^{2}\cos^{2}\phi$ [78] with

[TABLE]

and $\mu_{\pm}\equiv{\mu}(1\pm\beta)$ where $\beta\equiv\Gamma K\rho/\mu$ should be of the same order as $\alpha\equiv K\rho/\mu^{2}$ because the mobility $\Gamma\sim 1/\mu$. For conventional liquid crystals $\alpha,\,\beta\ll 1$.

When $R=0$, the extensile ($\sigma_{a}>0$) systems of interest here present a bend instability [see Eqs. (10) and (11)] with invasion speed $\sqrt{\sigma_{a}/\rho}$. For $v_{0}>0$ disturbances can outrun this invasive growth. The dimensionless combination $R$ describes this competition. Note that the contribution of $R$ vanishes for pure splay, Eq. (10) at $\phi=\pi/2$, so motility cannot stabilize contractile ($\sigma_{a}<0$) flocks in fluid.

II.2.1 Small-$q$ behaviour: the $O(q)$ and $O(q^{2})$ instabilities

Let us first examine the small-$q$ behaviour. Expanding Eqs. (10) and (11) up to order $q^{2}$ we then find

[TABLE]

for the splay-bend modes and

[TABLE]

for the twist-bend modes. Here $A(0)=A(\phi=0)=R-(1+\lambda)$. One note of caution: the small-$q$ expansion that led to Eq. (14) assumes $v_{0}q\cos\phi>q^{2}\mu/\rho$, which means that it does not apply for $\phi=\pi/2$, i.e., pure twist. It does however hold for any $\phi\in[0,\pi/2)$ but the closer $\phi$ is to $\pi/2$ the smaller $q$ must be for the result to apply.

Two of our main results now follow. If $R<1+\lambda$, Eq. (14) signals a bend instability with small-$q$ growth rate $\sim q$. This was discussed in the strictly apolar case $v_{0}=0$ in [10], and can be viewed as the small-$q$ extension of the Stokesian bend instability [10]. However, if $R>1+\lambda$, so that the $O(q)$ instability is averted, $0<1-(1+\lambda)/R<1$. If $R$ is not too large, this means the coefficient of $iq^{2}$ in Eqs. (13) and  (14) is positive, signalling a small-$q$ instability with diffusive growth. This $O(q^{2})$ instability exists for $R$ between $R_{1}=1+\lambda$ and

[TABLE]

For $R>R_{2}$ the flock is linearly stable. If $\beta\ll 1$ as in molecular systems, $R_{2}\gg R_{1}$, and the $O(q^{2})$ instability occupies a large range of $R$. In the $\beta=0$ limit the uniformly ordered flock is always linearly unstable, with small-$q$ growth rate $\sim q$ for $R<1+\lambda$ and $\sim q^{2}$ for $R>1+\lambda$. Fig. 1 summarizes the small-$q$ stability behaviour.

Note that the $O(q^{2})$ instability can be eliminated in the special case $\beta=1$, i.e., $\mu/\rho=\Gamma K$. Noting that $\Gamma$ should be roughly $1/\mu$, this condition implies $K=\mu^{2}/\rho$, an interesting condition that equates a Frank constant (which, recall, has units of force in three dimensions) to Purcell’s intrinsic force scale [79] $\mu^{2}/\rho$ for three-dimensional viscous fluids. As we remarked above, $\beta$ in molecular or colloidal systems is about $10^{-4}$ [80, 66], so requiring it to be order unity amounts to insisting that the swimmers have an exceptionally strong aligning interaction. This possibility cannot be ruled out a priori as alignment in living systems is likely to be active and behavioural, not a passive mechanical torque.

II.2.2 Large-$q$ dynamics and the Stokesian limit

Having established the general linearized behaviour of active extensile liquid crystals in the true hydrodynamic regime of small wavenumbers, we turn our attention to large wavenumbers. Two length scales are important here –

[TABLE]

below which viscosity overwhelms the inertial effects of self-advection and

[TABLE]

below which Frank elasticity dominates active stresses. For molecular or colloidal systems, for which, as we remarked earlier, $\alpha$ is exceedingly small [66, 80], ${\ell_{K}/\ell_{v}}=\sqrt{\alpha R}$ should be small too except in the unlikely condition of ultra-high self-advection speeds. The wavenumber range $\mbox{max}(\ell_{v}^{-1},\ell_{\sigma}^{-1})\ll q\ll\ell_{K}^{-1}$ should be substantial. Expanding Eqs. (10) and (11) for $q\gg\mbox{max}(\ell_{v}^{-1},\ell_{\sigma}^{-1})$, we find, to leading order in $\alpha$ and $\beta$, that the splay-bend mode that goes unstable at small $R$ has the form

[TABLE]

and the corresponding twist-bend mode has frequency

[TABLE]

where $A(\phi)$ is as defined in Eqs. (10) and (11). The Stokesian instability of active liquid crystals [10, 2, 11], with a single growth-rate scale $\sigma_{a}/\mu$, emerges from Eqs. (18), (19) if $\rho$ is set to zero. For $\sigma_{a}<0$, contractile or pusher suspensions, which is not the case we are focusing on in this work, Eq. (18) predicts an instability for $\phi>\pi/4$, which is splay-dominated. For $\sigma_{a}>0$, the extensile or pusher case, Eq. (18) predicts an instability for $\phi<\pi/4$, which is bend-dominated. More important, Eq. (19) predicts an bend instability for all directions other than pure twist $\phi=\pi/2$, although of course the value of $q$ below which the instability is seen approaches [math] as $\phi\to\pi/2$. In general, however, Eqs. (18) and (19) are not Stokesian expressions but short-wavelength limits of the linearized dynamics of a polar active suspension with inertia, which enters through $R$. We see in particular that the stability criteria in this large-$q$ regime are identical to those for the $O(q)$ mode at small $q$. Thus a twist-bend instability, with a growth rate $\sim\sigma_{a}/\mu$ for $\mbox{max}(\ell_{v}^{-1},\ell_{\sigma}^{-1})\ll q\ll\ell_{K}^{-1}$ takes place if $R<1+\lambda$. This establishes our claim that the $O(q)$ instability is the small-$q$ extension of the Stokesian instability [10] of active suspensions. The $O(q^{2})$ instability that intervenes at small $q$ as $R$ is increased does not reflect itself in the large-$q$ dynamics. Fig. 2 displays the growth or decay rates of the twist-bend mode as a function of wavenumber as $R$ is varied at $\beta=10^{-4}$.

It is important to keep in mind that the active stress $\sigma_{a}$ is a partial description of the mechanics of self-propulsion based on an estimate of the force-dipole concentration, and is not a priori determined by $v_{0}$. To take an extreme case, Stokesian swimmers with no force dipole exist, e.g., the pure quadrupole [79, 81, 82] Assuming a volume fraction of order unity, let us nonetheless try to estimate $R$ for typical swimmers of speed $v_{0}$ (the distinction between the speeds of self-propulsion and self-advection being unimportant for this discussion) and size $b$ (although we must remember that this size is notional in our coarse-grained description). For Reynolds number Re small at the scale of the individual organism it is plausible that $\sigma_{a}\sim\mu v_{0}/b$. In that case $R\equiv\rho v_{0}^{2}/2\sigma_{a}\sim\rho v_{0}b/\mu=\mbox{Re}\ll 1$, so we can replace Eqs. (18) and (19) by their Stokesian approximations. For swimmers at nonzero Reynolds number it is less obvious how to estimate $\sigma_{a}$. If we take it still to be a viscous stress then $R=\mbox{Re}$ continues to hold, so now $R$ dominates in Eqs. (18) and (19), or in Eqs. (10) and (11), guaranteeing stability. Even if $\sigma_{a}\sim\rho v_{0}^{2}$, $R\sim 1$ and it is plausible that the instability is averted [83].

Dominance of twist in the three-dimensional extensile instability – A noteworthy feature, to our knowledge not discussed in the literature, emerges in our three-dimensional analysis: there are two families of bend instability – mixed with splay as in Eqs. (13) and (18) and twist as in Eqs. (14) and (19). Interpolation with bend mitigates the instability in Eqs. (13) and (18), crossing over to stability for large enough $\phi$, but twist in Eqs. (14) and (19) has no such effect. The twist-bend instability Eqs. (14) and (19) should thus dominate, as it occurs for all $\phi$ except precisely $\pi/2$. This abundance of twisted unstable modes in Eqs. (14) and (19), independent of the roles of polarity and inertia, is doubtless the explanation of the numerical observations of Shendruk et al. [76] in their study of three-dimensional extensile active nematics.

We summarise this section by noting that, when inertia is taken into account, orientable active suspensions can have two types of linear instability at small wavenumber $q$, governed by the dimensionless parameter $R$ [Eq. (12)].The instability growth-rates are $O(q)$ for $R<R_{1}=1+\lambda$ where $\lambda$ is a flow-alignment parameter and $O(q^{2})$ for $R_{1}<R<R_{2}\sim R_{1}/\beta$ where $\beta$ is defined in Eqs. (10) and (11). Linearly stable behaviour is found for $R>R_{2}$. As $\beta\sim 10^{-4}$ in molecular systems, the $O(q^{2})$-unstable regime occupies a rather large range in parameter space. Indeed one could argue that the typical behaviour is that corresponding to the $\beta=0$ limit, in which the aligned state is always linearly unstable, either at $O(q)$ or at $O(q^{2})$. In §III we gain insight beyond this linear analysis through a detailed numerical study to discover the long-time fate of the system in these unstable regimes.

III Numerical studies of active hydrodynamics with inertia

In the following section we describe detailed numerical solutions of our equations, with an emphasis on the changes in behavior as $R$ is varied. We first verify the predictions of our linear stability analysis. Next, for extensile suspensions, we reveal an inertia-driven nonequilibrium phase transition from a disordered defect-turbulent state for $0<R<R_{1}$ to an ordered phase-turbulent state for $R_{1}\leq R<R_{2}$. We characterize these using the polar order parameter, that is, the macroscopic steady-state average of $\bm{p}$, correlation functions and energy spectra.

III.1 Direct Numerical Simulations (DNS)

We numerically integrate Eqs. (1) and (2) in square and cubic domains of volume $\mathcal{L}^{d}$ in dimensions $d=2$ and $3$. Spatial discretisation, with $N^{d}$ collocation points, is conducted by employing a pseudo-spectral method [84] for Eq. (1) and a fourth-order central finite-difference scheme for Eq. (2). For temporal integration we use a second-order Adams-Bashforth scheme [85]. Consistent with the linear stability analysis conducted earlier, we choose a uniform ordered state with transverse monochromatic perturbation as the initial condition, i.e. $\bm{u}=\bm{0}+A\hat{\bm{e}}_{\perp}\cos\bm{q}\cdot{\bm{r}}$, $\bm{p}=\bm{\hat{x}}+B\hat{\bm{e}}_{\perp}\cos\bm{q}\cdot{\bm{r}}$ where $\hat{{\bm{e}}}_{\perp}\equiv(\hat{\bm{y}}+\hat{\bm{z}})/\sqrt{2}$ is a unit vector in the plane perpendicular to the ordering direction, and we have made the arbitrary but acceptable choice $A=B=10^{-3}$.

We monitor the time-evolution of perturbations and, in the turbulent steady state, investigate the statistical properties of the velocity and the director fields. In Table 1, we summarize the parameters used in our DNS.

III.2 Initial growth of instabilities

We now present a comparison between the short-time growth obtained from the DNS with the analytical predictions of the linear stability analysis. The plot of the bend-twist dispersion curve given by Eq. (11) for $\phi=55^{\circ}$ is shown in Fig. 3. The black dots indicate the the initial temporal growth rate of perturbations obtained from our DNS, which shows excellent agreement with the analytical results. Furthermore, our simulations correctly capture the exponential and oscillatory characters of the growth for $R<R_{1}$ and $R_{1}<R<R_{2}$ respectively. Note that for $R<R_{1}$, the exponential growth rate of perturbations is much faster than the oscillatory kinematic contribution ${\rm Re}[\omega]=v_{0}q\cos\phi$. For $R_{1}<R<R_{2}$, $\rm{Re}[\omega]$ has contributions from both the kinematic and the inertial terms. Therefore, we observe an exponential growth of $|\bm{q_{\perp}}\times\bm{\delta p}_{\perp{\bm{q}}}|$ for $R<R_{1}$ [see Fig. (3A)], but oscillatory growth for $R_{1}<R<R_{2}$ [see Fig. (3B)].

III.3 A flocking phase transition

We now investigate the morphology and statistical properties of the orientation and flow emerging from the instabilities discussed above. Fig. (5) shows the typical flow structures observed in our DNS with increasing $R$ in the statistically steady state. For $0<R<R_{1}$, we observe hedgehog defects. The inter-defect spacing grows with increasing $R$. Unexpectedly, when $R$ increases past the first threshold $R_{1}$, a fluctuating but on average aligned state emerges. As we remarked in the Introduction, this is clear numerical evidence that $R=R_{1}$ marks a nonequilibrium phase transition from a statistically isotropic state to a flock or, in the terminology of spatiotemporal chaos, from defect turbulence to phase turbulence [86, 87, 51]. In the latter state long-wavelength statistical variation of the broken-symmetry variable is present but the amplitude of the order parameter is not destroyed by defects. We have not, however, measured the system-size dependence of the positive Lyapounov spectrum to establish spatiotemporal chaos quantitatively. We do not know the mechanism that serves to preserve macroscopic flocking order despite the $O(q^{2})$ instability. It appears that the growing amplitude of perturbations at small wavenumber $q$ triggers nonlinear effects which couples to large $q$ where the dynamics is stable. The behavior is reminiscent of that reported by Jayaprakash et al. [87] for the Kuramoto-Sivashinsky (KS) equation. The KS equation is a deterministic partial differential equation (PDE) with a negative diffusivity and hence a linear instability with growth rate $\propto q^{2}$ at small wavenumber $q$, peaking at a wavenumber $q_{*}$, stable behavior at large $q$ thanks to terms at higher order in $q$, and a nonlinearity that transfers weight from small to large $q$. Hayot et al. [87] carry out a numerical coarse-graining, i.e., a spatial low-pass filtering, on the two-dimensional KS equation to show that the effective equations of motion for the modes with $q<q_{*}$ are those of a stochastic PDE with a positive diffusivity. It is possible that such a mechanism is at work in our case, but to settle this issue will require a treatment analogous to that of [87] for our substantially more complicated equations. .

We now focus on the properties of the nonequilibrium phase transition. In Fig. 4(a), we plot the magnitude $|\langle{\bm{p}}\rangle|$ of the polar order parameter in the statistically steady state with increasing $R$, where angle brackets $\langle\cdot\rangle$ denote spatio-temporal averaging. For $R<R_{1}$, $|\langle{\bm{p}}\rangle|$ is consistent with zero. We observe an onset of polar order once $R$ increases beyond $R_{1}\equiv 1+\lambda$. Fig. 4(a) shows the order parameter for the largest system sizes studied; at large $R$ the values at half that system size are very similar. However, a detailed finite-size scaling analysis needs to be undertaken to find the correct scaling near the critical region [88].

In the defect-turbulence regime, we study the steady-state correlation function $C(r)=\langle{\bm{p}}({\bm{x}}+\bm{r})\cdot{\bm{p}}({\bm{x}})\rangle/\langle{{\bm{p}}(0)}^{2}\rangle$, where the angular brackets denote spatial averaging. We plot the correlation function $C(r)$ versus $r$ in Fig. 4(b) and evaluate the correlation length by fitting an exponential decay $\exp(-r/\xi)$ to the numerical data 111The value of $\xi$ obtained from the fit is comparable to the one obtained using the definition $\xi\equiv\int_{0}^{\mathcal{L}/2}C(r)dr$.. We see that the correlation functions for different values of $R<R_{1}$ fall on a single curve plotted against $r/\xi$. Moreover, from Fig. 4(c), $\xi$ grows and possibly diverges as $R\to R_{1}$; our limited data points are consistent with an exponent of unity. Further progress requires finite-size scaling studies and measurements of order-parameter correlations at asymptotically small wavenumber [88] for $R>R_{1}$ to test the nature of the ordered state.

III.4 Energy spectrum

A state of complex, correlated but disorderly flow is seen in a wide variety of suspensions of motile organisms and motorized biofilaments. It has been termed active turbulence and analyzed through the study of energy spectra as in conventional turbulence [90, 16, 34, 32, 29, 91, 17, 18, 26, 92, 93, 36]. However, these studies have all considered systems with negligible inertia. Here we examine numerically the spatial power spectral densities for the polar order parameter and the hydrodynamic velocity field, in the defect- and phase-turbulent regimes – the latter owing its existence to inertia. In keeping with typical turbulence studies, we use the shell-averaged energy spectra of the velocity and the order parameter

[TABLE]

where ${\bm{u}}_{\bm{m}}$ and ${\bm{p}}_{\bm{m}}$ are the Fourier coefficients of the velocity ${\bm{u}}$ and order parameter ${\bm{p}}$ fields.

Among the features of interest are Porod’s Law regimes corresponding to the fields of topological defects. In addition, for $R<R_{1}$ we find velocity correlations of Ornstein-Zernike form, with a correlation length much larger than that of the order parameter, whose origin we discuss below. For $R>R_{1}$, the phase-turbulent but ordered state, we present preliminary evidence of fluctuations of the broken-symmetry or Nambu-Goldstone mode. The behaviours of $E_{p}(q)$ and $E_{u}(q)$ for a range of values of $R$ are displayed in Figs. 6 and 7.

Energy spectra of the order parameter– We observe that for $R<R_{1}$, the spectrum peaks around $q\xi\sim 1$. At moderate values of $R$, $0.1<R<1$, because of exponential decay in orientational correlations, we expect $E_{p}(q)\sim 1/(1+\xi^{2}q^{2})$ for $q\xi>1$. On the other hand, for $R\ll 1$, because of defects, we expect a Porod’s scaling [94] $E_{p}(q)=q^{-3}$ in two dimensions and $E_{p}(q)=q^{-4}$ in three dimensions for $q\xi\gg 1$. Recent studies on dry active matter [75] using a scale-by-scale budget analysis revealed that, even in the presence of defects, the nonlinear transfer mechanisms could lead to a non-Porod scaling. Unfortunately, we do not have sufficient scaling range – due to modest grid-resolution at $R\ll R_{1}$ – to undertake such analysis. We emphasize that the quoted exponent values are empirically determined by conservatively selecting an appropriate dynamic range of wavenumbers away from the smallest ($\sim 1/L$) and the largest, viz., $q_{K}\equiv 2\pi/\ell_{K}$, beyond which Frank elasticity dominates.

In the phase-turbulent regime, $R_{1}<R<R_{2}$, we observe $E_{p}(q)\sim q^{-3}$ for $R$ close to $R_{1}$. As $R$ approaches $R_{2}$, the range of linearly unstable modes shrinks and is restricted to wave-numbers close to large scales [small $q$, see Fig. 6(b)]. For the linearly unstable modes we observe a weak $q$-dependence, whereas for wave-numbers outside the linearly unstable regime, the nonlinearities lead to a transfer of order-parameter fluctuations to small scales with a power-law spectrum $E_{p}(q)\sim q^{-3.8}$ [see Figs. 6(b) and (d)].

Energy spectra of the velocity – In Eq. (1) we expect the dominant balance to be between acceleration and activity as the Reynolds number obtained by comparing the advective and viscous terms, based on the root-mean-square hydrodynamic velocity, is small ($\mbox{Re}\equiv\rho u_{rms}\xi/\mu\leq 0.5$) 222We have verified that the results of our DNS do not change if the advective nonlinear term in Eq. (1) is absent.. We therefore expect for small $q$, $\omega u_{q}\sim\sigma_{a}q\sum_{k}{\bm{p}}_{k}{\bm{p}}_{q-k}$. If we assume that the dominant contribution to the convolution comes from terms with $|{\bm{k}}|=|{\bm{q}}-{\bm{k}}|=q$, i.e., on the same shell in Fourier space, we get $\omega u_{q}\sim\sigma_{a}qE_{p}(q)$. Using $\omega\sim v_{0}q$ [see Eq. (10)] we get, $E_{u}(q)\sim(\sigma_{a}/v_{0})^{2}E_{p}(q)^{2}$. The plot in Fig. 7(a) shows good agreement between $E_{u}(q)$ obtained from our DNS and the conjecture above for small $q$. For large $q>2\pi/\ell_{\sigma}$, we expect viscous dissipation to be dominant and therefore, similar to the dissipation range in hydrodynamic turbulence, we expect an exponential decay in the energy spectrum $E_{u}(q)\sim\exp(-ak^{\delta})$ [96, 97]. From our numerical simulations, we find $\delta=1$.

It is worth noting that although turbulence in an apolar active suspension is controlled by half-integer defects in 2D [17, 18] and disclination loops or line defects in 3D [98, 99], the flow energy spectrum $E_{u}(q)\sim q^{-3.5}-q^{-5}$ reported in those works on active nematics does not differ drastically from our observation $E_{u}(q)\sim q^{-4.8}-q^{-6}$ for defect turbulence in the $R<R_{1}$ regime in our polar system.

III.4.1 Energy spectra for $R\ll R_{1}$

For $R\ll R_{1}$, the inter-defect separation $\xi$ is comparable to $\ell_{\sigma}$; however, for length scales much larger than $\xi$, the system should in effect be an unsteady Stokes fluid with fluctuating stresses with short-ranged spatial correlation, and with a correlation time $\tau$. In such a scenario it is straightforward to show [100, 101] that the equal-time velocity correlator has Ornstein-Zernike form, so that $E_{u}(q)\sim q^{D-1}/[(q\ell_{\tau})^{2}+1]$ for $q\xi\ll 1$, where $\ell_{\tau}=\sqrt{\mu\tau/\rho}$ is the distance vorticity diffuses in a time $\tau$. Note: (a) although the analysis of ref. [100] contains this result, they do not emphasize the distinct roles of $\xi$ and $\ell_{\tau}$. (b) The power-law correlations discussed in the inertia-less treatment of [102, 93] amount to the $q\ell_{\tau}\ll 1$ regime of the above.

To investigate this regime, we perform high resolution simulations in two dimensions with large system size and small $R=0.01$ (run SPP11) to ensure ${\mathcal{L}}\gg\xi$. The plot in Fig. 8(a) shows that the $E_{p}(q)\sim q$ for $q\xi\ll 1$ indicating that order parameter fluctuations are uncorrelated. We, therefore, expect the active stresses to be spatially uncorrelated. Consistent with the arguments above, the plot in Fig. 8(b) shows that the kinetic energy spectrum follows $E_{u}(q)\sim q/[(q\ell_{\tau})^{c}+1]$ with $c\approx 2.3$ (obtained from a least-squares fit) close to the theoretically predicted value $c=2$. Note that the length scale $\ell_{\tau}$ is larger than the interdefect separation $\ell_{\tau}/\xi\approx 3$. As $R\to 0$, we expect $\xi\to 0$ and $\ell_{\tau}\gg\xi$; therefore, the peak of spectra in Fig. 8(b) would shift to very small $q\sim 1/\ell_{\tau}$. Thus our analysis naturally recovers and explains the recently observed $E_{u}(q)\sim q^{-1}$ scaling of active nematic turbulence in the Stokesian regime [93].

III.4.2 Evidence for the broken-symmetry mode?

We close our discussion of energy spectra with a speculation backed by qualitative numerical measurements. In any noisy ordered state in which a continuous symmetry has been spontaneously broken, the spatial power spectral density should contain information about the broken-symmetry modes, i.e., the components of the order parameter field perpendicular to the mean ordering direction, whose variance should diverge at small wavenumber [103, 56, 104]. This variance should be seen in the energy spectrum in the ordered but noisy phase-turbulent state we observe for $R_{1}<R<R_{2}$. We offer preliminary evidence for such fluctuations. Insets b1 and d1 to Fig. 6 show that for $R>R_{1}$ the contributions to the energy spectrum from components of $\bm{p}$ in a representative ($\perp$) direction perpendicular to the mean direction of ordering far outweigh those from components in the ordering direction, especially at small $q$. This is consistent with the expectations of an enhanced variance mentioned above. We check for consistency that the spectra for the disordered phase for $R<R_{1}$ (insets a1 and c1 to Fig. 6) show no such anisotropy.

A quantitative study of the spectrum of fluctuations at small wavenumber, to test whether the regime $R_{1}<R<R_{2}$ has the classic features of a broken-symmetry phase, makes high computational demands. The wavelength at which the linear instability growth rate is maximum can be viewed as the scale of energy injection, and therefore as the small-scale cutoff for a long-wavelength study of the ordered phase. At the same time, scales substantially smaller than this cutoff must be resolved so that the instability and the nonlinear effects leading to phase turbulence can operate. If such a simulation is realized, the approach of choice would be to emulate [87] and construct an effective stochastically-forced theory for the small-$q$ modes via numerical coarse-graining.

IV Summary and prospect

We have shown that extensile active polar liquid crystals and swimmer suspensions can outswim their viscous instability. Their fate is governed by a control parameter $R$, the ratio of the inertia of self-advection to the scale of active stress. Our stability analysis and numerical studies find evidence for a continuous flocking transition with a growing correlation length as $R$ increases past a threshold of order unity, from hedgehog-defect turbulence to a noisy but ordered phase-turbulent state. A quiescent, linearly stable ordered state sets in at larger $R$. These dramatic advances in the theory of flocks in fluid, whose instability [10] can now be seen as simply the Stokesian limit of a rich phase diagram, should stimulate a new wave of experiments on swimmers at nonzero Reynolds number. Important directions for the near future are studies of finite-size scaling and long-wavelength order-parameter correlations for $R>R_{1}$ to establish the nature of the ordered phase; the construction of an effective stochastic theory for the long-wavelength modes, as carried out [87] for the Kuramoto-Sivashinsky [105, 106] equation; the inclusion of active-particle concentration; the contribution of other polar terms to the dynamics; and the effect of added random forcing on our phase diagram. Meanwhile, we look forward to tests of our theory in experiments on collections of swimmers at small but nonzero Reynolds number [45], as well as particle-based numerical simulations featuring, for example, collections of “spherobots” [107, 108].

Acknowledgements.

S.R. was supported by a J. C. Bose Fellowship of the SERB (India) and by the Tata Education and Development Trust. R.C. acknowledges the support of the SERB (India). R.C., N.R., and P.P. acknowledge support from the Department of Atomic Energy (DAE), Government of India, under Project Identification No. RTI 4007. The simulations were performed using resources provided by TIFR, Hyderabad. We thank H Chaté for valuable comments on an early draft. R.C., N.R., P.P., and S.R. acknowledge the contributions of Aditi Simha, now sadly deceased, to the development and completion of this work.

Appendix A Stability analysis with concentration

We now present the linear stability analysis in presence of the concentration field. Linearizing Eqs. (1), (2) and (3) about the base state $(\bm{u}=0,\bm{p}=\hat{\bm{x}},c=1)$ we get:

[TABLE]

The dispersion relation for Eqs. 21, 22 and 23 are obtained using the same procedure as highlighted in Section II.2. The growth rate of the twist-bend modes is same as Eq. (11) because the terms containing concentration fluctuations in Eqs. 21 and 22 point in the direction of ${\bm{q}}_{\perp}$. The growth rate of the splay-bend modes is identical to Eq. (10) for $v_{1}=0$ because concentration fluctuations decouple from the orientation and velocity distortions. For $v_{1}\neq 0$ the splay-bend modes couple to the concentration fluctuations and are obtained by taking in-plane divergence ($\nabla_{\perp}\cdot$) on Eqs. 21 and 22. We compare the most unstable growth rate in Fig. 9 and show that the $O(q)$ and the $O(q^{2})$ behavior at remains unaltered small-q.

Appendix B Stability analysis including the neglected polar terms

We now explain the role of the additional polar terms on the stability thresholds. We first recall the complete Eqs. (1) and (2) with the polar terms:

[TABLE]

where ${\bm{\Sigma}}^{a}\equiv-\sigma_{a}\bm{p}\bm{p}-\gamma_{a}(\nabla\bm{p}+\nabla\bm{p}^{T})$,  ${\bm{\Sigma}}^{r}\equiv[(1-\lambda)\bm{p}\bm{h}-(\lambda+1)\bm{h}\bm{p}]/2-\ell(\nabla\bm{h}+\nabla\bm{h}^{T})$, and $\bm{h}=-\delta F/\delta{\bm{p}}$. These are sufficient for our current purpose because we have already established in Section II.2.1 that for extensile systems, which is our main focus, the dominant failure modes are the twist bend modes which are always decoupled from concentration fluctuations. The dispersion relation for these modes with these terms is:

[TABLE]

where $R\equiv\rho v_{0}^{2}/2\sigma_{a}$, $\mu_{\pm}={\mu}(1\pm\beta)$, $\beta=\Gamma K\rho/\mu$, $B(\phi)=(v_{0}\mu_{-}/\sigma_{a})\cos\phi$, $G(\phi)=-(\mu_{-}^{2}/2\rho\sigma_{a})+(K/2\sigma_{a})(1+\lambda\cos 2\phi)^{2}$, and $\bar{G}(\phi)=-(\mu_{-}^{2}/2\rho\sigma_{a})+(K/2\sigma_{a})(1+\lambda)^{2}\cos^{2}\phi$. Comparing (B) to the dispersion relation for $\gamma_{a},\ell=0$ [Eq. (11) of the main text]

[TABLE]

and working to order $q^{2}$ we see that the transition between $O(q)$ and $O(q^{2})$ instabilities, which is determined by $A(0)$, remains unchanged and that the relative shifts of the coefficients of $q^{2}$ in the mode frequencies are of order

[TABLE]

which resemble inverse capillary numbers given that $\gamma_{a}$ and $\ell\sigma_{a}$ have units of surface tension. The length $\ell$ controls polar couplings of orientation to flow in the passive theory and as such should be related, in a microscopic theory, to a geometrical measure of polarity of the constituent particles, such as a fore-aft size difference. $\gamma_{a}$, with units of force per unit length, governs polar active contributions to the stress tensor. It is then reasonable to expect that $\gamma_{a}\sim\sigma_{a}\ell$, and thus that $D_{1}$ and $D_{2}$ are similar in magnitude. For particles of size $b$ the viscosity-based estimate $\sigma_{a}\sim\mu v_{0}/b$ implies $D_{1},D_{2}\sim\ell/b$, while an inertia-based estimate $\sigma_{a}\sim\rho v_{0}^{2}$ yields $D_{1},D_{2}\sim\ell Re/b$ where $Re=\rho v_{0}b/\mu$ is the Reynolds number at the particle scale. With either estimate, $D_{1}$ and $D_{2}$ should typically be small because $\ell/b$ is the ratio of a fore-aft size difference to an overall size. Fig. 10 shows the results of linear stability analysis for nonzero $D_{1}=D_{2}$. The $O(q)$-unstable, $O(q^{2})$-unstable and linearly stable regimes exist over the entire range; for $D_{1},D_{2}\ll 1$ the linear stability analysis differs negligibly from that with $\gamma_{a},\ell=0$; and the $O(q^{2})$-unstable regime shrinks in extent as $D_{1},D_{2}$ grow to values of order unity. We therefore expect the qualitative features of the dynamical phase diagram as seen in our numerical studies to persist for nonzero $\gamma_{a}$ and $\ell$, but we have not carried out the corresponding direct numerical solutions of the hydrodynamic equations.