Minimal continuum theories of structure formation in dense active fluids

Self-sustained dynamical phases of living matter can exhibit remarkable similarities over a wide range of scales, from mesoscopic vortex structures in microbial suspensions and motility assays of biopolymers to turbulent large-scale instabilities in flocks of birds or schools of fish. Here, we argue that, in many cases, the phenomenology of such active states can be efficiently described in terms of fourth- and higher-order partial differential equations. Structural transitions in these models can be interpreted as Landau-type kinematic transitions in Fourier (wavenumber) space, suggesting that microscopically different biological systems can share universal long-wavelength features. This general idea is illustrated through numerical simulations for two classes of continuum models for incompressible active fluids: a Swift-Hohenberg-type scalar field theory, and a minimal vector model that extends the classical Toner-Tu theory and appears to be a promising candidate for the quantitive description of dense bacterial suspensions. We also discuss briefly how microscopic symmetry-breaking mechanisms can enter macroscopic continuum descriptions of collective microbial motion near surfaces and conclude by outlining future applications.


Introduction
Simple and complex life forms can exhibit remarkably similar collective behaviours over a wide range of length and time scales [1][2][3]. Well-known examples are flocking phenomena in swarms of birds [4] and self-sustained turbulent phases in schools of fish [5] that share several qualitative features with the meso-scale dynamics in bacterial suspensions [6][7][8][9] and films [10][11][12]. When studying such processes from a physicist's perspective, a main challenge consists in identifying generic models that capture the most essential aspects of their dynamics. To this end, it is often useful to regard biological systems that comprise a large number of elementary self-propelled units, such as motor proteins or swimming cells in suspensions, as active 'fluids' [2,3,13,14]. Unlike conventional liquids, which typically require external energy injection (stirring, shearing, shaking, etc) for the formation of large-scale patterns, active fluids are driven internally as their microscopic constituents are capable of transforming chemical into kinetic energy. The interplay between this intrinsic pumping and nonlinear elastic stresses due to physical or biological interactions facilitates the emergence of complex dynamical structures [1][2][3], whose systematic classification poses a formidable theoretical task.
Over the past two decades, intense efforts have been made to understand the phenomenology of microbial and other active fluids, but in spite of substantial progress it is still not entirely clear which of their characteristics are universal or system-specific [2,15,16], and which classes of dynamical equations are capable of providing adequate minimal descriptions. In recent years, a considerable number of continuum models for active systems have been proposed [2,[17][18][19][20][21][22][23][24], but most of them have yet to be tested against experiments [14,25]. Many of those theories focus on the couplings between two or more order-parameters (concentration, solvent velocity, orientation fields, etc) and typically involve a large number of parameters, thus making comparison with experimental data very difficult. To exploit current and future progress in experimental imaging and tracking techniques [5,[26][27][28][29], and to understand better the general ordering principles that govern active matter [1,2], it will be necessary to identify tractable minimal models that not only capture the essential instability mechanisms but also allow for quantitative comparison with experiments.
In this paper, we will analyse two such minimal continuum theories for active suspensions by focusing on generic structural properties and stressing formal analogies with classical phase transitions. Our approach is based on the hypothesis that dynamical transitions in many internally or externally driven systems, such as microbial [6-8, 14, 28, 30] or vibrated colloidal suspensions [31,32], can be phenomenologically modelled as Landau-type transitions in Fourier (wavenumber) space, which suggests that minimal hydrodynamic descriptions of active matter can be obtained in terms of higher-than-second-order partial differential equations (PDEs). Higher-order PDEs have been previously derived and studied for a wide range of nonlinear structure formation phenomena [31,[33][34][35][36], including Rayleigh-Benard convection [37], polymer and vesicle dynamics [38], quasi-crystal formation [39] and theories of ionic liquids [40]. However, to our knowledge, models of this type have rarely been considered in the context of microbial suspensions [41]. Therefore, one of our main objectives here is to draw attention to the possibility that one can obtain useful, testable continuum theories of bacterial and other active fluids, by restricting field variables to a minimal set of experimentally accessible order-parameters but admitting fourth-and higher-order spatial derivatives. Such theories can then be used as phenomenological models to obtain predictions for the behaviour of active fluids under shear [42] or in different confining geometries [43], thereby providing a conceptual basis for the interpretation of rheological measurements [44][45][46] and the design of optimized microfluidic devices for the control of microbial flow [47,48].
The basic idea is readily summarized as follows. It is well-known that the incompressible Navier-Stokes equation is capable of describing the dissipative flow dynamics v(t, x) of a wide range of conventional 'passive' fluids, regardless of their exact microscopic composition. A main reason for this is that these systems behave similarly at long-wavelengths (smallwavenumbers), so that the leading order viscous dissipation can be described by a term 0 v in the field equations, with partial information about the microstructure being retained in the viscosity coefficient 0 (throughout, = ∇ 2 denotes the Laplacian). Transforming to Fourier-space, the dissipative term yields a simple quadratic 'dispersion' relation ∼ 0 |k| 2 , which represents the dominant contribution in a systematic small-wavenumber expansion and leads to damping in the absence of external stimuli. By contrast, in active fluids, viscous dissipation competes with internal or external energy input and, in principle, one cannot exclude that higher-order contributions of the form 0 |k| 2 + 2 |k| 4 + · · · become relevant as well. In fact, they will certainly be needed to ensure stability if, due to the complex interplay of nonlinear interactions and energy input, the coefficient 0 should change its sign. In this context, it should be noted that, in the case of active fluids, the coefficients n and other transport coefficients [49][50][51] will depend on both the physical interaction parameters and the motility parameters of the microscopic fluid constituents. Formally, the inclusion of higherorder terms in the Fourier-space expansions is analogous to the well-known Landau-expansion of order-parameter potentials and, accordingly, sign-changes in the coefficients n can give rise to Landau-type kinematic phase transitions. When going back to position space, terms |k| 4 , |k| 6 , . . . will transform into higher-order spatial derivatives 2 , 3 , . . . . The inclusion of such terms 4 makes the theory successively more non-local. The physical origins of effectively non-local interactions can be manifold [52][53][54], ranging from global packing constraints to hydrodynamic and chemical interactions in biological systems. In microbial suspensions, such non-localities may arise naturally from active stresses that are generated by the swimming strokes of the organisms and transported through the fluid. It seems plausible that non-local stress contributions are also present in passive fluids, even though they are not dynamically relevant in this case since friction is dominated by the Laplacian viscosity term. By contrast, for active systems, recent studies of bacterial suspensions [14,55] have shown that a model with a negative coefficient 0 captures experimental observation like energy spectra and correlation functions in a quantitative manner. In order to be mathematically consistent, such models need to include higher-order derivative terms that provide the damping necessary to avoid unphysical short-wavelength singularities.
From the preceding considerations, it seems plausible that a systematic characterization of active fluids in terms of their asymptotic small-wavenumber expansions can help to distinguish specific from universal properties, thereby providing a basis for more systematic classification schemes similar to those for thermodynamic equilibrium phases in classical fluids or spin systems. Moreover, this analogy-driven approach promises analytically tractable models of active suspensions that are considerably simpler than many of the currently studied (potentially more accurate) multi-component theories [19][20][21], and will hopefully enable quantitative comparisons with experiments in the near future. Hence, the main goal of the present paper is to provide a coherent conceptual framework for the application [14,55] of higher-order structure formation theories to active fluids. Specifically, we will focus on theoretical aspects of fourth-order continuum models, starting with the simplest case, which is given by a Swift-Hohenberg-type scalar or pseudo-scalar field theory [33,37]. This model is used to illustrate how microscopic symmetry-breaking mechanisms [56] can enter macroscopic continuum descriptions of microbial motion near surfaces. As an example, we demonstrate how these ideas can help to understand the hexagonal sperm-vortex patterns in the experiments of Riedel et al [30], which to our knowledge have not yet been explained in terms of a simple continuum theory. Subsequently, we will generalize to non-scalar order-parameters by considering a minimal vector theory for incompressible active suspensions. The resulting flow model extends the seminal Toner-Tu theory [15,18] and, as two very recent experimental studies [14,55] have shown, it is a promising candidate for the quantitative description of highly concentrated bacterial fluids. Here, we focus on investigating previously unexplored properties of such vector theories, aiming to identify universal characteristics of their flow spectra. Throughout, we use results from two-dimensional (2D) continuum simulations to illustrate selected dynamical properties of the different models in more detail.

(Pseudo) scalar order-parameter theory
The minimal model considered in this section belongs to the class of generalized Swift-Hohenberg theories [31,37]. Our motivation for prepending a brief discussion of this well-known model here is two-fold: it is helpful to recall some of its basic properties before considering the generalization to vectorial order-parameters. This model is also useful for illustrating how microscopic symmetry-breaking mechanisms [56] can be incorporated into macroscopic descriptions of experimentally relevant microbial systems [30], as discussed in section 2.4 below.

Model equations
We consider the simplest isotropic fourth-order model for a non-conserved scalar or pseudoscalar order-parameter ψ(t, x), given by where ∂ t = ∂/∂t denotes the time derivative, and = ∇ 2 is the d-dimensional Laplacian. The force F is derived from a Landau-potential U (ψ) where c > 0 to ensure stability. We will assume throughout that the system is confined to a finite spatial domain ⊂ R d of volume adopting periodic boundary conditions in simulations. The derivative terms on the rhs of equation (1) can also be obtained by variational methods from a suitably defined energy functional (see appendix (A.1)). In the context of active suspensions, ψ could, for example, quantify local energy fluctuations, local alignment, phase differences, or vorticity. In this case, the transport coefficients (a, b, c, γ 1 , γ 2 ) in equations (1) and (2) will contain passive contributions due to steric or other physical interactions as well as active motilityrelated contributions. In general, it is very challenging to derive the exact functional dependence between macroscopic transport coefficients and microscopic interaction and motility parameters for active non-equilibrium systems, although considerable progress has been made in recent years [49][50][51]57]. With regard to practical applications, however, it is often sufficient to view transport coefficients as purely phenomenological parameters that can be determined by matching the solutions of continuum models, such as the one defined by equations (1) and (2), to experimental data [14,55]. This is analogous to treating the viscosity in the classical Navier-Stokes equations as a phenomenological fit parameter. The actual predictive strength of a continuum model lies in the fact that, once the parameter values have been determined for given a set-up, the theory can be used to obtain predictions for how the system should behave in different geometries or under changes of the boundary conditions (externally imposed shear, etc). In some cases, it may also be possible to deduce qualitative parameter dependences from physical or biological considerations. For instance, if ψ describes vorticity or local angular momentum in an isolated active fluid, say a bacterial suspension, then transitions from a > 0 to a < 0 or γ 0 > 0 to γ 0 < 0, which both lead to non-zero flow patterns, must be connected to the microscopic self-swimming speed v 0 of the bacteria. Assuming a linear relation, this suggests that, to leading order, a 0 = δ − αv 0 where δ > 0 is a passive damping contribution and αv 0 > 0 the active part, and similarly for γ 0 . It may be worthwhile to stress at this point that higherthan-second-order spatial derivatives can also be present in passive systems, but their effects on the dynamics will usually be small as long as γ 0 > 0. If, however, physical or biological mechanisms can cause γ 0 to become negative, then higher-order damping terms, such as the γ 2 -term in equation (1), cannot be neglected any longer as they are essential for ensuring stability at large wavenumbers. For completeness, one should also note that in the case of a conserved order-parameter field the field equations would either have to take the current-form ∂ t = −∇ · J( ) or, alternatively, one could implement conservation laws globally by means of Lagrange multipliers [38]. For example, for a dynamics similar to that of equation (1) and a simple global 'mass' constraint the Lagrange-multiplier approach yields the non-local equations of motions In the remainder of this section, however, we shall focus on the local dynamics defined by equations (1) and (2), since this well-known example will be a useful reference point for the discussion of the vector model in section 3.

Linear stability
The fixed points of equation (1) are determined by the zeros of the force F(ψ), corresponding to the minima of the potential U , yielding ψ 0 = 0 and Linearization of equation (1) Similarly, one finds for The unusual sign-convention in the exponential of the perturbation ansatz was so chosen as to emphasize the formal similarity of equations (5) and (6) with the quartic Landau potential (2), i.e. modes with σ < 0 are unstable. From equations (5) and (6), we see immediately that γ 2 > 0 is required to ensure smallwavelength stability of the theory and, furthermore, that non-trivial dynamics can be expected if a and/or γ 0 take negative values. In particular, all three fixed points can become simultaneously unstable if γ 0 < 0. The analogy with classical Landau-transitions is evident if we compare (5) and (6) with the order-parameter potential U in equation (2) for the symmetric case b = 0: changing the sign of γ 0 induces a dynamical transition (in Fourier space), which is formally similar to the standard 'configurational' second-order transition [2] in the vicinity of a = 0.

Numerical results in two dimension (2D)
We briefly illustrate the γ 0 -induced changes in the dynamics of the (pseudo-)scalar field ψ(t, x) through 2D numerical results. The discussion in this part merely serves as a reminder before considering symmetry-breaking in section 2.4.  (1) and (2). The unit time is defined by the damping time-scale t u = L 4 /γ 2 , where L is the length of the 2D simulation box, and the order-parameter is measured in units (1) and (2) in two space dimensions, we implemented a pseudospectral algorithm with periodic boundary conditions as commonly used in computational fluid dynamics [58]. The model equations are projected onto a Fourier space basis, and the remaining ordinary differential equations are solved numerically by an operator splitting method that computes the linear operator exactly [59]. The nonlinear terms were evaluated by applying the '2/3-rule' to suppress aliasing errors [60]. We simulated the model dynamics on 2D cubic grids with sizes ranging from 64 × 64 to 256 × 256 lattice points. The solver was written in Matlab, and its numerical stability was verified for a wide range of parameters and space-time discretizations. Rescaled dimensionless variables and parameters as adopted in the simulations are summarized in table 1. The rescaled time steps were typically of the order of t = 10 −1 . All simulations were initiated with isotropic, randomly chosen orderparameter values.

Structural transitions.
Results from the numerical simulations for the order-parameter field ψ(t, x) and two qualitatively different potentials U (ψ) are summarized in figure 1. In these simulation, the parameter γ 0 was varied between successive runs while keeping all other parameters fixed. To quantify changes in the quasi-stationary dynamics of ψ as a function of γ 0 , we measured the space-time averaged standard deviation σ 2 ψ = ψ 2 − ψ 2 (figures 1(a) and (d)). Regions with σ 2 ψ = 0 correspond to disordered structureless stationary states, whereas σ 2 ψ > 0 indicates the emergence of stationary or quasi-stationary dynamical structures. Singular points in the curve σ 2 (γ 0 ) signal qualitative changes in the order-parameter dynamics.
In the case of a mono-stable potential (a > 0), the quantifier σ 2 ψ undergoes a series of continuous transitions as γ 0 is lowered to negative values, see figure 1(a). Each of those transitions corresponds to an increase in the number of 'stripes' that are found to persist for long periods of time in the simulations (figures 1(b) and (c)). By contrast, in the case of bi-stable potentials (a < 0), the onset of pattern formation carries the signature of a first-order transition reflected by a sudden jump in σ 2 ψ ( figure 1(d)). However, while such singularities in σ 2 ψ share some formal similarities with macroscopic phase transitions, one could also argue that they merely signal a change in the typical number of excitable modes in the system. In fact, by , (c) Snapshots of the orderparameter field ψ at t = 500, scaled by the maximum value ψ m , for a monostable potential U (ψ) and homogeneous random initial conditions. After the first transition two stripes appear, and the number of stripes increases with the number of transitions. (e), (f) Snapshots of the order-parameter at t = 500 for a bi-stable potential. For γ 0 −(2π ) 2 γ 2 /L 2 , increasingly more complex quasi-stationary structures arise; see [31,61] for similar patterns in excited granular media and chemical reaction systems. viewing such elementary excitations as 'quasi-particles', the structural transitions in figures 1(a) and (d) appear to be more closely related to finite-systems singular points [62][63][64].
An estimate of the critical absolute value γ 0 for the first 'disorder-structure' transition can be obtained by dimensional analysis, or by equating the last two terms in equations (5) or (6), yielding γ c 0 ≈ −(2π ) 2 γ 2 /L 2 . For γ 0 γ c 0 , increasingly more complex quasi-stationary patterns may arise (figures 1(c) and (d)). Structures similar to those in figure 1 have been observed and widely studied [35] in granular media [31] and chemical systems [61]. In the next section, we shall demonstrate that certain aspects of collective microbial motion can be described within the same class of fourth-order PDEs.

Symmetry breaking
With regard to microbial suspensions, the minimal model (1) is useful for illustrating how microscopic symmetry-breaking mechanisms that affect the motion of individual organisms or cells [56,[65][66][67] can be implemented into macroscopic field equations. To demonstrate this, we interpret ψ as a vorticity-like 2D pseudo-scalar field that quantifies local angular momentum in a dense microbial suspension, assumed to be confined to a thin quasi-2D layer of fluid. If the confinement mechanism is top-bottom symmetric, as for example in a thin free-standing bacterial film [10], then one would expect that vortices of either handedness are equally likely. In this case, equation (1) must be invariant under ψ → −ψ, implying that U (ψ) = U (−ψ) and, therefore, b = 0 in equation (2). Intuitively, the transformation ψ → −ψ corresponds to a reflection of the observer position at the midplane of the film (watching the 2D layer from above versus watching it from below).
The situation can be rather different, however, if we consider the dynamics of microorganisms close to a liquid-solid interface, such as the motion of bacteria or sperms cells in the vicinity of a glass slide (figure 2). In this case, it is known that the trajectory of a swimming cell can exhibit a preferred handedness [56,[65][66][67]. For example, the bacteria Escherichia coli [56] and Caulobacter [65] have been observed to swim in circles when confined near to a solid surface. More precisely, due to an intrinsic chirality in their swimming apparatus, these organisms move on circular orbits in clockwise (anticlockwise) direction when viewed from inside the bulk fluid (glass surface). Qualitatively similar behaviour has also been reported for sea urchin sperm swimming close to solid surfaces [68].
Hence, for various types of swimming microorganisms, the presence of the near-by no-slip boundary breaks the reflection symmetry, ψ → −ψ. The simplest way of accounting for this in a macroscopic continuum model is to adapt the potential U (ψ) by permitting values b = 0 in equation (2). The result of a simulation with b > 0 is shown in figure 2(a). In contrast to the symmetric case b = 0 (compare figure 1(c)), an asymmetric potential favours the formation of stable hexagonal patterns ( figure 2(a))-such self-assembled hexagonal vortex lattices have indeed been observed experimentally by Riedel et al [30] for highly concentrated spermatozoa of sea urchins (Strongylocentrotus droebachiensis) near a glass surface ( figure 2(b)). Finally, to illustrate how the mean local angular momentum changes with asymmetry, we also show in figure 2(c) the average value ψ as function of parameter b. 5  (1) and (2) with b > 0, corresponding to a broken reflection symmetry ψ → −ψ. Blue regions correspond to clockwise motions. (b) Hexagonal vortex lattice formed spermatozoa of sea urchins (S. droebachiensis) near a glass surface; from [30] adapted and reprinted with permission from AAAS. At high densities, the spermatozoa assemble into vortices that rotate in clockwise direction (inset) when viewed from the bulk fluid. (c) Mean local angular momentum ψ as a function of the asymmetry parameter b. Arrows indicate the transitions between striped and hexagonal states. Red filled circle corresponds to the hexagonal state shown in (a).

Vector model for an incompressible active fluid
We now generalize the preceding considerations to identify a minimal vector-field model for dense microbial suspensions. Previously developed continuum theories [2,[17][18][19][20][22][23][24] of microbial fluids typically distinguish solvent concentration, bacterial density, solvent velocity, bacterial velocity, and various orientational order-parameter fields (polarization, Q-tensors, etc). Aiming to identify a minimal hydrodynamic model, we construct a simplified higherorder theory by focusing exclusively on the dynamics of the mean bacterial 6 velocity field v(t, x) and restricting ourselves to the incompressible limit. By construction, the resulting v-only theory, which is essentially a minimal Swift-Hohenberg-type [37] extension of the Toner-Tu model [17,18], may not be applicable to swarming or flocking regimes, where density fluctuations are dominant, but it can provide a useful basis for quantitative comparisons with experiments and simulations on highly concentrated active suspensions [14,55]. Another assumption implicit to the vector model below is that the energy input, required to maintain non-zero velocity patterns, is quasi-stationary. Relaxation of this assumption would imply the need for additional energy balance equations that account for spatial and temporal variations in the conversion of chemical into kinetic energy of motion. In other words, the v-only theory formulated below only applies to situations where concentrations of nutrients, oxygen, etc in a microbial suspension are approximately constant during the observation period.
In practice, v can be determined applying suitable coarse-graining procedures (PIV algorithms, local averaging, etc) to discrete experimental or numerical velocity data [14,69].

Model equations
Postulating incompressibility, which is a good approximation for very dense suspensions [14], 7 we assume that the dynamics of v is governed by the generalized Navier-Stokes equation The pressure p(t, x) is the Lagrange multiplier for the incompressibility constraint. Similar to the scalar case, equation (2) above, the (A, C)-terms in equation (8) Physically, the inclusion of a polar ordering potential accounts for the fact that microorganisms typically exhibit head-tail asymmetries that may favour polar alignment, as manifested in the 'bionematic' jets that form in bacterial suspensions [8,70]. For A > 0 and C > 0, the potential is mono-stable and the fluid is damped towards a disordered state with v = 0. By contrast, for A < 0, equation (9) describes a d-dimensional mexican-hat (sombrero) potential with fixedpoints |v| = √ −A/C corresponding to global polar order. However, the fact that polar ordering appears only locally but not globally in suspensions of swimming bacteria [7,8,70] suggests that other instability mechanisms must be at work [23]. To capture this mathematically, one must either introduce additional order parameters [2,17,18] or destabilize the theory by identifying a suitable phenomenological ansatz for the effective stresses [37].
Adopting the latter approach, we postulate that the components of the symmetric and traceless rate-of-strain E tensor are given by where is a d × d-dimensional mean-field approximation to the Q-tensor, representing active nematic stresses [23,49] due to swimming (δ i j is the Kronecker tensor). Although the S-term does not affect the linear stability of the model, general hydrodynamic arguments [22] imply that S < 0 for pusher-swimmers like E. coli [71] or Bacillus subtilis, whereas S > 0 for pullertype microswimmers such as Chlamydomonas algae [72]. The 0 -term in (10) is dictated by the requirement that the model contains the Navier-Stokes equations as a limit case, and the 2 -damping term is motivated by generic stability considerations, as recent experiments [14,55] suggest that 0 can become negative in dense bacterial suspensions. Inserting equations (10) and (11) into equation (8), and defining we obtain The standard Navier-Stokes equations for a passive fluid are recovered for S = A = C = 2 = 0 and 0 > 0. A variational formulation of the combined field equations (7) and (13) is given in appendix A.2. For 0 > 0 and 2 = 0, equation (13) reduces to an incompressible version of the classical Toner-Tu model [2,17,18]. It is, however, the combination of the two -terms with the non-variational convective derivative that turns out to be crucial for the formation of self-sustained quasi-chaotic flow patterns. The linear -terms are reminiscent of the higherorder spatial derivatives in the classical Swift-Hohenberg theory [37], see equations (1) and (13) with 0 < 0 and 2 > 0 yields a simple-if not the simplest-generic continuum description of turbulent meso-scale instabilities observed in dense bacterial suspensions [14]. More generally, equation (13) can provide a satisfactory phenomenological model whenever interaction terms in more complex field theories, that lead to instabilities in the v-field, can be effectively approximated by a fourth-order Taylor expansion in Fourier space. This is likely to be the case for a wide range of active systems. Phrased differently, the last two terms in equation (13) may be regarded as the Fourier-space analogue of the Toner-Tu driving terms, which correspond to a series expansion in terms of the order-parameter. Hence, similar to the higher-order gradient terms in the scalar theory from equation (1), the ( 0 , 2 )-terms in equation (13) describe intermediate-range interactions, and their role in Fourier-space is similar to that of the Landau potential in velocity space.

Linear stability analysis
To support the qualitative statements in the preceding paragraph, we now perform a stability analysis for the 2D case relevant to the simulations discussed below, assuming 0 < 0 and C > 0, 2 > 0.
The fixed points of equations (7) and (13) are given by the extrema of the quartic velocity potential U (v). For arbitrary values of A, equations (7) and (13) have a fixed point that corresponds to a disordered isotropic state (v, p) = (0, p 0 ) where p 0 is a constant pressure. For A < 0, an additional class of fixed points arises, corresponding to a manifold of globally ordered polar states (v, p) = (v 0 , p 0 ), where v 0 is constant vector with arbitrary orientation and fixed swimming speed |v 0 | = √ −A/C =: v 0 . Linearizing equations (7) and (13) for small velocity and pressure perturbations around the isotropic state, v = and p = p 0 + η with |η| | p 0 |, and considering perturbations of the form we find Multiplying the second equation by k and using the incompressibility condition implies that η = 0 and, therefore, Assuming 0 < 0 and 2 > 0, and provided that 4A < | 0 | 2 / 2 , we find an unstable band of modes with σ 0 (k) < 0 for k 2 − < |k| 2 < k 2 + , where For A < 0 the isotropic state is generally unstable with respect to long-wavelength (i.e. small-|k|) perturbations. We next perform a similar analysis for the polar state (v 0 , p 0 ), which is energetically preferred for A < 0 and corresponds to all active particles swimming in the same direction ('global order'). In this case, when considering small deviations it is useful to distinguish perturbations perpendicular and parallel to v 0 , by writing = || + ⊥ where v 0 · ⊥ = 0 and v 0 · || = v 0 || . Without loss of generality, we may choose v 0 to point along the x-axis, v 0 = v 0 e x . Adopting this convention, we have || = ( || , 0) and ⊥ = (0, ⊥ ), and to leading order Linearization for exponential perturbations of the form where with I = (δ i j ) denoting the identity matrix. Multiplying equation (23) with ik, and using the incompressibility condition (22), giveŝ Inserting this into equation (23) and defining is the orthogonal projector of k, we obtain σˆ = −M ⊥ˆ .
The eigenvalue spectrum of the matrix M ⊥ is given by The zero eigenvalues correspond to the Goldstone modes. The non-zero eigenvalues have eigenvectors (−k y , k x ), implying that, for 0 < 0, there will be a range of exponentially growing modes in the direction perpendicular to k. The right column shows the rescaled dimensionless parameters used in the simulations of the vector model from equation (13). The unit time is defined by the damping time-scale t u = L 4 / 2 , where L is the length of the 2D simulation box, and the order-parameter is measured in units Model parameter Rescaled dimensionless parameter Equations (17) and (28) predict that, when A < 0 and 0 < 0, isotropic and polar fixed points become simultaneously unstable, thereby signalling the existence of spatially inhomogeneous dynamic attractors. More generally, within the class of standard PDEs, the two -terms in equation (13) appear to provide the simplest 'linear way' of obtaining a v-only theory that exhibits non-trivial stationary dynamics. In principle, one could also try to model instabilities by combining odd or fractional powers of |k| in equations (17) and (28); this would be analogous to replacing the quartic Landau potential by a more general function of |v|. However, when considering eigenvalue spectra based on odd or non-integer powers of |k|, the underlying dynamical equations in position space would become fractional PDEs. Such fractional models could potentially be useful for describing active suspensions with long-range or other types of more complex interactions, but their analysis goes far beyond the scope of this paper.

Numerical results in 2D
We simulated the vector model, defined by equations (7) and (13), in two space-dimensions using an algorithm similar to that described in section 2.3. The primary difference compared with the simulations for the scalar model is an additional pressure correction subroutine that ensures the incompressibility of the flow (see [14] for details). Table 2 summarizes characteristic units and rescaled parameters as adopted in the computations. All simulations were initiated with random initial conditions, and the typical time discretization was t = 0.1 (in characteristic time units t u ).

Kinematic transitions.
We first study how a decrease of the 'viscosity' parameter 0 affects the stationary dynamics for mono-stable and mexican-hat (polar-ordering) potentials. To this end, simulations were performed with fixed potential functions at three different values of the pusher/puller parameter S, while varying 0 between successive runs. Changes in the quasi-stationary dynamics are quantified by measuring the space-time averaged variance figure 3. Similar to the scalar model, the first transition from an isotropic state with v ≡ 0 to a non-trivial stationary dynamics with σ 2 v > 0 is found to occur at 0 ≈ −(2π ) 2 2 /L 2 . To illustrate kinematic changes in the flow dynamics in more detail, we show quasistationary snapshots from simulations with a mexican-hat potential for different values ( 0 , S) in figure 4. In the special case S = 1, corresponding to a vanishing convective derivative in equation (13), stationary cubic vortex lattices form, with an increasing number of vortices as 0 is decreased (figures 4(a) and (b)). By contrast, for S = 1, nonlinear convective effects cause distortions in the vortex lattices. As a consequence, the dynamical system no longer approaches a time-independent stationary state but instead exhibits a complex non-equilibrium dynamics (figures 4(c) and (d)). For pushers (S < 0), the resulting turbulent flow patterns look remarkably similar to those observed in dense suspensions of B. subtilis [7,8,14,55,70]. 8

Flow spectra.
To resolve the structure of the flow fields obtained from equations (7) and (13), we calculated the energy spectrum E(k), formally defined by where k = |k|. By virtue of the Wiener-Khinchine theorem [73], E(k) can be estimated by Fourier-transformation of the equal-time two-point velocity correlation function, yielding in d dimensions Traditionally, spectral flow analysis has been an important tool in the investigation of classical turbulence phenomena [73]. Flow spectra for our numerical data are summarized in figure 5. The critical case S = 1 provides a useful reference point, since in this case the stationary dynamics becomes static, as illustrated in figures 4(a) and (b) for the sombrero potential. Accordingly, the spectra for S = 1 exhibit a sharp peak that reflects the typical vortex size in the stationary state (see green curves in figures 5(a) and (b)). In the case of the mono-stable potential, changing the pusher/puller parameter to values S = 1 affects the asymptotic slopes of the spectrum but leaves the position of the maximum practically constant (figure 5(a)). By contrast, for the polar-ordering potential, both slope and position of the maximum change as S is decreased or increased from unity ( figure 5(b)). The large-|S| pusher-spectra for the mexican-hat potential agree well with those measured in dense quasi-2D B. subtilis [14] suspensions.

Bionematic jets.
To illustrate qualitatively how the velocity potentials affect the stationary dynamics-and, hence, the spectral flow properties-we present in figure 6 snapshots from computations with an intermediate-size simulation volume and 0 −(2π) 2 2 /L 2 . For the mono-stable potential, we observe fairly homogeneous vortex structures (figures 6(a) and (b)) in agreement with the relatively sharp spectral peaks in figure 5(a). By contrast, in the case of the mexican-hat potential (figures 6(c) and (d)), extended jet-like regions form that look very similar to the 'zooming bionematic' phases in dense B. subtilis suspensions [7,8,70].

Nonlinear limit and 'universality'.
Interestingly, our numerical results suggest that the spectra approach universal functional forms in the nonlinear limit |λ 0 | = |S − 1| 1, even though the exact asymptotic behaviour depends on the type of the potential. Moreover, for pushers with λ 0 1, the spectra obtained from equation (13) agree well with those measured for dense quasi-2D B. subtilis suspensions [14], and a very similar small-k scaling was also In the strongly nonlinear regime, corresponding to |S| 1, the spectra of both pullers (S > 0) and pushers (S < 0) seem to approach 'universal' limit functions, the exact shape of which depends on the type of the potential. Simulation parameters: grid size 256 × 256, 0 = −1 in all simulations.
observed in 2D particle simulations of self-propelled rods [14]. These observations hint at some 'universality' in the spectral properties of dense active particle systems. Future investigations of this question will require larger simulations as well as systematic experimental investigations of larger systems, which allow to extract asymptotic scaling laws and other spectral details with higher accuracy. Generally, in the strongly nonlinear regime |λ 0 | = |1 − S| 1, the vector model predicts phenomenologically similar behaviour for both pullers (S > 0) and pushers (S < 0), see figure 6. However, while values of λ 0 1 appear to have been realized in dense suspensions of pusher-type B. subtilis bacteria [8,14,70], it seems unclear at present whether puller (e.g. algal) suspensions with S 1 can be achieved in experiments.

Conclusions and future challenges
Identifying 'universal' features of active fluids is one of the key challenges en route to a more systematic physical classification of biological matter. Here, we have argued that, in many cases, phenomenological aspects of dynamical non-equilibrium phases can be naturally described by minimal models that focus on a limited number of experimentally accessible order-parameters and are based on higher-than-second-order spatial derivatives. If true, this would imply that a variety of active fluids with different microscopic constituents can exhibit very similar longwavelength behaviour. More generally, adopting the view that kinematic transitions in living fluids reflect qualitative changes in small-wavenumber expansions, and thus may be interpreted as Landau-type transitions in Fourier space, can help to catalogue non-equilibrium systems according to their asymptotic behaviour at long wavelengths, similar to the classification of phase transitions in equilibrium thermodynamics.
Simplified continuum descriptions of active fluids remain practically relevant because many of the recently proposed multi-order-parameter theories feature a large number of unknown transport coefficients and, therefore, cannot be tested in detail with present data. The pusher flow field in (c) agrees qualitatively with experimentally observed flow fields for dense B. subtilis [7,8,14,70] suspensions.
While our above analysis has focused on the most basic scalar and vector models, the approach can be easily extended to higher tensorial order-parameter fields (e.g., Q-tensor descriptions [74] of active nematics). To test whether microscopically different active fluids do indeed share universal hydrodynamic long-wavelength characteristics will require combined analytical, numerical and experimental efforts. First attempts to apply the above ideas to very dense bacterial suspensions give promising results [14,55] but further quantitative studies will be needed to decide whether fourth and higher-order PDEs are capable of providing a sufficiently accurate phenomenological description of experimentally observed active phases. To expedite quantitative comparisons with experiments, it would be desirable to develop alternative computation schemes that will allow for faster simulations of higher-order PDE models. Promising candidates could be suitably adapted lattice-Boltzmann algorithms [21] that implement negative 'viscosities' [75].
Finally, the theoretical analysis of higher-order PDEs poses a number of future mathematical challenges, one of which being their derivation from underlying microscopic, multi-component or kinetic models through systematic projection methods [76]. The good agreement of the vector model with experimentally measured flow structures in quasi-2D B. subtilis suspensions [14,55] suggests that derivations of active continuum theories from microscopic models, which typically involve gradient expansions, should go beyond the most frequently considered second-order approximations. Another interesting question relates to the formulation of consistent boundary conditions at interfaces. This problem, which was circumvented in our simulations by considering toroidal domains, is also encountered in other higher-order structure formation and transport theories [31,35,40]. To identify physically reasonable boundary conditions for active fluids at solid interfaces is not only relevant from a mathematical perspective but also from the experimental point of view, as it will directly affect predictions for effective shear viscosities [77] and other measurable quantities. In particular, future simulations of the vector model in confined geometries can help to interpret rheological measurements in microbial suspensions [42,45,46], and they also promise insights into the effects [43] of surface structures and microfluidic channel design on active fluids.

A.2. Vector model
In component form, the vector model dynamics defined by equations (7) and (13) reads Using equation (A.2) with φ = ( p, v), these field equations can be written as (A.9) As evident from the representation (A.8), apart from the convective derivative, which can be rescaled by active hydrodynamic stresses via λ 0 , see equation (10), the vector model dynamics can be understood in terms of an optimization of the effective 'free-energy' functional (A.9).