Generalized thermodynamics of Motility-Induced Phase Separation: Phase equilibria, Laplace pressure, and change of ensembles

Motility-induced phase separation (MIPS) leads to cohesive active matter in the absence of cohesive forces. We present, extend and illustrate a recent generalized thermodynamic formalism which accounts for its binodal curve. Using this formalism, we identify both a generalized surface tension, that controls finite-size corrections to coexisting densities, and generalized forces, that can be used to construct new thermodynamic ensembles. Our framework is based on a nonequilibrium generalization of the Cahn-Hilliard equation and we discuss its application to active particles interacting either via quorum-sensing interactions or directly through pairwise forces.

we can construct a hydrodynamic description that fits within the framework of Section 1. The latter can then be used to predict quantitatively the phase diagram of QSAPs and its finite-size corrections.
In Section 3 we then turn to active particles with constant propulsion forces interacting via isotropic, repulsive pairwise forces (pairwise force active particles, or PFAPs) [3][4][5][6]. For these models, the slowdown triggering MIPS is due to collisions. Contrary to QSAPs, there is no method in the literature allowing to map the hydrodynamics of PFAPs onto the general framework of Section 1. Nevertheless, we show that we can still account for the phase equilibria of PFAPs following the ideas presented in Section 1.
Finally, we show in Section 4 how the generalized thermodynamic variables identified using our formalism play the role of generalized forces when changing ensembles. In particular, we show that using an externally imposed mechanical pressure, i.e., considering an isobaric ensemble, only leads to a Gibbs phase rule when mechanical and generalized pressures coincide.

General framework
We consider a continuum description of non-aligning active particles with isotropic interactions. The vectorial degrees of freedom corresponding to the particle orientations are then fast degrees of freedom and do not enter a hydrodynamic description. The sole hydrodynamic field is thus the conserved density ρ(r, t), obeyingρ = −∇·J. By symmetry, the current J vanishes in homogeneous phases. Its expansion in gradients of the density involves only odd terms under space reversal. At third order, we use: g[ρ] = g 0 (ρ) + g 1 [ρ] where g 1 = λ(ρ)(∇ρ) 2 − κ(ρ)∆ρ.
Note that for general κ(ρ) and λ(ρ), g[ρ] cannot be written as the derivative of a free energy. Eq. (1) is perhaps the simplest generalization of the Cahn-Hilliard equation out of equilibrium and has been argued to be relevant for the phase separation of active particles in the past [1,6,11,26,29]. For a non-constant M [ρ], it allows for circulating currents with non-zero curls. A generic third order expansion is formally equivalent to (1), at this order in the gradient expansion, using for instance M = 1 + ( β α − λ α )(∇ρ) 2 + ( ζ α + κ α )∆ρ and g 0 such that g 0 (ρ) = α(ρ), where the prime denotes a derivative with respect to ρ. Such choices, however, can lead to a change of sign or a divergence of M so that, in what follows, we restrict ourselves to dynamics of the form (1) with positive definite M . Such a restriction does not matter when considering fully-phase separated profiles in the macroscopic limit Maxwell equal-area construction on the pressure. In both panels, we used a double-well potential for illustrative purpose.
but was recently proved important when describing curved interfaces [30] where generic currents of the form (2) may lead to a richer phenomenology than that of Eq. (1). The spinodal region of a phase-separating system can easily be predicted from Eq. (1). A homogeneous profile of density ρ 0 is indeed linearly unstable whenever g 0 (ρ 0 ) < 0 and the sign of g 0 (ρ 0 ) hence defines the spinodal region.

Warm-up exercise: the equilibrium limit
Before deriving the binodal curve predicted by Eq. (1) in its most general form, it is illuminating to first review the corresponding equilibrium limit, i.e., the standard Cahn-Hilliard equation [27,28] which corresponds to 2λ + κ = 0 [26]. In this case, the dynamics (1) corresponds to a steepest descent in a free energy landscape F[ρ]: (1) is then the chemical potential, defined as the functional derivative of F with respect to ρ: where g 0 (ρ) = f (ρ) and The free energy functional F is extensive so that, in a macroscopic phase-separated system, the contribution of the interfaces is sub-dominant. The term 1 2 κ(ρ)(∇ρ) 2 in F can then be neglected and the phase equilibria can be determined from the bulk free energy density f (ρ): The coexisting densities ρ g and ρ in the gas and liquid phases are the one minimizing the free energy under the constraint that the average density ρ 0 is fixed. They are obtained through a common tangent construction on f (ρ) or, equivalently, as the densities satisfying the equalities of chemical potential f (ρ g ) = f (ρ ) =μ and pressure P (ρ g ) = P (ρ ) =P , with the pressure P defined as P (ρ) = ρf (ρ) − f (ρ). Alternatively, the coexisting densities can be constructed using a Maxwell equal area construction νg ν P (ν) −P dν = 0 (6) where ν ≡ 1/ρ is the volume per particle, ν g/ ≡ 1/ρ g/ . The two thermodynamic constructions are illustrated in Fig. 1. Note that, instead of relying on a free energy, the equality of pressures and chemical potentials between coexisting phases can also be derived directly from the dynamics (1). First the vanishing of the flux J in Eq. (1) immediately imposes a uniform chemical potential g, which is thus equal between coexisting phases: g 0 (ρ ) = g 0 (ρ g ). To derive the equality of pressures, we rewrite Eq. (1) asρ where σ is the stress tensor, whose expression in Cartesian coordinates is Note that, similar to g, σ is related to the free-energy functional through [31]: In fully phase-separated, flux-free steady states, one can get the equality of pressure between coexisting homogeneous phases from Eqs. (7)- (8). For finite systems, Eq. (7) can also be used to derive the finite-size corrections to the binodals due to Laplace pressure [28].

Generalized thermodynamic variables
For generic functions λ(ρ) and κ(ρ), which do not satisfy 2λ(ρ)+κ (ρ) = 0, the free energy structure breaks down because the gradient terms in g cannot be written as a functional derivative: A common tangent construction on a free energy density defined through f (ρ) = g 0 (ρ) then does not lead to the correct coexisting densities [11,26]. However, as we show below, g can be written as the functional derivative of a generalized free energy G with respect to a non-trivial new variable R, which depends on the functional forms of κ and λ. Although the dynamics (1) are a priori out of equilibrium, its steady states correspond to extrema of this generalized free energy and, as we show below, we recover the full structure of the equilibrium case described above. We now derive this mapping and show in Section 1.4 how it can be used to compute the binodals of Eq. (1) exactly. Finally, we turn to their finite-size corrections in Section 1.5.
To proceed, we consider the one-to-one mapping R(ρ) defined by where the derivatives are taken with respect to ρ. Direct inspection shows that g can now be written as a functional derivative with respect to R [26]: where we have defined a generalized free energy density φ(R) such that The dynamics of ρ is now written as the derivative of a generalized free energy functional: Note, however, that the structure of (15) differs from the equilibrium case (4) since the functional derivative is taken with respect to R instead of ρ. Nevertheless, the steady-state solutions of (15) correspond to extrema of G with respect to R ‡. Comparing Eq. (15) to the equilibrium case (3), we note that the former can be seen as driven by gradients of a generalized chemical potential g = δG δR . Similarly to the equilibrium case, we now show that the dynamics (15) can also be written so as to appear driven by the divergence of a generalized stress tensor. Specifically, the current J can be rewritten as with a tensor σ reading in Cartesian coordinates where we have defined § Once again, the generalized stress tensor can be deduced from the generalized free energy through In the following, we identify the diagonal coefficients of σ, the normal stresses, with generalized (potentially anisotropic) pressures. Again, we split h = −σ xx into a local function and an interfacial contribution: (20) ‡ The dynamics of R itself can be easily deduced asṘ = R ∇ · [M ∇ δG δR ]. Note that, in particular, R is not a conserved quantity. § Alternatively, h 0 can be obtained through h 0 = ρ R(ρ)g 0 (ρ)dρ, or, introducing υ = 1/R, through h 0 = − d(φυ) dυ .

GAS LIQUID G L
x ρ ρ g ρ x g x Figure 2. Schematic representation of the mean density field of a fully phase-separated system in 2d. We consider the density profile connecting gas and liquid phases along a horizontal cut so that the interface is oriented alongŷ (center). In the macroscopic limit, the interface is locally flat in the transverse directionŷ (left) and the problem simplifies into an effectively 1d domain wall computation for the density profile (right).
We emphasize here that σ and h need not have any connection to mechanics and momentum transfer.
Finally, we stress that the equilibrium case is easily recovered using 2λ + κ = 0: Eq. (11) then implies that R = ρ (up to multiplicative and additive constants that play no role in phase equilibria and can thus be discarded). All our generalized quantities then reduce to their equilibrium counterparts.
Before we turn to the derivation of the binodal curve, we first note that The spinodal region, defined as g 0 (ρ) < 0, thus corresponds to the region in which the generalized free energy density is concave, d 2 φ dR 2 < 0, provided R is chosen positive. Furthermore, from Eq. (18) one finds that h 0 (ρ) = Rg 0 (ρ) (22) so that the spinodal region can equivalently be defined from h 0 (ρ) < 0. Finally, we note that, contrary to the generalized free energy density φ which depends on λ and κ through R, the spinodal region is unaffected by the gradient terms in g. We now show how the above results directly yield the binodal curve of our generalized Cahn-Hilliard equation by considering fully phase-separated systems. We then discuss in Section 1.5 the corrections to the binodal curve for finite-size systems.

Phase coexistence in the large system size limit
A macroscopic droplet of, say, the dense phase has an infinite radius of curvature in the large system size limit, so that curvature effects are negligible. As in equilibrium, computing the coexisting densities reduces to studying a one-dimensional domain-wall profile perpendicular to the interface [28], whatever the original number of spatial dimensions. To do so, we consider a flat interface, orthogonal tox, between coexisting gas and liquid phases at densities ρ g and ρ (see Fig. 2).
For such a profile, any derivative with respect to a direction normal tox vanishes so that Eq. (16) directly implies that g and σ xx are constant. For coexisting homogeneous phases, this leads directly to and where R ,g ≡ R(ρ ,g ). These two constraints thus fully determine the coexisting densities and are equivalent to a common tangent construction on φ(R) since g 0 = dφ dR and h 0 = R dφ dR − φ. In stark contrast to equilibrium liquid-gas phase separation, the interfacial terms g 1 or h 1 affect the coexisting densities through the definition of R, Eq. (11), which depends on λ(ρ) and κ(ρ). Note that the sole knowledge of the dynamics in Eq. (1) allows us to determine the coexisting densities using the constructions above, without the need to solve for the full density profile at the interface. The common-tangent construction on φ leads to coexisting densities which are independent of the mean density ρ 0 . The lever rule for determining the phase volumes V and V g therefore still applies: ρ V + ρ g V g = ρ 0 V . Note that this lever rule applies to ρ and not to R since the latter is not a conserved quantity.
The Maxwell construction. As in equilibrium, the common tangent construction on φ is equivalent to a Maxwell construction on h 0 . We now derive the latter because it will be useful when considering PFAPs, and also since it provides a simpler numerical route to computing the binodal curve from the expression of h.
As we shall do for PFAPs, we start from a current given by Eq. (16) so that the flux free condition in a situation as depicted in Fig. 2 implies that the generalized pressure is constant, recalling that the curvature of the interface is negligible: Then, we introduce the generalized volume per particle and compute the integral where the spatial integral is computed along the direction normal to the interface. After some algebra, h 1 can be rewritten as This allows us, after some algebra, to show that h 1 ∂ x υ is a total derivative In turn, this leads to a generalized Maxwell construction on h 0 :

Finite size effects
Let us now consider what happens if one takes into account the finite curvature of the phaseseparated domains. Again, thanks to our mapping, the derivation below resembles closely the one done in equilibrium for the Cahn-Hilliard equation [28]. We consider a radial cut along the interface of a circular domain in 2D, as in Fig 2. By symmetry, the current J vanishes in steady state. Eq. (16) then immediately gives ∇g = 0 so that one still has an equality of generalized chemical potentials between the two phases: g 0 (ρ g ) = g 0 (ρ ). On the other hand, ∇ · σ = 0 does not lead to a uniform σ xx in this circular geometry, which highlights the different behaviors of the generalized chemical potential and the generalized stress tensor for finite systems.
To proceed, we integrate the radial component of ∇ · σ along the path depicted in Fig. 2. To highlight the spherical geometry, we parametrize this path as rr. Using the expression for the divergence of a tensor in spherical coordinates (polar in 2D) leads to Using the expression (17) of σ in this geometry then leads to where we have used that the isotropic terms in σ cancel and derivatives with respect to θ vanish by symmetry. When the width of the interface is small compared to the droplet radius r d , expanding r around r d and using that (∂ r ρ) 2 vanishes outside the interface leads to where we have introduced a generalized surface tension γ: Note that, as for h 0 and σ, γ need not have any mechanical interpretation for generic phaseseparating active matter systems. To leading order in 1/r d , γ can be computed across a flat interface (using a slab geometry as in Fig. 2). For an interface perpendicular to the x-axis, it then reads Finally, let us comment on the sign of γ which has recently attracted interest since it has been measured negative for PFAPs [32] (see Section 3.4 for a discussion of that case). Here, since κ need to be positive for stability reasons, we see from Eq. (34) that γ has the sign of R . Starting from the dynamics in Eq. (1), the sign of R is arbitrary, corresponding to an integration constant when solving Eq. (11). The generalized surface tension γ can then be either positive or negative, although with different expressions for the generalized pressure h 0 (ρ). On the other hand, starting from an expression for the stress tensor in Eq. (17), as will be the case for PFAPs in Section 3,  R is fixed by the expression for σ and can take either sign. Our framework thus supports both positive and negative γ.

Illustration of our general framework for a scalar active matter model
In this section we show on a particular example that our generalized thermodynamic construction predicts exactly the phase equilibria of our nonequilibrium Cahn-Hilliard equation (1) through Eq. (23). To this end, we numerically integrate this equation in 2d for the particular (and rather arbitrary) choice: To check the theory, we first numerically solve Eqs. (1) and (35) using a semi-spectral integration scheme (linear terms are computed in Fourier space, non-linear terms in real space) with Euler time stepping. For each value of r, we start from a phase-separated state with two arbitrarily chosen densities (ρ g = 1 and ρ = 5) and measure the coexisting densities once the system has relaxed to its steady state.
To compare with the theoretical predictions for the binodals, we first determine the function R(ρ) using Eq. (11), which (up to two unimportant integration constants) gives R(ρ) = log ρ. We then use either the common tangent construction on φ(R) or the Maxwell construction on h 0 (υ), shown in Fig. 3(b,c), as described in the previous section. The comparison with the coexistence densities measured in the simulations of Eq. (1) is shown in Fig. 3(a): the difference between theory and simulations is found to be smaller than 0.5% for every point, thus confirming that the dynamics does indeed yield the stationary state analyzed in Section 1.
These numerical results are obtained in systems where a straight band of liquid coexists with a dilute gas phase so that finite-size curvature effects are negligible. On the contrary, when finitesize liquid droplets coexist with a gaseous background, the coexisting densities differ from those  predicted by Eq. (23) due to the finite-size corrections discussed in Section 1.5. In this case, a jump of the generalized pressure through the interface is indeed measured numerically, and found to be given quantitatively by the generalized surface tension (33) (See Fig. 4a). Similarly, there are density shifts in each of the phases which scale as 1/r d as shown in Fig. 4b.

QSAPs.
We now turn to a microscopic model for QSAPs, for which we derive a hydrodynamic description and compare the predictions of our formalism with direct numerical simulations of the microscopic model. We consider particles labeled by i = 1 . . . N , moving at speed v along body-fixed directions u i which undergo both continuous rotational diffusion with diffusivity D r and complete randomization with tumbling rate α. The equations of motion are given by the Langevin dynamicṡ where η and ξ are delta-correlated Gaussian white noises of appropriate dimensionality. In addition to continuous angular diffusion, we have included in (36) a non-Gaussian noise accounting for tumbling events: the t i are Poisson distributed with a rate α and the δθ j 's are drawn from a uniform distribution between 0 and 2π. Each particle adapts its speed, v[ρ(r i + εu i )], to a local measurement of the density: with K(r) an isotropic coarse-graining kernel, andρ(r) = i δ(r − r i ) the microscopic particle density. Note that the local density is measured with an offset εu i which allows for anisotropic quorum sensing. This effect, which does not create alignment interactions, captures a slowdown of particles that would arise, for instance, due to a large density of particles in front of them. This can thus model, say, a visual quorum-sensing or steric hindrance. In a different context, anisotropic sensing has been shown to lead to a rich phenomenology for aligning active particles [33].

Hydrodynamic description of QSAPs.
Deriving hydrodynamic equations from microscopics is generally difficult, even in equilibrium [34]. For QSAPs we can follow the path of [1,24,35], taking a mean-field approximation of their fluctuating hydrodynamics. We first assume a smooth density field and a short-range anisotropy so that the velocity can be expanded as where ρ is evaluated at r i and 2 = r 2 K(r)dr. Following [24,35], the fluctuating hydrodynamics of QSAPs, derived in Appendix A, is then given by: with Λ a unit Gaussian white noise vector and The mean-field hydrodynamic equation of QSAPs is then Eq. (1) with the coefficients in Eq. (40). As mentioned earlier, the spinodal region is defined from the criterion g 0 (ρ) < 0, which leads here to a modification of the standard linear instability criterion for QSAPs [1]: To construct the phase diagram for a given choice of v(ρ), using the generalized thermodynamic procedure, we first solve for R(ρ) using Eq. (11) and from it obtain both φ(R) and h 0 (R). The binodals then follow via a common-tangent construction on φ(R) or, equivalently, by setting equal values of h 0 and g 0 in coexisting phases. Note that since 2λ + κ = 0 one has R = ρ. The phase diagram thus cannot be found by globally minimizing a free energy density f (ρ) defined from f (ρ) = g 0 (ρ) as discussed before [1,24]. Indeed, such a procedure correctly captures the equality of g 0 in both phases but predicts a common tangent construction on f which is violated. We now turn to describe the numerical simulations of microscopic models of QSAPs.

Comparison between theory and numerics.
In what follows we study models where the densityρ is computed according to Eq. (37) with the bell-shaped kernel  Here Θ is the Heaviside function, Z a normalization constant and we used an interaction radius of r 0 = 1. In addition we take the velocity to be This interpolates smoothly between a high velocity v 0 at low density (ρ ρ m ) and a low velocity v 1 at high density (ρ ρ m ). In addition to the 2d continuous space model described above, we also conducted simulations of QSAPs in 1d on lattice [2]. In this case, we consider run-and-tumble particles (RTPs): particle i has a direction of motion u i = ±1 which is flipped at rate α/2. It then jumps on the lattice site in direction u i with rate v[ρ(x i + εu i )]. Fig. 5a shows the phase diagrams predicted by our generalized thermodynamics and those measured in QSAP simulations for a symmetric sensing of the density (ε = 0). Overall, the agreement between predicted and measured binodals is excellent, in contrast to the common tangent construction on f (ρ). It is remarkable that, for QSAPs, we can quantitatively predict the phase diagram of a microscopic model without any fitting parameters, something rare even for equilibrium models. Fig. 5b shows the binodals measured in 1d simulations on lattice with ε = 0 together with the corresponding theoretical predictions. The dependence of the binodals on the asymmetry ε is apparent in both cases. It results from the explicit dependence of g 0 (ρ) on ε established in Eq. (40). This dependence probably explains why run-and-tumble particles hopping on lattices with excluded-volume interactions [2] are not well described by the coarse-grained theory proposed so far for QSAPs which did not account for any asymmetric sensing [1]. We can see that our theoretical predictions are more accurate for small ε, as expected from the derivation of the Our theoretical predictions for the phase diagram of QSAPs rely on two different approximations. First, we use a mean-field approximation to derive the specific expression (40) for g [ρ]. For our choice of v(ρ), MIPS occurs only at large densities so that this approximation works very well except in the small and numerically unresolved Ginzburg interval close to the critical point. Second, our general theory disregards higher order gradient terms in Eq. (1). This probably explains why the hydrodynamic description works best fairly close to the critical point, where interfaces are smoothest and the gradient expansion, Eq. (38), most accurate. The quantitative limitations of our gradient expansion highlights that gradient terms directly influence the coexisting densities through Eq. (11), unlike the equilibrium case.
In addition to giving quantitative predictions for the phase diagrams, our approach sheds light on the observed universality of the MIPS in QSAPs. For example, the phase diagram does not depend on the exact shape of the kernel K, which enters Eq. (40) through 2 which then cancels in the nonlinear transform R(ρ). Similarly, Fig. 5 also shows lattice simulations of QSAPs in 1d where complete phase separation is replaced by alternating domains (with densities given by the predicted binodal values). This confirms the equivalence of continuous (ABP) and discrete (RTP) angular relaxation dynamics for QSAPs [24,35]. Our results, however, also expose sensitivity to other microscopic parameters such as the fore-aft asymmetry ε which enters g 0 and therefore affects the binodals. This might explain the different collective behaviors seen in swarms of robots that adapt their speeds to the density sampled in either the forward or the backward direction [36].

Finite-size corrections.
Similar to the scalar active model of Section 1.6, we expect finite size corrections when a liquid droplet is formed in a finite system: for a droplet of radius r d , we expect to leading order in the droplet radius an effective pressure jump across the interface (32): Accordingly, one expect the finite-size corrections to the co-existing densities to decay as ∝ 1/r d . In Fig. 6a, we show that the measured binodals indeed converge towards their asymptotic values in a manner consistent with a 1/r d decay.
A quantitative check of (44) is difficult since, first, our derivation of h 0 (R(ρ)) is based on a number of approximations, and, second, our numerical measurements of the binodals are necessarily noisy. To proceed, we measure ρ g (r d ) and ρ (r d ) and construct ∆h 0 (r d ) = h 0 (ρ (r d )) − h 0 (ρ g (r d )). ∆h 0 (r d ) does not vanish exactly in the large system size limit, nor in a slab geometry in which we measure a small correction ∆h slab 0 /h 0 (ρ ) ≈ 0.1%. This systematic error can stem from several origins, from the gradient expansion to the mean-field approximation, through limitations in the numerical accuracy of the density measurement. Though very small, this error becomes comparable to the Laplace pressure jump for radii r d 20, highlighting the numerical challenges in measuring these finite-size effects. Nevertheless we show in Fig. 6b that ∆h 0 (r d ) − ∆h slab 0 converges to its asymptotic value consistently with a 1/r d decay. Furthermore, the prediction of Eq. (44) can be checked by measuring the prefactor γ ≡ r rg κ R (∂ r R) 2 dr of this decay in a slab geometry. The corresponding prediction is shown as a dashed line in Fig. 6b and agrees semi-quantitatively with our numerical results, without any fitting parameters.
To understand why we observe a quasi-quantitative agreement despite relatively small values of r d , it is useful to explicitly expand Eq. (31) as The first order correction to ∆h 0 = γ r d is thus given by Using that, from the definition (11), (R κ) = −2λR , the prefactor R κ can then be expanded around r = r d as Since λ = 0 for QSAPs, we are left with which is zero for a symmetric interface so that in that case For our choice of v(ρ), the density profile is indeed very close to a hyperbolic tangent (data not shown) and the lack of first order corrections for such profiles probably explains the semiquantitative agreement of our numerical results with the 1/r d behaviour.

PFAPs.
We now consider the case of self-propelled particles interacting via pairwise forces, which has attracted considerable interest over the past few years [3,4,8,17,25,29,37,38]. We define the model in Section 3.1 and construct its hydrodynamic description in Section 3.2. Contrary to QSAPs, there is no available method to derive accurate estimates of the coefficients λ(ρ) and κ(ρ) or to rule out the existence of other terms [30]. We discuss in Section 3.3 how we can nevertheless follow the path laid out using our generalized thermodynamic formalism to understand how coexistence densities are selected in PFAPs. Finally, finite-size effects are considered in Section 3.4.

Model.
We consider N self-propelled particles in two dimensions interacting via the repulsive, pairwise additive, Weeks-Chandler-Andersen potential: with an upper cut-off at r = 2 1/6 σ, beyond which V = 0. Here σ defines the particle diameter, determines the interaction strength, and r is the center-to-center separation between two particles. Particle i evolves in two dimensions, with periodic boundary conditions, according to the Langevin equations:ṙ Here, u i = (cos θ i , sin θ i ) indicates the direction of self-propulsion and η i , ξ i are unit Gaussian white noises. For simplicity, we only include continuous rotational diffusion but our results are expected to extend to run-and-tumble dynamics since these two types of orientational noise have been shown to lead to the same phase diagram [35]. The full phenomenology of this model requires scanning a three-parameter phase diagram, parametrized for instance by the Péclet number Pe = 3v 0 /(σD r ) , the packing fraction (π/4)σ 2 ρ and the potential stiffness µ /(v 0 σ). Here, we focus on the onset of MIPS as the Péclet number and the packing fraction are varied, disregarding the role of the potential stiffness [3][4][5]39]. In practice, we fix = 1, σ = 1, v 0 = 24, µ = 1 and vary D r and ρ. MIPS then occurs at high Historically, the Péclet number was defined as Pe = v 0 σ/D t with translational diffusion D t and a Brownian rotational diffusion D r = 3D t /σ 2 . It was latter realized that in simulations of PFAPs exhibiting MIPS, the translational diffusion has a negligible effect on the phase diagram and could be set to zero. This explains the factor 3 in the definition of Pe, although a dimensionless run length l r = v 0 /(σD r ) would seem more natural. enough densities when the run-length v 0 /D r is much larger than the particle size σ, namely when Pe exceeds a threshold value Pe c ≈ 50 [3][4][5][6].

Hydrodynamic description.
Following [17,40], we start from the exact Itō-Langevin equation for the microscopic density of particlesψ(r, θ) whereÎ (θ) (r, θ) = − dr µ∇V (|r −r|)ρ(r )ψ(r, θ),ρ(r) = ψ (r, θ)dθ is the fluctuating density, and η and ξ are unit-variance Gaussian white noises of appropriate dimensionality. Denoting averages over noise realizations by angular brackets we define ρ(r) = ρ(r) , m(r) = m(r) and Q(r) = Q (r) . Herem(r) = dθuψ(r, θ), is the orientation vector,Q(r) = dθ (u : u − 1/2)ψ(r, θ), is the nematic tensor, and 1 the identity matrix. Integrating Eq. (52) over θ and averaging over noise realizations, the dynamics of ρ(r, t) reads with The dynamics of m is then obtained similarly by multiplying Eq. (52) by u and integrating over θ. This yields, with an implied summation over repeated indices, where the last term is obtained by integration by parts and we have defined We stress that, so far, Eq. (52) and (55) are exact, although they are not closed since they feature Q and the microscopic correlators in I (0) and I (1) which depend on higher moments ofψ. As a first approximation, we use that, contrary to ρ(r, t), m is a fast mode decaying at a rate D r . On time scales much larger than D −1 r , one can thus assume that m α relaxes locally to The current in Eq. (53) is then given by Interestingly, Eq. (58) can be rewritten as the divergence of a stress tensor where we have followed Irving and Kirkwood [41] (and Ref. [31] in a similar context) and write We now turn to relate these results to the formalism derived previously.
Generalized pressure and equation of state.
The resulting dynamics for ρ, with the current given by Eq. (59), should be compared to the generalized Cahn-Hilliard equation of Section 1 with the current driven by the generalized stress tensor as in Eq. (16). We see that PFAPs correspond to the special case M/R = µ, the microscopic mobility. This has important consequences for the mechanical interpretation of σ. Indeed, one can see that imposing an external potential U on the particles leads to In a flux free steady state, J = 0 and Eq. (62) becomes a force balance. Integrating (62) from a point in the bulk to infinity shows the normal component of σ to be equal to the total force per unit area exerted on a boundary. Indeed, the normal component of σ exactly coincides in homogeneous phases with the equation of state (EOS) found previously for the mechanical pressure P of PFAPs [17]. Generalized and mechanical pressure thus coincide for PFAPs and we note, following Section 1.1 where we have defined, following earlier notation [17], the "active" contribution to the pressure P A and a "direct" passive-like part P D : Note that P A is sometimes also called "swim pressure" [12], even though neglecting the pressure of the surrounding fluid to describe the phase separation of actual swimmers is problematic. The value of the pressure in a homogeneous phase of density ρ 0 is then given by where P A 0 and P D 0 are the values taken by P A and P D in homogeneous disordered phases of density ρ 0 . This allows us to identify, in analogy with Eq. (20), Note that while the structure is similar to Eq. (20), there is no gradient expansion taken hereh 1 is exact, formally containing gradients of all orders. Its expression is given by where P contains the interfacial contributions to the active and the direct pressures. The terms in m x and Q xx are purely interfacial since they vanish in the (disordered) bulk phases. We now show that the phase equilibria in PFAPs can be understood using these results with the ideas of Section 1.

Phase equilibria in PFAPs.
One way forward would be to construct an explicit gradient expansion for h in terms of ρ and obtain closed expressions for h 0 and h 1 . This would then allow us to find R(ρ) and φ(R) analytically as was done for QSAPs in Section 2. Despite the extensive literature on PFAPs, such a gradient expansion has not yet been presented, but could be accomplished for instance by using a low-density virial approximation. Such a route would possibly lead to qualitative predictions for the phase diagram, but our goal here is to show that our formalism quantitatively accounts for the phase equilibria of PFAPs, and we thus do not want to rely on such approximations. We thus proceed differently, using an approach where we instead measure the gradient terms to quantitatively verify the validity of our formalism for PFAPs.
As with the other systems, we first consider the case of a macroscopically phase separated system, for which the liquid-gas interface is locally flat and perpendicular tox. As in Section 1.4, in a flux-free steady state, h =h is constant across the interface so that the pressure is equal in coexisting phases To construct the phase diagram, we need to complement this equality by a second constraint.
Since we do not have any closed expression for the interfacial terms h 1 , we cannot use a Maxwell construction in the (h 0 , υ = R −1 ) plane as was done in Section 1.4. Instead, we measure the violation of the equilibrium Maxwell construction in the (h 0 , ν ≡ ρ −1 ) plane, schematically depicted in Fig. 7, with ν = 1/ρ the free volume per particle: Here h 0 (ν) is the pressure-volume EOS,h is the pressure of coexisting phases and ∆A = 0 directly quantifies the violation of the Maxwell construction for PFAPs [17]. Given the value of ∆A, Eqs (68) and (69) are two independent constraints satisfied by ρ and ρ g . A fully predictive theory would thus evaluate ∆A analytically and then solve (68) and (69) to obtain the values of the binodals and the coexisting pressureh. Here, instead, we use a numerical measurement of ∆A to construct the phase diagram. Although less predictive than knowing h 0 and h 1 analytically, our method clearly illustrates how the violation of the Maxwell construction, due to the role played by the interfaces, selects the binodals.

Numerical strategy and results
To numerically construct the phase diagram, we first derive an approximation to the bulk equation of state h 0 (ρ). Then, we measure h(x) numerically via Eq. (63) from which we subtract h 0 (ρ(x)) to obtain h 1 , which is integrated to obtain the numerical value of ∆A. The right hand side of Eq. (69) is then held constant at this value, and the binodals are determined as the intersect between the EOS h 0 (ν) and a horizontal line of ordinateh whose value is adjusted until it satisfies Eq. (69). Note that, for the parameter range of interest here, the two contributions to h proportional to D t are negligible and we thus discard them hereafter.
(i) We first construct an analytical approximation for the pressure h 0 (ρ 0 ) by measuring the active and direct pressures from an ABP simulation in the homogeneous region (Pe < Pe c ) using Eqs. (64). Following the route proposed in [17]: we then apply scaling arguments to extrapolate the EOS into the two-phase region Pe > Pe c . Fig. 8 explains and verifies the proposed scaling in the low-Péclet region and shows the resulting EOS for P A 0 and P D 0 . Details about the numerical procedure (refined with respect to Ref. [17]) can be found in Appendix B.3.
(ii) The next step is to numerically determine h 1 using Eq. (67) and, through it, the value of ∆A. In numerical simulations of phase-separated systems in a slab geometry (see Fig.  9a), we thus measure the profiles ρ(x), P A (x), P D (x) and Q xx (x) across the interface (see Fig. 9b-c). Using the EOS P A 0 (ρ) and P D 0 (ρ) from (i) together with the measured density profile ρ(x) we obtain the gradient contributions to the active and direct pressures as P (Fig. 9d). Together with Q xx (x), this directly provides h 1 (x) and hence the value of ∆A in Eq. (69).
(iii) Using the equation of state h 0 (ν), we now adjusth in Eq. (69) until ∆A matches the value computed in step (ii) as shown in Fig. 10a. The resultingh and the corresponding two values of ν g and ν constitute our prediction for the pressure at coexistence and the binodals.
As seen in Fig. 10b, the predicted coexistence densities match very well the measured ones. We stress again that this is not a first principle prediction, since we do not use an analytic expression for the gradient terms, which thus have to be measured numerically. Nevertheless, the excellent agreement confirms the scenario proposed in Section 1 for MIPS: unlike in equilbrium, the interfacial contributions are essential in fixing the coexistence densities. Indeed, the equilibrium Maxwell construction (equivalent to taking ∆A = 0 in Eq. (69)) clearly fails to account for the phase diagram of PFAPs, as shown in Fig. 10b. Therefore, the interfacial contributions have to be accounted for, either by defining an effective density as in Sections 1, 1.6 and 2, or by quantifying the violation of the equilibrium constructions, as demonstrated here.
We finally note that the behavior of the interfacial terms P A 1 and P D 1 in Fig. 9d can be qualitatively understood as arising from the polarization of the gas-liquid interface. Since a particle at the interface is on average oriented towards the (denser) liquid phase, i.e., up the density gradient, it experiences a more efficient collisional slow-down than it would in an isotropic environment at the same local density. Since P A [ρ] is proportional to the effective swim-speed v, this will yield a lower P A than in the isotropic phase, and thus a negative P A 1 . Conversely, as P D (c) Profiles of the total pressure h(x) and its three non-negligible components P A , P D and v 2 0 Q xx /µD r (solid lines). The dashed lines correspond to the local contributions P A 0 (ρ(x)) and P D 0 (ρ(x)) that are predicted by the equation of state for a homogeneous system at density ρ(x). (d) The interfacial contributions to the pressure, entering h 1 in Eq. (67).
is proportional to the amount of repulsive particle contacts experienced by the particle, the same argument will lead to a positive interfacial contribution P D 1 to the direct pressure, confirming the observations in Fig. 9d. Since these two terms give the dominant contributions to h 1 , we thus conclude that, at the microscopic level, the phase coexistence densities in PFAPs is controlled by the polar ordering of particles at the gas-liquid interface.

Finite-size corrections.
We now turn to study the finite-size corrections to the phase equilibria of PFAPs. As previously, we consider a circular droplet of radius r d (see Fig. 2). Following Section 1.5, the pressure jump across the interface is given at leading order in 1/r d by where the surface tension γ is measured across a planar interface perpendicular tox. We follow the same route as for QSAPs and measure independently ∆h 0 and γ in numerical simulations to characterize the finite-size corrections to the phase coexistence.
To understand the different contributions to γ, we introduce the difference between the xx and yy components for each term in the stress tensor (60) (recall that the terms proportional to D t are negligible): The surface tension is then given by These three contributions and their sum, σ yy − σ xx , are plotted in Fig. 11, as measured across a flat (on average) interface; the resulting integral yields the estimate γ ≈ −140 in the units of our simulations. Interestingly, the contribution of the direct term δP D is completely negligible, in contrast to the equilibrium case in which the phase separation is due to attractive forces, which also determine the surface tension. Here, the main contributions stem from the anisotropy of the active pressure in the interface, as well as from the anisotropic nematic order of the particles in the interfacial region. We furthermore note that the resulting value of the surface tension is negative, which confirms the finding of Ref. [32] and can be rationalized following [42] by considering the escape angle of an active particle exiting a curved interface. We now evaluate the effective Laplace pressure ∆h 0 for curved droplets of different radii. Although this quantity is in principle directly measurable in simulations, it is numerically challenging due to the large fluctuations in the local pressure. We thus instead proceed similarly as for QSAPs, by first accurately measuring the coexisting densities in finite systems in which a liquid droplet of radius r d coexists with a vapor background. These are shown in Fig. 12a, showing that the liquid phase is effectively depleted for finite r d , hence confirming the heuristic argument given in [42]. The correction to the coexistence densities is again found to be compatible with a 1/r d decay. The pressure jump can then be computed using the equation of state and the measured densities as ∆h 0 = h 0 (ρ ) − h 0 (ρ g ), shown in Fig. 12b. To extract the leading order behavior in 1/r d , we fit ∆h 0 (r d ) with two parameters, using a function c 1 /r d + c 2 /r 2 d . The second-order term is necessary because the width of the liquid-vapor interface is large (≈ 40, see Fig. 9b) so that the assumption of large r d does not hold. The leading-order coefficient from the fit in Fig. 12 corresponds to a value of γ ≈ −230, to be compared to γ = −140 measured across the straight interface in Fig. 11. The sign and order of magnitude are thus correctly captured, in spite of the many approximations and numerical difficulties inherent in these measurements.
We stress that the procedure we detail above retains all the gradient terms entering σ through h 1 , and hence accounts for the negative value of γ. As explained before, we have not, however, carried out explicitly a gradient expansion of σ. Therefore, we do not know whether PFAPs can be quantitatively described by Eq. (16) and (17). We have shown that, in the formalism of Section 1, phase-separated solutions are compatible with a negative γ. However, the finite size corrections derived in Section 1.5 are constrained by the equality of generalized chemical potential in the two phases which imposes that the density correction ρ − ρ ∞ take the same sign in the two phases. This is at odds with the observation of Fig. 12(a), thereby suggesting that PFAPs are not fully described by our generalized Cahn-Hilliard equation. A promising suggestion is that the finite size effects of PFAPs are best described by a more general gradient expansion which would imply the analogue of a Laplace pressure jump for the chemical potential [30]. Figure 12. (a) Coexisting densities measured numerically as a function of the droplet radius for Pe = 100, normalized with the corresponding densities ρ ∞ in the slab geometry (i.e., r d = ∞). The solid lines indicate fits to the measured data using the function ρ/ρ ∞ = 1+c 1 /r d +c 2 /r 2 d , with c 1 and c 2 fitting parameters. Note that these measurements are very sensitive to the definition of coexisting densities, e.g. using the positions of the peaks of maximum probability in the distribution P (ρ) of local density ρ vs using the average of such peaks, so that we can expect at best semiquantitative agreement with theory. The droplet radii r d are estimated from the phase volumes obtained from the integral of the respective peaks in P (ρ). (b) The corresponding difference in coexistence pressure ∆h 0 , obtained from the densities in (a) using the numerical EOS, and normalized with the pressure h ∞ for a flat interface. Solid lines show fits to ∆h 0 /h ∞ = c 1 /r d + c 2 /r 2 d , where the fitting parameter c 1 h ∞ ≈ −230 is an estimation for the surface tension γ.

Change of ensembles
One powerful aspect of equilibrium thermodynamics is that it relates the physical states of a system under different environmental constraints. Beyond its engineering value, the existence of several ensembles provides useful theoretical tools to study phase transitions [43]. Similar developments for non-equilibrium systems have however proven difficult [44][45][46]. Interestingly, our formalism allows some progress.
We adapt our previous constant volume (isochoric) simulations to consider an isobaric (constant pressure) ensemble. PFAPs or QSAPs are now confined by mobile harmonic walls, Figure 13. PFAPs in the isobaric N, P w , Pe ensemble. (a): Snapshots from PFAP simulations with a mobile wall imposing a pressure P w at Pe = 100 during a slow upwards (left) and downwards (right) pressure ramp (for movies, see [47]). In the isobaric ensemble, the phase transition becomes discontinuous, in contrast to the phase coexistence observed in constant-volume simulations. subject to a constant force density P w which imposes a mechanical pressure P = P w (see Fig. 13a and movies in [47]). Since P = h 0 is a generalized thermodynamic variable for PFAPs, we expect, as in equilibrium, that the coexistence region of the isochoric (N, V, Pe) ensemble collapses onto a coexistence line in the isobaric (N, P, Pe) case, corresponding to the pressure at coexistence in the isochoric ensemble (see Fig 13b). Inposed-pressure loops carried out by slowly ramping up and down P w then lead to small hysteresis loops around the value of P w corresponding to coexistence. These loops would vanish in the large system size limit for quasi-static ramping of P w (see Fig 13c).
In contrast, for QSAPs the mechanical pressure P is unrelated to either of the generalized variables g 0 , h 0 . The same value of P w may thus lead to different states of the system depending on its history: the Gibbs phase rule does not apply for QSAPs in this ensemble. This translates into large hysteresis loops when slowly cycling P w , as shown in Fig 14. On a fundamental level, the different relationship between thermodynamical and mechanical observables can be related to the presence or absence of an effective momentum conservation in the steady state [48]. From a more practical point of view, this can be traced back to the fact that adding an external potential U to PFAPs gives a simple force balance equation in a flux-free  Figure 14. For QSAPs, the volume (or here the density at fixed particle number N = 150000) is not single-valued in the imposed mechanical pressure P w , leading to large hysteresis loops. Note that the mechanical pressures P w corresponding to liquid and gas binodals are different, as expected. Parameters: ρ m = 25, v 0 = 20, v 1 = 5, τ = 1, vertical size L y = 50.
steady state This makes the mechanical pressure a state variable for PFAPs while the more complicated relationship between g 0 , h 0 and U for QSAPs breaks this link [1]. This explains the different roles of pressure in these two systems when considering change of ensembles.

Conclusion
In this article, we have shown how to derive the phase equilibria of MIPS for a number of different systems. At the hydrodynamic scale, the simple gradient terms that drive Active Model B [11] out of equilibrium still allow for the construction of a generalized thermodynamics, which leads to the definition of generalized chemical potential, pressure and surface tension. Using this formalism, we account quantitatively for the binodal curve of fully-phase separated systems as well as for its finite-size corrections. For quorum-sensing active particles, we have shown how to build a hydrodynamic description that fits within our generalized thermodynamic framework, using a combination of a local meanfield approximation and a gradient expansion. Despite these approximations, our formalism accounts quantitatively for the phase diagram of QSAPs. For particles interacting via repulsive pairwise forces, no closed hydrodynamics description including the relevant gradient terms exist in the literature. We thus followed an alternative route and showed how the binodals are selected by an equality of mechanical pressure complemented by a violation of the equilibrium Maxwell construction due to interfacial contributions.
Our identification of the relevant intensive variables governing the phase equilibrium of MIPS is important to define thermodynamic ensembles, which we have illustrated by considering the isobaric ensemble for QSAPs and PFAPs. We hope that our approach will pave the way towards a more general definition of intensive thermodynamic parameters [44][45][46] for active systems. Building a thermodynamic of active matter would further improve our understanding and control of these intriguing systems and has become a central question in the field [1,[11][12][13]17,20,25,32,35,37,[49][50][51][52][53][54].

. Hydrodynamic equation
The Fokker-Planck equation (A.15) for ϕ is equivalent to an Itō-Langevin dynamics for the position of the QSAP. From there, one can derive the collective dynamics of N QSAPs using Itō calculus, as was done many times in simpler settings [1,35,55]. For simplicity, we consider here the case D t = 0. One then finds the coarse-grained N -body density of QSAPs to follow the stochastic dynamicsρ = ∇ · ρD(ρ)∇[ log ρv(ρ) + ε τ v(ρ) ] + 2ρD(ρ)η , (A.17) We can now expandρ(r) in gradients of the density field ; As hinted before [1,11,35], the non-locality of the density sampling results in a 'surface tension generating' term κ(ρ). Interestingly, the asymmetric sensing ε affects both the free energy density f (ρ) and the gradient terms κ(ρ).

Appendix B.1. Constant-volume simulations
Simulations in the isochoric (constant-volume) ensemble were carried out in rectangular boxes of size L x ×L y with periodic boundary conditions using a modified version of the LAMMPS molecular dynamics package [56]. For simulations in slab geometry at coexistence, we chose L x = 500, L y = 300, and N = 115,000 particles. In order to ensure a stable, flat (on average) interface spanning theŷ-direction, these simulations were initiated by first equilibrating the particles in a smaller box of size 300 × 300 with v 0 = 0, yielding a near-close-packed phase. After this initial equilibration, the box was expanded in thex-direction and the activity was turned on, after which the system relaxed towards a phase-separated steady state. The simulations were run for a time t = 1000, the data being collected over the second half of this time. We compute binodal densities by coarse-graining the local density using a weighting function w(r) ∝ exp[−r 2 cut /(r 2 cut −r 2 )], where r is the distance between the particle and the measuring point, and r cut is a cut-off distance. Histograms of the density then show two peaks that we identify as the coexisting densities.
In order to handle the relatively large fluctuations in the position of the interface, the density and pressure profiles measured on each timestep was translated to a common origin. This point was taken to be the point where the density (averaged over y) has fallen below ρ = 0.95.

Appendix B.2. Constant-pressure simulations
Here, we used a simulation box with L y = 100, N = 10000 and periodic boundary conditions in the vertical direction only. In the x-direction, the system was confined by two walls modeled by harmonic potentials: where k controls the stiffness of the walls (we take k = 5). The right wall is fixed at x = x R , while the position x L of the left wall obeys the deterministic overdamped dynamics: where x i is the abscissa of particle i, P w the pressure externally imposed on the wall and γ its mobility, taken to be γ = 0.1. The dynamics is integrated by Euler time stepping, with the same time step as for the particles. To compute the phase diagram in the isobaric (N, P w , Pe) ensemble, we ramp P w up and down very slowly. For each Pe, the system was first equilibrated at a starting pressure P w , which was then incremented or decremented in steps of ±1 every 10 8 time steps.