Polar patterns in active fluids

Abstract

We study the spatio-temporal dynamics of a model of polar active fluid in two dimensions. The system exhibits a transition from an isotropic to a polarized state as a function of density. The uniform polarized state is, however, unstable above a critical value of activity. Upon increasing activity, the active fluids displays increasingly complex patterns, including traveling bands, traveling vortices and chaotic behavior. The advection arising from the particles self-propulsion and unique to polar fluids yields qualitatively new behavior as compared to that obtained in active nematic, with traveling-wave structures. We show that the nonlinear hydrodynamic equations can be mapped onto a simplified diffusion-reaction-convection model, highlighting the connection between the complex dynamics of active system and that of excitable systems.

---

I. Introduction

Bacterial suspensions, in vitro mixtures of cytoskeletal filaments and motor proteins, and migrating epithelial cell layers are examples of active fluids composed of interacting units that consume energy and collectively generate motion and mechanical stresses. Due to their elongated shape, active particles can exhibit orientational order at high concentration and have been likened to “living liquid crystals”.¹ They have been modeled by continuum equations built by modifying the hydrodynamics of liquid crystals to include nonequilibrium terms that account for the activity of the system,^2–9 or derived from specific microscopic models.^10,11

The theoretical and experimental study of active fluids has revealed a wealth of emergent phenomena, including spontaneously flowing states,^7,12–14 unconventional rheological properties,^15–21 and new spatiotemporal patterns not seen in passive complex fluids.^22,23 Large-scale swirling patterns and “flocking” have recently been observed in motility assays consisting of a highly concentrated suspension of actin filaments propelled by myosin molecular motor proteins tethered to a plane.^24,25 Experiments in bacterial suspensions^26,27 and simulations of self-propelled hard rods^28–31 have also shown the formation of polar high-density clusters traveling trough a low density background. In a recent paper one of us and collaborators examined the spatio-temporal dynamics of active fluids with nematic symmetry in two dimensions, demonstrating the onset of spatially modulated relaxation oscillations and the close analogy between the dynamics of active fluids and that of excitable media.²²

In this paper we consider the spatio-temporal dynamics of an active suspension with polar symmetry. The model is appropriate to describe systems such as bacterial suspensions, where the particles are swimmers, hence head–tail asymmetric or polar. As a result, the active fluid can order in a macroscopic polar state, characterized by a nonzero value of a vector order parameter P, describing the mean polarization of the system. The order parameter can be written in terms of a magnitude P and a polar director p (with |p| = 1), denoting the direction of spontaneously broken symmetry in the ordered state. While in the nematic state the director n is a headless unit vector because the ordered state is invariant for n → —n, the polar state does not possess this symmetry and the polar director p is a true unit vector. As a result, the continuum equation for a polar active fluids contain terms that are forbidden by symmetry in the case of an active nematic. These terms describe swimming of the active particles relative to the bulk fluid and yield new self-advection contributions to the continuum equations.

Continuum descriptions of the collective dynamics of polar active particles have been discussed before in the literature in two cases: for particles moving on a frictional substrate (“dry” system) and for particles swimming in a fluid (active suspension). The hydrodynamics of the dry system is a continuum theory of the Vicsek model of flocking, consisting of a collection of self-propelled particles moving on a frictional substrate with noisy aligning rules. It was first introduced phenomenologically by Toner and Tu^2,3 and recently derived by explicit coarse graining of the microscopic dynamics.^32,33 In this case the only conserved variable is the density of active particles as momentum is not conserved. The theory is then formulated entirely in terms of density and polarization fields and the equations lack Galileian invariance due to the presence of the substrate. The hydrodynamic equations of active suspensions have been written down phenomenologically on the basis of symmetry considerations and also derived from specific microscopic models of a collection of interacting active particles. The equations used here were first obtained by Liverpool and Marchetti for a suspensions of cytoskeletal filaments with crosslinking motor proteins.¹¹ While in this case the filaments may be better described as a “mutually propelled”, rather than “self-propelled”, as the activity resulting in propulsion arises from the forces exchanged among filaments by the active cross-linkers, the form of the hydrodynamic equations does not depend on such microscopic details.

This work focuses on pattern formation in polar active suspensions. While the behavior of polar active films in a quasi-1D geometry was discussed before in ref. 14, here we present a systematic study of a two-dimensional system with periodic boundary conditions. New results include: the derivation of a “phase diagram” for the system that unifies the analysis of instabilities discussed previously in the literature in various contexts (Fig. 2), the prediction of the traveling vortices regime shown in Fig. 8, the derivation of a low dimensional model that captures the onset of the traveling waves regime, and the prediction that the transitions between the various regimes are discontinuous, allowing for the possibility of coexistence and hysteresis.

We begin by identifying the linear instabilities of the ordered state of polar active fluids, summarized in Fig. 2. These instabilities have been discussed before in the literature, either in the context of dry systems or of active suspensions. Here we highlight the connection between the behavior of the suspension and the dry limit. To go beyond the linear case, we then integrate numerically the nonlinear hydrodynamic equations for increasing values of activity. In a polar system there exist two classes of active terms. The first class, with strength controlled by a parameter we call α, is present in both systems of nematic and polar symmetry and describes active stresses induced by coupling of orientation and flow. The sign of α identifies active systems that generate contractile stresses (α > 0) versus those that generate tensile stresses (α < 0). The second class, with strength controlled by a parameter we call w, is unique to polar fluids and describes a variety of nonequilibrium advective contributions to the dynamics arising from the particles' self propulsion. Both parameters are ultimately controlled by the activity of the system, as embodied for instance by the rate of adenosine-triphosphate (ATP) consumption in motor-filament suspensions or by the forces exerted by swimmers on the surrounding fluid in a bacterial suspension. For this reason, and to reduce the number of independent parameters, in most of the following we take w = fα, with f a numerical factor of order one, and discuss the behavior of the system for increasing values of the activity α. Upon increasing α we observe a variety of increasingly complex patterns. A regime that is unique to polar systems consists of bands extending in the direction orthogonal to that of mean order (x direction) and traveling along x. Although traveling bands have been seen in dry systems, either via direct simulations of Vicsek models²⁹ or by numerical solution of the continuum theory,^32,34 there are important qualitative differences between the traveling bands in dry systems and in suspensions. In dry systems the bands appear to be solitary waves of ordered regions of high density and polarization traveling in a disordered background. In suspensions, in contrast, the traveling waves are mainly controlled by coupled oscillations in the flow and the direction of the polarization. They can be understood in terms of a simplified reaction-diffusion-advection model that highlights the crucial role played by active stresses in providing the energy input needed for setting up the pattern formation and by advective terms unique to polar systems in controlling the traveling nature of the pattern. As the activity is further increased, the traveling wave pattern is replaced by more complex spatial structures, including oscillating flows with pairs of traveling vortices, and eventually chaotic behavior.

Finally, spatial patterns in active polar fluids in two dimensions have also been discussed by Voituriez et al.³⁵ These authors include a coupling between splay and density fluctuations that is present in both equilibrium and active polar fluids, but neglect the convective couplings unique to active polar systems that are responsible for the oscillatory or traveling nature of some of the spatial patterns found in the present work. The modulated flowing phases identified in ref. 35 are characterized by a finite steady flow induced by activity, but do not represent traveling patterns. They are “equilibrium-like” in the sense that they have a corresponding analogue in equilibrium polar systems.³⁶ When polar convective terms are included the modulated patterns are replaced by traveling vortices.

In Sec. II we summarize the hydrodynamic equations of a polar active suspension. The homogeneous states of these equations and their linear stability are discussed in Sec. III, where we also show the connection between the instabilities of the suspension and those of a dry system. The results of the linear stability analysis are summarized in the phase diagram shown in Fig. 2. In Sec. IV we describe the spatio-temporal patterns obtained by numerical solution of the nonlinear hydrodynamic equations in a planar sample in two dimensions. These include traveling bands, traveling vortices and chaotic flow. We also show that the nonlinear equations can be mapped onto a simplified diffusion-reaction-convection model, highlighting the connection between the complex dynamics of active system and that of excitable systems. We conclude with a brief discussion.

II. Hydrodynamics of polar active suspensions

We consider a suspension of rod-like active particles of length [small script l] and mass M in a fluid. We assume that the dynamics of the system can be described by continuum fields consisting of conserved quantities: the concentration c of active particles, the total density ρ = ρ_solvent + Mc of the suspension, assumed constant, and the total momentum density g = ρv, with v the flow velocity. In addition, to describe the possibility of polar order we consider the dynamics of the polarization P. The concentration and polarization vector of the active particles can be written as¹¹where r_n(t) is the position of the center of mass of the n-th swimmer and [small nu, Greek, circumflex]_n(t) is a unit vector along the long axis of the swimmer. We consider here uniaxial swimmers, which are the simplest type of active polar particles.

Although the active suspension is a nonequilibrium system for which a free energy cannot be defined, it is convenient to write the hydrodynamic equations in terms of a “free energy” given by

with δc = c − c₀ the deviation of the concentration from its mean value, c₀. The first two terms on the right hand side of eqn (2) determine the homogeneous steady state of the system. The parameter a₂ ∼ c^* − c changes sign at a critical concentration c^*, while a₄ > 0 The third term describes the energy cost of spatially inhomogeneous deformations of the order parameters, with K a stiffness. For simplicity we have used a one-elastic constant approximation and equated the elastic constants for splay and bend deformations. Finally, the last term describes a coupling between concentration and splay deformation that is only allowed in a polar fluid. We have neglected polar couplings between splay and fluctuations in the magnitude of the order parameter ∼|P|²∇·P that yield advective-type terms and pressure corrections in the equation of motion for the order parameter.

The dynamics of the concentration of active particles is governed by a convection-diffusion equation

∂_tc + ∇·c(v + w₁P) = ∇·(D∇c), (3)

where D_ij = D₀δ_ij + D₁P_iP_j is the anisotropic diffusion tensor. The equation for the polarization is given by:

[∂_t + (v + w₂P)·∇]P_i + ω_ijP_j = λu_ijP_j + γ⁻¹h_i, (4)

where h = −δF/δP plays the role of the nematic molecular field, with:Also, γ is a friction and λ a reactive parameter that controls the coupling between orientation and flow. As in a nematic liquid crystal, order parameter variations are coupled to flow gradients, described by the antisymmetric and symmetric parts of the rate of strain tensor u_ij = (∂_iv_j + ∂_jv_i)/2 and ω_ij = (∂_iv_j − ∂_jv_i)/2.

Finally, w₁ and w₂ are characteristic velocities (taken independent of c), proportional to the self-propulsion speed of the active particles, and describe self-advection of the active particles relative to the fluid. These terms are unique to a polar active fluid. The evolution of the flow velocity is controlled by the Navier–Stokes equation

ρ(∂_t + v·∇)v_i = ∂_jσ_ij, (6)

with ∇·v = 0, as required by incompressibility. The stress tensor σ_ij can be written as

σ_ij = 2ηu_ij + σ^r_ij + σ^a_ij. (7)

The first term on the right hand side of eqn (7) is the dissipative part and we have used an isotropic-viscosity approximation. The second term is the reversible part, given by

with Π the pressure and [small lambda, Greek, macron] another reactive parameter. In a nonequilibrium system the coefficients of the various terms in the constitutive equation for the stress tensor are not in general related to those in the molecular field h, as is instead required in equilibrium. Here we have taken them to be the same. This simplifying assumption does not change the structure of the equation, but only serves to limit the number of independent parameters in the equations. It is a direct consequence of having written the reactive part of the fluxes in terms of the derivative of a free energy. Finally, the last term in eqn (7) is the active contribution. It can be written as the sum of two terms of different symmetry, σ^a_ij = σ^α_ij + σ^β_ij, with

σ^β_ij = βc(∂_iP_j + ∂_jP_i + δ_ij∇·P), (9b)

where σ^α_ij describes tensile/contractile active stresses present in both nematic and polar fluid, while σ^β_ij is unique to polar fluids and describes active stresses due to self-advection and treadmilling.³⁷ Note that α > 0 corresponds to contractile stresses as generated by pullers or by motor-filament systems, while α < 0 describes tensile stresses as generated by most swimming bacteria, such as E. coli.³⁸ In the following we will ignore the polar contribution to the active stress and simply incorporate the effect of polarity via the advective terms proportional to w_i. We will further assume w₁ = w₂ = 2w̃₃/γ = w. We make the system dimensionless by scaling all lengths with the length [small script l] of the active particles, scaling time with the relaxation time of the director τ_p = [small script l]γ/K and scaling stress using the elastic stress σ = K/[small script l]².

III. Homogeneous steady states and their stability

There are two homogeneous, steady state solutions of the hydrodynamic equations, both with constant density c₀. For c₀ < c^*, or a₂₀ ≡ a₂(c₀) > 0, the homogeneous solution is the isotropic state, with P = 0. For c₀ > c^*, or a₂₀ < 0, rotational symmetry is spontaneously broken and there is an ordered polarized state, with . In the numerical work described below we have used the following form for the coefficients a₂ and a₄:

a₂ = c^* − c, a₄ = c^* + c.

This choice ensures that for c₀ → c^*₊ and P₀ ∼ 1 for c₀ ≫ c^*.

Next, we examine the linear stability of the homogeneous states against spatially inhomogeneous fluctuations in the continuum fields. The isotropic state is always linearly stable and will not be discussed further. In order to ensure the incompressibility condition, it is convenient to rewrite the Navier–Stokes equation in terms of vorticity and stream function ψ by expressing:

v_x = ∂_yψ, v_y = −∂_xψ.

The incompressibility condition is then automatically satisfied. The vorticity field is given by

ω = 2ω_xy = ∂_xv_y − ∂_yv_x, (10)

and the stream-function ψ is related to the vorticity ω through a Poisson equation:

∇²ψ = −ω. (11)

The Navier–Stokes equation can be written in terms of ω and ψ as

∂_tω = ηΔω + ∂²_xτ_yx + ∂_xyτ_yy − ∂_yxτ_xx − ∂²_yτ_xy, (12)

where we have separated out the viscous part of the stress tensor and defined τ_ij as

To analyze the stability of the polarized state, we choose the x axis along the direction of order and compactify the notation by introducing a vector φ = {c, P_x, P_y, ω}. The homogeneous polarized state is characterized by φ⁽⁰⁾ = {c₀, P₀, 0, 0}. We then introduce the fluctuations from the homogeneous values by letting

φ(x, t) = φ⁽⁰⁾ + εφ⁽¹⁾(x, t)

and consider the dynamics of fluctuations for ε ≪ 1. To enforce periodic boundary conditions on a square domain, we look for solutions of the form
where q_n = 2πn/L and q_m = 2πm/L. The Fourier components of the stream-function are related to those of the vorticity by:
With this choice the linearized hydrodynamic equations reduce to a set of coupled linear ordinary differential equations for the Fourier modes φ_nm, given by

∂_tφ_nm = A_nmφ_nm. (14)

The matrix A_nm can be expressed in the following block-form

The explicit expression for the matrix A_nm is given in Appendix A. In order for the homogeneous state φ₀ to be stable, the real part of the eigenvalues of the matrix A_nm must be negative. It is useful to examine the behavior of the eigenvalues for (n, m) modes (1, 0) and (0, 1), corresponding to wavevectors longitudinal and transverse to the direction of mean order, respectively. In both cases (δc, P_x) fluctuations decouple from (P_y, ω) fluctuations and the dynamics can be examined by solving two independent 2 × 2 eigenvalue problems.

We first consider the (1, 0) or longitudinal modes. The matrix b₁₀ = 0, hence A₁₀ is block diagonal. The coupled dynamics of c and P_x fluctuations (note that δP_x describes fluctuations in the magnitude of the order parameter) is controlled by the eigenvalues of the 2 × 2 matrix a₁₀, which does not depend on the active parameter α, but only on the advective velocity w. The (10) eigenvalues are given by

with:

C = 2|c₀ − c^*| + (K + D)q²,

q = 2π/L and D₀ = D₁ = D. The fluctuations of concentration and order parameter magnitude in the longitudinal (1, 0) exhibit unstable growth for arbitrarily small values of w as the mean field order-disorder transition is approached from above, i.e., c₀ → c^*+. The real part of the eigenvalues of the matrix a₁₀ is shown in Fig. 1. This instability does not couple to the flow and is precisely the instability obtained in ref. 32, 34 for dry systems consisting of self-propelled particles with aligning interactions moving on a frictional substrate. It occurs for both tensile and contractile systems, but it is unique to polar active fluids. The instability only occurs in a narrow range of densities above the putative mean-field continuous transition point c^* between polarized and isotropic states, as shown in the stability phase diagram of Fig. 2. It has been argued that it may be associated with the renormalization of the order of the transition, which has been shown to be discontinuous in recent simulations of Vicsek-type models.²⁹ We note that the magnitude of the order parameter P_x is not a hydrodynamic variable in the strict sense as its decay rate has the finite value −2|a₂₀| at large wavelengths. At the transition, however, a₂₀ = 0 and fluctuations in P_x become long-lived, driving the instability. In the absence of noise, the numerical solution of continuum equations obtained by coarse-graining the Vicsek model has yielded in this region of the parameters a spatially inhomogeneous state consisting of solitary waves or bands of ordered regions traveling in a disordered background.³² We refer to the region of parameters where this linear instability occurs as longitudinal traveling Vicsek-type waves (V) regime.

Fig. 1 Real part of the eigenvalues s_± of the matrix a₁₀ displaying the coupled instability of (δc, δP_x) for w = 0.5 and c₀ = 0.01.

Fig. 2 Phase-diagram in the plane (c₀ − c^*, α) for w = α and q = 0.2, λ = 0.5, η = D₀ = D₁ = 1. The labeled regions correspond to the homogeneous steady state (HSS), the spontaneous flow regime (SF), the traveling waves regime (TW) and the region of longitudinal traveling waves of Vicsek-type (V).

The coupled dynamics of fluctuations in the transverse component P_y of the order parameter, describing director fluctuations, and vorticity ω yields additional instabilities, both in the (1, 0) and (0, 1) directions, controlled by the d_nm matrix. The first unstable modes are the longitudinal mode (1, 0) and the transverse mode (0, 1). These modes becomes unstable for negative and positive values of the combination α(1 − λ), respectively. These instabilities arise from the coupling of director deformations and flow and occur even for w = 0, corresponding to systems with nematic symmetry. For α(1 − λ) > 0, transverse modes (0, 1) go unstable, corresponding to coupled splay deformations of the director and fluctuations of the v_y velocity. For λ < 1, the instability occurs above a positive critical value of α given by

More precisely this instability occurs for tumbling (|λ| < 1) or disk-like (λ < −1) pullers for α > α₀₁ > 0 and for rod-like (λ > 1) pushers for α < α₀₁ < 0. The interplay between the parameter λ controlling the flow alignment properties of an orientationally ordered suspension and the activity α in controlling the instabilities of active fluids was discussed in ref. 14 and we refer to that work for details. The critical value α₀₁ of transverse or bend instability does not depend on w and is identical to the value obtained for active nematic, with the replacement P₀ → S, where S is the magnitude of the nematic order parameter. This instability has been obtained before in the literature⁴ and is commonly referred to as the generic instability. For α(1 + λ) < 0, corresponding to tensile active stresses as generated by pushers, such as E. coli, longitudinal modes (1, 0) go unstable, corresponding to coupled bend deformations of the director and fluctuations in the v_x velocity. This instability occurs for λ > −1 and α < α₁₀ < 0 and for λ < −1 and α > α₁₀ > 0, with

Notice that the critical value α₁₀ depends on w (a finite w in this case actually suppresses the instability). When w = 0, this coincides with the generic instability discussed in [4] for pushers and active nematic. For finite w, however, the modes that goes unstable is a propagating wave, suggesting that the system may support traveling waves solutions in this region of parameter. This region is labelled at TW in Fig. 2. Finally, the splay/bend instabilities obtained here are simply the 2D generalization of the instability to a spontaneously flowing state discussed in ref. 12 and 14 for nematic and polar active fluids, respectively, in a quasi-one-dimensional strip geometry. Their origin lies in the long-range nature of the hydrodynamic interactions between active swimming particles, as shown in ref. 39. Some additional detail on the eigenvalues of the d_nm matrix can be found in Appendix B.

A. Relation to dry systems

For completeness we compare the behavior of the active suspension considered here to that of active particles on a frictional substrate, referred to so far as dry systems. In this case the only conserved field is the concentration of active particles and the continuum equations are written solely in terms of concentration and polarization fields. The momentum of the system is not conserved and there is no equation for the flow velocity.⁴⁰ The low density longitudinal instability of the ordered state obtained near the mean-field order order-disorder transition (c₀ → c^*+) is unique to polar active fluids with aligning interactions, corresponding for instance to coarse-grained versions of the Vicsek model, and does not couple to the flow velocity. It is present in both dry systems and suspensions. It was discussed for the case polar particles on a substrate in ref. 32 and 34.

In contrast, the splay and bend instabilities obtained at high density are due curvature-induced currents arising from the coupling of director deformation and flow.

One can obtain the limit of dry systems from our general equations for a suspension by adding a drag −ζv to the right hand side of eqn (6), describing the coupling frictional coupling to a substrate. In the limit of large friction ζ, one can then neglect inertial terms and write the solution of eqn (6) to lowest order in the gradients as

where the active stress σ^a_ij is given in eqn (9a) and we must satisfy ∇·v = 0. If eqn (19) is used to eliminate v from the equations for concentration and polarization, one then recovers both the splay and bend instabilities. Note that when eqn (19) is used to eliminate the flow velocity in favor of polarization and density the terms responsible for the splay and bend instabilities are of order (∇²P³). These terms have been neglected in the continuum models of polar fluids on a substrate discussed in the literature, which explain why such instabilities were not obtained to linear order.⁴¹ When considering the full nonlinear equations similar effects are, however, provided by other nonlinearities.

IV. Nonlinear spatio-temporal patterns

In this section we report some results obtained from a numerical solution of eqn (3), (4) and (6). To explore different regimes at increasing level of activity we set w = 0.1α and use α as a variable parameter. Eqn (3), (4) and (6) are then integrated using a vorticity/stream-function finite difference scheme on a collocated grid of lattice spacing Δx = Δy = 0.078. The time integration is performed via a fourth order Runge–Kutta method with time step Δt = 10⁻³. As mentioned in Sec. III, the vorticity/stream-function method requires one to solve a Poisson equation at each time step in order to calculate the components of the flow velocity. This was efficiently preformed with a V-cycle multigrid algorithm. As initial configurations we consider a homogeneous system whose director field was aligned along the x axis subject to a small random perturbation in density and orientation. Thus c = c₀ + ε, θ = ε, and v_x = v_y = 0, where ε is a random number of zero mean and variance 〈ε²〉 = 10⁻². The equations where then integrated from t = 0 to t = 10³, corresponding to 10⁶ time steps. For all the numerical calculations described in this section we used η = c^* = D₀ = D₁ = 1, λ = 0.1 and L = 10, in the units described in Sec. II. Upon varying the activity parameter α five different regimes have been observed: (1) homogeneous stationary state; (2) spontaneous steady-flow; (3) oscillating flow with orthogonal traveling bands; (4) oscillating flow with pairs of traveling vortices and (5) chaotic flow. These regimes are described below in more detail.

A. Homogeneous state and spontaneous steady flow

For small values of |α|, the system quickly relaxes to a homogeneous state with the polarization vector aligned along the x direction. The concentration is equal to the initial concentration c₀ and the polar order parameter equates its equilibrium value . Upon raising α above the critical value (17), the homogeneous state becomes unstable to a non-homogeneous flowing state (top of Fig. 3). The structure of this state is the same of active nematic fluids and consists of two bands flowing in opposite directions.²² The direction of the streamlines is dictated by the initial conditions which, in this case, favor a flow in the x direction. This solution has been intensively discussed in the literature of active matter (see for instance ref. 12–14) and will not be commented here any further.

Fig. 3 The velocity field (left) and the polarization direction (right) superimposed to a density plot of the concentration and the polarization magnitude for α = 1.5 (top) and α = −2.1 (bottom). The colors indicate regions of large (red) and small (green) density and large (red) and small (blue) polarization. For positive values of α the homogeneous state undergoes a splay instability that leads to the formation of two stationary bands. For negative α a bending instability results in the appearance of vertical bands traveling at speed v₀ ∼ w ∼ α (hence from right to left).

B. Traveling waves

For negative values of α larger in magnitude than the critical value (18) (α ≈ −1.19 for our choice of parameters) traveling waves form in the system. As shown in Fig. 3 the waves consist of a shear-flow along the y direction, that translates at constant speed along x: v_x = 0 and v_y = v_y(x − v₀t) with v₀ the wave velocity. The concentration c and the polarization field P follows the flow, thus forming traveling bands across the x direction (Fig. 4 and 5). The polarization direction p, on the other hand, oscillate while keeping the mean orientation orthogonal to the bands. The frequency and the amplitude of the traveling waves are shown in Fig. 6 as a function of α.

The physical origin of the traveling waves in this model of active polar fluids relies on the fact that the polarization field is simultaneously a broken-symmetry variable, a mechanical stress and a physical velocity. Being a broken-symmetry variable, the polarization relaxes toward a non-zero value in the characteristic time scale τ_p (see Sec. II). Correspondingly, shear stress σ^a_xy ∼ αP_xP_y is injected in the system at rate τ⁻¹_a ∼ α. When τ_a ≈ τ_p this residual stress is accommodated in the passive structures leading to an elastic distortion and flow as discussed in Sec. III. Unlike nematic fluids, here the existence of a direction of broken-symmetry results in a directed motion of particles with velocity v^a ∼ wP, hence the pattern produced by the interplay between active stress, liquid crystalline elasticity and flow is advected along P resulting in the traveling waves described above. Note that the pattern generated by the splay instability (top of Fig. 3), is not accompanied by the formation of traveling bands since the system is translational invariant along the polarization direction, hence there is no advective transport: i.e.P·∇ = P_x∂_x = 0.

To elucidate the properties of this solution, it is useful to reduce the hydrodynamic equations to a tractable form. In the following we will see that under drastic but reasonable simplifications, the hydrodynamic equations presented in Sec. II map into a reaction-diffusion-advection system that admits traveling waves solutions. Diffusion-reaction systems have played a fundamental role in the study of biological waves and pattern formation, especially in the context of morphogenesis since the seminal 1952 paper by Turing.⁴²

As a starting point we notice that at the onset of the traveling bands regime, the variations in c and P_x are only slight with respect to those in ω and P_y (see Fig. 4 and 5). This suggests that the dynamics of the traveling bands is dominated by the latter two fields, as also indicated by the linear stability analysis. Thus we approximate c and P_x as constants and focus our attention on the equations for ω and P_y:

[∂_t + (v + wP)·∇]P_y = λ(u_yxP_x + u_yyP_y) − ω_yxP_x + h_y.

Fig. 4 The quantities c, v_y, P_x and P_y along the x axis in the traveling waves regime for α = −2.1. The yellow region indicate the band visible in the bottom panel of Fig. 3. The slant of the yellow region indicate the direction of motion of the band.

Fig. 5 The quantities c, ω, P_x and P_y at the point x = y = L/2 as a function of time in the traveling waves regime for α = −2.1.

Next, we notice that v_x = 0 and that there is no dependence on the y coordinate. Incompressibility ensures then: u_xx = u_yy = 0. Morover u_xy = u_yx = ω_xy = −ω_yx = ω/2. The previous expression and eqn (12) simplify then to:

∂_tω = η∂²_xω + ∂²_xτ_yx. (20b)

Fig. 6 The frequency ν of the waves and the magnitude of P_y as a function of α in the traveling waves regime.

Without loss of generality, we take a₄ = 1 and a₂ + a₄P²_x = −a. Now, in order to understand the origin of the traveling bands, one does not need to consider the full stress tensor τ_yx, but only the active contribution τ_yx = αP_xP_y, which is responsible for injecting shear stress inside the system at the rate α. Thus, taking for simplicity w = α and transforming α → αP_x we arrive to the following set of equations:

Ṗ + αP′ = (a − P²)P + P′′ + bω, (21a)

[small omega, Greek, dot above] = ηω′′ + αP′′, (21b)

where the dot and the prime indicate respectively the time and space derivative, b = P_x(1 + λ)/2 and we have renamed P = P_y for shortness. Eqn (21) have the classical form of a diffusion-reaction-advection system. The cubic nonlinearity is the trademark of excitable dynamical systems and, as discussed in ref. 22, gives the active fluid the behavior of an excitable medium. The advection term αP′, on the other hand, is an exclusive feature of active polar systems and is ultimately responsible for the existence of traveling waves.

For small values of α the solution of eqn (21) is the homogeneous state . A small perturbation of the homogeneous state of the form δP(t) = P_k(t)e^ikx and δω(t) = ω_k(t)e^ikx evolves according to the equations:

Thus, upon lowering α below the critical value α_c ∼ −2aη/b (for small k) the homogeneous state becomes unstable and the solution consists of a wave train traveling along the negative x direction at constant velocity. In order to describe the behavior of the traveling waves at the onset of the transition, we can construct an approximate solution of eqn (21) using the method of harmonic balance. We then look for a solution of (21) of the following form:

with the bar indicating the complex conjugation. Substituting this in (21) and taking only the lowest wave numbers yields the following system of ordinary differential equations for the amplitudes:

These equations admits solutions of the form:

Without loss of generality we can take ϕ₁ = 0 and ϕ₂ = ϕ. Then, substituting this form in (22) we find, after some manipulations, the following set of algebraic equations for the quantities ν, ϕ, |P_k| and |ω_k|:

3|P_k|³ + (k² − a)|P_k| − b|ω_k|cos ϕ = 0, (23a)

αk²|P_k| + (ηk²cos ϕ − νsin ϕ)|ω_k| = 0, (23b)

(kα + ν)|P_k| − b|ω_k|sin ϕ = 0, (23c)

νcos ϕ + ηk²sin ϕ = 0. (23d)

Solving this equations to leading order in α gives:

The larger spatial mode k = 2π/L thus evolves in time through traveling waves of the form:

Fig. 7 shows a plot of the amplitudes |P_k| and |ω_k| for the largest wavelength corresponding to k = 2π/L in eqn (24) together with those obtained from a numerical solution of the simplified eqn (21).

Fig. 7 The amplitudes of the traveling waves resulting from the numerical solution of the simplified model (21) together with the predictions of the harmonic balance analysis (24) for the largest spatial mode k = 2π/L (solid lines). At α ≈ −2aη/b the system undergoes a subcritical Hopf bifurcation with a bistability region in the range −2aη/b < α < 0. The red squares and blue/green circles correspond respectively to the stable homogeneous steady state and the stable large-amplitude limit cycle.

The combined linear stability and harmonic balance analysis show that the Hopf bifurcation that leads to the formation of traveling waves is subcritical with bistable region in the range −2aη/b < α < 0. As the activity parameter α is lowered below zero, the system goes through a subcritical bifurcation where there is a stable stationary state and a stable large-amplitude limit cycle with an unstable periodic orbit acting as the separatrix between the two states. The existence of bistability at small negative values of α could have significant effects in the collective motion of active polar suspensions such as the presence of hysteresis when the amount of activity is cycled across the bifurcation point.

C. Traveling vortices and chaos

If we further decrease α, the traveling waves described in the previous section become unstable. In the new regime, the concentration is distributed in the form of kidney-shaped clusters that travels along the positive x direction at constant speed (Fig. 8 right). The flow field, in turn, consists of a pair of oppositely spinning vortices whose centers also move at constant speed following the concentration clusters as well as two saddle stagnation points due to the toroidal topology of the domain (Fig. 8 left). The order parameter appears concentrated along “valleys” whose curved conformation partially resembles the banded structure of the previous solution. Once again the relative variation in polarization is larger than that in concentration so that the value of the order parameter at the bottom of the valleys is almost 50% lower than along the ridges.

Fig. 8 The velocity field (left) and the polarization direction (right) superimposed to a density plot of the concentration and the polarization magnitude at two different times for α = −5. The flow is characterized by two vortices and two stagnation saddle points traveling across the system from left to right.

The formation of the vortices is associated with a break-down of the translational invariance along the longitudinal direction of the bands. At the onset of vortex formation the vorticity function for the spatial mode of largest wave-length is given heuristically by:

ω_k(x, t) ≈ ω_⊥cos k(x − αP_xt) − ω_∥cos ky (25)

whose corresponding velocity field consists of two vortices and two saddles traveling along the x direction with velocity αP_x. From the numerical data shown in Fig. 9 we see that P_x is everywhere negative, hence the vortices move toward the positive x direction. In contrast with the traveling waves, where the polarization vector oscillate around its initial orientation (x̂ in this case), the formation of traveling vortices is anticipated by an inversion of the polarization vector (i.e.p ∼ −x̂). The origin of this inversion is likely due to the short time scale at which the active stress is initially injected in the system, which makes the polarization field undergo dramatic distortions before it is able to settle on a limit cycle.

Upon further decreasing α, the traveling vortices become unstable and the system undergoes a transition to a chaotic regime (Fig. 10). The polarization valleys now continuously form and disappear in unpredictable fashion. The flow still exhibits larger vortices, which however move chaotically across the sample with variable direction and magnitude. Interestingly, the polarization direction is continuously distorted and does not posses regions of partial alignment. This differs from the spatio-temporal chaos found active nematics²² where the chaotic flow is organized in grains of uniform orientation separated by chaotically moving grain boundaries. The nonlinear spatiotemporal patterns described in this section are summarized in the phase-diagram of Fig. 11. The Vicsek-type traveling waves shown in the linear phase diagram of Fig. 2 are also obtained from the numerical solution of the nonlinear equations, but at a much lower density that shown in Fig. 11.

Fig. 9 The quantities c, ω, P_x and P_y at the point x = y = L/2 as a function of time in the traveling vortices for α = −5.

Fig. 10 The velocity field (left) and the polarization direction (right) superimposed to a density plot of the concentration and the polarization magnitude for α = −6. The flow consists of large vortices that move chaotically across the system with variable direction and magnitude.

Fig. 11 Phase diagram of the nonlinear spatiotemporal patterns: homogeneous steady state (HSS), traveling waves (TW), traveling vortices (TV) and chaos (CH). The dots have been obtained from the numerical integration of eqn (3), (4) and (6) for λ = 0.1, η = D₀ = D₁ = 1, L = 10 and w = 0.1α. The solid blue line separating the homogeneous steady state from the traveling waves regime is given by eqn (18). The dashed line separating the traveling waves and the traveling vortices regime has been inferred from the data. The color gradient between the traveling vortices and the chaotic regime indicates a fuzzy boundary line.

V. Discussion

We have analyzed the complex spatio-temporal patterns of an active polar fluid in two dimensions. The continuum model considered is appropriate to describe systems such as bacterial suspensions, consisting of self-propelled particles in a solvent. The self-propulsion of the active units yields a number of advective terms in the hydrodynamic equations that are unique to polar fluids. These advective terms are responsible for the traveling nature of the large scale structures, which include traveling stripes oriented normal to the direction of propagation and traveling vortices. These structures are remarkably stable in a wide region of parameters and resemble the patterns seen recently in actin motility assays at high actin density, suggesting that systems may indeed have an underlying polar symmetry. We stress that the transitions between regimes characterized by different stable large-scale structures can only be obtained by allowing for spatial variations of the order parameter. A novel feature of the present work is that the transitions between these various regimes appear to be discontinuous, allowing for the possibility of coexistence and hysteresis. Discontinuous transitions between disordered and ordered states have also been seen in dry active systems via numerical studies of the Vicsek model.²⁹ Understanding the nature and order of the transition in both suspensions and active particles on substrates will, however, require further investigation.

The present analysis is limited to two dimensional systems and may be applicable to situations like motility assays where the system is quasi-2D. One should, however, be cautious about generalizing the results found in 2D to 3D systems, as new patterns may arise one the polarization field is allowed to rotate out of the plane. Very recent work by Fielding et al.²³ has indeed shown that the behavior of a sheared active nematic in two dimensions can be quite different from that in one dimension.

Our work, together with recent work on pattern formation inactive nematic, suggests a strong analogy between the spatio-temporal dynamics of active fluids and the behavior of excitable systems. In both cases the complex behavior obtained from the nonlinear hydrodynamic equations can be reproduced qualitatively by simplified lower dimensional models that resemble reaction-diffusion or reaction-diffusion-advection equations familiar from the study of excitable systems. This mapping may provide a useful framework for understanding and classifying active systems.

Appendix A

In this appendix we give the explicit expressions for the elements of the matrix A_nm given in eqn (15). The various block elements are given bywith:

C²²_nm = 2c₀P₀[q_mq_nλ + α(q²_m + q_nq_m − q²_n)] + q_mq_nλP₀[(q²_m + q²_n) − 2c^*]

and:
where a₂₀ = a₂(c₀) = c^* − c₀ < 0, D₀ = D₁ = D and the prime denotes a derivative with respect to concentration, i.e., a^′_2,4 = (∂a_2.4/∂c)_c=c0, which gives −(a^′₂ + a^′₄P²₀) = 2c^*/(c₀ + c^*).

Appendix B

Analytical expression for the values of α associated with these instabilities can be found by diagonalizing the matrix d_nm. For both the longitudinal (10) and the transverse (01) directions, it is the coupled modes of P_y and ω that go unstable. For α > 0 the modes go unstable along the 01 direction in this case P_y corresponds to splay deformations and the vorticity describes fluctuations in the v_x velocity. For α < 0 the modes go unstable along the 10 direction in this case P_y corresponds to bend deformations and the vorticity describes fluctuations in the v_y velocity. Both these instabilities exist also when w = 0. The matrices d₁₀ and d₀₁ are given by:
and:
with q = 2π/L. The corresponding eigenvalues are:

Acknowledgements

MCM was supported by the NSF on grants DMR-0806511 and DMR-1004789. LG was supported by the Harvard-NSF MRSEC, the Harvard-Kavli Nano-Bio Science and Technology Center and the Wyss Institute. We thank Aparna Baskaran and Paolo Paoletti for illuminating discussions.

References

1. H. Gruler, U. Dewald and M. Eberhardt, Eur. Phys. J. B, 1999, 11, 187.
2. J. Toner and Y. Tu, Phys. Rev. Lett., 1995, 75, 4326.
3. J. Toner, Y. Tu and S. Ramaswamy, Ann. Phys., 2005, 318, 170.
4. R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101.
5. S. A. Edwards and J. M. Yeomans, Europhys. Lett., 2009, 85, 18008.
6. K. Kruse, J. F. Joanny, F. Jülicher, J. Prost and K. Sekimoto, Phys. Rev. Lett., 2004, 92, 078101.
7. K. Kruse, J. F. Joanny, F. Jülicher, J. Prost and K. Sekimoto, Eur. Phys. J. E, 2005, 16, 516.
8. D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2007, 99, 058102.
9. D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2008, 100, 178103.
10. T. B. Liverpool and M. C. Marchetti, Phys. Rev. Lett., 2003, 90, 138102.
11. T. B. Liverpool and M. C. Marchetti, in Cell Motility, ed. P. Lenz, Springer, New York, 2007.
12. R. Voituriez, J. F. Joanny and J. Prost, Europhys. Lett., 2005, 70, 118102.
13. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E, 2007, 76, 031921.
14. L. Giomi, M. C. Marchetti and T. B. Liverpool, Phys. Rev. Lett., 2008, 101, 198101.
15. Y. Hatwalne, S. Ramaswamy, M. Rao and R. A. Simha, Phys. Rev. Lett., 2004, 92, 118101.
16. T. B. Liverpool and M. C. Marchetti, Phys. Rev. Lett., 2006, 97, 268101.
17. M. E. Cates, S. M. Fielding, D. Marenduzzo, E. Orlandini and J. M. Yeomans, Phys. Rev. Lett., 2008, 101, 068102.
18. L. Giomi, T. B. Liverpool and M. C. Marchetti, Phys. Rev. E, 2010, 81, 051908.
19. D. Saintillan, Phys. Rev. E, 2010, 81, 056307.
20. S. Rafaï, L. Jibitu and P. Peyla, Phys. Rev. Lett., 2010, 104, 098102.
21. A. Sokolov and I. S. Aranson, Phys. Rev. Lett., 2009, 103, 148101.
22. L. Giomi, L. Mahadevan, B. Chakraborty and M. F. Hagan, Phys. Rev. Lett., 2011, 106, 218101.
23. S. M. Fielding, D. Marenduzzo and M. E. Cates, Phys. Rev. E, 2011, 83, 041910.
24. V. Schaller, C. Weber, C. Semmrich, E. Frey and A. R. Bausch, Nature, 2010, 467, 73.
25. T. Butt, T. Mufti, A. Humayun, P. B. Rosenthal, S. Khan, S. Khan and J. E. Molloy, J. Biol. Chem., 2010, 285, 4964.
26. D. Kaiser, Nat. Rev. Microbiol., 2003, 1, 45.
27. D. R. Zusman, A. E. Scott, Z. Yang and J. R. Kirby, Nat. Rev. Microbiol., 2007, 5, 862.
28. F. Peruani, A. Deutsch and M. Bär, Phys. Rev. E, 2006, 74, 030904(R).
29. H. Chaté, F. Ginelli, G. Grégoire and F. Raynaud, Phys. Rev. E, 2008, 77, 046113.
30. F. Ginelli, F. Peruani, M. Bär and H. Chaté, Phys. Rev. Lett., 2010, 104, 184502.
31. Y. Yang, V. Marceau and G. Gompper, Phys. Rev. E, 2010, 82, 031904.
32. E. Bertin, M. Droz and G. Grégoire, J. Phys. A: Math. Theor., 2009, 42, 445001.
33. T. Ihle, Phys. Rev. E, 2011, 83, 030901(R).
34. S. Mishra, A. Baskaran and M. C. Marchetti, Phys. Rev. E, 2010, 81, 061816.
35. R. Voituriez, J. F. Joanny and J. Prost, Phys. Rev. Lett., 2006, 96, 028102.
36. D. Blankschtein and R. M. Hornreich, Phys. Rev. B, 1985, 32, 3214.
37. The parameters β and α have different dimensions from those used in ref. 11 as here we have incorporated an additional factor of [small script l]² in their definition.
38. We use the sign convention introduced in ref. 10, which is opposite to that used on much later work in acto-myosin systems (e.g., ref. 7).
39. A. Baskaran and M. C. Marchetti, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 15567.
40. The hydrodynamics of a dry polar active systems was first considered by Toner and Tu in ref. 2. The velocity used by these authors as the order parameter is simply proportional to our polarization P and should not be confused with the flow velocity introduced in the present work.
41. A splay instability is obtained in ref. 34 by a different mechanism.
42. M. Turing, Philos. Trans. R. Soc. London, Ser. B, 1952, 237, 37.

---