Thermodynamically consistent flocking: from discontinuous to continuous transitions

We introduce a family of lattice-gas models of flocking, whose thermodynamically consistent dynamics admits a proper equilibrium limit at vanishing self-propulsion. These models are amenable to an exact coarse-graining which allows us to study their hydrodynamic behavior analytically. We show that the equilibrium limit here belongs to the universality class of Model C, and that it generically exhibits tricritical behavior. Self-propulsion has a non-perturbative effect on the phase diagram, yielding novel phase behaviors depending on the type of aligning interactions. For aligning interaction that increase monotonically with the density, the tricritical point diverges to infinite density reproducing the standard scenario of a discontinuous flocking transition accompanied by traveling bands. In contrast, for models where the aligning interaction is non-monotonic in density, the system can exhibit either (the nonequilibrium counterpart of) an azeotropic point, associated with a continuous flocking transition, or a state with counterpropagating bands.


Introduction
Flocking is a prominent phenomenon in nonequilibrium statistical mechanics [1,2], comprising the collective motion of a large group of aligning self-propelled agents.It appears in systems ranging from bird flocks [3] to human crowds [4] and synthetic self-propelled colloids [5].Various theoretical models lead to a qualitatively similar flocking behavior, which suggests possible universality.A pioneering theoretical account of flocking is due to Vicsek et al. [6], who proposed an agent-based model almost three decades ago, followed by Toner and Tu [7] who offered a hydrodynamic description.Since then, these models and their variants have played a leading role in the study of active matter [8].
In the Vicsek model, locally aligning spins move with fixed speed and individually fluctuating orientation, which drives the dynamics out of equilibrium.Strikingly, under local alignment, this model undergoes an ordering transition even in two dimensions, which would be precluded in equilibrium by the Mermin-Wagner theorem [9,10].When increasing either density or alignment, the system transitions from a homogeneous disordered gas to a state of homogeneous collective motion with long-ranged order.Originally thought of as a continuous transition [6], it took almost another decade to establish flocking as a discontinuous transition [11]: in between the two homogeneous phases, there is a region in parameter space where dense ordered bands propagate against a dilute disordered background.Unlike equilibrium liquid-gas transitions, the critical point is at infinite density.This feature has been attributed to the fact that the dense and dilute phases have different symmetries, so that they cannot be connected via a continuous transformation [12].
A recent lattice-based model, known as the active Ising model (AIM), has led to a series of analytical predictions for the phase diagram of the flocking transition [13,14].In this simple extension of the Ising model, motile spins with discrete symmetry can only self-propel along a single axis.Importantly, this model is amenable to coarse-graining, yielding an exact hydrodynamic description [12,15,16].Although the phenomenology of AIM is clearly distinct from any equilibrium analogue, it is tempting to try and rationalize how AIM departs from an equilibrium reference model, with proper reversible dynamics, as was done for repulsive active particles [17,18].Interestingly, AIM does not reduce to equilibrium at vanishing self-propulsion, as one might expect.This inconsistency reveals that AIM is actually not compatible with the theory of stochastic thermodynamics [19], which assumes that the system is in contact with equilibrium reservoirs and driven by external forces (so that it naturally relaxes to equilibrium in the absence of drive).The same holds for other flocking models with continuous spin symmetry [6], or with topological interactions that involve nonreciprocal alignment [20,21], which also lead to irreversibility even at vanishing self-propulsion.
It is then natural to examine flocking models where the microscopic dynamics is amended to entail an equilibrium limit in the absence of self-propulsion, which we refer to as thermodynamic consistency.Several questions remain open: Would novel qualitative behavior emerge from thermodynamic consistency?To which universality class does the equilibrium limit of such flocking models belong?Inspired by AIM, we explore these issues in the context of flocking systems with a discrete symmetry.To make analytical progress, we choose lattice models which, like AIM, admit an analytically exact and tractable hydrodynamic description.
This work is part of a broader current effort in recasting various active systems in a thermodynamically consistent framework such as in phase field crystal models [22], active fluids [23,24] , and reaction-diffusion systems [25].Although our models are permanently out of thermodynamic equilibrium, we find that their behavior can be usefully analyzed in terms of the underlying free energy that controls their equilibrium limit at vanishing self-propulsion.Indeed, the free energy plays an informative role even outside the regime of linear response to the nonequilibrium drive.As noted previously [26,13,12], the standard scenario of a discontinuous flocking transition stems from a density modulated coupling of the alignment interaction.In our models, we reveal a strikingly different phenomenology between the cases of monotonic and non-monotonic density modulations.For models where ordering increases monotonically with density, the equilibrium limit exhibits a single tricritical point [16], and introducing activity is a non-perturbative effect: the tricritical point diverges to infinite density, and the flocking transition is discontinuous all across phase space.In contrast, for models where ordering is non-monotonic with density, the equilibrium limit exhibits a pair of tricritical points.At finite self-propulsion, such points can either collide into a single azeotropic point, yielding a continuous flocking transition, or instead a new collective state emerges with counterpropagating bands.
This paper is organized as follows.In Sec. 2, we set up the general definition of our thermodynamically consistent flocking models.We propose a microscopic dynamics with a generic type of aligning interactions, and derive the corresponding hydrodynamic description through systematic coarse-graining.Then, we examine in Secs. 3 and 4 the phase diagrams for monotonic and non-monotonic dependence of the aligning strength on density, respectively.Overall, our results illustrate how various choices of density-dependent alignment affect the flocking transition, yielding unprecedented phenomenology in the flocking context.Yet, some of this phenomenology echoes that of equilibrium systems, so that our thermodynamically consistent description allows one to anticipate and classify the topologies of nonequilibrium phase diagrams in terms of their equilibrium counterparts.

Thermodynamic consistency: From microscopics to hydrodynamics
In this Section, we introduce a class of lattice dynamics for active spins that obey thermodynamic consistency.We take locally aligning interactions with mesoscopic range to capture a transition to collective motion.Such a form of interactions is amenable to exact coarse-graining, which allows us to derive the corresponding hydrodynamic equations for some relevant fields.Importantly, our derivation is given in a functional form with an arbitrary dependence of the alignment strength on density.

Microscopic dynamics of aligning active spins
Flocks with discrete symmetry, like in AIM [12], are defined by two species of particles (+ and −) which hop with a bias to the right or left (respectively) on a periodic lattice of L ≫ 1 sites, see Fig. 1.We consider for simplicity a one dimensional lattice, yet our fluctuating hydrodynamics is trivially extended to higher dimensions.The dynamics does not exclude multiple occupancy, so that each site i can have any number of particles.The total number of particles in the system is N = ρ 0 L, where ρ 0 is the mean density.
The dynamical update rules consist of diffusive hops biased by self-propulsion, and spin flips + ↔ −.Differently from previous models [12,16], at zero self-propulsion we make both these rules derive from the same Hamiltonian H via a detailed balance constraint [27,19].We leave H unspecified at this stage.Self-propulsion is then added as a weak bias of hops in the direction of the particle's spin, yielding the following update rules (see Fig. 1): (i) Site hoping: Any particle jumps to a neighboring site with rate D 0 e the jump is in the direction of its spin, and D 0 e − β∆H 2 otherwise.Here, ∆H is the energy difference between the configurations of the system before and after the jump.
Here, D 0 is the bare diffusion constant, λ is the self-propulsion strength, and γ sets the tumble rate.The scaling of the rates with L ensures that in the hydrodynamic limit (L → ∞ at fixed ρ 0 ) all processes occur on diffusive time scales [15,28].Indeed, the time it takes for particles to traverse a macroscopic system of size L, either via diffusive motion or on account of self-propulsion, scales as L 2 which is also the time scale for tumbling events.This L-dependent scaling of the rates is standard for other lattice-based diffusive models [15,28].The main difference with AIM [13,12,26] is that our site hopping now accounts for the energy difference ∆H.This feature renders our models reversible at vanishing self-propulsion, thus setting a proper equilibrium limit, as a consequence of thermodynamic consistency [27,19].
We now restrict H such that the interaction range spans a mesoscopic scale ∆x = L δ , with 0 < δ < 1.This scale is sub-extensive in system size, yet contains a large number of particles 1 ≪ ∆x ≪ L, and is the scale appearing in the coarsegraining procedure detailed below.Within this choice of scaling each particle interacts with a vanishing fraction of the system of the order of ∼ ∆x/L ≪ 1, yet it still enables a phase transition in one dimension even in the equilibrium limit.Indeed, in defiance of the usual Landau argument (which states that having only a finite energetic cost of domain walls in one dimension destroys long range order [29]) let's consider the free energy cost of the interface between phases.The configurational entropy of a single domain wall scale as log L (as there are L sites to center it upon them).Now since each particle across the domain wall interacts with about ∼ ∆x other particles, then the energtic cost of a domain wall would scale at least as ∆x ∼ L δ , which is much larger then the entropic log L contribution (compare this with a local interaction range where the energetic cost is O (1)).Thus, the large O L δ free energy cost will make domain walls thermodynamically unfavorable and keep the ordered state stable.
Note that it was recently found [30] that the flocking transition does not survive in the one-dimensional AIM.Instead it is replaced by either direction-reversing ordered aggregates or static asters composed of opposing clusters of left and right movers that block each other.These arise by virtue of the purely local interaction range of the AIM and do not appear in our model as a result of our choice of a mesoscopic interaction range.
In practice, interactions are given in terms of some local average occupancies for the density and 'magnetization' variables, respectively η ρ i = η + i + η − i and η m i = η + i − η − i , where η + i and η − i are the number of + and − particles at site i.The local averages around site i are given by ρi We then choose the aligning Hamiltonian to take the form For f (ρ i ) = 1, this is the sum of L fully connected Ising Hamiltonians within some restricted range ∆x, where the coupling is scaled with the total number of interacting spins, see also [12].Indeed, denote by s k i ∈ {+1, −1} with k = 1, . . ., η ρ i the spin variable of the k'th particle on site i.Then the Ising Hamiltonian with all-to-all interactions within the range ∆x of site i reads The factor 1/2 avoids double counting.The normalization 2∆xρ i ensures the interactions remain extensive with the total number of interacting spins, as is common for the fully-connected Ising Hamiltonian.Then summing the right hand side (3) over all disjoint compartments of size 2∆x is equivalent to the sum (2) with f (ρ i ) = 1 which runs over all lattice points (up to sub-leading corrections in L).Then the density-dependent modulation f (ρ i ), which plays a key role in the following, tunes the alignment strength.For AIM, a similar modulation appears at a coarse-grained level, as an effective renormalization due to the hydrodynamic fluctuations [26].In our model, since the microscopic dynamics already contains a mesoscale interaction range, we choose to directly introduce such a modulation at this level.Various forms of f can be motivated on account of different types of close neighbor interactions [31,32,33,34].
In the following, we consider the hydrodynamic limit of this model.As a preliminary step, note that for a single isolated particle, there are two important length scales: ξ p = λL/γ is the average displacement of particles between tumbles due to selfpropulsion, and ξ d = L D 0 /γ is the typical diffusive displacement, over the same time.These distances are measured in units of the lattice spacing and they are both O(L), because of the L-dependent hopping rates of the model.

Hydrodynamic equations for density and magnetization
We follow a standard coarse-graining procedure, inspired by similar lattice models with diffusive scaling [35], which has already been deployed in some active lattice gas models [36,37].We start with a spatial coarse-graining over the diffusive displacement ξ d .I.e., we define the macroscopic spatial coordinate x = i/ξ d ∈ [0, ℓ s ] where ℓ s = L/ξ d = γ/D 0 is the system size in the hydrodynamic representation (measured in units of ξ d , as in [15]).The reason for this choice is that we will later analyse phase-separated profiles: their interfacial widths are of order 1 in this hydrodynamic representation and it will be convenient to take ℓ s ≫ 1 so that the system is large compared to this interfacial width.For consistency with these choices we make a diffusive scaling of time, setting t = γ t/L 2 with t the time variable of the microscopic system.
The local density in the vicinity of point x is defined by the mesoscopic coarsegraining (1) as: Since ρ is conserved and m is not, the dynamics of these fields must take the general form where J ρ , J m are the conservative fluxes and K is the net rate of change of the ρ − density due to spin flipping.As detailed in Appendix A, to leading order in system size we find the following expressions The Péclet number Pe = ξ p /ξ d is the standard measure of activity in similar active systems with diffusive scaling [38].Here, is the free-energy functional of the equilibrium (λ = 0) system with the corresponding free energy density given by the function while C and M are mobility coefficients given by The expressions (6-7) are exact as long as Pe ̸ = 0. Yet in order to describe correctly the equilibrium (Pe = 0) hydrodynamics for non convex F, one has to retain surface tension terms that are sub-leading . Still, as is usual in the thermodynamics of phase separation, the bulk term F suffices for determining the binodal (coexisting) densities ‡.
The dynamics (5-6) has a thermodynamically consistent structure.Indeed, for Pe = 0, all terms derive from the free energy F, which then serves as a Lyapunov function, as discussed in Appendix A. As shown in what follows, this equilibrium limit belongs to the universality class of Model C [16,39].Note that, as written, the dynamics (5-6) neglects some sub-leading noise terms of order L −1/2 ; these can be readily deduced from the detailed-balance condition (and are included by definition in Model C).
The general structure of (5-6) should hold rather broadly for all thermodynamically consistent flocking models with discrete symmetry and under diffusive scaling, including models with purely local interaction range (albeit with a renormalized free energy depending on dimensionality).The conservative fluxes in (6) are linear in the free-energy derivatives, as is common when expanding close to equilibrium.That these expressions hold also far away from equilibrium is a consequence of the L-dependent scaling of the microscopic rates, and the diffusive scaling of space and time; any higher order spatial gradients will be sub-leading in system size L −1 .This reasoning does not apply for the flipping term K, since it does not involve spatial variations of the fields.Indeed, this term goes beyond linear order in the free energy derivatives even at the hydrodynamic level [40].Activity (Pe > 0) leads to the extra terms Pe m and Pe ρ in the conservative fluxes (6) but otherwise does not interfere with the remaining bare equilibrium terms.Again, this is a consequence of the diffusive scaling, where self-propulsion enters as a weak bias in the microscopic dynamics.
Our thermodynamically consistent dynamics in (5)(6) is easily related to the nonconsistent one of AIM [13,12,26].Our flipping term K is the same as in the hydrodynamic equations of AIM, under appropriate matching of f .The conservative ‡ The profiles presented in the following correspond to finite Pe ̸ = 0 and were computed without the need for the sub-leading interfacial tension terms.
fluxes J ρ and J m are different from those of AIM, since they stem from a different choice of hopping rules.Nevertheless, one can reproduce the corresponding terms in AIM by sending the temperature that enters these terms to infinity (namely, taking β = 0 in the free energy that enters the conservative fluxes J ρ and J m ), as discussed in Appendix A. Indeed, denote by F 0 the free energy density evaluated at β = 0 and write F 0 for the corresponding free energy functional.Then the hydrodynamics of the AIM are also given by (5,6), except that one replaces F → F 0 in the expressions for J while retaining the full F in K. To this extent, the source of residual irreversibly in the zero self-propulsion limit of AIM can be regarded as a mismatch between two temperatures controlling respectively the conservative fluxes (J ρ , J m ) and the flipping term K.

Phase diagram for monotonic alignment strength
Various forms of the density-dependent alignment strength f (ρ), entering in the microscopic alignment (2) and in the hydrodynamic free energy (7), can be motivated to reflect different types of close neighbor interactions [31,32,33,34].We first explore the phase behavior for a monotonic f (ρ) by considering the form appearing previously in AIM [13,12,26]: This choice is convenient mathematically, although on physical grounds one might instead favor choices that are positive everywhere, and do not diverge at ρ → 0. Yet this does not change the phenomenology reported below.Other monotonic forms have been recently proposed in extensions of AIM yielding a similar phase behavior [16].In some of these models, the alignment term is not purely quadratic in the magnetization as assumed here (see (2)).Yet, its expansion close to criticality reduces to this form, and the topology of the phase diagram actually follows from the behavior close to criticality.The resulting phase diagram, given in terms of T = 1/β and ρ 0 , is composed of three phases shown in Fig. 2: (i) a homogeneous disordered phase (H; m = 0), (ii) a homogeneous ordered phase (H; m ̸ = 0), which becomes a collective motion (C.M.) for non-vanishing self-propulsion, and (iii) a phase-separated state (P.S.) between different discrete symmetries (m = 0 and m ̸ = 0), which becomes a traveling band (T.B.) for nonvanishing self-propulsion.In this Section, we study the topology of this phase diagram and derive its phase boundaries.We address first the equilibrium limit (Pe = 0), and then the nonequilibrium regime (Pe > 0).

The equilibrium limit (Pe = 0): Model C dynamics
As pointed out in [16], the equilibrium, zero-propulsion limit of thermodynamically consistent flocks with discrete symmetry must generically lie within the universality class of Model C [39].In our model, this should remain true for any non-trivial choice of f (ρ), while allowing different global phase behaviors even at equilibrium.Model C describes an equilibrium dynamics which couples a conserved and a non-conserved scalar   (10) with r = 1.Curves mark the different branches of solutions for the magnetization of homogeneous states (14).Dashed blue and red lines mark instabilities of the two branches.Both originate at the density value ρ c given by the curve (13).The instability of the magnetized branch (dashed red) terminates at the density value given by the density spinodal (15).The region of densities enclosed between these two values is linearly unstable for both branches.Red and blue dots mark the binodal densities of the phase-separated state (P.S.), which is the globally stable state for any mean density between these points (including densities for which all homogeneous states are unstable).The black dot marks a critical point.
order parameter, and the canonical Landau form close to criticality is given by [41,42] where higher order terms and terms purely linear in the density are omitted.Expanding our free energy in (7) around a homogeneous neutral state (ρ = ρ 0 , m = 0), one can derive the expressions of the different coefficients in (11), which are determined by the choice of f (ρ) §.Importantly, the coupling coefficient e in ( 11) reads which vanishes for a constant f (ρ); in this special case, our model falls outside the Model C universality class.(As discussed in Sec. 4, this also makes the nonequilibrium Pe > 0 regime different from the usual flocking scenario.)As detailed in [42], the topology of the phase diagram is set by the values of the coefficients in (11).As we have the exact form of the free energy at our disposal (7), we follow the procedure presented in [42] to map the phase diagram across the whole phase space, even beyond criticality.We find the stable phases by minimizing the free energy (7) under the constraint of conserved density (1/ℓ s ) ℓs 0 ρdx = ρ 0 .These belong to one of three possibilities: the homogeneous neutral state (ρ = ρ 0 , m = 0) (marked in Fig. 2 by H.m = 0), the homogeneous magnetized state (ρ = ρ 0 , m = m 0 ̸ = 0) (marked by H.m ̸ = 0) or a phase separated state (P.S.) with spatially varying magnetization and density.In all of the stable states, since magnetization is not conserved, its value is locked to the local value of the density via the minimization of F with respect to the magnetization.A non-trivial minimizer m 0 ̸ = 0 then emerges at ∂ mm F| m=0 = 0 which defines the curve Generally, the relation (13) only establishes linear instability of the neutral state (m = 0).In regions that also admit phase separation, as discussed below, this curve will then mark a spinodal instability which we will term the magnetization spinodal, see dashed blue curve in Fig. 2(a).Only in the absence of phase separation, this curve marks a second-order transition between disordered and ordered states, see solid black curve in Fig. 2(a).In the ordered phase, the value m 0 (ρ, T ) is found from ∂ m F = 0, which gives the magnetization in an implicit form: where the second equality defines the dependence of z on ρ and T .A state with magnetization ( 14) minimizes the free energy only for homogeneous profiles.To find phase separation, one has to follow the usual procedure of a common-tangent construction over the single-variable function F [ρ, m 0 (ρ, T )] [42].Non-convexity of this function marks linear instability against spatial density variations, and is equivalent to non-convexity of the two-variable free energy F: where |Hess (F)| is the Hessian determinant of F (ρ, m).An explicit expression for (15) is presented in Appendix B.
The relation (15) defines a spinodal curve that is complementary to the magnetization spinodal (13), see dashed red line Fig. 2(a), which we will term the density spinodal (since instability here also involves density modulations).The two spinodal curves meet at a tricritical point, as detailed in Appendix B, whose location obeys for the density ρ tri The region enclosed between the two spinodals ( 13) and ( 15) defines the region of spinodal decomposition.This is shown by the bifurcation diagram Fig. 3 which marks instabilities along a constant temperature cross section of Fig. 2(a).At T > T tri (Fig. 3 (a)) the neutral branch m = 0 loses stability at the density given by ( 13) while the entire branch of magnetized states m 0 ̸ = 0, prescribed by the nontrivial solution to ( 14), is stable.In contrast, at T < T tri (Fig. 3 (b)), the magnetized branch has an unstable segment which terminates at the spinodal density (15).Thus, any homogeneous state with density within the region enclosed between the two spinodals ( 13) and ( 15) is linearly unstable.Any monotonic f (ρ) that vanishes at a finite density while saturating at infinity will have at least one root of ( 16), yielding a phase diagram like that in Fig. 2(a).At T > T c (ρ tri ), the curve (13) marks a line of critical points separating the homogeneous disordered (m = 0) and ordered (m ̸ = 0) phases.At T < T c (ρ tri ) the transition to the ordered phase becomes discontinuous with a region of coexistence between an ordered liquid and a disordered gas.Within the miscibility gap, the relation ( 13) marks the magnetization spinodal, and the complementary density spinodal given by ( 15) lies within the symmetry-broken phase.The location of the binodal curves follows from a common tangent construction, as detailed in Appendix B. Note that, although coexisting bulk phases are found within this procedure, properly resolving the interfaces would necessitate retaining sub-leading gradient terms, which quantify surface tension in the free energy.Note also that Model C can exhibit additional phase diagram topologies, such as critical points and double critical end points [42], under appropriate tuning of f (ρ).These are absent for the curves (10), but could arise for other sufficiently elaborate choices, including strictly monotonic curves.

The nonequilibrium regime (Pe > 0): Discontinuous flocking transition
The phase diagram for the monotonic f (ρ) obeying (10) at finite Pe = 0.5 is shown in Fig. 2(b).The tricritical point shifts to infinite density, turning the second-order transition line into a first order coexistence region (miscibility gap) along its entire length.The same phase behavior is observed for AIM and its recent variants [13,16].Within the miscibility gap, the phase separated state becomes a traveling band (T.B.), as shown in Fig. 2(c).Indeed, at non-vanishing self-propulsion, any magnetized bulk phase attains a non-vanishing density flux given by Pe m, according to (6).Then, simple flux balance across the interface implies that it must propagate with velocity given by where ∆ρ and ∆m are the density and magnetization difference between phases.Moreover, we have found numerically that the propagation velocity is bounded below as |V | ≥ Pe, for all models that we have examined.
In the absence of a free energy minimization principle, there is no simple recipe to establish the binodals in the nonequilibrium case.
Throughout this paper, the nonequilibrium phase diagrams were found by numerically integrating the hydrodynamics (5-6) using a finite difference scheme for space and time.Then binodals are identified as the plateaus of the traveling density profiles as shown in Fig. 2(c).For more details see Appendix D.
To establish more broadly the generic topology of the phase diagram in Fig. 2(b), we analyze the linear instability of the dynamics.First, the magnetization spinodal line (13) remains unchanged even at finite Pe > 0, since it describes instabilities against homogeneous variations.As the active terms in the dynamics (5-6) only couple to spatial flux derivatives, they do not affect this instability.In contrast, the complementary density spinodal (15) is changed.As we show in Appendix C the instability in ( 15) is now modified into the following criterion where the mobility M is given in (8), and the magnetization m is to be evaluated at m = m 0 (ρ, T ) in (14).
Based on other flocking models, it was argued that the flocking transition should generically be discontinuous everywhere [2].In our model, this holds true if the conditions in ( 13) and ( 18) predict a finite spinodal gap at any Pe > 0. It is useful to notice that, at the magnetized state m = m 0 (ρ, T ) in ( 14), we get ∂ ρm F/∂ mm F = −∂m 0 /∂ρ, so that the density spinodal (18) becomes Moreover, close to the magnetization spinodal (13), we have Given that (20) diverges close to the magnetization spinodal (13), a finite neighborhood of the magnetization spinodal remains unstable at arbitrary small Pe no matter where one is on the phase diagram.Therefore, the non-perturbative effect of activity on the phase diagram can be attributed to the singularity of the order disorder transition of the equilibrium limit, as encoded in m 0 .Remarkably, however, this argument fails for any non-monotonic f (ρ) at its maximum f ′ (ρ 0 ) = 0.As discussed next in Sec. 4, the case of non-monotonic f actually leads to novel phase behaviors and a diagram very different from Fig. 2(b).

Phase diagram for non-monotonic alignment strength
Non-monotonic alignment strength functions f (ρ) are likely to arise physically in systems where alignment is restricted by visual obstruction, as was recently proposed to model natural bird flocking [43].Non-monotonicity opens the door to novel phase behaviors, not encountered for the monotonic case in Sec. 3 above.In equilibrium (Pe = 0), when f (ρ) features a single maximum, one typically finds two roots for ( 16) marking a pair of tricritical points.This scenario leads to two distinct miscibility gaps, as shown in Fig. 4(a).The miscibility gap to the left of the maximum has a disordered low-density phase coexisting with an ordered higher-density phase.The miscibility gap to the right of the maximum has a disordered high-density phase coexisting with an ordered lowerdensity phase.
In the rest of this Section, we consider how self-propulsion (Pe > 0) affects this equilibrium picture.The regime of non-monotonic f is more elaborate compared with the monotonic case in Sec. 3. In particular, we distinguish two qualitatively different phase behaviors, depending on the curvature of f (ρ) close to its maximum.For illustration purposes, we use particular examples in a specific family of non monotoniccurves: with α > 0. In what follows, we show that the topology of the phase diagram is controlled by local properties of f close to the maximum, so that the exponential cutoff at large densities is immaterial for our purposes.

Azeotropic point in the nonequilibrium case (Pe > 0)
For mildly curved functions f (ρ), we find numerically that the tricritical points collide at the maximum of f at any finite Pe, see Fig. 4(b).As for the case of monotonic f (ρ) in Sec. 3, the effect of self-propulsion on the phase diagram is non-perturbative: the equilibrium critical line expands almost everywhere into a finite miscibility gap, except at a single point which remains critical.Thus, in contrast with the case of monotonic f (ρ), there does remain a critical point on the phase diagram at finite Pe.Therefore, a non-monotonic f (ρ) can allow for a continuous transition from homogeneous disorder (H., m = 0) to collective motion (C.M.).
For equilibrium systems showing a phase-diagram topology equivalent to Fig. 4(b), the critical point is called an azeotropic point [44], so that we also adopt the same term here.This point marks the pairwise meeting of two branches of first order lines, with miscibility gaps that display distinct behavior.The left gap features the usual traveling bands (T.B.) as in 2(c).The right gap features reversed bands (R.B.), shown in Fig. 4(c).These bands propagate in the direction opposite to the bulk magnetization.Indeed, in the corresponding miscibility gap, the coexistence is between a magnetized lower-density phase and a paramagnetic high-density phase.Then, according to the flux balance condition (17), it follows that the phase boundaries propagate in the reverse direction from usual.
Our scenario shows that the flocking transition need not necessarily be discontinuous, at variance with a widely held view among those studying flocking [2], although it remains an exceptional case.In order to rule out some continuous transition scenarios, note that the disordered (H., m = 0) and the ordered phases (H., m ̸ = 0 at Pe = 0 or C.M. at Pe ̸ = 0) entail different symmetries.Thus, it is not possible to pass between these pure phases without crossing a (symmetry-breaking) phase transition of some kind.This argument precludes the standard scenario of liquid-liquid phase separation, where a single pair of binodals terminate at a critical point.Indeed, one sees from the equilibrium phase diagrams in Figs.2(a) and 4(a) that the binodals terminate at tricritical points, and any path between the pure phases must cross a phase transition line.In the nonequilibrium phase diagrams of Figs.2(b) and 4(b), one observes discontinuous transitions as well as an azeotropic point, all of which are consistent with the relevant symmetry principles.

5.
Phase diagrams for (a) the equilibrium case (Pe = 0) and (b) the nonequilibrium regime (Pe = 1), with non-monotonic f given by ( 21) with α = 1/2.colored dashed and solid curves are spinodals and binodals, respectively.The black solid curve is a critical line.The black dots mark tricritical points.The black dashed line marks the maximal magnetization curve (24).A C.B. phase emerges in the nonequilibrium regime.(c) The counterpropagating bands profiles for Pe = 1, ρ 0 ≃ 1.8, and T ≃ 0.17.

Counterpropagating bands in the nonequilibrium case (Pe > 0)
When f (ρ) is sufficiently curved, two tricritical points in the equilibrium limit are closer to the maximum of f (ρ), see Fig. 5(a).Here, we observe a different behavior at finite Pe compared with Sec.4.1.Instead of the azeotropic point reported in Fig. 4(b), there now appears a region with counterpropagating bands (C.B.), as shown in Fig. 5(b).This happens when the two magnetised binodals (red solid lines) merge, at a finite temperature T C.B. .For temperatures between T C.B. and T max = max [f (ρ)] , the system enters the C.B. phase (see Fig. 5(c)), a strongly fluctuating dynamical steady state whose properties can be understood in terms of the partial restoration of a broken symmetry, as we discuss next.
Recall that within the symmetry-broken T.B. and R.B. phases, the steady state consists of a single magnetized band, which travels through the system.The magnetization may be either positive or negative: this spontaneous breaking of spinreversal symmetry is familiar from the Ising model.In turn, this symmetry breaking drives a particle current, which causes the band to move either left or right.This additionally breaks the spatial reflection symmetry of the model.Due to the periodic boundary conditions of the system, the resulting steady state is also time-periodic, with period ℓ s /V .
In the C.B. phase, both left-and right-moving bands are present at the same time, as shown in Fig. 6.Since they propagate in opposite directions, the two bands regularly meet each other.It turns out that they interact weakly (we recall that the number of particles on each site is unbounded, in contrast with exclusion processes).The resulting situation is most simply represented in terms of the density of + and − particles As shown in Figs.7(a-b), the space-time dependence of these profiles is wellapproximated by ρ + ≃ ρ (x − V t) and ρ − ≃ ρ (−x − V t).In this case, the density ρ(x, t) is still time-periodic, with behaviour invariant under spatial reflection, so that it is analogous to a standing-wave solution, while the T.B. and R.B. are analogous to travelling-wave solutions.These analogues of standing waves and travelling waves are strongly anharmonic, as expected since the governing equations are non-linear.In fact, ρ is close to a top-hat profile, as it was in the T.B./R.B. phases: there are extended bulk regions at densities (ρ 1 , ρ 2 ), separated by interfaces whose widths are O (1) ≪ ℓ s .Plotting the density profile at generic times, as shown in Fig. 6, one typically finds up to four distinct regions whose densities are the following.Outside of either band the density is that of the dilute binodal ρ g = 2ρ 1 .Inside both bands, the density is that of the dense binodal ρ l = 2ρ 2 .Finally, inside one band and outside the other, the density (a-b) The ± counterpropagating bands profiles are perfectly antisymmetric: ρ + ≃ ρ (x − V t) and ρ − ≃ ρ (−x − V t) coincide when reflected and superimposed.(c) The periodic evolution of the phase volumes of the intermediate and liquid phases ν i and ν l respectively, for parameter values in Fig. 6.The complementary gas phase volume (not shown) is simply is that of the intermediate binodal Hence the intermediate binodal is always halfway between the extreme binodals, while their difference equals twice the magnetization m 0 : An important consequence of ( 23) is that the propagation velocity (17) attains its limiting minimal value |V | = Pe.The magnetization in (23) marks the bulk magnetization of the intermediate binodal m 0 = 0 (ρ i T as given explicitly by (14).Therefore, both binodals ρ g,l in ( 23) are fully determined by the value of the intermediate binodal ρ i .In practice, we find empirically that the binodal ρ i coincides with the maximal magnetization curve ρ max , defined by ∂ ρ m 0 (ρ max , T ) = 0, as shown in dashed black line in Fig. 5(b).We get an implicit representation of ρ max (T ) by differentiating ( 14) with respect to ρ, yielding Overall, Equations (23-24) specify all three binodals of the C.B. phase in closed analytic form.
We now address the phase volumes of the bulk phases, defined as the fraction of space that each phase occupies, and the corresponding highly dynamical nature of the C.B. phase state.Since the profiles of ρ + and ρ − are symmetric, each one of them carries exactly half of the total mass Lρ 0 /2.Besides, given that ± bands counter-propagate at constant speed, their overlap region varies periodically with time.We report in Fig. 7(c) the corresponding dynamics of the phase volumes for a mean density below the intermediate binodal (ρ 0 < ρ i ).The maximal phase volume of the intermediate phase ν max i is set by the usual lever rule between the intermediate and gas binodals: Then, the maximal phase volume for the liquid phase occurs at perfect overlap of the ± profiles and is exactly ν max i /2, as shown in Fig. 7(c).Interestingly, one can formulate a necessary condition for the emergence of the C.B. phase state.Indeed, the liquid and gas binodals must lie outside the spinodals ρ s g and ρ s l : The condition (25) must hold across the C.B. phase.In particular, it should hold in the vicinity of the maximum of f , where the C.B. phase first emerges.Approximating f (ρ) by a parabola in this region, the spinodal densities ρ s g,l defined in (13) are symmetrically positioned around the maximum ρ * : Performing a similar expansion for the maximal magnetization curve in (24) yields Substituting the expression ( 27) in ( 23), and recalling that ρ i = ρ max in the C.B. phase, we find that the liquid and gas binodals are also symmetrically positioned around the maximum: Finally, the condition (25) for existence of the C.B. phase is equivalent to the ordering ∆ b > ∆ s of the binodal and spinodal gaps.Using the expressions for these gaps, respectively in ( 26) and (28), we arrive at the following condition for the curvature near the maximum: Remarkably, the condition ( 29) is independent of the Péclet number, which corroborates that the C.B. phase appears as a singular perturbation of the underlying equilibrium phase diagram at any finite activity, no matter how small.

Discussion
Our results show that the density-dependence of the aligning interactions can play an unexpectedly major role in flocking models, yielding a richer flocking phenomenology than any previously reported.These results are found by choosing a thermodynamically consistent approach which entails a proper equilibrium limit, pertaining to the universality class of Model C [39,16], at vanishing self-propulsion.Such rich phenomenology should not be restricted to thermodynamically consistent models provided an appropriate choice of density dependent interactions.However, we find that thermodynamically consistent models are best suited for anticipating and understanding these phenomena.In particular, the spinodal condition (13) which is given in terms of the equilibrium free energy, helps rationalize why a tricritical point in the equilbirium limit (Fig. 2 (a)) 'runs away' to infinite density in the presence of self propulsion (Fig. 2 (b)), recovering the standard scenario of discontinuous flocking [2].Moreover, for non-monotonic alignment f , one has that the pair of tricritical points that exist in the equilibrium limit (Fig. 4 (a)) must collide at the maximum of f in the presence of self propulsion, giving way to more exotic phase topologies which could involve either an azeotropic point (Fig. 4 (b)) or a region of counterpropagating bands phase (Fig. 5 (b)).Remarkably, in contrast with standard flocking models [2], the flocking transition is continuous in the presence of an azeotropic point, while remaining discontinuous elsewhere on the phase diagram.
Our models employ a diffusive scaling of the microscopic rates to allow for an analytically exact derivation of a hydrodynamic description.
However, the phenomenology is expected to persist for models without such scaling.Indeed, although the hydrodynamic description of the AIM employs a similar diffusive scaling, microscopic simulations show that the predicted flocking transition phenomenology holds without such scaling, see [12].The mesoscopic interaction range ∆x used here enables these phenomena to appear in our simple one-dimensional setting.We expect that in higher dimensions much of the phenomenology would survive also for locally interacting models.
Given the rich topology of these phase diagrams, it would be interesting to examine how the entropy production rate (EPR) varies across the various phase transitions that we have unveiled.In thermal systems [27,19], EPR measures both the breakdown of time-reversal symmetry and the amount of energy dissipated by nonequilibrium currents.For flocking systems, irreversibility measures have been provided in both microscopic models [45,46,47,48] and field theories [49,50].Yet, these measures should coincide with a quantification of dissipation only for thermodynamically consistent models [51,52,53].In that respect, our approach opens the unprecedented opportunity to properly study how the dissipation varies across the flocking transition.For instance, it could lead to a spatial resolution of dissipation, thus uncovering the profile of energy dissipation associated with travelling bands.are functions of the local fields to leading order: rate Importantly, the leading order terms are symmetric hops with rate D 0 .These dominate over both the asymmetric O (L −1 ) jump rate, and over the slower O (L −2 ) tumbling rate defined in Sec. 2. This timescale separation guarantees that the measure is controlled by the simple symmetric hop rules corresponding to noninteracting random walkers [28].
It is then straightforward to write down the coupled hydrodynamics of the ρ + and ρ − fields The hydrodynamic fluxes are found by coarse-graining, using the hydrodynamic coordinates x = ℓ s i/L and t = γ t/L 2 , with t the time variable of the microscopic dynamics: and, following the same steps, the tumbling rate reads Note that, when sending β → 0 in (A.10), one recovers the fluxes of AIM [13].Finally, one obtains the expressions (6-7) in terms of ρ and m with a few steps of algebra.Interestingly, in the equilibrium limit (Pe = 0), the free energy F serves as a Lyapunov function: where we have used that M > 0 and C is positive semi-definite.The case dF/dt = 0 is achieved in steady state.In the presence of noise, this case corresponds to the system eventually relaxing towards the minimum δF/δρ = δF/δm = 0. Notice, however, that the well-posedness of this hydrodynamic description must break in the equilibrium limit Pe = 0 for any non-convex free energy F (7).It is cured by retaining sub-leading interfacial terms O (L −2 |∇ρ| 2 , L −2 |∇m| 2 ).

Appendix B. Spinodals and binodals in the equilibrium limit
As mentioned in Sec.3.1, the density spinodal (15), which is complementary to the magnetization one (13), marks the boundary of the convex region of the single-variable function F [ρ, m 0 (ρ, T )] [42].Indeed, this boundary is given by where we used the fact that, in the ordered state m 0 (ρ 0 , T ), we have Differentiating (B.2) with respect to ρ yields ∂m 0 (ρ, T ) ∂ρ Plugging the relation (B.3) into (B.1), it follows that finding the density spinodal is equivalent to locating the onset of local non-convexity of the two-variable free energy: The density spinodal can be expressed in an implicit parametric form, since all partial derivatives entering (B.4) are explicit functions of ρ and m.Indeed, we can use the parametric representation ( 14) The binodals are found by the common-tangent construction over the singlevariable free energy F [ρ, m 0 (ρ 0 , T )], which can result in several possible phase-diagram topologies [42].The binodals correspond to the coexistence between dense and dilute phases.For monotonic f (ρ) in Sec. 3, it is always the dense and the dilute phases which are respectively ordered (m = m 0 (ρ, T )) and disordered (m = 0).Yet, one can have the opposite case for non-monotonic f (ρ), as explained in Sec. 4. For illustration purposes, considering here the case of monotonic f (ρ), the condition for common tangent is given by F [ρ l , m 0 (ρ l , T )] = F (ρ g , m = 0) + ∂ ρ F (ρ g , m = 0) (ρ l − ρ g ) , ∂ ρ F [ρ l , m 0 (ρ l , )] = ∂ ρ F (ρ g , m = 0) , (B.7) which explicitly reads (B.8) where we have used the parametric representation z = m 0 (ρ l , T )/ρ l .Substituting the temperature relation (B.5) into (B.8),we get which finally gives the two binodals in an implicit parametric representation for ρ l (z), ρ g (z) and T = T (z).The tricritical point is defined as the meeting point of the two binodals (B.9), which is also the meeting point of the two spinodals (B.6) and (13).This occurs at vanishing z = m 0 /ρ 0 .Indeed, expanding either one of these expressions at small x leads to (16).

Appendix C. Linear instability in the nonequilibrium case
To derive the expression of the density spinodal in (18), we examine the linear stability of the dynamics (5-6) around the homogeneous solutions ρ (x, t) = ρ 0 + δρ (x, t) and m (x, t) = m 0 (ρ 0 , T )+δm (x, t), where m 0 is given by (14).Using the Fourier convention for an arbitrary field X: (C.9) The condition of vanishing real part for the eigenvalue (C.9) yields the ordered spinodal (18).noted in Appendix A, the equilibrium hydrodynamics becomes ill-defined for any nonconvex free energy F since this free energy is missing the standard interfacial energy terms and the width of the interfacial profile will vanish.The self propulsion terms in (6) cure this and assign a finite interface width, but this decreases with decreasing temperature and is hard to resolve numerically.Luckily, due to the same reasoning, close to vanishing temperatures the binodals of the non-equilibrium phase diagram Pe > 0 must approach those of the equilibrium Pe = 0 limit.As explained in Appendix B, the latter are determined analytically without the need of solving the hydrodynamics (5-6), which for the equilibrium case are ill-defined.We use these equilibrium asymptotics to extend the binodal curves of the non-equilibrium phase diagrams all the way to zero temperature (the numerically computed curves at small temperatures are found to match them smoothly).

Figure 1 .
Figure 1.The active flocking model with thermodynamic consistency.Each particle can be in either one of two states (+ and −) which determine the direction of biased diffusion (respectively to the right and left).The aligning Hamiltonian H constrains both the change of particle states and their diffusive hops.

Figure 2 .
Figure 2. Phase diagrams for (a) the equilibrium case (Pe = 0) and (b) the nonequilibrium regime (Pe = 0.5), with monotonic f given by (10) with r = 1.The colored dashed and solid curves are spinodals and binodals, respectively.The black solid curve is a critical line.The black dot marks a tricritical point.(c) Traveling band (T.B.) profiles for Pe = 0.5, ρ 0 ≃ 3.62, and T ≃ 0.71.

Figure 3 .
Figure 3. Equilibrium (Pe = 0) bifurcation diagrams for (a) T > T tri and (b) T < T tri , with monotonic f given by(10) with r = 1.Curves mark the different branches of solutions for the magnetization of homogeneous states(14).Dashed blue and red lines mark instabilities of the two branches.Both originate at the density value ρ c given by the curve(13).The instability of the magnetized branch (dashed red) terminates at the density value given by the density spinodal(15).The region of densities enclosed between these two values is linearly unstable for both branches.Red and blue dots mark the binodal densities of the phase-separated state (P.S.), which is the globally stable state for any mean density between these points (including densities for which all homogeneous states are unstable).The black dot marks a critical point.

Figure 6 .
Figure 6.Time evolution of the counterpropagating bands for Pe = 1, ρ 0 ≃ 1.8, and T ≃ 0.17.(Top) The counterpropagating bands of the ρ ± fields, in cyan and magenta, barely scatter upon encounter.(Middle) The corresponding density ρ = ρ + + ρ − and magnetization m = ρ + − ρ − display an oscillating dynamical state with volume fraction periodically varying in time.(Bottom) the corresponding space time plot for the density profile.
Figure 7.(a-b) The ± counterpropagating bands profiles are perfectly antisymmetric: ρ + ≃ ρ (x − V t) and ρ − ≃ ρ (−x − V t) coincide when reflected and superimposed.(c) The periodic evolution of the phase volumes of the intermediate and liquid phases ν i and ν l respectively, for parameter values in Fig.6.The complementary gas phase volume (not shown) is simply ν g = 1 − ν i − ν l .
which gives z = z(ρ 0 ) in an implicit form, from which follows the curve T (ρ 0 ) for the density spinodal.This procedure allows us to determine the corresponding curves in Figs.2(a), 4(a) and 5(a).