Exact hydrodynamics and onset of phase separation for an active exclusion process

We consider a lattice model of active matter with exclusion and derive its hydrodynamic description exactly. The hydrodynamic limit leads to an integro-differential equation for the density of particles with a given orientation. Volume exclusion results in nonlinear mobility dependent on spatial density. Such models of active matter can support motility-induced phase separation, which occurs despite the absence of attractive interactions. We study the onset of phase separation with linear stability analysis and numerical simulations.


Introduction
Active matter systems are composed of interacting agents that consume energy from their environment.They exhibit a rich variety of collective phenomena including flocking [1], lane formation [2,3], clustering [4][5][6] and pattern formation [7].Active matter systems can be found in many natural and synthetic contexts [1], ranging from biological systems such as bacteria [8,9], microtubules [10] and animal groups [11], to physical and chemical systems such as colloids [4] and self-propelled rods [7].Several mathematical models have been developed using individual or agent-based approaches to capture the essential features of these systems.These models can help us understand the emergence of complex behaviours from simple rules at the microscopic level.
An important class of active matter models is based on self-propelled particles.In some of these models, the propulsion directions of nearby particles tend to align with each other, which naturally leads to flocking behaviour [1].Without such alignment, one often observes a remarkable phenomenon of spontaneous condensation [12], occurring in the absence of attractive forces between particles.This behaviour, called Motility-Induced Phase Separation (MIPS), arises from the interplay between particle crowding and velocity reduction.When particles cluster together due to their self-propulsion, they also reduce their speed due to collisions and interactions.This creates a positive feedback loop that enhances the clustering and leads to a separation between a dense liquid-like phase and a dilute gas-like phase.MIPS has been extensively reviewed in [12].It is a striking example of how active matter can display emergent collective behaviours that are different from those of passive particles.
Numerical simulations have proven extremely valuable for understanding the emergence of MIPS and other complex phenomena in active matter models [5,[13][14][15][16][17]. Analytical results can complement such studies by providing exact and general insights into this behaviour, but they are challenging to obtain, due to the non-equilibrium and nonlinear nature of active matter dynamics.One possible approach is to analyse coarse-grained models that capture the essential features of active matter systems at a macroscopic level, for example through density and velocity fields.Several phenomenological coarse-grained models have been proposed in the literature [12,16,[18][19][20][21][22].Such models have yielded important insights, but deriving them from microscopic models relies on various assumptions and approximations, that are hard to justify.
A more rigorous and systematic approach is to derive coarse-grained descriptions directly from models of interacting particles [23,24].This can be achieved via the theory of hydrodynamic limits, in which macroscopic evolution equations are derived directly from microscopic stochastic rules.It is most naturally applied to lattice gas models, where particles hop randomly on a discrete lattice according to some local rules.By taking a limit where the lattice spacing tends to zero, one can obtain continuum equations that exactly describe the large-scale behaviour of the system.This method has been successfully applied to various classes of lattice gas models [24][25][26] and can provide valuable insights into the physics of active matter [23].
Following the recent interest in such lattice models of active systems [16,23,27], this paper considers an Active Lattice Gas (ALG) where particles interact via exclusion, with diffusive dynamics for their orientations.We show that the macroscopic limit of ALG can be described exactly by an integro-differential equation.To explore the onset of MIPS, we analyse the stability of the homogeneous stationary state, in the parameter space of occupied volume fraction and Péclet number, through linear stability analysis and numerical simulations of both the microscopic and macroscopic models.

Previous results
A variety of active lattice gas models have been studied in recent years [16,23,27].In these models, each particle is equipped with an orientation variable, which determines the direction of its self-propulsion.In the active lattice gas model of [23], the orientation can point in a finite number of drift directions, aligned with the lattice axes.In the terminology of [25], that model is of gradient type (see Definition 5.2), and the macroscopic model can be derived exactly and explicitly.This enables a derivation of the phase diagram, including the limit of stability of the homogeneous state, and the density profiles for systems with inhomogeneous steady states.The macroscopic description of a more complex alignment model, with continuous particle orientation, was proven in [27]: the latter model is of non-gradient type, which makes the analysis more challenging [26,28].
In this work, we consider a different two-dimensional model where particles' orientations are unit vectors that undergo diffusive motion on a torus, similar to off-lattice models of active Brownian particles [5,6].This model was proposed in [16], where an approximate macroscopic description was derived within a mean-field approximation.The resulting behaviour was compared with other coarse-grained off-lattice models, including a characterisation of the linear stability around a homogeneous state and the onset of MIPS.They found that short-range interactions were important for observing MIPS because a similar model with long-range interactions did not exhibit phase separation.In the following, we will characterise the exact hydrodynamic limit for this model without relying on the mean-field assumption.Using this description, we revisit the stability of the homogeneous state and the onset of MIPS.We verify the accuracy of this macroscopic description via a comparison with numerical simulations of the underlying microscopic models.

Structure of the article
The structure of the paper is as follows.In Section 2, we define the microscopic model, present its hydrodynamic limit, and summarise our main results for the stability of the homogeneous state.The linear stability analysis is detailed in Section 3. In Section 4, we present numerical computations of both the macroscopic and microscopic models.We compare the linear stability computation with the numerical simulations and demonstrate excellent agreement between simulations of the microscopic model and numerical solutions of the macroscopic equation.The mathematical proof of the hydrodynamic limit is very technical and follows closely the derivation of [27]: we sketch the main arguments of the proof in Section 5, focusing on the steps where our derivation differs from that of [27].
2 Model definition and summary of main results

Microscopic model: an Active Lattice Gas
The model that we consider is an Active Lattice Gas (ALG).It is defined on a twodimensional N × N periodic lattice We define the occupation configuration η ∈ {0, 1} T 2 N .Each site z of T 2 N is either occupied by a particle with orientation θ ∈ S = [0, 2π), with periodic boundary conditions, (η z = 1, θ z = θ) or empty (η z = 0, θ z = 0 by convention).The site configuration at z is ηz := (η z , θ z ).The set of all configurations is given by Each point z ∈ T 2 N corresponds to a macroscopic point z N ∈ T 2 in the unit torus, on which the macroscopic profile will later be defined.
The initial configuration is chosen in a state of local equilibrium, close to a smooth profile ζ : T 2 × S → R + .More precisely, each site z ∈ T 2 N is initially occupied by a particle with probability ζ(z/N ) := S ζ(z/N, θ)dθ.(For this construction to make sense, we assume that ζ is such that ζ(x) ∈ [0, 1] for all x ∈ T 2 ).The angle θ z is then sampled according to the probability distribution ζ(z/N, θ)dθ/ζ(z/N ).We denote by µ ⋆ N, ζ this initial distribution for our ALG.Informally, the dynamics of particles in the ALG can be described as follows 1.Each particle attempts to perform a nearest neighbour random walk, weakly biased in the direction of their orientation.In particular, a jump in the u direction is attempted at rate , where e θ = (cos(θ), sin(θ)) T is the drift direction, D E is the spatial diffusion constant and v 0 is the self-propulsion speed.2. If the target site of a jump is occupied, the jump is aborted; if the site is otherwise empty, the jump is executed.This is called the exclusion rule.3.Each particle's orientation diffuses according to a Brownian motion on S with diffusion constant D O .
The N -dependence of the particles' hop rates corresponds to a parabolic scaling [25] and enables the exact characterisation of the hydrodynamic limit [27].It also ensures that for N → ∞, a single particle's motion converges to an Active Brownian particle with diffusion constant D E , speed v 0 and angular diffusion D O .The physical interpretation of this scaling for the hydrodynamic limit of interacting particles is discussed in Sec.2.3, below.
To formalise the model, we identify the generator of this stochastic process, which can be decomposed into three parts, according to the N -dependence of the transition rates.This generator is whose elements we now define.The generator is defined on The nearest-neighbour simple symmetric exclusion process generator L S is defined as Similarly, the asymmetric exclusion operator L A is: On its own, L A is not a Markov generator because it has negative jump rates.However, when added to the symmetric part, N 2 L S + N L A becomes well-defined for N large enough.All numerical simulations in this work satisfy this condition.In the expressions above, ηz,z ′ denotes the configuration where the occupation variables ηz and ηz Finally, L O is the generator for the diffusion of the orientation of each particle Together, the initial distribution µ ⋆ N, ζ and the Markov generator L N uniquely define a Markov process (η(t)) t≥0 , whose macroscopic behaviour we want to characterise.

Hydrodynamic Limit of the ALG
We now introduce the exact hydrodynamic limit of the orientation density for the ALG, whose derivation is outlined in Section 5.
The macroscopic equation of the ALG depends on the self-diffusion coefficient of a tagged particle in a simple symmetric exclusion process (SSEP).(The same dependence appears in the models studied in [27].)Definition 2.1 (Self-diffusion coefficient).In Z d , consider an infinite SSEP and a tracer particle placed initially at the origin.We place a particle at each site in Z d \ {0} independently with probability ρ.The position of the tracer at time t is X t ∈ Z d and we denote by Q ρ the distribution of the resulting process.Then the self-diffusion coefficient is Spohn [29] obtained a variational formula for the self-diffusion coefficient.The function d s was shown to be Lipschitz continuous by Varadhan in [30] for d ≥ 3 and was later shown to be of class C ∞ for d ≥ 2 by Landim, Olla and Varadhan in [31].Various approximations to d s are discussed in Section 2.3, below.
The hydrodynamic equation describes the evolution of the local macroscopic density of particles with angle θ, which we denote by f (t, x, θ).It can be interpreted formally as where ε is the side length of a small macroscopic box, in which we measure particle density, where 1 ≫ ε ≫ 1/N ≫ dθ.
We will sometimes change our notation to f t (x, θ) or f θ (t, x) to emphasise the dependency on one of the variables (see e.g.Definition 5.1 and remark 5.1 below).We also define the polarisation p and the particle density ρ, by [Recall that e θ = (cos(θ), sin(θ)) T is the drift direction of a particle with orientation θ.]It is also convenient to introduce two coefficients D and s as With these definitions, the macroscopic evolution equation for the density field f of the ALG defined in (10) can be derived rigorously from the microscopic model as In the simpler models studied in [23], particle jumps are only biased in two directions.However, a straightforward generalisation to a continuum of directions would also lead to integro-differential equations similar to (13), but in a simpler form due to the gradient nature of the models considered there.Similar integro-differential equations have been studied in the context of flocking models (see e.g.[32] and references therein), but unlike in our case, where the evolution equation can be derived rigorously from the microscopic model, deriving such evolution equations typically relies on mean-field assumptions, meaning particle's interactions are averaged over a neighbourhood.
Note that to make this hydrodynamic limit as easy to read as possible, we stated it here in the informal form (13) without regard for the well-posedness of the equation and without proper justification for the convergence result (10).The rigorous statement of the ALG's hydrodynamic limit involves the process's empirical measure π N t , defined by its integral against test functions H ⟨π N t , H⟩ := We postpone the full statement of the hydrodynamic limit to Section 5, where we give a precise definition of what we mean by solutions of (13), and state in Theorem 5.1 that the empirical measure converges in probability to the local density field (10), namely The proof of the hydrodynamic limit for the ALG (Theorem 5.1) involves severe mathematical challenges [27], which mostly arise because the model is of non-gradient type (See Definition 5.2).However, the ALG considered here is very close to those analysed in [27]: the only significant difference is that the orientation dynamics is a diffusion on S instead of a Markov jump process.While this means that the hydrodynamic limit proofs of [27] cannot be used verbatim in this case, the proof of Theorem 5.1 is a straightforward generalisation of that work.The main steps of the proof are described in Sec. 5, and the aspects that differ from [27] are explained.(The diffusive orientational dynamics simplifies some aspects of the proof, see Remark 5.3) Whilst the hydrodynamic limit is an extension of [27], numerical and analytical analysis of ( 13) is hindered by the density dependence of the diffusion constant: this enters via d s (ρ), for which no explicit formula is available.As noted in [27] the proof of the hydrodynamic limit can easily be extended to higher dimensions d ≥ 2.
To understand the physical structure behind the hydrodynamic equation ( 13), we observe that if the particles are passive, v 0 = 0, the microscopic model is a reversible Markov chain and its invariant measure is a product over the lattice sites.In this case, we explain in Remark 5.2 that its free energy is We also define a parametric distribution which plays the role of a mobility where f (θ) := f (•, •, θ) and δ θ (dθ ′ ) is the Dirac distribution at θ on S, see Remark 5.1 for a further discussion.Then the hydrodynamic equation ( 13) can be factorised as In the passive case v 0 = 0, this is a gradient-flow for S, as one would expect for the hydrodynamic limit of a microscopic model with time reversal symmetry [33].
In general, the objects in square brackets are thermodynamic forces that drive the evolution of the model, note that the active self-propulsion is not the gradient of any free energy, consistent with the fact that the microscopic model has a non-equilibrium steady state.This form for the hydrodynamic equation is discussed further in Remark 5.2.

Properties of the hydrodynamic equation
The macroscopic evolution described by ( 13) depends on five parameters.Three of them appear in (13) itself: the spatial diffusion constant D E , self-propulsion speed v 0 , and angular diffusion constant D O .There are also two further parameters: the domain size, and the occupied volume fraction of the active particles.
The lattice spacing (or 'particle size') in the microscopic model has already been taken to zero in a domain of fixed size.The parameters of the model encode two additional length scales where ℓ p (resp.ℓ D ) represents the typical distance travelled by a particle because of its drift (resp.because of its symmetric motion) before its angle reorients significantly.
For the model of [23], the widths of the interfaces separating the dense and dilute phases of MIPS are proportional to ℓ D .For a physical interpretation of the hydrodynamic limit, we emphasize once again that the N -dependence of the generator L N ensures that ℓ p and ℓ D are of order unity as N → ∞, while the range over which particles interact by exclusion is tending to zero.In physics, it is more common to take ℓ p as a fixed multiple of the particle size [5,6,20], on the grounds that both length scales are intrinsic properties of individual particles.In this case, various approximate descriptions are available for the (fluctuating) hydrodynamic behaviour [15,17,20], but the hydrodynamic limit cannot be characterised exactly.Hence, our focus is on the generator L N , for which exact computations are possible.
Overall, a natural set of dimensionless parameters that is sufficient to capture the behaviour of (13) consists of which are the occupied volume fraction, the Péclet number and the diffusive length scale respectively.In the following, the domain size is set to unity, rescaling time as t ′ = t/D E .In dimensionless parameters, (13) becomes The analysis of ( 13) is hindered by the density dependence of the diffusion constant: this enters via d s (ρ), for which no explicit formula is available.However, previous work has led to very accurate numerical approximations.Our approach here is to approximate d s by a polynomial derived in [34], which is asymptotically exact at both low and high density.This approximation -defined in (34) below -is used for all numerical computations in the following.An alternative asymptotically accurate formula was derived in [35]; other approximations are available as mean-field-type approximations [36,37], variational bounds [38,39] or low-rank tensor approximations [40,41].Our polynomial approximation is employed for its accuracy at high and low density and ease of use.Other approximations, such as mean-field [16], lead to similar qualitative behaviour, but the quantitative predictions for the onset of MIPS are less accurate.
With this numerical approximation of d s in hand, we have analysed MIPS in our ALG in two different ways.In Section 3, we study the linear stability of the homogeneous solution to (13).A linear instability at this level signals that MIPS is taking place.A detailed discussion is given below, but we summarise the behaviour in Fig. 1, which shows the onset of linear stability as a function of Pe and volume fraction ϕ.This is compared with the behaviour of numerical solutions of (13).For low densities, we observe that the system is always linearly stable.For higher densities, the onset of instability in the numerical solutions coincides with the boundary of linear stability.
In addition to these results for stability, Section 4 presents numerical solutions of the macroscopic equation (13), which is compared with particle-based stochastic simulations of the microscopic model.The agreement is excellent, including in the MIPS regime: this is consistent with the rigorous derivation of the hydrodynamic limit, combined with the numerically accurate approximation for d s .

Linear Stability
We consider the linearisation of equation ( 13) around the homogeneous steady state f * ≡ ϕ 2π .We insert f = f * + δ f , and ρ = ϕ + δ ρ, for δ ≪ 1. Furthermore e θ dθ = 0 and therefore p = p[f ] = δ p := p[δ f ].The resulting linearised problem is: The linear stability analysis considers solutions of the form f ∝ e λt , in which case (22) becomes an eigenproblem.Since the domain is periodic, the eigenfunctions take the general form In the spatial variable, the equation is closed at each frequency, and therefore we can consider each spatial frequency ω separately.To find the eigenvalue with the largest real part, we seek the lowest non-trivial frequency in x, so it is sufficient to consider ω = (ω, 0) = (2π, 0), which corresponds to a wave along one axis.Furthermore, taking a wave in the x 1 direction, equation ( 22) immediately yields B k = 0. Hence it is sufficient to consider Then (22) becomes where and For k ≥ 2, (25c) are the difference equations derived from the Mathieu Equation [42].
In that case, (25c) applies for all k, including k = 0, 1.In our case, there are additional terms at k = 0 and k = 1 due to the terms ρ and p in (22), which are integrals of f .The simplicity of the underlying integrals means that these terms only appear for k = 0, 1, which means that we are still able to close the equations, as follows.Equation (25c) has two linearly independent solutions [43], one diverges and the other behaves as as k → ∞.Hence a convergent solution to the full problem ( 25), (A k ) k≥0 ∈ L 2 , only exists when λ takes specific values, called the characteristic numbers.For numerical computation, it is useful to rewrite the recurrence relation as an eigenvalue problem of a countable tridiagonal matrix whose W is defined in (26).For each eigenvalue of ( 28), there is a corresponding eigenvalue λ n of the corresponding truncated n × n matrix: the resulting sequence (λ n ) n∈N behaves for large n as (cf.[44]) Combining with (27) . Thus we can approximate the eigenvalue with the largest real part λ max by λ max n for large n. Figure 2 plots Re(λ max n ) = 0 for n = 2, 4, 6, 40, displaying the rapid convergence to Re(λ max ) = 0, the boundary between stability and instability of the linearised problem (22).Remark 3.1 (Sharp-interface limit).In the limit ℓ ≪ 1, the inhomogeneous steady states of (13) consist of uniform regions separated by narrow interfaces.These are the phase-separated states that occur in MIPS, where it is conventionally assumed that the interfacial widths are much smaller than the sizes of the domains of the individual phases.For a fixed ω, D E and ϕ, let Pe * (ℓ) be the critical Peclet such that λ max = 0.In the limit ℓ → 0 assume that Pe * (ℓ) → Pe * (0).Then ℓ −k A k (ℓ) → c k , and in particular (25) implies det(W ) → 0. Setting lim ℓ→0 det(W ) = 0 yields the limiting curve which has asymptotes at ϕ = 1 and when Pe Fig. 2 Plots of Re(λ max n ) = 0 as determined from (29).As n → ∞, these curves converge to the boundary between stability and instability of the linearised problem (22).The case n = 6 is indistinguishable from n = 40 due to the rapid convergence of the scheme.The area above the limiting curve corresponds to the linear problem having a positive eigenvalue.We take the representative parameter value ℓ = 0.5.
positive in some neighbourhood of ϕ = 0. Therefore, for volume fractions below some critical level, the homogeneous steady state is stable for all Pe.Remark 3.2 (Run and tumble particles).In the case of run and tumble particles, the diffusive dynamics of particle orientation is replaced by uniformly random reassignment at a given rate D O .In (13) As explored in [45], run and tumble dynamics result in similar phenomenology to diffusive orientations.However, in the linear stability calculations above resulting in increased dependence on higher-order Fourier modes.
Several comments should be borne in mind when interpreting the results of this linear stability analysis.For systems described by ordinary differential equations (with sufficient regularity), the boundary of linear and nonlinear stability are the same [46,Sec. 9.3].They are also equivalent in quasilinear partial differential evolution equations where terms containing the highest order spatial derivatives are linear so that the regularising effect of the linear terms outweighs the nonlinear terms [47,Sec. 4.3], allowing the solution to be written explicitly using a semigroup generator.However, this does not apply to (13) because of the nonlinear terms involving secondorder derivatives.We expect agreement between linear and nonlinear stability, but a rigorous proof is beyond the scope of this paper.Similarly, nonlinear terms can cause difficulty in rigorously proving exponential convergence to equilibrium for linearly stable cross-diffusion systems [48].
For the physical setting of MIPS, the limit of linear stability of the homogeneous state is known as the spinodal (30).In the unstable regime, Equation ( 13) has a stable inhomogeneous steady-state solution, which corresponds to phase separation.(In fact, there is a whole family of such solutions related by spatial translations.)However, the physics of phase separation is more complicated than this, because one generally expects a parameter regime where (13) supports homogeneous and inhomogeneous steady solutions, which are both stable.The boundary of this region is known as the binodal [12].In this regime, the microscopic model exhibits metastability [49]: it supports fluctuations which can lead to spontaneous transitions between the metastable states.The hydrodynamic description does not capture these; they would require an analysis of large deviations of the empirical measure [24,25,50,51], which is beyond the scope of this work.

Numerical PDE scheme
We use a first-order finite-volume scheme to obtain numerical solutions to (13).We first rewrite the equation in the form where M x1 , M x2 , M θ are scalar mobilities and U = (U x1 , U x2 , U θ ) is the velocity vector.The mobilities are defined by M x1 = M x2 = f d s , M θ = f and the velocities are given by where . We note that (31) is not a Wasserstein gradient flow as the velocity U cannot be written as the derivative of an entropy (unless v 0 = 0).
The function d s is not known explicitly, so we approximate it by the polynomial approximation derived in [34] which is asymptotically exact to the first order at both low and high density: where α = π/2 − 1.
We discretize the phase space Υ the time interval [0, T ] is discretized by t = ∆t, 2∆t, . . ., T .Define the cell averages We use the finite-scheme where F is the upwind flux and similarly for F x2 and F θ .We discretize the integrals where e θ = (cos(θ), sin(θ)) T .The velocities are approximated by centred differences or centred averages, and similarly for U x2 , U θ .Finally, we discretize the system of ODE (36) by the forward Euler method with an adaptive time stepping condition satisfying where a x1 = max{|U x1 i,j,k |}, a x2 = max{|U x2 i,j,k |}, and a θ = max{|U θ i,j,k |}.In [52], they derive a CFL condition for a nonlocal model of active particles and show that it leads to a positivity-preserving numerical scheme.A key difference is that their scheme is second-order in phase space as they use a linear density reconstruction at the interfaces that preserves positivity.In contrast, we follow [16] and use the values at the centre of the cells instead.In our numerical tests (in which ∆t ≤ 10 −5 ), we observe (40) to be sufficient to preserve positivity.
Figure 3 shows two examples of instability from an initial random perturbation The particle density ρ is displayed in the first row, Fig. 4 Evolution of the L2 distance (42) of the particle density f from uniform, with initial condition f 0 given by (42).For parameters D E = 1.0, ϕ = 0.92, ℓ = 0.5 with stable case, Pe = 8.0, the unclassified case Pe = 9.0 and unstable case Pe = 10.0.
and the polarisation p is displayed in the second row.In both cases, the system phase separates into high and low-density regions.We observe that larger spatial gradients occur in regions with significant polarization.Roughly speaking, the polarisation tends to lie parallel with the density gradient, so that the diffusive and advective terms balance in (13).

Nonlinear instability
To analyse the stability of the homogeneous solution to (13), we solve the macroscopic model ( 13) using the finite volume scheme described in Subsection 4.1 with an initial perturbation around the homogeneous state.In particular, we take the initial condition to be the eigenfunction corresponding to the eigenvalue with the largest positive real part, as given by (24).That is where Re(x) represents the real part of x.Having obtained a suitable numerical solution, we study the growth of the perturbation over time using a discretized L 2 norm: In particular Pe = 10.0 is classified as unstable as it surpasses 2δ = 2 × 10 −4 , indicated with a dashed line.The case Pe = 8.0 is classified as stable as it is decreasing and below 2δ, whereas Pe = 9.0 is unclassified because it is below 2δ but increasing at the final time t.
To obtain Figure 1, we iterated this process and performed a parameter sweep for ϕ = 0.44, 0.46, . . ., 1.0 and Pe = 0, 1, . . .40.For reference, we also plot the limit of linear stability (28).The numerically computed instability matches the limit of linear stability across all densities.

Comparison of macroscopic and microscopic systems
In this section, we compare the results obtained from the macroscopic equation with stochastic simulations of the microscopic models.For comparison with the macroscopic case, we define the local particle density in the microscopic model: and the local polarisation is We initialise the model according to Section 2.1, with N 2 sites and a uniform initial profile ( ζ ≡ ϕ 2π ).The system dynamics combines both a jump and diffusion process.A fixed timestep ∆t is used to discretize the diffusion process.Due to the separation of scales, the timestep can be chosen such that N −2 D E ≪ ∆t ≪ 1/D O , we take ∆t = 0.001.For each time increment ∆t, we run a Gillespie algorithm, with fixed θ z , until the time elapsed δt surpasses ∆t.The Gillespie algorithm simulates a finite state space Markov chain exactly [53] .At each time step we first sample the next time interval according to an exponential distribution.We then sample a jump, weighted by each jumps probability of occurring first.We then execute the jump and update the corresponding jump rates.Once the time elapsed δt surpasses ∆t, we update each θ x , sampling from N (θ x , 2D O δt) mod 2π.
Figure 5 shows two examples of the system at t = 0.0, 0.1, 1.0 for N = 32, ϕ = 0.7 and Pe = 10, 15.The particles are represented by triangles of size 1  N pointing towards e θx with colours corresponding to their orientations.Each cell is also shaded in proportion to the local particle density, with ε = 1 16 .Figure 6 shows two examples of the system at t = 0.0, 0.5, 1.0 for N = 128 and ϕ = 0.5, Pe = 30 and ϕ = 0.7, Pe = 12.These results show that MIPS is present for both Pe = 12, 30.The particles in each system self-organise into a dilute region with almost zero density and a dense region with particles in close packing.We also observe a strong particle alignment in the boundary between the two regions.
We emphasise that while the initialisation protocol of the system is homogeneous and invariant under spatial translations, dense and dilute regions have developed.This is a spontaneous breaking of translational symmetry in the occupation.This effect can be quantified by defining [54,55]  Estimates of Φ(η t ), averaged over 10 samples, are shown in Figure 7 for N = 64, 128 and Pe = 6, 8, 10, 12.They indicate that the system is macroscopically homogeneous for Pe = 6, 8, but the symmetry is spontaneously broken for Pe = 10, 12.The stability of the macroscopic model displayed in Figure 1 is in good agreement with Figure 7, for both which Pe = 6, 8 are stable and Pe = 10, 12 are unstable.
We close this Section with a quantitative comparison between the steady states of the microscopic and macroscopic models.To this end, we take parameters ϕ = 0.5, Pe = 30, ℓ = 0.5.For the microscopic model, we compute the local density ρ εN , with ε = 1  16 , for N = 32, 64, 128.After waiting for the system to reach equilibrium, we averaged over times t = 2.0, 2.1, . . ., 4.0.Combining data from six such realisations, we plotted histograms of the local density, shown in Fig. 8. Since the homogeneous steady state is unstable, the system phase separates, and two peaks emerge; the low-density peak corresponds to the dilute ("gaseous") phase; the dense phase is a "liquid".The local maxima of the histogram are the modal densities of the two phases.
To compare this microscopic result with the macroscopic description, we define a corresponding local density which we compute from our numerical solutions of the hydrodynamic equation, after initialisation with a random perturbation around the uniform steady state, and averaged over times t = 2.0, 2.1 . . ., 4.0.Histograms of this local density are also shown in Figure 8.For increasing N , the histograms computed from the microscopic model converge towards the result of the macroscopic model, as predicted by Theorem 5.1.

Hydrodynamic Limit
In this section, we give a precise statement of the hydrodynamic limit formally described in (13).To arrive at this macroscopic description of our ALG, we define its empirical measure.The hydrodynamic limit is then obtained as a characterisation of this (random) measure in the limit N → ∞.Let M(T 2 × S) denote a non-negative measure on the continuous configuration space endowed with the weak topology and denote the space of right continuous and left limited trajectories on M(T 2 × S).Each trajectory of the process η = (η(t)) t∈[0,T ] admits a natural image in M [0,T ] through its empirical measure where δ x,θ ∈ M(T 2 × S) is Dirac measure at (x, θ).Note that π N t can equivalently be characterised by its integral against a test function H, see (14).We endow M [0,T ] with Skorohod's metric, and the set P(M [0,T ] ) of probability measures on M [0,T ] with the weak topology.We now define as the distribution of the trajectory π N = (π N t ) t∈[0,T ] of the empirical measure of our process η = (η(t)) t∈[0,T ] .In deriving the hydrodynamic limit of this model, we explain that this trajectory π N t (dx, dθ) converges as N → ∞, in probability, to a deterministic trajectory π t (dx, dθ) = f (x, θ, t)dxdθ, whose density f represents the local density of particles with angle θ, and is a solution of a hydrodynamic equation (13).Throughout this section, through abuse of notation, we will use a subscript to denote the time argument, f t (x, θ) = f (t, x, θ).Definition 5.1 (Weak solution to ( 13)).We call a measure-valued trajectory (π t ) t∈[0,T ] ∈ M [0,T ] a weak solution to (13), with initial condition ζ(x, θ), if the following are satisfied 2. For all t ∈ [0, T ], the measure π t is absolutely continuous w.r.t the Lebesgue measure on T 2 × S, that is, there exists a density profile f : )), and is therefore spatially differentiable.

For all functions H
where ⟨f, g⟩ = T 2 ×S f (x, θ)g(x, θ)dxdθ.We can now state precisely what it means for (13) to be the hydrodynamic limit of our system.Theorem 5.1 (Hydrodynamic Limit).The sequence (Q N ) N ∈N (49) is weakly relatively compact, and any of its limit points Q * are concentrated on trajectories π ∈ M [0,T ] which are weak solutions of (13) in the sense of Definition 5.1.
where we defined for any smooth function for the microscopic approximation of ∂ xi H t .The currents j s , j a are defined as the angular distributions which represents the symmetric current of particles with angle in [θ, θ + dθ] going through the edge (z, z + e i ), and which represents the asymmetric current going through (z, z + e i ), and is the strength of the driving asymmetry in direction i for a particle with angle θ.
Note that because of the angular diffusion, the last term in (53) appears directly as a function of the empirical measure π N s .Further note the extra factor N multiplying the symmetric current j s , which we must balance by a discrete gradient.To prove the hydrodynamic limit, one needs to close Equation ( 53) by replacing the two remaining microscopic quantities, j s and j a , by functions of the empirical measure, which are their macroscopic counterparts.This will allow us to take the limit N → ∞ in (51).Although j s and j a are distributions in θ, we will, by abuse of language, refer to them as local functions of the configuration, since they only depend on the configuration through a finite number of sites.Our goal is to replace local functions g z (η, dθ) at site z by their average value, which is a function Φ(ρ εN x , ρεN x (dθ)) of the mesoscopic state of the process, where is the empirical mesoscopic density around site x = z/N ∈ T 2 , and is the empirical distribution of the angles in a mesoscopic box around position x = z/N ∈ T 2 in the macroscopic box.
This substitution, of local functions by their average value, typically stated as a Replacement Lemma, holds as N → ∞ and then ε → 0. One can see, for fixed ε and letting N → ∞, that both ρ εN x and ρεN x can be well approximated by functions of the empirical measure: given a site z = N x, for example, one easily checks that where where

Grand-canonical states and local equilibrium
A crucial element in the Replacement Lemma, which amounts to a local law of large numbers for the process, is to identify the local distribution at any given point z of the configuration, namely the probability of observing at time t a local configuration {η z ′ , z ′ ∈ Λ(z)} in a large microscopic box Λ(z) around z.Because the symmetric displacement of particles happens at a very fast rate N 2 , non-conserved quantities relax quickly to equilibrium, and one can expect that the local distribution of the process will be closely approximated by the stationary states of L S , parametrised by its conserved quantities.
The only conserved quantities for L S are the total number of particles and their angles.Hence we fix a distribution ρ on S with total mass ρ := S ρ(dθ) ∈ [0, 1].One then defines the grand-canonical states µ ⋆ Λ, ρ in a finite box Λ ⊂ Z 2 as a product distribution where each site of Λ is occupied by a particle independently with probability ρ, and particles' angles are chosen independently with distribution ρ/ρ : We denote by E ⋆ Λ, ρ the corresponding expectation.
For the reasons listed above, we expect that the distribution of our Markov process at site z and macroscopic time t is locally close, as N → ∞ and then ε → 0, to µ ⋆ Λ, ρεN z .More precisely, given a local function g of the configuration, we expect that for any Λ containing the support of g, (1). (64) This property, called local equilibrium, is a crucial argument in the proof of hydrodynamic limits.
Fix a density profile f : T 2 × S → R + , satisfying ρ(x) := S f (x, θ)dθ ∈ [0, 1], ∀x ∈ T 2 .We associate to f a local equilibrium distribution, defined analogously to (63) This means that under µ ⋆ N,f each site z is occupied w.p. ρ(z/N ), and if it is occupied, the angle of the particle at site z is distributed as f (z/N, θ)dθ/ρ(z/N ); this is similar to our prescription for initialisation of the ALG.Given a reference density ρ ⋆ ∈ [0, 1], we simply denote by ν ρ ⋆ := µ ⋆ N,f ρ ⋆ the equilibrium distribution associated with the uniform profile f ρ ⋆ (x, θ) ≡ ρ ⋆ 2π , in which each site is occupied with constant probability ρ, and particle's angles are chosen uniformly in S.
which respectively represent the cross-diffusion coefficient and the cross-mobility coefficient.We now introduce which can be interpreted as cross-compressibility.Above, Cov f (•, •) represents the covariance of the two variables with respect to the equilibrium state (65) on Λ = Z 2 .Since the latter is a product distribution, all contributions in (67) vanish, except for z = 0, which gives the second identity in (67).Given a function f (t, x, θ) and a distribution µ on S, define and the vector e : θ → e θ .The hydrodynamic equation (13) straightforwardly rewrites Then, straightforward computations show that the Einstein relation holds, namely Remark 5.2 (Entropy and gradient-flow structure).We now define the relative entropy Assuming that f is smooth, it is straightforward to show, by the law of large numbers, that the relative entropy, once suitably rescaled in N , is related to the free energy (16) as for a suitable constant C 1 = C 1 (ρ ⋆ ).This relative entropy can be interpreted as the large deviation rate functional for the empirical measure π N in the ALG's stationary state ν ρ ⋆ .The functional derivative of S is given by which yields the factorised form of the hydrodynamic equation given in (18).

Non-gradient replacement Lemma
Note that the gradient ∂ xi,N H s in (53) comes from the difference of particle currents between current going in (from z − e i to z) and currents going out (from z to z + e i ) at any site z, and balances out a factor N in the symmetric and asymmetric parts of the generator L N (cf.Equation ( 3)).By "balancing out", we mean that the extra factors N rescaling the exclusion processes' generator (both symmetric and asymmetric parts) get absorbed into the discrete derivative of the test function so that the limiting contribution is of order 1.One main difficulty in deriving the hydrodynamic limit of the ALG comes from the non-gradient nature of the generator, meaning that the symmetric part L S of L N does not readily act as a local Laplacian.More explicitly, the definition of gradient type is given in Remark 2.4 [25] as Definition 5.2 (gradient-type).Let E be a subset of N, corresponding to particle occupation number.A translation invariant nearest neighbour particle system η(t) on E T d N with generator L N is said to be gradient type if there exist cylinder functions g i,n and finite range functions p i,n such that the current can be written where τ is the translation operator and In our case, recalling that particles have orientations, this entails that the symmetric current j s z,z+ei (55) cannot be written in form of Definition 5.2 for any local function g of the configuration.Note that the gradient condition (76) can be interpreted as a microscopic Fick's law and bears no relation to the gradientflow form explored in Remark 5.2.Without a gradient decomposition (76), one cannot perform a second summation by parts on the test function, which would allow us to balance out the second factor N appearing in the symmetric current.
As a result, we need to prove that N j s z,z+ei can be replaced by a discrete gradient, which is also a function of the empirical measure.This replacement is much more complicated than for gradient systems, for which local equilibrium (64) is enough, because instead of replacing a function of order 1 by its average, we need to do so for a function of order N , meaning that the correction term in (64) no longer vanishes.In other words, non-gradient local equilibrium states are distorted so that we cannot replace j s by its average under the grand-canonical distribution, and lower order correlations to local equilibrium need to be taken into account in (64).Furthermore, although the asymmetric current j a is expected to be replaced by its local equilibrium average, lower-order corrections to local equilibrium (64) need to be considered.This is a manifestation of Einstein's relation (cf.Remark 5.1 above).xIn particular, the asymmetric contribution in the hydrodynamic limit appearing in ( 13) is not simply the gradient of E ⋆ Λ,f dθ (j a 0,ei ).More precisely, we have the following result: Proposition 5.1 (Non-gradient Replacement Lemma).We introduce as well as where x = z/N ∈ T 2 .Then for any H ∈ C 1 ([0, T ] × T 2 × S), and i ∈ {1, 2}, lim sup ε→0 lim sup Note that both Y ε,N z,i (η, dθ) and Z ε,N z,i (η, dθ) above are differences between microscopic currents and their mesoscopic averages, and that the latter are functions of the empirical measure according to (60) and (61).Further note that we recovered in Y ε,N z,i (η, dθ) microscopic gradients which will ultimately balance out the last factor N by summation by parts in (53).The proof of Proposition 5.1 is very technical, we do not reproduce it here and refer the interested reader to [27, p. 128, Cor.6.7.3 and Lem.6.7.4] for a detailed implementation in the same setting.

Proof of the hydrodynamic limit
We are now in a position to close Equation (53) and take the limit N → ∞.Recall that according to ( 52) and ( 53), vanishes in probability as N → ∞.Thanks to Proposition 5.1, the total current N j s z,z+ei + j a z,z+ei in the last term can be replaced by where x = z/N ∈ T 2 .To consider the limit N → ∞, consider the distribution Q N of the empirical measure's trajectory (π N t ) t≤T , defined in (49).We have the following result.Proposition 5.2 (Relative compactness and regularity of the density).The sequence (Q N ) N ∈N is relatively compact, and any of its limit points Q * is concentrated on trajectories that are 1.absolutely continuous with respect to the Lebesgue measure on T 2 × S, i.e. such that there exists a function f such that π t (dx, dθ) = f t (x, θ)dxdθ.

such that the local density
Note that this result says nothing, a priori, about the spatial regularity of f itself, so that in the limit N → ∞, there is no reason for ρεN x+ei/N − ρεN x to converge to a welldefined quantity (in this case, ∂ xi f (x, θ)dθ.This is not a problem however, because ρ admits spatial derivatives, and the self diffusion coefficient d s is C ∞ [31], so that the discrete gradient ρεN x+ei/N − ρεN x can be transferred onto d s (ρ εN x ) and to the test function, as in Definition 5.1.
In the limit N → ∞, according to the previous proposition, for z = N x we can replace ρ εN x (η(t)) by 1 ε 2 ⟨π N , 1 [x−ε,x+ε] ⟩, which in turn can be replaced in the limit ε → 0 by ρ t (x).Similarly, ρεN x (η(t), dθ) can be replaced as N → ∞ then ε → 0 by f t (x, θ)dθ.Finally, spatial averages can be replaced by integrals and discrete spatial derivatives by their continuous counterparts.This, together with the fact that (81) vanishes in probability and by replacement (82), proves that any limit point Q * of (Q N ) is concentrated on trajectories π satisfying i) and ii) above, and such that ⟨f t , H t ⟩ − ⟨f 0 , H 0 ⟩ + ⟨f t , [∂ vanishes in Q * -probability.Note we abuse our previous notation, in that we also denote by ⟨•, •⟩ the inner product in L 2 (T 2 × S).This proves Theorem 5.1.Remark 5.3.In Proposition 5.2 above, we slightly expanded on results in [27].In [27] it was not shown that π t (dx, dθ) = f t (x, θ)dxdθ, but rather that there exists for any t, x an orientation distribution ft (x, dθ) on S such that π t (dx, dθ) = ft (x, dθ)dx.
Here we chose, to make it more readable, to write statement i) of Proposition 5.2 this way, assuming that π t is, at any x, absolutely continuous with respect to the Lebesgue measure on S.However, this is not a problem because our orientation dynamics are given by angular diffusion, and the absolute continuity of π t on S can be straightforwardly shown in two ways.
1. Either by showing directly on the microscopic system that, since particle's angles diffuse independently, the distribution of angles ρεN x in a mesoscopic box of size εN must be, for any positive time, absolutely continuous w.r.t. the Lebesgue measure for any t, x. 2. Or by proving the statement at a macroscopic level, by first showing a weakened hydrodynamic limit taking the form ft (x, dθ) in the spirit of [27], and then using the maximum principle on the orientation diffusion to show that any such weak solution ft (x, dθ) is actually continuous with respect to the Lebesgue measure and can therefore be written ft (x, dθ) = f t (x, θ)dθ Note that this was not true for the ALG studied in [27] because the dynamics considered there was a jump dynamics with alignment, which had no such instantaneous mixing properties on particles' orientations.The fact that ft (x, dθ) is absolutely continuous with respect to the Lebesgue measure on S is, therefore, a direct consequence of the angular diffusion dynamics we chose for particle's orientations in our ALG.

Discussion
In this paper, we considered a lattice model for active matter and derived its exact hydrodynamic limit.This model exhibits a phenomenon of motility-induced phase separation (MIPS), which results from the interplay between self-propulsion and particle crowding.When self-propelled particles form clusters, they tend to slow down due to collisions and interactions.This creates a positive feedback loop that enhances the clustering and leads to a phase separation between a dense (liquid-like) phase and a dilute (gas-like) phase.
We obtain an exact macroscopic description of the active lattice gas model in terms of a PDE and use it to characterise the large-scale behaviours that emerge from the microscopic dynamics.Note that the hydrodynamic limit derived in Section 5 is exact for any finite time in the limit N → ∞.However, as the numerical simulations show, large but finite-sized systems are subject to fluctuations around the hydrodynamic limit.Macroscopic fluctuation theory (or Large deviation theory) is required to account for such effects [24,50,51].These terms are needed to understand metastable solutions, but their rigorous derivation remains an open problem.
Our analysis of MIPS consists of two approaches: linear stability analysis of the PDE model and numerical simulations with perturbed initial conditions.The latter approach allows us to observe the emergence and evolution of different patterns from the homogeneous state.Unlike previous studies [12,23] that obtained analytical solutions for the stationary states, we resort to numerical methods due to the complexity of our model, which involves nonlinear coefficients and arbitrary orientations.
While most previous work on hydrodynamic limits focuses on gradient systems [23,56], we do expect non-gradient systems to be ubiquitous, for example, in multicomponent systems with simple exclusion [34,57] but also elsewhere, such as in systems with general exclusion [58].Among such systems, the one we have here is well-chosen because we can handle non-explicit coefficients stemming from the non-gradient aspect via our knowledge of the self-diffusion coefficient d s (ρ).
Our results demonstrate the complexity and richness of active matter systems and highlight several promising avenues for future research.For example, the phenomenology of MIPS usually includes metastable behaviour between the spinodal and binodal [12], in which case (13) will support both homogeneous and inhomogeneous steady solutions and transitions between these states would require an analysis of large deviations.Phase separation occurs in the limit of fast rotational dynamics that cause the interfacial widths to shrink to zero.This limit leads to a macroscopic description in terms of the density ρ alone (instead of the orientation-dependence density f ): it would be useful to have a rigorous mathematical formulation of this limit.

)Fig. 1
Fig.1Illustration of the linear stability boundary(24) for the homogeneous solution of (13) (solid line), compared with numerical solutions of the equation.Results are shown for the illustrative parameter value ℓ = 0.5 The coloured points indicate whether the numerical results indicate stable or unstable behaviour.The comparison between analytical and numerical results is discussed in the main text.

φ 5 Fig. 3
Fig. 3 Example of two-dimensional patterns of the macroscopic model (13) corresponding to ϕ = 0.5, Pe = 30.0and ϕ = 0.7, Pe = 12.0 for D E = 1.0, ℓ = 0.5, Nx 1 = Nx 2 = 128, N θ = 64 starting from a random uniform perturbation of the uniform steady state.Time advances from left to right.The first and third rows correspond to the density ρ(x, t), and the second and fourth rows show stream plots of the polarisation p(x, t).

) 5 Fig. 5
Fig. 5 Snapshots of an output of the microscopic model at T = 0.0, 0.1, 1.0 The particles are represented by triangles of size 1N pointing towards e θx with colours corresponding to their orientations.Each cell is also shaded in proportion to the local particle density.