DOI:
10.1039/C1SM06077E
(Paper)
Soft Matter, 2012,
8, 129-139
Polar patterns in active fluids
Received
8th June 2011
, Accepted 16th September 2011
First published on 31st October 2011
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”.1 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 suspensions26,27 and simulations of self-propelled hard rods28–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.22
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 Tu2,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.11 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 models29 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.35 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.36 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
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 as11 |  | (1a) |
|  | (1b) |
where rn(t) is the position of the center of mass of the n-th swimmer and
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
|  | (2) |
with
δc =
c −
c0 the deviation of the concentration from its mean value,
c0. The first two terms on the right hand side of
eqn (2) determine the homogeneous steady state of the system. The parameter
a2 ∼
c* −
c changes sign at a critical concentration
c*, while
a4 > 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|
2∇·
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 + w1P) = ∇·(D∇c), | (3) |
where
Dij =
D0δij +
D1PiPj is the anisotropic diffusion tensor. The equation for the polarization is given by:
| [∂t + (v + w2P)·∇]Pi + ωijPj = λuijPj + γ−1hi, | (4) |
where
h = −
δF/
δP plays the role of the nematic molecular field, with:
|  | (5) |
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
uij = (∂
ivj + ∂
jvi)/2 and
ωij = (∂
ivj − ∂
jvi)/2.
Finally, w1 and w2 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·∇)vi = ∂jσij, | (6) |
with
∇·
v = 0, as required by incompressibility. The stress tensor
σij can be written as
| σij = 2ηuij + σrij + σaij. | (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
|  | (8) |
with Π the pressure and
![[small lambda, Greek, macron]](https://www.rsc.org/images/entities/char_e0cc.gif)
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,
σaij =
σαij +
σβij, with
|  | (9a) |
| σβij = βc(∂iPj + ∂jPi + δ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.
37 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.
38 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
wi. We will further assume
w1 =
w2 = 2
3/
γ =
w. We make the system dimensionless by scaling all lengths with the length
![[small script l]](https://www.rsc.org/images/entities/char_e146.gif)
of the active particles, scaling time with the relaxation time of the director
τp =
γ/
K and scaling stress using the elastic stress
σ =
K/
2.
III. Homogeneous steady states and their stability
There are two homogeneous, steady state solutions of the hydrodynamic equations, both with constant density c0. For c0 < c*, or a20 ≡ a2(c0) > 0, the homogeneous solution is the isotropic state, with P = 0. For c0 > c*, or a20 < 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 a2 and a4:
a2 = c* − c, a4 = c* + c. |
This choice ensures that
for c0 → c*+ and P0 ∼ 1 for c0 ≫ 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:
The incompressibility condition is then automatically satisfied. The vorticity field is given by
| ω = 2ωxy = ∂xvy − ∂yvx, | (10) |
and the stream-function
ψ is related to the vorticity
ω through a Poisson equation:
The Navier–Stokes equation can be written in terms of ω and ψ as
| ∂tω = ηΔω + ∂2xτyx + ∂xyτyy − ∂yxτxx − ∂2yτxy, | (12) |
where we have separated out the viscous part of the stress tensor and defined
τij as
|  | (13) |
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, Px, Py, ω}. The homogeneous polarized state is characterized by φ(0) = {c0, P0, 0, 0}. We then introduce the fluctuations from the homogeneous values by letting
φ(x, t) = φ(0) + εφ(1)(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
qn = 2π
n/
L and
qm = 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
The matrix Anm can be expressed in the following block-form
|  | (15) |
The explicit expression for the matrix Anm is given in Appendix A. In order for the homogeneous state φ0 to be stable, the real part of the eigenvalues of the matrix Anm 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, Px) fluctuations decouple from (Py, ω) 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 b10 = 0, hence A10 is block diagonal. The coupled dynamics of c and Px fluctuations (note that δPx describes fluctuations in the magnitude of the order parameter) is controlled by the eigenvalues of the 2 × 2 matrix a10, which does not depend on the active parameter α, but only on the advective velocity w. The (10) eigenvalues are given by
|  | (16) |
with:
C = 2|c0 − c*| + (K + D)q2, |
q = 2π/
L and
D0 =
D1 =
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.,
c0 →
c*+. The real part of the eigenvalues of the matrix
a10 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.
29 We note that the magnitude of the order parameter
Px is not a hydrodynamic variable in the strict sense as its decay rate has the finite value −2|
a20| at large wavelengths. At the transition, however,
a20 = 0 and fluctuations in
Px 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.
32 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 a10 displaying the coupled instability of (δc, δPx) for w = 0.5 and c0 = 0.01. | |
 |
| Fig. 2
Phase-diagram in the plane (c0 − c*, α) for w = α and q = 0.2, λ = 0.5, η = D0 = D1 = 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 Py of the order parameter, describing director fluctuations, and vorticity ω yields additional instabilities, both in the (1, 0) and (0, 1) directions, controlled by the dnm 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 vy velocity. For λ < 1, the instability occurs above a positive critical value of α given by
|  | (17) |
More precisely this instability occurs for tumbling (|λ| < 1) or disk-like (λ < −1) pullers for α > α01 > 0 and for rod-like (λ > 1) pushers for α < α01 < 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 α01 of transverse or bend instability does not depend on w and is identical to the value obtained for active nematic, with the replacement P0 → S, where S is the magnitude of the nematic order parameter. This instability has been obtained before in the literature4 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 vx velocity. This instability occurs for λ > −1 and α < α10 < 0 and for λ < −1 and α > α10 > 0, with
|  | (18) |
Notice that the critical value α10 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 dnm 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.40 The low density longitudinal instability of the ordered state obtained near the mean-field order order-disorder transition (c0 → 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
|  | (19) |
where the active stress
σaij 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 (∇
2P3). 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.
41 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−3. 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 = c0 + ε, θ = ε,
and vx = vy = 0, where ε is a random number of zero mean and variance 〈ε2〉 = 10−2. The equations where then integrated from t = 0 to t = 103, corresponding to 106 time steps. For all the numerical calculations described in this section we used η = c* = D0 = D1 = 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 c0 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.22 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 v0 ∼ 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: vx = 0 and vy = vy(x − v0t) with v0 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 σaxy ∼ αPxPy is injected in the system at rate τ−1a ∼ α. 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 va ∼ 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·∇ = Px∂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.42
As a starting point we notice that at the onset of the traveling bands regime, the variations in c and Px are only slight with respect to those in ω and Py (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 Px as constants and focus our attention on the equations for ω and Py:
[∂t + (v + wP)·∇]Py = λ(uyxPx + uyyPy) − ωyxPx + hy. |
 |
| Fig. 4 The quantities c, vy, Px and Py 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, ω, Px and Py at the point x = y = L/2 as a function of time in the traveling waves regime for α = −2.1. | |
Next, we notice that vx = 0 and that there is no dependence on the y coordinate. Incompressibility ensures then: uxx = uyy = 0. Morover uxy = uyx = ωxy = −ωyx = ω/2. The previous expression and eqn (12) simplify then to:
|  | (20a) |
| ∂tω = η∂2xω + ∂2xτyx. | (20b) |
 |
| Fig. 6 The frequency ν of the waves and the magnitude of Py as a function of α in the traveling waves regime. | |
Without loss of generality, we take a4 = 1 and a2 + a4P2x = −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 = αPxPy, which is responsible for injecting shear stress inside the system at the rate α. Thus, taking for simplicity w = α and transforming α → αPx we arrive to the following set of equations:
| Ṗ + αP′ = (a − P2)P + P′′ + bω, | (21a) |
| = ηω′′ + αP′′, | (21b) |
where the dot and the prime indicate respectively the time and space derivative,
b =
Px(1 + λ)/2 and we have renamed
P =
Py 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) = Pk(t)eikx and δω(t) = ωk(t)eikx 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:
|  | (22a) |
|  | (22b) |
|  | (22c) |
|  | (22d) |
These equations admits solutions of the form:
Without loss of generality we can take ϕ1 = 0 and ϕ2 = ϕ. Then, substituting this form in (22) we find, after some manipulations, the following set of algebraic equations for the quantities ν, ϕ, |Pk| and |ωk|:
| 3|Pk|3 + (k2 − a)|Pk| − b|ωk|cos ϕ = 0, | (23a) |
| αk2|Pk| + (ηk2cos ϕ − νsin ϕ)|ωk| = 0, | (23b) |
| (kα + ν)|Pk| − b|ωk|sin ϕ = 0, | (23c) |
| νcos ϕ + ηk2sin ϕ = 0. | (23d) |
Solving this equations to leading order in α gives:
|  | (24) |
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 |Pk| 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 − αPxt) − ω∥cos ky | (25) |
whose corresponding velocity field consists of two vortices and two saddles traveling along the
x direction with velocity
αPx. From the numerical data shown in
Fig. 9 we see that
Px 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 with combining circumflex]](https://www.rsc.org/images/entities/b_char_0078_0302.gif)
in this case), the formation of traveling vortices is anticipated by an inversion of the polarization vector (
i.e.p ∼ −
![[x with combining circumflex]](https://www.rsc.org/images/entities/b_char_0078_0302.gif)
). 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 nematics22 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, ω, Px and Py 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, η = D0 = D1 = 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.29 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.23 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 Anm given in eqn (15). The various block elements are given by |  | (A1) |
with:
C22nm = 2c0P0[qmqnλ + α(q2m + qnqm − q2n)] + qmqnλP0[(q2m + q2n) − 2c*] |
and:
where a20 = a2(c0) = c* − c0 < 0, D0 = D1 = D and the prime denotes a derivative with respect to concentration, i.e., a′2,4 = (∂a2.4/∂c)c=c0, which gives −(a′2 + a′4P20) = 2c*/(c0 + c*).
Appendix B
Analytical expression for the values of α associated with these instabilities can be found by diagonalizing the matrix dnm. For both the longitudinal (10) and the transverse (01) directions, it is the coupled modes of Py and ω that go unstable. For α > 0 the modes go unstable along the 01 direction in this case Py corresponds to splay deformations and the vorticity describes fluctuations in the vx velocity. For α < 0 the modes go unstable along the 10 direction in this case Py corresponds to bend deformations and the vorticity describes fluctuations in the vy velocity. Both these instabilities exist also when w = 0. The matrices d10 and d01 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
- H. Gruler, U. Dewald and M. Eberhardt, Eur. Phys. J. B, 1999, 11, 187 CrossRef CAS.
- J. Toner and Y. Tu, Phys. Rev. Lett., 1995, 75, 4326 CrossRef CAS.
- J. Toner, Y. Tu and S. Ramaswamy, Ann. Phys., 2005, 318, 170 CrossRef CAS.
- R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101 CrossRef.
- S. A. Edwards and J. M. Yeomans, Europhys. Lett., 2009, 85, 18008 CrossRef.
- K. Kruse, J. F. Joanny, F. Jülicher, J. Prost and K. Sekimoto, Phys. Rev. Lett., 2004, 92, 078101 CrossRef CAS.
- K. Kruse, J. F. Joanny, F. Jülicher, J. Prost and K. Sekimoto, Eur. Phys. J. E, 2005, 16, 516 CrossRef.
- D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2007, 99, 058102 CrossRef.
- D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2008, 100, 178103 CrossRef.
- T. B. Liverpool and M. C. Marchetti, Phys. Rev. Lett., 2003, 90, 138102 CrossRef.
-
T. B. Liverpool and M. C. Marchetti, in Cell Motility, ed. P. Lenz, Springer, New York, 2007 Search PubMed.
- R. Voituriez, J. F. Joanny and J. Prost, Europhys. Lett., 2005, 70, 118102 CrossRef.
- D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E, 2007, 76, 031921 CrossRef CAS.
- L. Giomi, M. C. Marchetti and T. B. Liverpool, Phys. Rev. Lett., 2008, 101, 198101 CrossRef.
- Y. Hatwalne, S. Ramaswamy, M. Rao and R. A. Simha, Phys. Rev. Lett., 2004, 92, 118101 CrossRef.
- T. B. Liverpool and M. C. Marchetti, Phys. Rev. Lett., 2006, 97, 268101 CrossRef CAS.
- M. E. Cates, S. M. Fielding, D. Marenduzzo, E. Orlandini and J. M. Yeomans, Phys. Rev. Lett., 2008, 101, 068102 CrossRef CAS.
- L. Giomi, T. B. Liverpool and M. C. Marchetti, Phys. Rev. E, 2010, 81, 051908 CrossRef.
- D. Saintillan, Phys. Rev. E, 2010, 81, 056307 CrossRef.
- S. Rafaï, L. Jibitu and P. Peyla, Phys. Rev. Lett., 2010, 104, 098102 CrossRef.
- A. Sokolov and I. S. Aranson, Phys. Rev. Lett., 2009, 103, 148101 CrossRef.
- L. Giomi, L. Mahadevan, B. Chakraborty and M. F. Hagan, Phys. Rev. Lett., 2011, 106, 218101 CrossRef CAS.
- S. M. Fielding, D. Marenduzzo and M. E. Cates, Phys. Rev. E, 2011, 83, 041910 CrossRef CAS.
- V. Schaller, C. Weber, C. Semmrich, E. Frey and A. R. Bausch, Nature, 2010, 467, 73 CrossRef CAS.
- T. Butt, T. Mufti, A. Humayun, P. B. Rosenthal, S. Khan, S. Khan and J. E. Molloy, J. Biol. Chem., 2010, 285, 4964 CrossRef CAS.
- D. Kaiser, Nat. Rev. Microbiol., 2003, 1, 45 CrossRef CAS.
- D. R. Zusman, A. E. Scott, Z. Yang and J. R. Kirby, Nat. Rev. Microbiol., 2007, 5, 862 CrossRef CAS.
- F. Peruani, A. Deutsch and M. Bär, Phys. Rev.
E, 2006, 74, 030904(R) CrossRef.
- H. Chaté, F. Ginelli, G. Grégoire and F. Raynaud, Phys. Rev. E, 2008, 77, 046113 CrossRef.
- F. Ginelli, F. Peruani, M. Bär and H. Chaté, Phys. Rev. Lett., 2010, 104, 184502 CrossRef.
- Y. Yang, V. Marceau and G. Gompper, Phys. Rev. E, 2010, 82, 031904 CrossRef.
- E. Bertin, M. Droz and G. Grégoire, J. Phys. A: Math. Theor., 2009, 42, 445001 CrossRef.
- T. Ihle, Phys. Rev. E, 2011, 83, 030901(R) CrossRef.
- S. Mishra, A. Baskaran and M. C. Marchetti, Phys. Rev. E, 2010, 81, 061816 Search PubMed.
- R. Voituriez, J. F. Joanny and J. Prost, Phys. Rev. Lett., 2006, 96, 028102 CrossRef CAS.
- D. Blankschtein and R. M. Hornreich, Phys. Rev. B, 1985, 32, 3214 CrossRef CAS.
- The parameters β and α have different dimensions from those used in ref. 11 as here we have incorporated an additional factor of
2 in their definition.
- 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).
- A. Baskaran and M. C. Marchetti, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 15567 CrossRef CAS.
- 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.
- A splay instability is obtained in ref. 34 by a different mechanism.
- M. Turing, Philos. Trans. R. Soc. London, Ser. B, 1952, 237, 37 CrossRef.
|
This journal is © The Royal Society of Chemistry 2012 |
Click here to see how this site uses Cookies. View our privacy policy here.