Alexandre P Solon

Jean-Baptiste Caussin

Denis Bartolo

Hugues Chaté

Julien Tailleur

Pattern formation in flocking models: A hydrodynamic description

Keywords: number(s): 05.40.-a, 87.23.Cc, 05.65.+b, 47.27.ed

We study in detail the hydrodynamic theories describing the transition to collective motion in polar active matter, exemplified by the Vicsek and active Ising models. Using a simple phenomenological theory, we show the existence of an infinity of propagative solutions, describing both phase and microphase separation, that we fully characterize. We also show that the same results hold specifically in the hydrodynamic equations derived in the literature for the active Ising model and for a simplified version of the Vicsek model. We then study numerically the linear stability of these solutions. We show that stable ones constitute only a small fraction of them, which, however, includes all existing types. We further argue that, in practice, a coarsening mechanism leads towards phase-separated solutions. Finally, we construct the phase diagrams of the hydrodynamic equations proposed to qualitatively describe the Vicsek and active Ising models and connect our results to the phenomenology of the corresponding microscopic models.

I. INTRODUCTION

Collective motion is the ability of large groups of motile agents to move coherently on scales much larger than their individual sizes. It is encountered at all scales in nature, from macroscopic animal groups, such as bird flocks, fish schools, or sheep herds, down to the cellular scale, where the collective migration of cells [1] or bacteria [2] is commonly observed. At the subcellular level, in vitro motility assays of actin filaments [3] or microtubules [4] have shown the spontaneous emergence of large vortices. Collective motion is also observed in ensembles of human-made motile particles such as shaken polar grains [5], colloidal rollers [6], self-propelled droplets [7], or assemblies of polymers and molecular motors [3,4,8]. Despite the differences in their propulsion and interaction mechanisms, these seemingly very different systems share common macroscopic behaviors that can be captured by minimal physical models. Of particular interest is the emergence of directed collective motion, which was first addressed in this context in a seminal work by Vicsek and coworkers [9]. The Vicsek model consists of point particles moving at constant speed and aligning imperfectly with the direction of motion of their neighbors. When the error on the alignment interaction is decreased, or the density of particles is increased, a genuine phase transition from a disordered to a symmetry-broken state is observed. This flocking transition gives rise to an emergent ordered phase, with true long-range polar order even in two dimensions (2D), where all the particles propel on average along the same direction. Toner and Tu showed analytically, using a phenomenological fluctuating hydrodynamic description, how this ordered state, which would be forbidden by the Mermin-Wagner theorem at equilibrium [10], is stabilized by self-propulsion [11]. The transition to collective motion in the Vicsek model has a richer phenomenology than originally thought. As first pointed out numerically in Ref. [12], at the onset of collective motion, translational symmetry is broken as well. In periodic simulation boxes, high-density ordered bands of particles move coherently through a low-density disordered background. The transition between these bands and the homogeneous disordered profile is discontinuous, with metastability and hysteresis loops. These spatial patterns and the first-order nature of the transition can be encompassed in a wider framework, which describes the emergence to collective motion as a liquid-gas phase separation [13,14]. The traveling bands result from the phase coexistence between a disordered gas and an ordered polarized liquid. This framework captures many of the characteristics of the transition, from the scaling of the order parameter to the shape of the phase diagram. This phase-separation picture is robust to the very details of the propulsion and interaction mechanisms. More specifically, it has also been quantitatively demonstrated in the active Ising model [13] in which particles can diffuse in a 2D space but self-propel, and align, only along one axis. However, the specifics of the emergent spatial patterns and the type of phase separation depend on the symmetry of the orientational degrees of freedom. While the active Ising model model shows a bulk phase separation, the Vicsek model is akin to an active XY model and is associated with a microphase separation where the coherently moving polar patterns self-organize into smectic structures [14] (see Fig. 1).

In this paper, building on the two prototypical models that are the Vicsek model and the active Ising model, we provide a comprehensive description of the emergent patterns found at the onset of the flocking transition from a hydrodynamic perspective. We first recall the definitions and phenomenologies of these two models in Sec. II. In Sec. III, we provide a simplified hydrodynamic description of the flocking models. In line with Refs. [15,16], we show that these models support nonlinear propagative solutions whose shape is described using a mapping onto the trajectories of pointlike particles in one-dimensional potentials. Finding such solutions thus reduces to a classical mechanics problem with one degree of freedom. For given values of all the hydrodynamic coefficients, and hence of all underlying microscopic parameters, we find an infinity of solutions, describing both phase and microphase separations, that we fully characterize. We then show that the same results hold specifically for the hydrodynamic equations explicitly derived for the active Ising model [13] and for a simplified version of the Vicsek model [15]. Next, we investigate the linear stability of these solutions as solutions of the hydrodynamic equations in Sec. VI and their coarsening dynamics in Sec. VII. Finally, we provide full phase diagrams constructed from the hydrodynamic model in Sec. VIII. We close by discussing the similarities and differences with the phenomenology of the agent-based models and conjecture on the role of the hydrodynamic noise in the selection of the band patterns.

II. PHENOMENOLOGY OF MICROSCOPIC MODELS

Let us first briefly recall the phenomenology of the Vicsek and active Ising models. They are both based on the same two ingredients: Self-propulsion and a local alignment rule. The major differences between the two models are thus the symmetries of the alignment interaction and of the direction of motion.

A. Vicsek model

In the Vicsek model [9], N pointlike particles, labeled by an index i, move at constant speed v 0 on a rectangular plane with periodic boundary conditions. At each discrete time step t = 1, the headings θ i of all particles are updated in parallel according to

θ i (t + 1) = θ j (t) j ∈N i + η ξ t i , (1)
where N i is the disk of unit radius around particle i, ξ t i a random angle drawn uniformly in [-π,π], and η sets the level of noise, playing a role akin to that of a temperature in a ferromagnetic XY model. Then particles hop along their new headings:

r i (t + 1) = r i (t) + v 0 e t+1 i
, where e t+1 i is the unit vector pointing in direction given by θ i (t + 1).

B. Active Ising model

In the active Ising model [13], particles carry a spin ±1 and move on a 2D lattice with periodic boundary conditions. Their dynamics depend on the sign of their spin: A particle with spin s jumps to the site on its right at rate D(1 + sε) and to the site on its left at rate D(1sε), where 0 ε 1 measures the bias on the diffusion. On average, +1 particles thus self-propel to the right and -1 particles to the left at a mean velocity v 0 ≡ 2Dε. Both types of particles diffuse symetrically at rate D in the vertical direction.

The alignment interaction is purely local. On a site i, a particle flips its spin s at rate

W i (s → -s) = exp - s T m i ρ i , ( 2
)
where T is a temperature and m i and ρ i are the magnetization and number of particles on site i. (An arbitrary number of particles is allowed on each site since there is no excluded volume interaction.)

C. A liquid-gas phase transition

The phase diagrams in the temperature-density (or noise-density) ensemble are shown for both models in Fig. 2, highlighting their similarity. At high temperature (or noise) or low density both systems are in a homogeneous disordered gas state. At low temperature (or noise) and high density they are homogeneous and ordered; in these liquid phases, all particles move in average in the same direction. In the central region of the phase diagram, inhomogeneous profiles are observed, with liquid domains moving in a disordered gaseous background.

The phase transitions of both models have all the features of a liquid-gas transition, exhibiting metastability and hysteresis close to the transition lines [12][13][14]. The main difference between the two models lies in the coexistence region: In the active Ising model, the particles phase separate in a gaseous background and an ordered liquid band, both of macroscopic sizes [13]. The coexisting densities depend only on temperature and bias but not on the average density; in the coexistence region, increasing the density at fixed T ,ε thus results in larger and larger liquid domains whose density remains constant, as shown in Fig. 1. Conversely, in the Vicsek model, the system forms arrays of ordered bands arranged periodically in space which have a finite width along their direction of motion: A micro-phase separation occurs [14]. As shown in Fig. 1, increasing the density at constant noise, the number of bands increases but their shape does not change [14].

Three types of propagating patterns can thus be observed at phase coexistence, all shown in Fig. 1: (i) localized compact excitations, (ii) Smectic microphases, and (iii) phase-separated polar liquid domains. In the vicinity of the left coexistence line, collective motion emerges in the form of localized compact excitations in both models [17]. At higher density, phaseseparated domains are found in the active Ising model and periodic "smectic" bands in the Vicsek model. Understanding the emergence of these three types of solutions will be the focus of the rest of the paper.

III. HYDRODYNAMIC EQUATIONS

A lot of attention has been given in the literature to hydrodynamic equations of flocking models. Two different approaches have been followed, starting from phenomenological equations [11,16,18] or deriving explicitly coarse-grained equations from a microscopic model [13,15,[19][20][21]. All these equations describe the dynamics of a conserved density field ρ( r,t) coupled to a nonconserved magnetization field, the latter being a vector m( r,t) for continuous rotational symmetries, as in the Vicsek model, or a scalar m( r,t), for discrete symmetries, as in the active Ising model.

We first introduce in Sec. III A two sets of hydrodynamic equations derived by coarse-graining microscopic models which will be discussed in this paper. Then we turn in Sec. III B to a simpler set of phenomenological hydrodynamic equations on which we will establish our general results in Sec. IV.

A. Coarse-grained hydrodynamic descriptions

We first consider the equations proposed by Bertin et al. to describe a simplified version of the Vicsek model [15], in which one solely considers binary collisions between the particles. One can then use, assuming molecular chaos, a Boltzmann equation formalism to arrive at the following hydrodynamic equations for the density field and a vectorial magnetization field [22],

∂ρ ∂t = -v 0 ∇ • m, (3) ∂ m ∂t + γ ( m • ∇) m = ν∇ 2 m - v 0 2 ∇ρ + κ 2 ∇(| m| 2 ) -κ( ∇ • m) m + (μ -ζ | m 2 |) m. ( 4
)
The mass-conservation equation (3) simply describes the advection of the density by the magnetization field. Equation (4) can be seen as a Navier-Stokes equation complemented by a Ginzburg-Landau term (μζ | m| 2 ) m, stemming from some underlying alignment mechanism and leading to the emergence of a spontaneous magnetization. Because particles are self-propelled in a given frame of reference, these equations break Galilean invariance so one can have γ = 1 and κ = 0 unlike, e.g., in the Navier-Stokes equation. In Eqs. ( 3) and ( 4), to which we refer to as "Vicsek hydrodynamic equations" hereafter, all the coefficients γ , ν, κ, μ, and ζ depend on the local density; see Ref. [15] for their exact expression.

The second set of equations, which we refer to as "Ising hydrodynamic equations" in the following, has been derived to describe the large-scale phenomenology of the active Ising model [13]. In this case, the dynamics of the density field and the scalar magnetization-corresponding to the Ising symmetry-are given by

∂ρ ∂t = D ρ -v 0 ∂ x m, ( 5
)
∂m ∂t = D m -v 0 ∂ x ρ + 2 β -1 - r ρ m -α m 3 ρ 2 , ( 6
)
where β = 1/T , α and r are positive coefficients depending on β only, and v 0 = 2Dε. The advection term v 0 ∇ • m of Eq. ( 3) is here replaced by a partial derivative v 0 ∂ x m because, in the active Ising model, the density is advected by the magnetization only in the x direction.

B. Phenomenological hydrodynamic equations with constant coefficients

Coarse-grained hydrodynamic equations derived from microscopic models have the advantage of expressing the macroscopic transport coefficients in terms of microscopic quantities (noise, self-propulsion speed, etc.). However, these possibly complicated relations may not be relevant to understand the qualitative behavior of the models. Thus, before discussing the Vicsek and Ising hydrodynamic equations in Sec. V, we first study in detail, in Sec. IV, a simpler model,

∂ t ρ = -v 0 ∇ • m, ( 7
)
∂ t m + ξ ( m • ∇) m = D∇ 2 m -λ ∇ρ + a 2 m -a 4 | m| 2 m, ( 8
)
where the transport coefficients v 0 , ξ , D, λ, and a 4 are constant. In the following, we refer to these equations as the "phenomenological hydrodynamic equations". This simplified model is very similar to that first introduced by Toner and Tu from symmetry considerations [11]. However, unlike the original Toner and Tu model, we keep an explicit density dependence in a 2 : a 2 (ρ) = ρϕ g , which is essential to account for inhomogeneous profiles [13,15,18]. The stability criteria of the homogeneous solutions [ρ(r,t) = ρ 0 , m(r,t) = m 0 ] of Eqs. (7) and ( 8) are readily computed:

(a) For ρ 0 < ϕ g (a 2 (ρ 0 ) < 0) only the disordered solution (ρ 0 , | m 0 | = 0) exists and is stable.

(b) For ρ 0 > ϕ g (a 2 (ρ 0 ) > 0) the disordered solution becomes unstable and ordered solutions [ρ 0 , | m 0 | = (ρ 0ϕ g )/a 4 ] appear.

(c) The ordered solutions are linearly stable only when ρ 0 > ϕ = ϕ g + 1 4a 4 v 0 +2λ . Thus, in the range ρ 0 ∈ [ϕ g ,ϕ ], homogeneous solutions are linearly unstable. In the language of the liquid-gas transition, ϕ g and ϕ are the gas and liquid spinodals, between which the homogeneous phases are linearly unstable and spinodal decomposition takes place. In the next section we address the existence of heterogenous ordered excitations propagating through stable disordered (gaseous) backgrounds. This analysis will make it possible both to identify all the possible flocking patterns and to further understand the first-order nature of the flocking transition.

IV. PROPAGATIVE SOLUTIONS

Let us now establish and classify all the inhomogeneous propagating solutions of Eqs. (7) and (8). In order to do so, we first recast this problem into a dynamical system framework in Sec. IV A. We then show in Sec. IV B that three types of propagating solutions exist with different celerities c and densities of the gaseous background ρ g . Sections IV C, IV D, and IV E are dedicated to a detailed study of how these solutions depend on c and ρ g . Section IV F shows how, once the average density is fixed, we are left with a one-parameter family of solutions. Last, Sec. IV G is devoted to cases where the inhomogeneous profiles can be studied analytically.

A. Newton mapping

Following Ref. [15], we look for inhomogeneous polar excitations invariant along, say, the y direction and which propagate and/or relax solely along the x direction. We thus assume m y = 0 and reduce Eqs. (7) and (8) to:

∂ t ρ = -v 0 ∂ x m, ( 9
)
∂ t m + ξm∂ x m = D∂ 2 x m -λ∂ x ρ + (ρ -ϕ g )m -a 4 m 3 , ( 10
) where we wrote m = m x to ease the notation. We look for solutions propagating steadily at a speed c. Introducing the position z = xct in the frame moving at c: ρ(x,t) = ρ(z), m(x,t) = m(z), we obtain

c ρ -v 0 ṁ = 0, ( 11
)
D m + (c -ξm) ṁ -λ ρ + (ρ -ϕ g )m -a 4 m 3 = 0, ( 12
)
where the dots denote derivation with respect to z. Solving Eq. (11) gives ρ(z) = ρ g + v 0 c m(z). If ρ(z) is localized in space, then ρ g has a simple meaning. Since ρ(z) = ρ g when 7) and (8). Center: Magnetization profile in the comoving frame z = xct or, equivalently, trajectory m(z) of a point particle in the spurious time z. Right: Phase portrait corresponding to the trajectory m(z). Parameters:

D = v 0 = λ = ξ = a 4 = ϕ g = 1.
m(z) = 0, the integration constant ρ g is the density in the gaseous phase surrounding the localized polar excitation. We can then insert the expression of ρ in Eq. ( 12) and obtain the second-order ordinary differential equation

D m + c - λv 0 c -ξm ṁ -(ϕ g -ρ g )m + v 0 c m 2 -a 4 m 3 = 0. (13)
Following Refs. [15,16], we now provide a mechanical interpretation of Eq. ( 13) through the well-known Newton mapping. Rewriting Eq. ( 13) as:

D m = -f (m) ṁ - dH dm , ( 14
)
H (m) = -(ϕ g -ρ g ) m 2 2 + v 0 3c m 3 - a 4 4 m 4 , ( 15
)
f (m) = c - λv 0 c -ξm, ( 16
)
it is clear that this equation corresponds to the mechanical equation of motion of a point particle. The position of the particle is m, z is the time variable, D is the mass of the particle, H (m) is an energy potential, and f (m) is a position-dependent friction. The trajectory m(z) of this fictive particle exactly corresponds to the shape of the propagative excitations of our hydrodynamic model in the frame moving at a speed c (see Fig. 3). We shall stress that for a given hydrodynamic model, Eq. ( 14) is parametrized by the two unknown parameters c and ρ g which a priori can take any value. Each pair (c, ρ g ) gives different potential H and friction f and hence different trajectories m(z). We now turn to the study of these trajectories and of the corresponding admissible values for the celerity c and the gas density ρ g .

B. Three possible propagating patterns

The original problem of finding all the inhomogeneous propagative solutions m(x,t), ρ(x,t) of the hydrodynamic equations is now reduced to finding all the pairs (c, ρ g ) for which the corresponding trajectories m(z) exist. Mass conservation, Eq. ( 9), imposes the boundary condition m(z = -∞) = m(z = +∞) so we are looking for solutions of ( 14) which are closed in the (m, ṁ) plane. An example of propagative solutions m(x,t), ρ(x,t) together with the corresponding trajectory m(z) and its phase portraits is shown in Fig. 3.

0 m H(m) 0 m - m + m H(m) 0 m - m + m H(m)
To put a first constraint on ρ g and c, let us rule out the potentials which cannot give such physical solutions. The extrema of H , solutions of H (m) = 0, are located at m = 0 and m = m ± with

m ± = v 0 2a 4 c 1 ± 1 - 4a 4 (ϕ g -ρ g )c 2 v 2 0 . ( 17
)
We can already discard the case where H (m) has two complex roots since H then has a single maximum at m = 0, and all trajectories wander to m = ±∞ (see Fig. 4, left). This leads to a first condition on c, ρ g :

(ϕ g -ρ g )c 2 < a 4 v 2 0 . (S1)
Without loss of generality we can assume that c > 0 and look only for solutions with m 0. This rules out the (c, ρ g ) values for which m -< 0 and m + > 0 which give oscillations between negative and positive values of m (see Fig. 4, central panel). At the hydrodynamic equation level, such solutions would indeed correspond to different parts of the profiles moving in opposite direction. The corresponding condition

ρ g < ϕ g (S2) imposes 0 < m -< m + .
The potential H then has two maxima, at m = 0 and m = m + , and one minimum, at m = m -. The typical shape of potential which gives admissible solutions is shown in Fig. 4 along with examples of potentials ruled out by conditions (S1) and (S2).

From the admissible shape of the potential H , we can now list all possible trajectories m(z) and the corresponding fields m(x,t), ρ(x,t):

(i) Limit cycles, whose corresponding magnetization profiles are periodic bands, as shown in the first row of Fig. 5.

(ii) Homoclinic orbits, that start infinitely close to a maximum of H , hence spending an arbitrary large time there, before crossing twice the potential well in a finite time to

f > 0 f < 0 ×10 -0.2 0.0 0.2 0.4 0.6 -4 -2 0 2 m H(m) f > 0 f < 0 ×10 -0.2 0.0 0.2 0.4 0.6 -4 -2 0 2 m H(m) f > 0 f < 0 ×10 FIG. 5. ( Color
D = v 0 = λ = ξ = a 4 = ϕ g = 1.
finally return to the same maximum of H at z = ∞. These trajectories correspond to isolated solitonic band profiles, as shown in the second row of Fig. 5.

(iii) Heteroclinic orbits that spend an arbitrary large time close to a first maximum of H , cross the potential well in a finite time, and spend an arbitrary large time close to the second maximum of H before returning to the first maximum. These trajectories correspond to phase-separated profiles. The arbitrary waiting times at the two maxima of H then reflect the arbitrary sizes of two phase-separated domains (see the third row of Fig. 5).

A third condition on ρ g ,c arises from the nonlinear friction term. Following the classical mechanics analogy, we define an energy function

E = 1 2 D ṁ2 + H. ( 18
)
Multiplying the equation of motion ( 13) by ṁ, we get

dE dz = -f (m) ṁ2 . ( 19
)
Energy is injected when f (m) < 0 and dissipated when f (m) > 0. On a closed trajectory, the friction f must thus change sign. Since f is a decreasing function of m, this imposes f (0) > 0 for trajectories with m(z) > 0 or, equivalently,

c > λv 0 . (S3)
The conditions (S1), (S2), and (S3) thus provide loose bounds on the subspace of the (c,ρ g ) plane which contains the three types of trajectories m(z) described above. These trajectories correspond to the three types of inhomogeneous profiles seen in the microscopic models [23]. Before studying the stability and coarsening of these propagative solutions, we first need to understand precisely how they are organised in the (c, ρ g ) plane. In order to do so, we first analyze the phase portrait of the dynamics (14). We then study how the phase portrait evolves when ρ g and c are varied.

C. Stability of the fixed points

The structure of the phase portrait is most easily captured by locating the fixed points of ( 14) and studying their stability. We first rewrite (14) as a system of two first-order differential equations:

d dz m ṁ = ṁ -f (m) D ṁ -H (m) D . ( 20
)
The fixed points are the solutions satisfying ṁ = 0 and H (m) = 0, i.e., the constant solutions extremizing H . As explained before, because of the condition (S2), the extrema of H at m = 0,m -,m + are such that 0 < m -< m + , so 0 and m + are two maxima and m -is a minimum of H . Linearizing the dynamics around one of the fixed points, we define m = m 0 + δm with m 0 = 0, m -,m + , so ṁ = δ m and

d dz δm δ m = 0 1 -H (m 0 )/D -f (m 0 )/D δm δ m . ( 21
)
The stability of the fixed points is given by the eigenvalues λ 1,2 of the 2 × 2 matrix which read

λ 1,2 (m 0 ) = -f (m 0 ) 2D ± f (m 0 ) 2D 2 - H (m 0 ) D . ( 22
)
(i) At the maxima m 0 = 0 and m 0 = m + of H , H (m 0 ) is negative and the two eigenvalues are thus real with opposite signs. These fixed points are saddle points with one unstable direction (λ 1 > 0) and one stable direction (λ 2 < 0).

(ii) At the minimum m 0 = m -of H , H (m 0 ) is positive and the real part of the two eigenvalues have the same sign, given by -f (m -). The fixed point is stable when f (m -) > 0 and unstable when f (m -) < 0. Physically, when the friction of the fictive particle is negative around m = m -, small perturbations are amplified, driving the trajectory away from the fixed point. Conversely, a positive friction damps any initial perturbation, leading to trajectories converging towards m -. When c and ρ g are such that f (m -) = 0, λ 1,2 are complex conjugate imaginary numbers. A Hopf bifurcation takes place, leading to the apparition of a limit cycle.

At the onset of a Hopf bifurcation, a limit cycle appears around the fixed point whose stability changes. In the following sections we elucidate how the interplay between the saddlepoint and the Hopf dynamics rules the nonlinear dynamics of the fictive particle and hence the polar-band shape.

D. Hopf bifurcation

Let us first provide a comprehensive characterization of the Hopf bifurcation. It happens when the real part of λ 1,2 vanishes, i.e., when

f (m -) = c - λv 0 c -ξm -(c,ρ g ) = 0, ( 23
)
where m -, which depends on both c and ρ g , is given by Eq. (17). Equation ( 23) is satisfied on the line

ρ H g (c) = ϕ g + (-c 2 + v 0 λ)(-a 4 c 2 + a 4 v 0 λ + v 0 ξ ) c 2 ξ 2 , ( 24
)
which we call the Hopf transition line.

Following standard text books in bifurcation theory [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF], the type of Hopf bifurcation (super-or subcritical) is given by the sign of

= ξ 2 16ωD 3 H m -,ρ H g ∂m - ∂ρ g ρ H g , ( 25
)
where ω = H (m -,ρ H g )/D > 0 is the imaginary part of the eigenvalues at the bifurcation point. Moreover,

∂m - ∂ρ g = -1 v 2 0 c 2 -4a 4 (ϕ g -ρ g ) (26)
is always negative because of condition (S1). The sign of is thus given by the sign of H (m -,ρ H g ), which changes at c = c * with Two different scenarios occur depending on whether c is larger or smaller than c * .

c * = √ v 0 (3a 4 λ + ξ ) √ 3a 4 . ( 27
)
(a) When c < c * , the Hopf bifurcation is supercritical ( > 0). The system branches from a stable fixed point for ρ g > ρ H g (case I, Fig. 6) to a stable limit cycle surrounding an unstable fixed point for ρ g < ρ H g (case II, Fig. 6). (b) When c > c * , the Hopf bifurcation is subcritical ( < 0). The system branches from an unstable fixed point when ρ g < ρ H g (case IV, Fig. 6) to an unstable limit cycle surrounding a stable fixed point when ρ g > ρ H g (case III, Fig. 6).

The organization of these four typical cases in the (c,ρ g ) plane is illustrated in Fig. 6. We thus see that, when c < c * , limit cycles exist for ρ g smaller than ρ H g , whereas when c > c * , they exist for ρ g larger than ρ H g . The Hopf bifurcation line is thus a boundary of the domain of existence of periodic propagative solutions of the hydrodynamic equations. Let us now consider what happens when we explore the c,ρ g plane further away from the Hopf bifurcation line.

E. Structure of the (c, ρ g ) solution space

So far, we have shown that three different types of trajectories m(z) (periodic, homoclinic, and heteroclinic) can be found by varying the values of c, ρ g . The subspace where these physical solutions can be found was first bounded by the conditions (S1), (S2), and (S3). In the previous section, we further found that the Hopf transition line ρ H g (c) given by 12 in the (c, ρ g ) space for c < c * . Top: Three phase portraits along the line, from the vicinity of the Hopf bifurcation line (A) until the homoclinic is found (C). Bottom-left: Size of the limit cycle m defined in Eq. ( 28). The limit cycle disappears when m is large enough that the orbit reaches m = 0, where the trajectory is homoclinic. Bottom-right: Average density ρ 0 of the solutions. Parameters:

D = v 0 = λ = ξ = a 4 = ϕ g = 1.
Eq. ( 24) is the upper boundary for the admissible values of ρ g when c < c * and the lower boundary when c > c * .

To explore the remaining (c, ρ g ) space, we numerically integrated the dynamical system (20) using a Runge-Kutta scheme of order 4. Starting from different initial conditions, one easily finds the basins of attraction of the different solutions. To locate unstable fixed points and limit cycles, we integrated the dynamics backward in time since they are attractors when z → -∞. As c and ρ g vary, so do the shapes and sizes of the limit cycles. To quantify these variations, we measured the "amplitude" of a cycle, defined as the difference between the two extrema of m(z),

m ≡ max z [m(z)] -min z [m(z)]. (28)
We systematically vary ρ g at fixed c, first focusing on the case c < c * where the Hopf bifurcation is supercritical. Decreasing ρ g , a stable limit cycle of vanishing amplitude appears at ρ g = ρ H g [Fig. 7,panel A]. The amplitude of the cycle then increases as ρ g decreases [Fig. 7,panel B] until it hits the fixed point at m = 0 where the limit cycle becomes an homoclinic trajectory [Fig. 7,panel C]. For even lower ρ g the particle escapes to m = -∞. The variation of the cycle amplitude with ρ g shown in Fig. 7 can be qualitatively explained. When ρ g decreases, the distance m -m f increases, where m f = 1 ξ (cλv 0 /c) is the value of m where the friction changes sign, i.e., f (m f ) = 0. More energy is thus injected in the system and, to dissipate this energy, the trajectory need to go closer to m = 0.

A symmetric behavior is observed when c > c * for the subcritical Hopf bifurcation. Increasing ρ g , an unstable limit cycle of vanishing amplitude appears at ρ g = ρ H g . The amplitude of the cycle then increases with ρ g until the trajectory hits the point m = m + where we have an (unstable) homoclinic solution that starts from m = m + as shown in

D = v 0 = λ = ξ = a 4 = ϕ g = 1.
Fig. 8. Physically, when increasing ρ g , m fm -increases so the friction around the stable fixed point at m = m -becomes larger and thus its basin of attraction (whose boundary is the unstable limit cycle, see Fig. 6) becomes larger.

All in all, the central results of this section is that all the admissible solutions lie in a band delimited by the Hopf bifurcation line ρ H g (c) and a line where the homoclinic trajectories are found, as shown in Fig. 9. Inside this band there exists stable nondegenerate limit cycles, corresponding to periodic propagating profiles. The unique heteroclinic trajectory is located exactly at c = c * where the Hopf bifurcation changes from supercritical to subcritical. We thus observe a twoparameter family of periodic solutions, a line of homoclinic trajectories, and a unique heteroclinic trajectory. Going back to the original pattern formation problem, they correspond to a two-parameter family of micro-phase-separated profiles, a line of isolated solitonic bands, and a unique phase-separated state where a macroscopic polar liquid domain cruises through a disordered gas.

F. Working at fixed average density

In the microscopic models and the original hydrodynamic equations the average density ρ 0 is a conserved quantity fixed by the initial condition. On the contrary, when considering the trajectories of the fictive particle m(z), ρ 0 is not a priori fixed and varies between the different solutions. To compute the mean density on a trajectory m(z) we simply average ρ(z) = ρ g + v 0 m(z)/c over time.

As shown in Fig. 7 (bottom-right), we find that at fixed c < c * , ρ 0 decreases when ρ g decreases. It ranges from

ρ 0 = ρ g + v 0 c m -when ρ g = ρ H
g to ρ 0 = ρ g at the homoclinic trajectory where the portion of the trajectory with m(z) = 0 becomes negligibly small. Note that, at the heteroclinic trajectory, ρ 0 can take a large range of values. Since the size of the gas and liquid domains are arbitrary, the average density can take any value in [ρ h g ,ρ h ], where ρ h g and ρ h are the densities in the gas and liquid domains, respectively.

Fixing ρ 0 adds a constraint that selects a line of solutions in the (c, ρ g ) space, as shown in Fig. 10 (left). For all ρ 0 ∈ [ρ h g ,ρ h ] these lines end at the heteroclinic trajectory. We also observe that, at fixed ρ 0 , the closer the trajectories are to the heteroclinic solution, the larger their amplitude (see Fig. 10,right). This means that along a line ρ 0 = cst, the higher the amplitude the faster the band excitations propagate. This point will turn out

D = v 0 = λ = ξ = a 4 = ϕ g = 1.
to be crucial when discussing the coarsening dynamics at the hydrodynamic level in Sec. VII.

Until now, we have shown that three different types of possible trajectories m(z) exist, which correspond to all the propagative solutions observed in the microscopic models of flying spins. We have further identified the subset of values of the propagation speed c and the gas density ρ g for which these solutions exist. We can now turn to the study of their dynamical stability at the hydrodynamic equation level. However, we first discuss analytically in the next section the shape of inhomogeneous solutions.

G. Exact solution for the heterocline

There are no general analytic solutions for the propagating inhomogeneous profiles. However, progress is possible for some limiting cases. In the following we show that a complete solution for the heterocline-its position in the (c, ρ g ) plane and its shape-can be determined exactly [13]. In Appendix A, we then show that, although exact solutions are not available, progress regarding the shape of the homoclinic solutions is achievable in the small D limit.

To compute the shape of the heterocline, let us start from the ansatz

m 1,2 (z) = m 2 [1 + tanh (k 1,2 (z -z 1,2 ))]. (29)
Each of m 1 (z) and m 2 (z) describe an interface centered around z = z 1,2 between a disordered phase with m = 0 and an ordered phase with m = m . The complete heteroclinic trajectory then consists of two fronts glued together: An ascending front m 1 (z) with k 1 > 0 and a descending front m 2 (z) with k 2 < 0, with z 2 z 1 (see Fig. 11); being part of the same profile, the two fronts share the same celerity c and density ρ g .

Moreover, we know that m must be located at the second maximum of H so

m = m + = v 0 2a 4 c 1 + 1 - 4a 4 (ϕ g -ρ g )c 2 v 2 0 . ( 30
)
Plugging m 1 and m 2 in Eq. ( 13) and replacing m by its expression, one obtains for each of the fronts

g(c,ρ g ,k 1,2 ) + h(c,ρ g ,k 1,2 ) tanh (k 1,2 (z -z 1,2 )) = 0, (31)
D = v 0 = λ = ξ = a 4 = ϕ g = 1.
where g and h are complicated functions that we omit for conciseness. Equation (31) can be true only if g and h vanish independently for both k 1 and k 2 .

We can express k 1 and k 2 as functions of c and ρ g by linearizing the ansatz (29

) around m = 0. When k 1,2 (z - z 1,2 ) → -∞, one has m 1,2 ∼ exp[2k 1,2 (z -z 1,
2 )] so we can identify k 1,2 with the two eigenvalues of the linear stability analysis Eq. ( 22). The ascending front is associated with the unstable direction k 1 = λ 1 /2 and the descending front with the stable direction

k 2 = λ 2 /2.
Replacing k 1,2 by their values in Eq. ( 31), we have four equations for the two unknowns c and ρ g . After some algebra, one obtains a unique solution (c h , ρ h g ) with

c h = c * = √ v 0 (3a 4 λ + ξ ) √ 3a 4 , ( 32
)
ρ h g = ϕ g - 2v 0 9a 4 λ + 3ξ . ( 33
)
This gives us the magnetization m and the the fronts steepness k 1,2 as

m = 2v 0 √ 3a 4 v 0 (3a 4 λ + ξ ) , ( 34
)
k 1 = v 0 (8a 4 D + ξ 2 ) - √ v 0 ξ 4D √ 3a 4 (3a 4 λ + ξ ) , ( 35
)
k 2 = -v 0 (8a 4 D + ξ 2 ) - √ v 0 ξ 4D √ 3a 4 (3a 4 λ + ξ ) . ( 36
)
In Fig. 11, we show that this solution matches exactly the heteroclinic orbit found by numerical integration of Eq. (20).

V. BACK TO THE VICSEK-LIKE AND THE ACTIVE ISING MODELS: NONLINEAR SOLUTIONS OF THE HYDRODYNAMIC EQUATIONS

In Secs. III and IV we consider phenomenological hydrodynamic equations and assumed the simplest possible dependencies of their coefficients with density. Here we extend our study to the more realistic hydrodynamic equations presented in Sec. III A. We first consider in Sec. V A the Vicsek hydrodynamic equations before turning to the Ising hydrody-namic equations in Sec. V B. For sake of completeness, we also consider in Appendix B a more general case where the potential H appearing in the dynamical system for m(z) is not a polynomial in m. While not directly relevant for the hydrodynamic equations studied in this paper for the Vicsek and active Ising model, such Hamiltonians cannot be ruled out and may arise, for instance, from a density dependence of the coefficient a 4 in the phenomenological hydrodynamic equations (8).

A. Vicsek-like equations

Let us consider the hydrodynamic equations ( 3) and ( 4) introduced by Bertin and coworkers to describe a simplified version of the Vicsek model [15]. These equations have the same structure as the phenomenological equations ( 7) and ( 8) with two additional gradient terms κ 2 ∇(| m| 2 ) and -κ( ∇ • m) m. Importantly, all the coefficients γ , ν, κ, μ, and ζ depend on the density.

We follow the same approach as before. Looking for propagative solutions invariant in the transverse direction we set m y = 0, write m x = m, and go to the comoving frame z = xct, to obtain

c ρ -v 0 ṁ = 0, ( 37
)
ν m + c - v 2 0 2c ṁ -γ m ṁ - v 0 2 ρ + μm -ζ m 3 = 0. ( 38
)
Note that after setting m y = 0 the two κ gradient terms of Eq. ( 4) cancel each other. As before, Eq. (37) directly yields

ρ(z) = ρ g + v 0 c m(z). ( 39
)
As we show in Appendix C, Eq. (38) can also be greatly simplified using the explicit density dependence of its coefficients. Introducing γ = γ /ν, ζ = ζ /ν, and writing μ = μ 1 ρμ 2 and ν -1 = ν 1 ρ + ν 2 , one obtains the second-order ordinary differential equation

m + (α -ξm) ṁ -a 2 m + a 3 m 2 -a 4 m 3 = 0, (40)
where the coefficients are all function of c and

ρ g α = c - v 2 0 2c (ν 1 ρ g + ν 2 ) ξ = v 3 0 2c 2 ν 1 + γ a 2 = μ 2 ν 2 + (μ 2 ν 1 -μ 1 ν 2 )ρ g -μ 1 ν 1 ρ 2 g a 3 = v 0 c (2ρ g μ 1 ν 1 + μ 1 ν 2 -μ 2 ν 1 )
a 4 = ζ - v 2 0 c 2 μ 1 ν 1 . (41)
Interestingly, Eq. (40) has exactly the same form as Eq. ( 13)the dynamical system stemming from the phenomenological hydrodynamic equations studied in Sec. IV-although with coefficients depending in a more complicated way on ρ g and c. Their propagative solutions will thus be the same, up to a slightly different organization in the c,ρ g plane. However, an important difference between the Vicsek and phenomenological hydrodynamic equations is the scaling of the magnetization with density in the ordered phase. [Eq. ( 40)]. The color code and different phases are the same as for the a 4 = cst case (Fig. 9). We obtain the same three types of inhomogeneous profiles: periodic, homoclinic, and heteroclinic. σ = 0.7.

The Vicsek hydrodynamic equations ( 3) and ( 4) of Bertin and coworkers indeed predict that the homogeneous ordered solution satisfies

| m 0 | ρ 0 = 1 ρ 0 μ ζ ---→ ρ 0 →∞ μ 1 ν 1 ζ , ( 42
)
which is consistent with what is observed in microscopic models such as the Vicsek model. On the contrary, because the coefficient a 4 in Eq. ( 8) does not depend on density, the phenomenological equations ( 7) and ( 8) would yield m 0 ρ 0 → 0 as ρ 0 increases. In this region of parameter space, which is not the main focus of this paper, the Vicsek hydrodynamic equations are thus more faithful to the phenomenology of microscopic models studied in Sec. II than the phenomenological equations ( 7) and ( 8). As we show in Appendix B, however, these phenomenological equations can recover a scaling akin to that of the Vicsek model if the coefficient a 4 of Eq. ( 8) is allowed to depend on density. This leads to a slightly more complicated dynamical system [16] that we study in the appendix for completeness.

Coming back to the dynamical system (40), following the same method as that introduced in Sec. IV, one can derive analytical expressions for the Hopf bifurcation line ρ H g (c) and the speed c * where the bifurcation becomes subcritical. We show in Fig. 12 the phase diagram in the (c, ρ g ) plane and examples for the three types of inhomogeneous trajectories. Again an exact solution for the heteroclinic trajectory can be derived and (the dashed lines in Fig. 12) is found at speed c h = c * . As expected, there is no qualitative difference with the simpler case studied in Sec. III.

B. Active Ising model equations

The active-Ising hydrodynamic equations [13] seem to differ a priori since they contain none of the nonlinear gradient terms found in the phenomenological and Vicsek hydrodynamic equations. However, these terms are in fact generated by the diffusion term in the dynamical equation (5) for the density field [13], so we will once again recover the same types of propagative solutions.

Looking for stationary profiles in the comoving frame z = xct, Eq. ( 5) and ( 6) reduce to

D ρ + c ρ -v 0 ṁ = 0, ( 43
)
D m + c ṁ -v 0 ρ + 2 β -1 - r ρ m -α m 3 ρ 2 = 0. ( 44
)
Equation ( 43) can be solved by expanding ρ in gradients of m.

Introducing the ansatz

ρ(z) = ρ g + ∞ k=0 α k d k m dz k (45)
into Eq. ( 43), and solving order by order, we get the following recursion relation:

a 0 = v 0 c a k+1 = - D c a k , (46)
from which we obtain

ρ(z) = ρ g + v 0 c ∞ k=0 - D c k d k m dz k . ( 47
)
In the following we retain only the first-order terms in gradient: 48) To simplify Eq. ( 44) we linearize around the density ϕ g = r/(β -1) where the disordered profile becomes linearly unstable. As shown in Ref. [13], this is a good approximation when β 1 because ϕ g → ∞ while ρϕ g remains finite for inhomogeneous solutions. We then obtain

ρ(z) = ρ g + v 0 c m(z) - Dv 0 c 2 ṁ(z) + D 2 v 0 c 3 m(z). (
D m + c ṁ -v 0 ρ + 2r ϕ 2 g (ρ -ϕ g )m -α m 3 ϕ 2 g = 0. ( 49
)
Finally, inserting Eq. ( 48) in the previous equation we obtain a second-order ordinary differential equation with exactly the same terms as those obtained from the phenomenological and Vicsek hydrodynamic equations, Eq. ( 13) and (40):

D m + c - v 2 0 c -ξm ṁ -a 2 m + a 3 m 2 -a 4 m 3 = 0, ( 50
) with D = D 1 + v 2 0 c 2 ξ = 4rDv 0 c 2 + v 2 0 ϕ 2 g , ( 51
)
a 2 = 2r(ϕ g -ρ g ) ϕ 2 g a 3 = 2rv 0 cρ 2 g a 4 = α ρ 2 g . ( 52
)
This equation has again the same qualitative behavior as in the phenomenological theory and the same three types of inhomogeneous solutions are found with the same organization in the (c, ρ g ) parameter space.

VI. LINEAR STABILITY OF THE PROPAGATIVE SOLUTIONS IN THE 1D HYDRODYNAMIC EQUATIONS

In Secs. IV and V, we have found and classified all the propagative solutions of different hydrodynamic equations. We have shown that in all cases three types of such solutions exist: periodic patterns of finite-size bands, solitary band solutions, and phase-separated solutions. These solutions were found as stable limit cycles, homoclinic, and heteroclinic orbits m(z) of the reduced dynamical system (13). This study does not tell us anything about the local and a fortiori global stability of these solutions as solutions m(xct) and ρ(xct) of the original hydrodynamic partial differential equations. Indeed, a stable orbit of the reduced dynamical system can very well be unstable to spatiotemporal perturbations when reexpressed as an inhomogeneous propagative solution of the hydrodynamic equations. For example, consider the fixed point at m = m -which is stable for ρ g > ρ H g , in regions I and III in Fig. 9. Without performing any calculation, we know that the corresponding homogeneous solution (ρ 0 = ρ g + v 0 c m -, m 0 = m -) is unstable in the hydrodynamic equations because it lies inside the spinodal region where no homogeneous solution is linearly stable.

In this section we study the linear (local) stability of these solutions at the hydrodynamic level. Although, it is a well-defined linear problem, determining the stability of inhomogeneous solutions of nonlinear partial differential equations cannot be done analytically even in the (rare) cases where these solutions are known analytically. A direct numerical study is possible (for an example, see, e.g., Ref. [START_REF] Chaté | [END_REF]) but is rather tedious and all the more so as we deal with 2D hydrodynamic equations. In the following, we study mostly the linear stability of the hydrodynamic equations reduced to one dimension (that of propagation), using a simple numerical procedure explained below. For an account of some fully 2D preliminary investigation, see the end of this section.

A. Numerical procedure

We investigated numerically the stability of the propagative solutions of the phenomenological hydrodynamic equations [Eqs. (7) and ( 8)] reduced to one space dimension, that of propagation. To do so, we select a solution of the corresponding classical mechanics problem (13) and use it as initial condition of the numerical integration of Eqs. (7) and (8).

The numerical integration is done using a semispectral algorithm, the linear terms being computed in Fourier space, the nonlinear ones in real space and a semi-implicit time stepping. This method, where the fields ρ and m are represented by their N first Fourier modes (with N large enough that the simulation has converged), is well suited to simulate systems with diffusion terms and periodic boundary conditions.

Preparing the initial condition of the hydrodynamic equations always brings in discretization errors. Indeed, because of the periodic boundary conditions in the hydrodynamic equations, we need to select a portion of the solution m(z) which is a multiple of the period. This is done with an error of x m t = 0 t = Ts FIG. 13. (Color online) Stability of homoclinic orbits in the hydrodynamic equations. a c is the difference in the amplitude of the solution measured between t = 0 and t = T s . The solution is declared to be linearly stable if δm(T s ), defined in Eq. ( 53) is smaller than a threshold that we take at 10

-3 . v 0 = λ = ξ = a 4 = ϕ g = 1, c = 1.125 (right).
order dx, the space discretization step used in the numerical integration of ( 7) and ( 8).

Accordingly, we observe a rapid relaxation at short times due to the discretization errors. Subsequently, we find that the original solution is either stable at long times or is quickly destabilized and converges to another solution. To analyze systematically the stability of the propagative solutions, we defined a quantitative criterion for the stability: We choose a time T s = 2000 much larger than the relaxation time of the initial perturbation (which happens in a time ∼100) but not too large to test solely the linear stability regardless of a possible long-time coarsening dynamics that could be induced by numerical noise. We then measure the amplitude of the solution | m|(t) defined in Eq. (28) as a function of time. If

δm(T s ) ≡ | m(T s ) -m(0)| (53)
is smaller than 10 -3 , then the solution is said to be stable and unstable otherwise (see Fig. 13). This protocol does not give exact answers to the question of linear stability, since in particular the small but finite initial perturbations may take the initial condition out of the basin of attraction of a (stable) solution. But the results presented below are relatively robust to changing our numerical resolution and the conditions used to decide stability, and we are thus confident that they represent well the "true" subset of linear stable solutions. The precise criterion does not matter much for the results since we find an abrupt transition from stable to unstable trajectories (see Fig. 13) which is visible on a variety of observables (norm of the fields, period, max and min values, etc.).

B. Results

Figure 14 contains one of the central results of our paper: Only a very small subset of the propagative solutions are stable at the hydrodynamic level. However, this subset includes the three possible types of trajectories: periodic bands, solitonic bands, and phase-separated profiles. The linearly stable solutions are all found in the region of the (c, ρ g ) plane close to the heteroclinic trajectory and close to the line of homoclinic orbits. Examples of stable and unstable solutions are shown in Fig. 15.

To understand why only the vicinity of the heteroclinic and homoclinic solutions is stable, one can argue that the solutions must have a large-enough amplitude to be dynamically stable. Indeed, the periodic solutions oscillate around m = m -, which lies inside the spinodal region of the hydrodynamic equations, where no homogeneous solutions is stable. Small-amplitude oscillations around this point should thus also by dynamically unstable and only large-enough amplitude cycles, found near the homoclinic line and the heteroclinic trajectories, are stable.

Note that the region where stable solutions are found has a rather rough boundary in Fig. 14. This is most probably an artifact due to the initial discretization error which is not controlled and varies from one propagative solution to another. Close to the threshold of linear instability, this can easily make a (linearly) stable trajectory unstable.

We report finally on fully 2D simulations performed on rectangular boxes of width L y = 100, with small noise added to each grid point on a 1D solution extended trivially along y. These yielded essentially no change with respect to the results presented above. While this encourages us to believe 7) and (8). The initial conditions are stable limit cycles of the dynamical system (13) for m(z). Top row: Unstable solution with c = 1.135, ρ g = 0.8351. Bottom row: Stable solutions with c = 1.14, ρ g = 0.8341. System size 340 × 100, 7) and ( 8) from three different initial conditions. Top: Two homoclinic trajectories glued together. Middle: Disordered initial condition. Bottom: Ordered initial condition. System size 500 × 100 (top and bottom), 1000 × 100 (middle). The profiles are averaged along the y direction in which the system is invariant.

dt = 0.1, dx = 0.4. v 0 = λ = ξ = a 4 = ϕ g = 1.
v 0 = λ = ξ = a 4 = ϕ g = 1, dt = 0.1, dx = 0.1.
that no unstable mode has components along y, and thus that the subset determined above corresponds to the linearly stable solutions of the full 2D hydrodynamic equations, we remain cautious. As a matter of fact, recent results obtained in the case of hydrodynamic equations for active nematics have revealed the existence of (very) long wavelength instability along the homogeneous direction of band solutions [26]. A similar investigation is left for future studies.

VII. COARSENING IN THE HYDRODYNAMIC EQUATIONS

The two-dimensional subset of propagative solutions that are linearly stable still contains an infinity of solutions, including smectic patterns, solitary bands, and phase-separated profiles. We can thus wonder whether a single solution is selected in the hydrodynamic equation starting from a random initial condition.

One can obtain some insight into this question by studying the lines of propagative solutions having a given average density ρ 0 , shown in Fig. 10. Indeed, ρ 0 is conserved in the hydrodynamic equations and any stable propagative solutions has to lie on such a line. As discussed earlier, the larger the amplitude of a solution, the faster the propagation. We can thus expect that if several traveling bands coexist, they will encounter and merge because of their different speeds. This will in turn increase the sizes of the surviving objects and hence their speed. This mechanism would naturally lead the system toward the heteroclinic solution, i.e., to the phase-separated state.

We checked this scenario by simulating the phenomenological hydrodynamic equations ( 7) and ( 8) with three different initial conditions, as shown in Fig. 16: (i) We build an initial condition made of two homoclines glued together, both solutions of the differential equation (13) for different speed c and inside the stability domain of Fig. 14.

To avoid discontinuities we interpolate smoothly between the gas densities of the two solutions using a hyperbolic tangent function. We then observe that the two travelling bands get closer until, when close enough (though not in contact), the smaller one evaporates, its mass being transferred to the second band. The final solution is hence a single larger isolated band.

(ii) Starting from a random ordered solution with constant ρ and a magnetization m fluctuating around m 0 = 0, several bands form at short times that all go in the same direction. The system then coarsens because of the speed differences between the liquid droplets until only "phase-separated domains," that all have the same speed, remain. In practice, because the band speeds can be very close, a final state with a single phaseseparated profile may not be reached within the time scales of our simulations and a precise study of this coarsening regime is beyond our numerical capacities.

(iii) Starting from a disordered initial condition with m(x) fluctuating around 0 and a density inside the spinodal region, many domains of positive and negative magnetization form. These objects then encounter and merge, yielding a rapid coalescence process. Because of the periodic boundary conditions, this process typically yields a single phase-separated state for the system sizes we considered. For larger systems, the coalescence process could result in a larger number of bands propagating in the same direction. We should then observe the same type of coarsening as in the second case discussed above.

At the level of hydrodynamic equations, the fact that larger ordered domains travel faster leads to a natural coarsening towards the phase-separated states. This coarsening relies both on a coalescence process and, when bands traveling in the same direction are close enough, on a ripening during which a smaller band evaporates and is "swallowed" by its larger neighbor.

Note that we studied the stability of propagative solutions and their coarsening process using the phenomenological hydrodynamic equations ( 7) and (8). While the precise results will depend on the set of hydrodynamic equations under study, our simulations of the Vicsek and Ising hydrodynamic equations, (3) and (4) and then ( 5) and ( 6), did not suggest any qualitative difference. Their comprehensive investigation, which is beyond the scope of this paper, would nevertheless be interesting, especially to quantifiy the role played by the nonlinearities in the vectorial hydrodynamic equations (4).

VIII. PHASE DIAGRAM OF COARSE-GRAINED HYDRODYNAMIC EQUATIONS

We can now construct the full phase diagram for the Ising and Vicsek hydrodynamic equations in the temperaturedensity ensemble (or noise-density for Vicsek), Fig. 17. It is composed of two spinodal lines, which are the limit of linear stability of the homogeneous disordered and ordered profiles, and two binodal lines outside which inhomogeneous solutions cannot be observed.

A. Spinodal lines

The spinodal lines can be determined analytically from a linear stability analysis, as was previously done for the Vicsek model [15] and the active Ising model [13]. The lower spinodal line ϕ g is simply the density at which the coefficient of the term linear in m in the hydrodynamic equations changes sign (in our phenomenological equation, when a 2 = 0). It reads for the Ising and Vicsek hydrodynamic equations

ϕ I g (T ) = r/(β -1), ( 54
)
ϕ V g (σ ) = μ 2 μ 1 = π (1 -e -σ 2 /2 ) 4(e -σ 2 /2 -2/3) , ( 55
)
where β = 1/T . The higher spinodal ϕ can be determined exactly for the Ising hydrodynamic equations

ϕ I (T ) = ϕ g v 0 α v 2 0 κ + 8D(β -1) 2 + v 2 0 κ + 8Dα(β -1) 2v 2 0 κ + 8Dα(β -1) , ( 56
)
where κ = 2 + α -2β and α = β 2 (1β/3). For the Vicsek hydrodynamic equation, the exact determination of ϕ is much more cumbersome [15]. In Fig. 17 we show the line ϕ V (σ ) computed numerically by simulating the Vicsek hydrodynamic equations at different densities in the homogeneous ordered state and looking at the growth of a small perturbation.

B. Binodal lines

The binodal lines ρ g and ρ are defined as the minimum and maximum global densities beyond which inhomogeneous propagative profiles cannot be observed in simulations of the hydrodynamic equations. As explained in Sec. IV F, at the heteroclinic trajectory the size of the liquid and gas domains are arbitrary so phase-separated solutions can have any density in the range [ρ h g ,ρ h ]. We find that, for all parameters we tested, ρ h = ρ , i.e., no other stable solution has a larger density than the liquid domain of the heteroclinic solution.

The situation is more subtle for the lower binodal ρ g . Depending on the external parameters, the line of homoclinic trajectories solution of the dynamical system is not always monotonous as a function of c. For example, in the Vicsek hydrodynamic equation it is not monotonous at low noise, as shown in Fig. 18. In this case the minimum density accessible to propagative solutions is the minimum of the homoclinic line. This solution need not be stable in the hydrodynamic equations so one should repeat the stability analysis done in Sec. VI for each value of the noise to determine exactly the binodal line. For simplicity, the lower binodal of the Vicsek hydrodynamic equation shown in Fig. 17 is the gas density read from the phase-separated profile, which is true at high noise and a good approximation of ρ g at lower noise values.

For the Vicsek hydrodynamic equations, the coexisting densities of the heteroclinic trajectory are known exactly, whereas in the Ising case they can be determined analytically only after linearizing Eq. (44) around ρ = ϕ g . The binodal lines in the phase diagram of the Ising hydrodynamic equations shown in Fig. 17 are thus determined by integrating numerically the hydrodynamic equations and measuring the density of the liquid and gas domain of a phase-separated solution, whereas we plot the analytical solution in the Vicsek case.

IX. CONCLUSION

Before summarizing our results, let us discuss how the study of the hydrodynamic equations presented in this article compares with the phenomenology of the microscopic models. The phase diagrams of the Ising and Vicsek hydrodynamic equations shown in Fig. 17 are qualitatively similar. They are also consistent with the phase diagrams of the microscopic models shown in Fig. 2. The hydrodynamic equations thus provide a picture which is consistent with the liquid-gas framework discussed in Sec. II C. For instance, the asymptote at finite noise (or temperature) can again be seen as the simplest way of forcing the system to cross a transition line to go from its gas to liquid phase.

The comparison between the phase-separated regions of the microscopic models and hydrodynamic equations is, however, more subtle. We have indeed shown, using the phenomenological hydrodynamic equations ( 7) and (8), that the coarsening leads, at the hydrodynamic level, to the phase-separated solution. This is true in the active Ising model but not in the Vicsek model, which only shows micro-phase separation and thus a coarsening leading to a periodic solution (see Fig. 1). This difference between microscopic models and hydrodynamic equations is, however, not surprising since it was recently shown that fluctuations are essential to account for the selection of micro-phase-separated profiles in the Vicsek model [14].

More precisely, phase-separated and micro-phaseseparated solutions are both linearly stable in the hydrodynamic equations for vectorial order parameters. As noise is added to these equations, though, large bands are destabilized and break in periodic array of finite-size bands, in agreement with what is seen in microscopic simuations of the Vicsek model [14]. We can now tentatively put together these results to propose a coarsening process that would lead to the microphase-separated states seen in the Vicsek model. Starting from a profile with many coexisting bands, the coarsening at fixed ρ 0 would tend to lead to the formation of larger and larger bands. This coarsening would, however, be arrested by the fact that the fluctuations in the Vicsek model set a maximal size beyond which bands are (nonlinearly) unstable. How this size is selected, however, remains to be determined.

All in all, we have shown in this paper that the hydrodynamic equations describing polar flocking models generically share the same propagative solutions. We found three types of solutions that are all observed in microscopic models: periodic orbits that describe microphase separation, homoclinic orbits describing solitonic objects, and heteroclinic trajectories that correspond to phase separation. Only a subset is stable in the hydrodynamic equations, but it still contains the three types of solutions.

The coarsening in the hydrodynamic equations favors the fastest solutions which are also the largest patterns. It thus leads naturally to the phase-separated solution which controls the phase diagram of the hydrodynamic equations. The same behavior is observed in the microscopic active Ising model, whereas one can understand the microphase separation of the Vicsek model as an arrested coarsening, the largest pattern being nonlinearly unstable to the (giant) fluctuations of the liquid phase. such that F (m) = f (m). Our dynamical system can then be recast as

d dz m Z = 1 D [Z -F (m)] -H (m) . ( A2
)
Let us consider a large-amplitude orbits that start close to m = 0. An example of such orbits and the corresponding phase portrait for the (m, Z) variables is shown in Fig. 19. When D is small, ṁ relaxes quickly to zero so Z relaxes to the parabola Z = F (m) in a time ∼D. Following the trajectory in Fig. 19, starting from point A at m ≈ 0, the trajectory between A and B is above the parabola F (m). The distance with the parabola first increases when m < m -, which implies H (m) < 0 and Ż > 0, and decreases afterwards when m > m -. The distance with the parabola stays of order D, set by the relaxation time of ṁ. When Z = f (m) at point B, m relaxes to m = 0 in a time ∼D. This gives profiles with sharp leading fronts and long exponential tails, which are indeed consistent with the profiles seen at small temperature in the Vicsek model and its putative hydrodynamic description [12,15,19].

This picture is consistent with the fact that at leading order in D the eigenvalues at point (m = 0, ṁ = 0) read

λ 1 = ϕ g -ρ g c -λv 0 /c λ 2 = - c -λv 0 /c D , ( A3
)
so we have a slow unstable direction λ 1 and a fast stable direction λ 2 . These two eigenvalues indeed control a large part of the trajectory, as shown in the right panel of Fig. 19. On the contrary we observe in the microscopic Vicsek and active Ising models that, at large densities, | m 0 |/ρ 0 reaches a constant P 0 1 (set by the noise intensity in the microscopic models). This can be achieved by assuming that a 4 = (ρP 2 0 ) -1 as in Ref. [16] so

| m 0 | ρ = P 0 ρ -ϕ g ρ ---→ ρ→∞ P 0 . (B2)
Repeating the same treatment as before, we look for 1D propagative solutions with speed c which are the solutions of

D m + c - λv 0 c -ξm ṁ -(ϕ g -ρ g )m + v 0 c m 2 - m 3 P 2 0 ρ g + v 0 c m = 0, ( B3
)
which can be written in the same form as Eq. ( 14) with a different potential

D m = -f (m) ṁ -H (m), ( B4
) 6). Bottom row: Phase portraits of the solutions in phase V for the same c = 1.2421 and ρ g = 0.81228 (left) and ρ g = 0.81248 (right). When ρ g increases, the size of the stable cycle decreases while the size of the unstable cycle increases, until they collide as indicated on the phase diagram. v 0 = λ = ξ = P 0 = ϕ g = 1.

H (m) = - c 3 ρ 2 g P 2 0 v 3 0 m + c 2 ρ g 2P 2 0 v 2 0 - (ϕ g -ρ g ) 2 m 2 - c 3P 2 0 v - v 0 3c m 3 + c 4
One observes the same three types of trajectories already shown in Fig. 5: periodic, homoclinic, and heteroclinic. The phase diagram (ρ g , c), shown in Fig. 20 is also similar to the case a 4 = cst except in a very small region close to the heteroclinic trajectory. We find again a two-parameter family of periodic orbits and a line of homoclinic solutions which terminates at a unique heteroclinic trajectory.

As before, we can compute analytically the line ρ H g (c) where the Hopf bifurcation takes place and the speed c * such that the bifurcation is supercritical for c < c * and subcritical for c > c * . The difference with the previous case is that the heteroclinic trajectory, for which we do not have an exact solution anymore, is not found at c = c * but at c h > c * . As a consequence, we observe a new phase (shown as number V in Fig. 20) in which two limit cycles are found. A large stable cycle surrounds a smaller unstable limit cycle which itself encapsulates the stable fixed point m = m -. As seen from the phase portraits in Fig. 20, when increasing ρ g the size of the stable cycle decreases, whereas the size of the unstable one increases. Thus the upper limit of existence of this phase (the magenta line in the phase diagram) is the moment where the cycles collide and annihilate.

FIG. 1 .

 1 FIG. 1. (Color online) Top: Micro-phase separation in the Vicsek model. η = 0.4, v 0 = 0.5, ρ 0 = 0.83, 1.05, 1.93. Bottom: Phase separation in the active Ising model. D = 1, ε = 0.9, β = 1.9, ρ 0 = 1.5, 2.35, 4.7. System sizes 800 × 100. High-density bands propagate as indicated by the red arrows on the left snapshots.

FIG. 2 .

 2 FIG. 2. (Color online) Phase diagrams of the microscopic models. The red (upper) and blue (lower) lines delimit the domain of existence of (micro-) phase-separated profiles. The black horizontal lines and squares indicate the position of the snapshots shown in Fig. 1. v 0 = 0.5 for the Vicsek model, D = 1, and ε = 0.9 for the active Ising model.

FIG. 3 .

 3 FIG. 3. (Color online) Left: Density and magnetization profiles of a propagative solution of the hydrodynamic Eqs. (7) and (8). Center: Magnetization profile in the comoving frame z = xct or, equivalently, trajectory m(z) of a point particle in the spurious time z. Right: Phase portrait corresponding to the trajectory m(z). Parameters: D = v 0 = λ = ξ = a 4 = ϕ g = 1.

FIG. 4 .

 4 FIG. 4. (Color online) The green potential can give rise to physical (positive, nonexploding) solutions while the red ones are ruled out by our conditions (S1) (left) and (S2) (center).

  FIG. 5. (Color online) Example of the three types of trajectories. From left to right: Magnetization and density profiles, phase portrait, and potential H . Top row: Periodic trajectory, ρ g = 0.835, c = 1.14. Center row: Homoclinic trajectory, ρ g = 0.83412, c = 1.14. Bottom row: Heteroclinic trajectory, ρ g = 0.83333, c = 1.1547. Phase portrait: Crosses indicate saddle points at m = 0 and m = m + . Squares indicate stable fixed points at m = m -. Potentials: The blue dashed line indicates where the friction changes sign. The red portion of the potential is the one visited by the trajectory. Parameters: D = v 0 = λ = ξ = a 4 = ϕ g = 1.

  FIG. 6. (Color online)The four types of phase portrait (m, ṁ) obtained in our system when changing ρ g and c. The line ρ g = ρ H g is where the Hopf bifurcation takes place. The bifurcation is supercritical for c < c * and subcritical for c > c * . The plain (respectively, open) black squares denote stable (respectively, unstable) fixed points at m = m -. The plain (respectively, dashed) black lines denote stable (respectively, unstable) limit cycles. The crosses denote the saddle points at m = 0 and m = m + . The initial condition of each trajectory is marked by a magenta disk and the direction of "time" z indicated by a magenta arrow.

FIG. 7 .

 7 FIG.7. (Color online) Line at fixed c = 1.12 in the (c, ρ g ) space for c < c * . Top: Three phase portraits along the line, from the vicinity of the Hopf bifurcation line (A) until the homoclinic is found (C). Bottom-left: Size of the limit cycle m defined in Eq. (28). The limit cycle disappears when m is large enough that the orbit reaches m = 0, where the trajectory is homoclinic. Bottom-right: Average density ρ 0 of the solutions. Parameters: D = v 0 = λ = ξ = a 4 = ϕ g = 1.

FIG. 8 .

 8 FIG. 8. (Color online) Unstable homoclinic trajectory starting from m = m + for c = 1.16. Parameters: D = v 0 = λ = ξ = a 4 = ϕ g = 1.

FIG. 9 .

 9 FIG. 9. (Color online) Space of all solutions. A two-parameter family of periodic orbits is found inside the band delimited by the Hopf bifurcation and the homoclinic trajectories. The constraints (S2) and (S3) are indicated by the black dashed lines. The constraint (S1) lies out of the range of this plot. Bottom: Zoom of the plot above around the point c = c * where the Hopf bifurcation changes from supercritical to subcritical. This is also where the unique heteroclinic trajectory is found. The roman numbers refer to Fig. 6, indicating the type of phase portrait found in each region. Parameters: D = v 0 = λ = ξ = a 4 = ϕ g = 1.

FIG. 10 .

 10 FIG. 10. (Color online) Top: Lines of solutions having a fixed average density ρ 0 in the space of all solutions. Bottom: The amplitude m of the cycles, defined in Eq. (28), along the lines ρ 0 = cst increases when c increases. Parameters: D = v 0 = λ = ξ = a 4 = ϕ g = 1.

FIG. 11 .

 11 FIG. 11. (Color online) Comparison between the exact solution for the heterocline (dashed lines) and the result from numerical integration of Eq. (20) (blue line). Parameters: D = v 0 = λ = ξ = a 4 = ϕ g = 1.

FIG. 12

 12 FIG. 12. (Color online) Top left: Phase diagram for the propagative solutions of the hydrodynamic equations of Bertin et al.[Eq. (40)]. The color code and different phases are the same as for the a 4 = cst case (Fig.9). We obtain the same three types of inhomogeneous profiles: periodic, homoclinic, and heteroclinic. σ = 0.7.

FIG. 14 .

 14 FIG. 14. (Color online) The blue region is the subset of admissible propagative solutions which are linearly stable. None of the solutions for c > c * is stable in the hydrodynamic equations. v 0 = λ = ξ = a 4 = ϕ g = 1. We use system sizes L x = 300 or more (adapted to fit the solution) with resolution dx = 0.5 and dt = 0.1.

FIG. 15

 15 FIG. 15. (Color online) Evolution of two propagative solutionsin the hydrodynamic equations(7) and(8). The initial conditions are stable limit cycles of the dynamical system(13) for m(z). Top row: Unstable solution with c = 1.135, ρ g = 0.8351. Bottom row: Stable solutions with c = 1.14, ρ g = 0.8341. System size 340 × 100, dt = 0.1, dx = 0.4. v 0 = λ = ξ = a 4 = ϕ g = 1.

FIG. 16

 16 FIG. 16. (Color online) Coarsening in the phenomenological hydrodynamic equations (7) and (8) from three different initial conditions. Top: Two homoclinic trajectories glued together. Middle: Disordered initial condition. Bottom: Ordered initial condition. System size 500 × 100 (top and bottom), 1000 × 100 (middle). The profiles are averaged along the y direction in which the system is invariant. v 0 = λ = ξ = a 4 = ϕ g = 1, dt = 0.1, dx = 0.1.

FIG

  FIG. 17. (Color online) Phase diagram of the Visek (left) and Ising (right) hydrodynamic equations. The spinodal lines ϕ g and ϕ are known exactly except ϕ in the Vicsek hydrodynamic equation which we computed numerically by looking at the stability of systems of size 50 × 50. The binodals ρ g and ρ are the coexisting densities of phase-separated (heteroclinic) solutions as explained in the text. The dashed line indicate the asymptote above which only disordered solutions exist. Insets are close-ups on the high-density regions with a logarithmic scale on the x axis. Parameters: v 0 = D = 1 = r = 1 (Ising), v 0 = 1 (Vicsek).

7 FIG. 18 .

 718 FIG.18. (Color online) Propagative solutions of the Vicsek hydrodynamic equations for σ = 0.2 and 0.7. The solution with the lower mean density ρ min is found at the minimum of the line of homoclinic solutions. At low noise, it does not coincide with the heteroclinic solution (blue squares). v 0 = 1.

1 FIG. 19

 119 FIG. 19. (Color online) Left: Phase portrait for the variables m and Z = D ṁ + F (m), where F (m) is the red parabola. Right: The corresponding trajectory. The dashed lines are the linear approximation around m = 0. The inset in semilog scale shows the exponential tail close to m = 0. Parameters: D = 0.01, v 0 = λ = ξ = a 4 = ϕ g = 1.

  In the hydrodynamic equation(8) the homogeneous ordered solutions have a magnetization | m 0 | = √ a 2 /a 4 . Thus, in the simplified case studied in Sec. III where a 2 = ρϕ g and a 4 = cst, one observes that| m 0 | ρ = ρϕ g a 4 ρ 2 ---→ ρ→∞ 0. (B1)

FIG. 20

 20 FIG. 20. (Color online) Top row:Phase diagram for the propagative solutions of Eq. (B4). The right plot is a zoom on the region c * < c < c h . Phases I to IV are the same as in the case a 4 = cst (see Fig.6). Bottom row: Phase portraits of the solutions in phase V for the same c = 1.2421 and ρ g = 0.81228 (left) and ρ g = 0.81248 (right). When ρ g increases, the size of the stable cycle decreases while the size of the unstable cycle increases, until they collide as indicated on the phase diagram. v 0 = λ = ξ = P 0 = ϕ g = 1.

ACKNOWLEDGMENTS

We thank T. Dauxois, A. Peshkov, J. Toner, and V. Vitelli for many illuminating discussions. This project was supported through Agence Nationale de la Recherche (ANR) projects BACTTERNS and MiTra.

APPENDIX A: LIMIT OF SMALL D

While we do not have an analytic solution for the homoclinic profiles, lots of insight on their physics can be gained by studying the limit of small diffusion coefficient D. This is most easily done by introducing the auxiliary variable Z = D ṁ + F (m), where

APPENDIX C: PROPAGATIVE SOLUTIONS OF THE VICSEK HYDRODYNAMIC EQUATIONS

Here we show how the study of propagative solutions in the Vicsek hydrodynamic equations leads to the same dynamical system as for the phenomenological equations, but with coefficients depending in a more complicated way on ρ g and c.

We start from the hydrodynamic equations ( 3) and (4) in the comoving frame z = xct:

The coefficients in Eq. (C2) are greatly simplified by dividing the equation by ν, thus getting

where γ = γ /ν and ζ = ζ /ν depend only on the noise σ on the alignment interaction [15] γ = 8 π

while μ and ν depend also on the density

One can solve Eq. (C1) to get ρ(z) = ρ g + v 0 c m(z) as before and use it in Eq. (C3) to obtain the second-order ordinary differential equation

where the coefficients are all function of c and ρ g ,

where μ 1,2 are defined by μ = μ 1 ρμ 2 and ν 1,2 by ν -1 = ν 1 ρ + ν 2 .