Cluster phases and bubbly phase separation in active fluids: Reversal of the Ostwald process

It is known that purely repulsive self-propelled colloids can undergo bulk liquid-vapor phase separation. In experiments and large scale simulations, however, more complex steady states are also seen, comprising a dynamic population of dense clusters in a sea of vapor, or dilute bubbles in a liquid. Here we show that these microphase-separated states should emerge generically in active matter, without any need to invoke system-specific details. We give a coarse-grained description of them, and predict transitions between regimes of bulk phase separation and microphase separation. We achieve these results by extending the $\phi^4$ field theory of passive phase separation to allow for all local currents that break detailed balance at leading order in the gradient expansion. These local active currents, whose form we show to emerge from coarse-graining of microscopic models, include a mixture of irrotational and rotational contributions, and can be viewed as arising from an effective nonlocal chemical potential. Such contributions influence, and in some parameter ranges reverse, the classical Ostwald process that would normally drive bulk phase separation to completion.

It is known that purely repulsive self-propelled colloids can undergo bulk liquid-vapor phase separation. In experiments and large scale simulations, however, more complex steady states are also seen, comprising a dynamic population of dense clusters in a sea of vapor, or dilute bubbles in a liquid. Here we show that these microphase-separated states should emerge generically in active matter, without any need to invoke system-specific details. We give a coarse-grained description of them, and predict transitions between regimes of bulk phase separation and microphase separation. We achieve these results by extending the φ 4 field theory of passive phase separation to allow for all local currents that break detailed balance at leading order in the gradient expansion. These local active currents, whose form we show to emerge from coarse-graining of microscopic models, include a mixture of irrotational and rotational contributions, and can be viewed as arising from an effective nonlocal chemical potential. Such contributions influence, and in some parameter ranges reverse, the classical Ostwald process that would normally drive bulk phase separation to completion. Active colloidal fluids are non-equilibrium systems in which individual particles continually consume fuel in order to self-propel [1,2]. Time-reversal symmetry is broken locally, leading to features impossible in thermal equilibrium. One of these is motility-induced phase separation (MIPS) where an assembly of repulsive but active particles phase separates into bulk dense and dilute regions [3][4][5]. The MIPS phenomenology has been confirmed in simulations [5][6][7][8][9] and explored experimentally [10,11]. Although a far-from-equilibrium phenomenon, MIPS was first understood via an approximate mapping onto an effective equilibrium system undergoing standard liquid-vapor separation [3,4,12]. This has led to speculation that for MIPS, time reversal symmetry is effectively restored at the macroscopic level in steady state [3,9,[13][14][15][16][17].
However, simulations of repulsive active Brownian particles (ABPs) reveal that during coarsening, the macroscopic droplets of dense liquid arising via MIPS host a population of mesoscopic vapor bubbles that are continuously created in the bulk, coarsen, and are ejected into the exterior vapor [7]. This suggests that the steady state involves coexistence between a boiling liquid and excess vapor which we call bubbly phase separation. The lifecycle of the mesoscopic bubbles, if persistent even when coarsening is completed, would confirm that time reversal symmetry remains manifestly broken in this steady state. While not fully established for ABPs simulations, we present results below for continuum models that definitely show this behavior. The 'boiling liquid' is itself a microphase-separated state of vapor bubbles (or 'bubble phase') surrounded by dense liquid.
Puzzlingly, these observations on active Brownian particles are almost dual to the experimental findings on synthetic self-propelled particles [10,18,19], which can instead form dense living clusters. In these experiments, self-propelled colloids indeed aggregate but the clusters show highly dynamic particle exchange and, importantly, they never coarsen to reach the system size but saturate at a mesoscopic scale. (Such clusters can also show crystalline order, which we ignore here.) The 'boiling liquid' referred to above is akin to a cluster phase in which the dense and dilute regions have been interchanged. In the literature, several mechanisms for the emergence of cluster phases have been suggested, based on long range interactions (whether hydrodynamic [19][20][21] or chemotactic [22,23]), or on other, system-specific, effects [9,[24][25][26][27]. Cluster phases are also found in models of phase separation coupled to birth-and-death processes, where the number density of one chemical species is not conserved [28][29][30].
In this paper, we show that the observed occurrence of cluster phases and of bubbly phase separation do not require system-specific explanations, but already emerge from a completely generic continuum model of active phase separation. This differs from equilibrium phase separation models by allowing for broken time-reversal symmetry (TRS) at mesoscopic scales. At leading order in an expansion in gradients of the density, our model admits two distinct contributions to the particle current, only one of which was considered in previous studies of active phase separation [3,[31][32][33][34][35]. The additional current term was considered in [17], but its role in phase equilibria not studied there; here we show it to be crucial to predicting microphase separation. The unified description that we present, in which clusters and bubbles are dual to one another (see Section I), not only accounts for the existence of cluster phases (either alone or coexisting with an excess dense fluid), and bubbly phase separation, but also predicts the existence of a strict bubble phase in which the whole system is filled by a dense liquid that supports a fluctuating, but not coarsening, population of mesoscopic vapor bubbles. This phase was so far not reported in either experiments or particle-based simula- tions of active matter. Snapshots of the various steady states predicted by the continuum model considered in this paper, are given in Fig. 1.
In equilibrium systems subject to particle diffusion without momentum conservation, bulk phase separation is driven by the Ostwald process [36,37]. Here, local equilibrium requires the Laplace pressure jump across an interface to increase with its curvature. The chemical potential at a curved interface is thereby shifted, creating a compositional excess outside each droplet that decreases in magnitude with increasing droplet size. The local excess sets boundary conditions on a quasi-static solution of the diffusion equation, which slowly drives material down the density gradients from smaller to larger droplets. Small droplets evaporate, and large ones grow; the mean droplet size increases forever, ultimately leading to complete phase separation. Exactly similar arguments apply to the growth of vapor bubbles on the opposite side of the phase diagram.
A major achievement of the present work is to show analytically how the Ostwald process can become reversed in active fluids. This unexpected result is supported by numerical findings that demonstrate microphaseseparated steady states in the parameter regimes where the reverse Ostwald process is predicted. Interestingly, our model's parameter values can be chosen to stabilize either vapor bubbles dispersed in a homogeneous dense liquid, or liquid clusters dispersed in a vapor, but not both at once. This might help explain why computer simulations of hard-core active Brownian particles report only bubbles, and not cluster phases [5][6][7][8] with the reverse being true for experiments on autophoretic Janus colloids [10,18,19]. The microscopic interactions are clearly different in the two systems (and much more complex for the experiments) so perhaps they represent opposing halves of the parameter space of our coarsegrained model.
Our starting point is Model B, the φ 4 theory of diffusive phase ordering of a conserved order parameter φ with time-reversal symmetry [38]. Crucially though, and differently to previous attempts in this direction [3,[31][32][33][34][35], we consider the field theory obtained by adding all terms that break time-reversal symmetry to leading order in ∇ and φ, while retaining mass conservation, φ-independent mobility, and hence additive noise. (The latter is a crucial simplification for analytic work.) Following [39], we call the ensuing theory Active Model B+ (AMB+).
Like Model B [38], AMB+ is naturally introduced on the basis of symmetries, conservation laws, and the gradient expansion. However, we will show that a closely related continuum description -the main difference being that mobility is not constant and hence that noise is multiplicative -can be obtained from explicit coarse graining of microscopic models of self-propelled particles. This strongly suggests that the TRS-breaking terms in AMB+ are not forbidden by any microscopic mechanism that we might have overlooked, and hence that its prediction of reverse Ostwald regimes are generic, in the same way that Model B is generic for equilibrium binary fluid phase separation [36]. Linking microscopic to continuum field-theoretic models also paves the way to future investigations of cluster phases and bubbly phase separation in terms of microscopic parameters, and hence to the possibility of better controlling these new phases of matter experimentally. It might also suggest where to look experimentally for the strict bubble phase mentioned above.
The paper is organised as follows. In Section I we introduce AMB+ on the basis of symmetry arguments. In Section II we analyse the model at mean-field level, looking at stationary phase-separated solutions with spherical symmetry, and deriving coexisting densities as a function of the droplet radius. Such analysis is preliminary to the main analytical results of the paper, derived in Section III. There, we construct the mean-field phase diagram of AMB+ by analysing the Ostwald process. In Section IV, we finally switch on the noise in the field theory and, by means of direct numerical simulations, we show that microphase separations arise in direct correspondence to the regions where Ostwald process is reversed. Statistical properties of the microphase separated states are also analysed. In Section V, we show that the active currents responsible of microphase separation in AMB+ can be derived from explicit coarse-graining of a model of self-propelled particles. We conclude in Section VI with a discussion of future perspectives. Appendices contain details for those readers who are interested in the technical aspects of this work.

I. ACTIVE MODEL B+
An equilibrium system undergoing diffusive fluidfluid phase separation is governed on continuum scales by Model B in the classification of Hohenberg and Halperin [38]. Model B describes the dynamics of a scalar field φ(r, t) which, for the liquid-vapor transition, is a linear transform of the particle density ρ(r, t), chosen such that φ = 0 at the mean-field critical point density ρ c of the model [40]. The dynamics of φ follows a mass conservation law: where Λ(r, t) is a Gaussian white noise with zero mean and unit variance. M [φ], which is positive definite, is the mobility and D is the temperature. Model B has an equilibrium structure, meaning that the deterministic current J is proportional to the negative gradient of an equilibrium chemical potential µ eq , which is itself derived from a free energy functional, µ eq = δF/δφ. This is taken of square-gradient, φ 4 form: where f (φ) is the bulk free energy density and K > 0. By choosing φ to vanish at ρ = ρ c , we have set f (0) = f (0) = 0 so that there is no cubic term in f (φ). This is sufficient reason to work with φ rather than the particle density ρ. Note also that any linear term in f (φ) would contribute a constant to µ eq and therefore have no effect. Thus the local part of the free energy f is even to quartic order as written, with the critical point at a = 0 and phase separation for a < 0. (These statements hold without loss of generality, even if the microscopic interactions are asymmetric.) The neglect of terms in f beyond quartic, which would arise for example from expanding the ideal gas entropy term ρ ln ρ, is strictly justified only when the order parameter φ remains small. The coexisting mean-field binodal densities ±φ b obey φ b = −a/b; these should therefore map under the linear transform onto ρ = ρ c ± ∆ρ b with ∆ρ b ρ c . We assume this holds here, but also assume that the binodals lie outside the Ginzburg interval [40] so that we can avoid consideration of critical phenomena. (These are the topic of a separate paper [41].) In the language of Bray [36], we are thus interested in phase separation in the neighborhood of the zero-temperature dynamical fixed point, where noise terms might be significant for kinetics (for instance in allowing nucleation) but do not dominate the steady-state statistics, at least in the equilibrium case of Model B. The universality of this fixed point means that qualitative predictions, such as powerlaw scalings in time for the mean droplet size (R ∼ t 1/3 ), remain valid even beyond the regime of small φ to which the model at first appears restricted [36,40].
With the above proviso of avoiding the critical region, we can now rescale φ to set the bulk binodals of Model B at ±φ b = ±1. Alternatively we may achieve this by choosing which will apply from now on. The equilibrium interfacial tension then obeys σ eq = (8KA/9) 1/2 , and the interfacial width ξ eq = (K/2A) 1/2 [36]. Finally, we note that a mobility functional M [φ] appears both as a factor in the deterministic current J and as a multiplier of the noise Λ. This dual role is required by TRS, but it also emerges from explicit coarse-graining of active particle models that do not have TRS [3]. In this paper, as is normal in the Model B literature, we shall assume the mobility M to be φ-independent, and set it to unity hereafter. The technical advantages of avoiding multiplicative noise are substantial. This simplification is generally considered to be harmless because of the universality of the zero-temperature fixed point, although this can be broken by singular choices of M [φ] [42].
Detailed balance is broken in active matter, so that time-reversal symmetry, under which equilibrium Model B is built, should not longer be imposed on the continuum model at the level of (1,2). To describe active phase separation, we need to extend Model B to include terms that break TRS. Allowing that the noise remains nonmultiplicative, the mobility constant, and that there are no external fields coupled to gradients of φ (so that rotational symmetry is broken spontaneously if at all), the first nonlinear terms with such a property arise as contributions to ∂ t φ at order O(∇ 4 φ 2 ). This is because odd powers of ∇ are eliminated by isotropy and terms ∇ 2 φ n can be absorbed into the passive local free energy f (φ). Note that any term that can be absorbed by redefining F cannot break TRS.
The field theory we study, which will be called Active Model B+ (AMB+), is then obtained from Model B by adding all terms to order O(∇ 4 φ 2 ), while retaining constant M = 1 and D. It can be written Here arise two TRS breaking terms, in λ and ζ. Additionally, within F, the coefficient of the square-gradient term in (3) is modified to this order as K → K(φ) = K+2K 1 φ. AMB+ is the most general field theory describing the diffusive dynamics of a single scalar field to order O(∇ 4 φ 2 ), once the choice to keep constant mobility and additive noise has been made. (We must of course retain the term arising from the quartic part of the local free energy which is of order O(∇ 2 φ 3 ).) Any other term at order O(∇ 4 φ 2 ) can be written as a linear combination of λ, ζ and K 1 terms. The way TRS is broken in AMB+ is by destroying the free energy structure that underlies the equilibrium dynamics of Model B without altering the relation between mobility and noise. As shown below (and, for the λ term only, in [32,33]) this is what happens when one explicitly coarse-grains interacting active particle models. A quite different route would be to break the fluctuationdissipation relation between noise and mobility while leaving the free energy structure intact [43]. This alternative is interesting, but thus far we do not see a similar microphysical motivation for it.
There is a fundamental difference between K 1 , which merely modifies σ eq and does not break detailed balance, and the λ, ζ terms, which locally break TRS. The λ term can be thought of as locally defining a non-equilibrium chemical potential because, although λ|∇φ| 2 cannot be written as a derivative of any free energy functional, the resulting current remains of gradient form: The case of λ = 0 and ζ = 0 is Active Model B (AMB), derived and studied previously [32,33]. The qualitative physics of AMB is very similar to the one of passive Model B; in particular there is no microphase separation. Instead, the effect of λ = 0 is to change quantitatively the coexisting vapor and liquid densities with little or no qualitative change in coarsening behavior [33]. The full AMB+ model (5,6), and in particular the effect of the ζ term, remains unstudied until now in the context of phase equilibria and phase ordering kinetics. As we noted in earlier work [17], the ζ term is distinctive as it allows ∇ ∧ J = 0 so that, even in steady state, circulating real-space currents are possible. These are prohibited by detailed balance in all passive systems. They are separately prohibited in AMB (ζ = 0), despite broken TRS, by the gradient form of the current in that case.
At first sight one would not expect circulating current to play a role in phase equilibria or phase ordering because the rotational part of J is (by definition) divergence-free, and hence does not enter the equation of motion (5) for φ. Importantly though, when ζ = 0, although via Helmholtz decomposition one can always define a nonequilibrium chemical potential µ such that This is because Helmholtz decomposition and the gradient expansion do not commute. Nonlocality of µ[φ] will play a central role in the calculations below.
Before presenting those calculations, we note various further points. (i) Since it respects time reversal symmetry, our K 1 term should be relatively unimportant; we focus hereafter on K 1 = 0 but confirm this expectation in Appendix D. (ii) Reasons for choosing a symmetric local potential f (φ) = f (−φ) in Model B were given above and hold to quartic order, creating an invariance under (φ, λ, ζ) → −(φ, λ, ζ). This duality reduces by half the analytic and numerical calculations required; we exploit it extensively. But it is not fundamental, and we would expect no qualitative changes when using a nonsymmetric local free energy f . (iii) Although we have introduced AMB+ on the basis of conservation laws, symmetry arguments and the gradient expansion, we show in Section V that, alongside λ [32], the ζ term naturally emerges by explicit coarse-graining of particle-based models. (iv) In [17], a field theory containing all terms up to order O(∇ 4 φ 3 ) was defined (but not studied); our AMB+ model is a special case of it (with, in the notation of that paper, In the following, we study AMB+ with analytical tools and simulations in d = 2; our numerical methods are detailed in Appendix A. Unless otherwise specified, simulations are performed with A = 1/4 and K = 1. Analytical results are valid in any dimension d unless otherwise specified.

II. PHASE EQUILIBRIA AND THE PSEUDO-TENSION
We start by considering the mean-field limit of AMB+ (D = 0). In this Section, we compute the binodals φ 1 (of the vapor phase) and φ 2 (of the dense phase), defined as the coexisting densities for a flat interface at steady state. Then, we consider phase-separated stationary solutions with spherical symmetry and compute the corrections to the binodals φ ± (R) as a function of the droplet radius R. Our analysis generalises to AMB+ tools developed for AMB and other models whose currents take gradient form in [33][34][35]. It lays the foundations for the main analytic results of this paper, presented in Section III, where we calculate how λ and ζ affect the Ostwald process. Consider a flat interface. Clearly (without noise) the density profile depends only on the normal coordinate x and the problem is effectively one-dimensional. It is then easy to see that the ζ term and the λ term are interchangeable since ∇(∇φ) 2 = 2(∇ 2 φ)∇φ. Hence φ 1,2 only depend on the combination of activity parameters ζ−2λ. More explicitly, in d = 1, the deterministic current in (6) where (as previously) We look for the stationary solution φ(x) subject to boundary conditions φ(−∞) = φ 1 and φ(+∞) = φ 2 with a smooth interface around x = 0. Setting J = 0 implies that the chemical potential µ equals a constant valueμ, fixed by boundary conditions. Then Equation (8) alone is not enough to determine the binodals φ 1 and φ 2 . In the passive limit (λ = ζ = 0), a second condition is obtained by demanding equal thermodynamic pressure in the two bulk phases: P eq (φ 1 ) = P eq (φ 2 ), where the equilibrium bulk pressure is P eq (φ) = φf (φ) − f (φ). The same condition can also be derived, without reference to thermodynamics, by multiplying (8) by φ and integrating across the interface. A formal generalisation of the above result to AMB (and hence AMB+ so long as the interface is flat) was developed in [34]. This amounts to introducing a pseudodensity ψ(φ) and a pseudo-potential g(φ) satisfying: such that ψ → φ and g(φ) → f (φ) in the passive limit (λ → 0 and ζ → 0). Indeed, with the above definitions, multiplying (8) by ∂ x ψ and integrating across the interface gives where the bulk pseudo-pressure P (φ) is defined in terms of ψ(φ) and g(φ) as Thus the stationary conditions for finding the binodals φ 1,2 is to equate the bulk chemical potentials (8) and the pseudo-pressures (10,11) on both sides of the interface. Explicit results are easily obtained numerically because the solutions to (9) can be found explicitly. Details are given in Appendix B and the results shown in Fig. 2. Notice that, as first discovered in [33] for AMB, only in the passive limit (ζ = λ = 0) are the binodals at φ = ±1, which are the minima of the free energy density f (φ).
Despite its role in phase equilibria, P has no direct link to the mechanical pressure [33,34]. Indeed it appears so far to have no direct physical interpretation; the same applies to our other pseudo-quantities ψ, g and (as introduced below) σ. Their role instead is to exactly convert the new mathematics of active phase separation into forms that closely resemble the passive case, allowing familiar lines of argument to be applied here (albeit to give unfamiliar outcomes). This is a significant achievement since those lines of argument took many years to perfect [36], and without the transformation to pseudo-variables we could be waiting equally long for answers here. We will see below that our main results indeed match the passive ones under this transformation and also are more easily read in transformed variables than in the original ones.
We next generalize the same approach to curved interfaces, considering a dense droplet in a dilute environment (see inset of Fig. 3(a)). The case of a dilute bubble in a dense environment does not need to be analysed separately, as we recall that AMB+ is symmetric under the duality relation (φ, λ, ζ) → −(φ, λ, ζ).
Let us consider a stationary state given by a spherically symmetric droplet of radius R with bulk density φ + (R) inside the droplet and φ − (R) outside. Although metastable in the large system limit when the noise is switched on (D = 0) this configuration is, once φ ± are chosen correctly, a true stationary state of the mean-field dynamics. Denoting as r = |r| the radial coordinate from the centre of the droplet, we are interested in the stationary density profile φ(r) subject to boundary conditions We expect a smooth interfacial profile around r = R and are mainly interested in situations where R greatly exceeds the interfacial width. It will be crucial for what follows that, for d > 1, the ζ term in the current (6) can be written via Helmholtz decomposition as: Here the effective chemical potential related to J ζ is with ∇ −2 (|r − r |) the Green function of the Laplacian.
Observing that only the gradient part of J ζ affects the dynamics of φ, we can forget about A in the following. Crucially however, this construction comes at the price of a chemical potential µ ζ [φ] that is non-local. The full nonequilibrium chemical potential µ, defined via ∇ 2 µ ≡ −∇ · J, thus contains contributions from the passive (δF/δφ) and active terms (λ and ζ): which, using spherical symmetry, can be rewritten as Here φ = ∂ r φ. At steady state µ is a constant which we shall call µ I (R). Using the boundary conditions φ(0) = φ + and φ(+∞) = φ − , we have where To go from (17) to (18) we have assumed a sharp interface limit (R ξ eq ) so that we can approximate 1/r as 1/R. Equation (16) states that the full chemical potential µ inside the droplet is equal to that outside so that no current flows across the interface at steady state. However, the equilibrium part of the bulk chemical potential µ eq (φ) = f (φ) has a jump ∆ as we cross the interface. This jump is present to cancel one in the µ ζ term whose nonlocal form (13) resembles a Coulomb integral with a curvature-dependent charge density on the interface. Such a charge configuration would create an interfacial step discontinuity in the electrostatic potential and the same happens here for µ ζ in the sharp interface limit.
As we did for a flat interface, we proceed to find the second condition for φ + and φ − . This is done by multiplying (16) by ∂ r ψ, integrating across the interface, and again assuming that R is much larger than the interfacial width. We get where σ is a pseudo-tension defined by defined in terms of two constants that depends on the full shape of the interface: The pseudo-pressure balance (19) is formally very similar to what one has in equilibrium, where the Laplace pressure jump crossing the interface is ∆P eq = (d−1)σeq R , where σ eq = K ∞ 0 φ 2 (y)dy is the equilibrium interfacial tension. Equation (19) is therefore exactly what one would obtain from the classical equilibrium results upon replacing the pressure and the surface tension with their pseudo counterparts. Moreover, in the passive limit (ζ → 0 and λ → 0), our results reduce to the equilibrium ones, and in particular σ → σ eq . A related analysis of the pseudo pressure applicable to AMB (but not AMB+) was presented in [35].
The coexisting densities for a spherical droplet of radius R in a dilute environment are given by equations (16) and (19). However, no explicit values of φ + (R) and φ − (R) can be found unless S 0 and S 1 are known, which depends on the precise shape of the interface. In Appendix C, we explain how to obtain the shape of the interface perturbatively in 1/R. Doing so, and imposing the conditions (16) and (19), we compute perturbatively the values of φ ± up to corrections of order O(1/R 2 ). These results, which become exact in the sharp interface limit, are shown in Fig. 3 (lines) along with results from mean-field numerical simulation of AMB+. There is good agreement between the two, given the relatively modest values of R simulated.
We are finally able to derive the most important result of this Section. In equilibrium (λ = ζ = 0), the interfacial tension σ eq is always positive. However, in AMB+, σ can become negative; from (20) this happens when This is the criterion for σ < 0 for a dense liquid droplet in a vapor environment. The analogous condition for a vapor bubble surrounded by the dense liquid then follows from the duality transformation (φ, λ, ζ) → −(φ, λ, ζ). The phase diagram is thereby divided in three regions, as represented in Fig. 4(a), as follows. In region A, σ > 0 regardless of the sign of the interfacial curvature. Note that AMB (ζ = 0) lies entirely within this region. In region B, σ < 0 for liquid droplets in a vapor environment but σ > 0 for vapor bubbles in a liquid environment. In region C, conversely, σ is negative for vapor bubbles and positive for liquid droplets.
We have shown that the pseudo-tension σ determines the pseudo-pressure jump at a curved interface in the same way that, in equilibrium systems, the usual tension σ eq determines the Laplace pressure jump. However, this does not mean that σ has any direct connection to the mechanical tension as defined for active systems e.g. in [44]. Crucially, the presence of a negative pseudotension does not cause interfaces to become unstable: the coexisting densities φ ± (R) remain well defined. The effect of negative σ is instead to change how these quantities depend on R. In equilibrium systems, both φ ± are decreasing functions of R. This remains true in AMB+ for σ > 0, as shown for a particular choice of λ and ζ Ostwald ripening is reversed where σ < 0. In region A, σ > 0, and we have normal/forward Ostwald ripening for both bubbles and droplets. In region B (resp. C) Ostwald ripening is reversed for dense droplets dispersed in a dilute environment (resp. dilute bubbles in a dense environment). Points are results from simulations of two droplets (red/square for reverse and blue/circle for forward Ostwald ripening). (b) Evolution of two droplets with initial radii R1 > R2 for (ζ, λ) = (−1, 0.5) (red), corresponding to region A, and (−4, −1) (blue), corresponding to region B. Solid lines are direct numerical simulations of (5,6) with D = 0 and dashed lines are predictions from our theory (see Appendix E). We observe forward Ostwald ripening in the first case and reverse Ostwald ripening for the second, as predicted.
within region A, in Fig. 3. Instead, when λ and ζ belong to region B, while φ + (R) is still a decreasing function of R, φ − (R) now increases with R.

III. REVERSE OSTWALD RIPENING: MEAN-FIELD PHASE DIAGRAM
In this Section we show how negative pseudo-tension affects the Ostwald process, which is the kinetic pathway to bulk liquid-vapor phase separation in equilibrium as described by Model B. As in the previous Section, we focus on region B of Fig. 4(a); the behavior in region C then follows by the duality relation.
A dramatic change in behavior on passing from region A to region B in Fig. 4(a) is natural because, as explained in Section II and in Fig. 3, this changes the sign of ∂ R φ − (R). In region A, just as in equilibrium models, σ > 0 giving normal or forward Ostwald ripening in which the density outside a small droplet is higher than that outside a bigger droplet and thus the current flows from small droplets to big ones. In contrast, in region B, σ < 0 and the density outside a small droplet is lower than that outside a bigger droplet; the particle current now flows from big to small droplets. This is the regime of reverse Ostwald ripening in which small droplets grow at the expense of larger droplets, which shrink.
We now make the above reasoning precise by generalizing to AMB+ the standard arguments for Ostwald ripening in equilibrium models [36,37]. We work in d = 3 dimensions; for d = 2 there are correction terms to the following arguments involving the logarithm of the total phase volume occupied by droplets, but these do not change our conclusions (see, e.g., [45]).
We consider a spherical droplet of radius R and denote as r = |r| the distance from the centre of this droplet. On a rapid time scale the densities just inside and outside the droplet relax to their quasi-stationary values, φ + and φ − , which differ from the binodal densities φ 2 and φ 1 respectively by terms of order 1/R. Now we assume the droplet to be surrounded by a sea of distant droplets, such that the mean density of the vapor phase at r = ∞ is φ s = φ 1 + (t), where is the supersaturation. We wish to calculate the time evolution of the radius R(t) of our central droplet; this depends on (t). As is standard [36] our argument assumes a quasi-static diffusion so that ∇ 2 µ(r) = 0 at any time. This is to be solved with boundary conditions µ(∞) = f (φ s ) and µ(R) = µ I (R). Once the resulting solution µ(r) is known, mass conservation allows us to give an equation for the temporal evolution of the droplet radius R(t).
Rewriting (19) with φ ± = φ 2,1 + O(1/R), we find that the chemical potential at the interface of the droplet is: where ∆g = g(φ 2 ) − g(φ 1 ), ∆ψ = ψ(φ 2 ) − ψ(φ 1 ), and σ is the pseudo-tension defined in (20). Note that for λ = ζ = 0, one has g = f , ψ = φ and σ = σ eq , thereby recovering the classical result for equilibrium models [36]. Next we solve ∇ 2 µ = 0 subject to boundary conditions µ(r = R) = µ I (R) and µ(r = ∞) = f (φ s ) to obtain the quasi-static nonequilibrium chemical potential: The inward radial current at the interface of the droplet J I now obeys −J I = −∇µ R = (µ I (R) − f (φ s ))r/R. Finally using the conservation of mass: we obtain the time evolution of the droplet radius R(t): where R s is a fixed-point radius, β is a rate parameter, and ∆φ = φ 2 − φ 1 . The explicit expressions for R s is: and likewise for β: In equilibrium models, the surface tension σ → σ eq is always positive and thus β and R s are also positive. R s is then an unstable fixed-point radius: droplets smaller than R s shrink and droplets larger than R s grow. This conventional Ostwald process is maintained, for droplets and bubbles alike, throughout region A of the phase diagram of AMB+, as shown in Fig. 4(a), where σ > 0.
In AMB+, however, the pseudo-tension σ is negative when condition (23) is met, which implies that the rate β is also negative. To determine the sign of R s , we Taylor expand Eq. (28) for small to obtain: where we have used the fact that φ 1,2 satisfy (8) and (10). Since the binodal φ 1 is always outside the spinodal, f (φ 1 ) > 0 is always positive. We also note that when σ (hence β) is negative, the supersaturation , which is set by the average of φ − − φ 1 across the population of distant droplets, is necessarily also negative; see Appendix E below, and also compare the red and blue curve in Fig. 3(a)). We thus conclude that, for small in regime B (σ < 0), R s > 0 and β < 0. R s is now a stable fixed point radius: droplets of smaller than R s grow and those larger than R s shrink. This is the regime of reverse Ostwald ripening for liquid droplets dispersed in a continuous vapor phase. Regime C in the mean field phase diagram ( Fig. 4(a)) gives another region where Ostwald ripening is reversed, but now for vapor bubbles dispersed in a dense liquid environment. This can be derived explicitly as above by considering a bubble of radius R with density φ 1 inside and φ 2 − at infinity. Alternatively, it follows directly from the duality relation of AMB+. Note that there is no choice of parameters for which both droplets and bubbles are subject to the reverse Ostwald mechanism. In consequence, as confirmed numerically below, we expect microphase separation to occur in AMB+ on one side of the phase diagram or the other but not both.
To test our analytic mean-field predictions, we performed numerical simulations in the mean-field limit (D = 0) for several values of λ and ζ, using as initial condition two droplets with different initial radii R 1 (0) > R 2 (0). Fig. 4(b) shows the typical time evolution of the radii R 1 (t) and R 2 (t) (solid lines) for parameter values in regions A and B in Fig. 4(a). In region A we confirmed that the larger droplet grows at the expense of the smaller one whereas in region B the big droplet shrinks until R 1,2 reach the same value. For all the other values of λ and ζ that we also tried, our numerical results for two droplets likewise confirmed the analytic predictions for whether the normal or the reverse Ostwald process should be seen (see points in Fig. 4(a)).
As shown in Appendix E, the above calculations can be extended to track explicitly the evolution in time of the two radii R 1,2 (t) in a binary droplet system. (We have not seen a calculation of this kind in the literature even for equilibrium models. To achieve this one needs to assume that the second droplet acts as a distant boundary condition on the first (and vice versa), while retaining radial symmetry. This effective-medium approximation was used above (as is standard [36]) for the multi-droplet case, but it is clearly less valid when only two droplets are present. An additional source of error is that the calculation assumes the two droplets are separated by a distance large compared to their radii which is not really met in simulations due to computational limitations. Given these shortcomings, one does not expect quantitative agreement, but as shown in Fig. 4(b), these predictions are qualitatively in accord with our two-droplet simulations.

IV. SIMULATIONS WITH FINITE NOISE
Our mean-field analysis show that the reverse Ostwald process can arise generically in active phase separation, without the need for specific interaction mechanisms or fine-tuning of parameters. It is then natural to speculate that this leads to microphase-separated stationary states, at least when dynamical noise is present so that the steady-state is independent of initial conditions. Indeed, when noise is present, each time a droplet/bubble is nucleated, it will converge to roughly the same size as the others. To test this hypothesis, we have performed full stochastic simulations of AMB+ with finite noise D = 0.2 on a 192 × 192 lattice (see Appendix A for numerical details). We consider two sets of activity parameters: (ζ = 1, λ = −0.5) and (ζ = 4, λ = 1) corresponding to regions A and C of the mean-field phase diagram. These parameters have the same ζ − 2λ and thus the same binodals: φ 2 0.86 and φ 1 −1.10. Fig. 5(a) shows the steady state phase diagram as a function of global mean density φ 0 = 1 V dr φ(r, t). Outside the binodals, φ 0 < φ 1 or φ 0 > φ 2 , the uniform phase is always stable. In region A (top row), the phase diagram resembles that of passive Model B, with bulk phase separation into dense liquid and dilute vapor phases, except that the binodals φ 1,2 are now shifted. This is the same phenomenology as described previously for AMB [33], whose parameters (λ = 0, ζ = 0) indeed lie within region A.
The phenomenology in region C (bottom row) -or equivalently by duality in region B -is much more interesting. Besides the low-density and high-density uniform phases when φ 0 is outside the binodals, for φ1+φ2 2 φ 0 φ 2 in region C we observe initial coarsening of dilute bubbles. At long times, however, the coarsening arrests and the number of bubbles fluctuates around some average (Supp. Movie 1). This represents nonequilibrium microphase separation in the form of a homogeneous bubble phase, which becomes a cluster phase ( Fig. 1(b)) when we flip into region B by the duality transformation (ζ, λ, φ 0 ) → −(ζ, λ, φ 0 ). We have checked that this homogeneous microphase separation is the true steady steady within the above range of φ 0 , by establishing that it is reached not only from a uniform initial state but also one at bulk phase coexistence (Supp. Movie 2).
At global densities φ 0 in the range φ 1 φ 0 φ1+φ2 2 within region C, the system creates a globally inhomogeneous steady state resembling conventional bulk phase separation, except that within the dense phase, vapor bubbles are continuously formed and then expelled to the exterior bulk vapor (Supp. Movie 3 and 4). This 'bubbly phase separation' looks very similar dynamically to what was reported in simulations of active Brownian particles with hard-core repulsions [7]. Given what we have learned about the global phase diagram, it can best be viewed as a bulk coexistence between a homogeneous bubble phase and a conventional vapor. By duality we find a range of densities within region B in which there is bulk coexistence between a cluster phase and dense liquid. We have seen no reports of this in either experiments or particle-based simulations on active particles; but since neither has explored parameter space fully, we predict that such states will be found in future. The dynamical properties of the homogeneous bubbles phase (and its cluster-phase dual) are intriguing. In steady state, bubbles are created by nucleation and destroyed mainly by coalescence with others at low noise (Supp. Movie 1). They can, however, also disappear by shrinking especially at high noise (Supp. Movie 5). The steady state average bubble sizeR is then expected to depend strongly on noise, set by a balance between the noise-induced nucleation, the mean-field dynamics inducing reversed Ostwald process, the diffusion (which leads to coalescence), and the possible shrinking at high noise. Interestingly, see Fig. 5(c), we found thatR increases when decreasing the noise strength D, which is a nontrivial statistical property of this phase.  Fig. 4(a), and ζ = 4, λ = 1 corresponding to region C. In A, the steady state phase diagram resembles equilibrium Model B with full phase separation into dense (yellow) and dilute (blue) phases. In C, besides the low and high density uniform phases, we observe either bubbly phase separation or a homogeneous bubble phase. Since AMB+ is symmetric under (φ, λ, ζ) → −(φ, λ, ζ), the phase diagram in region B is obtained by exchanging dense with dilute regions. (b) Probability distribution of bubble radius for parameters ζ = 4, λ = 1 as a function of time, starting from an initial uniform density φ(r, t = 0) = φ0 = 0.6. The system reaches a steady state at t 32000. (c) Average radius of bubblesR and number density of bubblesn in the homogenous bubble phase as a function of D. Parameters: ζ = 4, λ = 1 and φ0 = 0.6.
We also measured the probability distribution of the bubble radius P (R), see Fig. 5(b), as a function of time when starting from a uniform initial state. We found that P (R) converges to a well defined shape at long times. We leave for the future any attempt to predictR and P (R) in the uniform bubble/cluster phase. In particular, it would be interesting to generalize the Lifshitz-Slyozov theory [36], using the mapping to pseudo-pressure and pseudo-tension, to predict P (R). However, the fact that R depends on the noise level suggests that this theory, which is strictly mean-field in character, will not suffice when σ < 0. Since the noiseless two-droplet problem then has a stable solution at equal droplet size, one possibility is that, for any given initial conditions, a noiseless theory at Lifshitz-Slyozov level lands the system somewhere on a manifold of static solutions with strictly monodisperse droplets. So long as they all have the same size, then at effective medium level, such droplets can be in stable equilibrium with each other and may have no reason to evolve further. If so, the selection of a mean droplet size would have to depend on a perturbation to this meanfield scenario, such as dynamical noise.R could then be set by a balance between noise and deterministic terms, as speculated in the previous paragraph.

V. EXPLICIT COARSE GRAINING FROM MICROSCOPIC DYNAMICS
In Section I, we introduced AMB+ on the basis of topdown arguments based on conservation laws and a gradient expansion. Here we show how the form of the TRSbreaking terms in the force density J/M in (6), crucial for the phenomenology found above, emerge from explicit coarse graining of a microscopic dynamics. More precisely, we obtain by coarse-graining a field theory closely related to AMB+, the main differences being that the mobility M is density-dependent, and the local free energy is not a polynomial. We do not consider these differences important, because they also arise for Model B, where constant M and the polynomial form of f (φ) are not fundamental properties inherited from some microscopic model, but strategic choices made to simplify the field theory into one of φ 4 form with additive noise only. AMB+ makes the same choices, so our task here is only to check that the ζ and λ terms can emerge as expected from microscopic models.
We start from a model of active Ornstein-Uhlenbeck particles (AOUPs) [9,14,46] where each overdamped particle i experiences first an active force v 0 u i , where u i is exponentially correlated unit noise with persistence time τ , and second a two-body interaction with other particles. In a similar spirit to liquid-state theory [47], we treat separately the hard-core part of the pair interaction from its softer tail. Following [3,4,32] we argue that the hard-core interactions effectively coarse-grain into an effective propulsion speed v(r, [ρ]) that depends functionally on the coarse grained density ρ of active particles.
Differently to previous approaches [34,48], we keep explicitly the soft pair interaction via a residual two body potential U = U (|r i − r j |). Thus, with unit damping, the particles' positions r i and directions u i obey: where η η η i and ξ ξ ξ i are independent Gaussian random variables with zero mean and unit variance. D t is a passive translational diffusion constant. An objection to this treatment of AOUPs is that replacing the hard-core interactions by v[ρ] sacrifices the equation of state between pressure and density that was present for the original model [49,50]. However, this refers to mechanical pressure which plays no part in the analyses of this paper. In any case, (31)(32) can equally well be viewed as a microscopic dynamical model in its own right. At this level it describes motile micro-organisms whose exponentially autocorrelated (AOUP-like) self-propulsion has a root-mean-square swim speed v(r, [ρ]) that depends on local density via a chemical (quorum-sensing) interaction rather than a mechanical one [34]. There is in addition a soft pair interaction U between the motile particles.
We sketch here the coarse-graining procedure, focusing on the effect of U , and refer to Appendix G for more details. Technically, the novelty of our approach stands in first taking the small τ limit, and then writing an equation of motion for the empirical density ρ d (r, t) = i δ(r−r i ) directly, instead of first writing one for ψ d (r, u, t) = i δ(r−r i )δ(u−u i ) and then passing to the hydrodynamic limit to obtain the evolution of ρ d as is more often done [34]. Our method is much more compact and gives us the possibility of treating U = 0, which would pose difficulties via the standard route. Still, for U = 0, we get exactly the same results as already known in the literature [3,34].
It is well known from the theory of slow-fast dynamical systems [51,52] that, in the limit of small τ , u i can be eliminated from (31,32) leading tȯ where, importantly, the equatily has to be interpreted in the Stratonovich sense, which is the meaning of = s . This equation can now be transformed to Itô form, givinġ We then write the Dean equation [53], which is the equation formally solved by ρ d and, finally, we replace ρ d by its smooth limit ρ. Note that only on this replacement is information lost and the coarse graining really done [39,54]. Such a procedure can be rigorously shown to give correctly the statistics of ρ d at least in the case where the two body interaction U is weak and the density large (the so-called Kac limit) [39,55]. More generally, passing from ρ d to ρ remains a standard procedure in the theoretical physics literature even far from this limit [56][57][58], although its proper justification is then an outstanding mathematical question [55].
The resulting dynamics is written as (Appendix G): where The effective force densityJ cg sets the deterministic current J cg via J cg = M [ρ]J cg . The force density due to quorum sensing and to the two-body potential add lin-earlyJ whereJ qs = −∇µ qs (38) Here denotes the convolution product and we have explicitly allowed that v is a functional of ρ. The subscript qs stands for quorum sensing (see above).
It is now clear that, if U = 0, the effective force densitỹ J cg becomes the gradient of a nonequilibrium chemical potentialJ cg =J qs = −∇µ qs , irrespective of the functional dependence of v on ρ and of whether an equilibrium structure is present (that is, whether µ qs = δF qs /δρ for some functional F qs [ρ]). This special case has been widely considered in the literature [3,33,34]. Explicitly, expanding v[ρ] at lowest order in gradients and retaining only terms allowed by rotational symmetry, we write v[ρ] = v(ρ) + β(ρ)|∇ρ| 2 + γ(ρ)∇ 2 ρ. Note that β = 0 iff the active particles can detect local gradients directly; the γ term instead describes particles whose speed depends on the average density in a small neighborhood [4]. Expanding (39), we obtain where F qs = dr f qs (ρ) + K qs + 2K qs 1 ρ 2 The expressions of f qs , K qs , K qs 1 and λ qs in terms of microscopic parameters are given in Appendix G and agree with those previously obtained in the literature.
The above coarse-graining at U = 0 represents a pure quorum-sensing model. Importantly, this gives only one of the two TRS-breaking terms in (6), namely λ. It does not give us the ζ term that distinguishes AMB+ from AMB and is responsible for all the new physics discussed in this paper. But we can now see immediately how the ζ term appears inJ cg in (37) as soon as the soft potential U is nonzero: it does so becauseJ r is, according to (40), not the gradient of any scalar function. More explicitly, performing a gradient expansion, we obtaiñ where the explicitly derived ζ parameter is Here v 0 , v 1 and γ 0 are defined by v(ρ) = v 0 + v 1 ρ + o(ρ) and γ(ρ) = γ 0 + O(ρ). Finally, the constants B and C only depend on the two body potential U where Ω is the surface of the unit sphere. The effective chemical potential µ r in (43) has the same form as (41) with the subscript qs replaced by r and the explicit dependence of f r , K r , K r 1 and λ r given in Appendix G. These results from coarse graining now pass to those of AMB+ via (i) the linear transformation from ρ to φ [32]; (ii) replacement of the local part of f cg (ρ) by an even quartic polynomial in φ; and (iii) suppression of the density dependence in the mobility M , allowing this to be set to unity so that J cg =J cg and the noise becomes additive. As already stressed at the beginning of this Section, these three steps exactly mirror those used when passing from coarse-grained microscopic models for diffusive phase separation in passive systems to their canonical incarnation as passive Model B [36].
In summary, our explicit coarse graining shows that both of the TRS-breaking terms considered in AMB+ (5,6) emerge naturally from at least one well-motivated model for the microscopic dynamics of interacting active particles. A more careful study of that coarse-grained model in its own right (without the three additional simplifications listed above) is postponed to future works.

VI. CONCLUSION
We have presented results for Active Model B+ (AMB+), a field theory that, at leading order in a gradient expansion, fully generalises the canonical model of passive phase separation (Model B) to break timereversal symmetry. The model has two separate activity parameters, λ and ζ. When both are small, Ostwald ripening operates normally and bulk phase separation is retained. It is also retained at large activity when λ and ζ have opposite signs. However, when both are large and positive, one finds a regime of reverse Ostwald ripening for vapor bubbles in a dense liquid; for both large and negative, one finds this for dense droplets in a vapor.
This unexpected reversal of the Ostwald process is captured mathematically by a negative pseudo-tension which plays the role of the equilibrium interfacial tension for the purposes of phase equilibria. This is however not a mechanical tension and indeed the interfaces between dense and dilute phases remain stable even when it is negative. The reversal of Ostwald ripening by active processes is interesting, because in passive systems Ostwald ripening is often difficult to arrest, let alone reverse. This is a main obstacle to the stabilisation of emulsions, where arrest can be achieved by trapping particles inside droplets that are insoluble in the surrounding fluid [59][60][61]. Ostwald ripening can however also be arrested by continuous chemical reactions [29,30] or by continuous shearing [62], both of which maintain the system away from equilibrium, as do our active current terms.
Our numerical studies confirm that the outcome of the reverse Ostwald process in AMB+ is a microphase separation, establishing this as a generic feature of active matter, without the neccessity of long-range interactions caused by hydrodynamics [19][20][21], chemotaxis [22,23] or other system-specific causes [9,[24][25][26][27]. The resulting microphase separated bubble and cluster phases should be seen as new (but intimately related) phases of active matter, whose detailed properties require further investigation in the future. We predict that each of these phases can either fill space in its own right, or coexist with a bulk phase comprising excess vapor (in the case of bubbles) or liquid (in the case of clusters). In the former case, within the microphase-separated region, bubbles constantly form, migrate to the boundary, and pop into the vapor: we have called this 'bubbly phase separation'. A similar phenomenon, with dense and dilute regions interchanged, is predicted at coexistence between a liquid cluster phase and excess liquid. Indeed, within AMB+ there is duality relation that maps every aspect of bubble behavior onto a corresponding one for clusters.
As stated in the introduction, cluster phases have been observed in experiments on active colloids [10,18,19], whereas bubbly phase separation has been observed in simulations of Active Brownian Particles with hard core repulsions [7]. This suggests a conjecture that, upon coarse graining to continuum level, these two systems map into opposite regions of the (λ, ζ) parameter space of AMB+. So far, we are not aware of any observationeither experimental or computational -of either a homogeneous bubble phase, or of coexistence between a homogeneous dense phase and the cluster phase. We hope our prediction of their generic existence will motivate further investigations in this direction.
The statistical properties of the microphase-separated states are intriguing and also merit further study. In particular, we observed that the average bubble/cluster size increases when decreasing the noise level in the field theory. We so far have only a speculative explanation for this, involving a steady-state balance between deterministic and noisy terms in the evolution equations.
Also of interest is our finding of a transition line separating normal phase separation from bubbly phase separation. It would be good to identify specific microscopic models where this transition line can be studied; this might require independent control of the noise level and the mobility or interaction parameters. (There may hence be some active systems, perhaps including hardcore active Brownian particles, whose phase separation is always bubbly.) This transition signifies the presence of two distinct regimes, one where time reversal symmetry is almost restored at mesoscopic spatial and temporal scales, as speculated in the literature for over a decade [3,9,[13][14][15][16][17], and a second one where it remains strongly broken (in at least one of the coexisting phases). It would be natural to start such an investigation from the minimal model we were able to coarse-grain explicitly in Section V, where two-body forces are added to a quorum sensing model.
Meanwhile our coarse graining of that model confirms that nothing prevents the emergence from microscopics of the non-gradient currents governed by the activity parameter ζ at continuum level. These currents lead directly to reverse Ostwald ripening and microphase separation via a nonlocal contribution to the effective nonequilibrium chemical potential µ (as defined by ∇ · J = −∇ 2 µ). Our linking of microscopic to macroscopic models should help pave the way to future investigations of cluster phases and bubbly phase separation in terms of microscopic parameters, leading to better control of these nonequilibrium states of organization in simulations and, eventually, in experiments.
In this paper we have mainly worked at mean-field level and therefore not addressed the effects of activity near the liquid-vapor critical point, where fluctuations are dominant even in the passive limit. A preliminary numerical investigation of the critical dynamics of ABPs has recently appeared [63] suggesting, albeit inconclusively, that the universality class is altered from the passive case. A Renormalization Group study of AMB+ in the critical region, to be published elsewhere [41], suggests that the combined action of ζ and λ can indeed strongly change the behavior, but only beyond a finite activity threshold. 740269. MEC is funded by the Royal Society.
Author contributions. ET and CN contributed equally to this work.
All numerical simulations are performed at dimension d = 2 in a square box of size L×L with periodic boundary conditions on all sides. We choose A = 0.25, K = 1 and fix K 1 = 0 since this term does not change the observed phenomenology (see Appendix D). We discretise time as t = n∆t and space as x = i∆x and y = j∆y where n, i, j are integers. We choose ∆t = 0.001 and ∆x = ∆y = 1; our numerical simulations are converged (with respect to averaged quantities such as coarsening rates) at these values of ∆t and ∆x.
The dynamics of the field φ n ij then follows Itô integration in time: where Γ n αij is now a Gaussian random variable with zero mean and variance Γ m αij Γ n βkl = δ mn δ αβ δ ik δ jl . Now to compute the spatial derivatives, we use finite difference at order O(∆x 8 ): and similarly for ∂ y φ ij . For second derivatives such as ∂ 2 x φ ij , we apply (A2) twice [17]. This is to ensure that the equilibrium limit of (A1) satisfies detailed balance exactly on the lattice. limit (λ = ζ = 0) is given by Then, finding the binodals amounts to find those φ 1,2 such that the two conditions on bulk chemical potentials (8) and on the pseudo-pressures (10) are satisfied. This is a simple numerical problem and the results are reported in Fig. (2). As can be seen [33], the binodals are at φ = ±1 only in the passive limit (ζ = λ = 0), which are the minima of the free energy density f (φ).
Appendix C: Interface shape and perturbative corrections to binodals for curved interfaces We give here details about finite size corrections to binodals for spherically symmetric stationary solutions of dense droplets in a dilute environment in the meanfield limit (D = 0), i.e. the computation of φ ± (R) as a function of R.
As emphasised in the main text, the conditions (16) and (19) are not sufficient to find φ ± (R). Indeed they depend implicitly, via S 0 , S 1 , on the shape of the interface. Here we discuss the structure of the interface of a spherical droplet, and show how to compute w(φ) = (∂ r φ) 2 perturbatively in 1/R. The knowledge of w(φ) allows us to compute S 0 = exp ζ−2λ K φ + φ− φ+ w(x)dx and K x dx (see equations (21)(22)). Combined with the stationary conditions (16) and (19), this fixes φ + (R) and φ − (R) perturbatively to order O(1/R 2 ).
We consider a spherical droplet with density profile φ(r). As all along the paper we consider a dense droplet in a dilute environment. In (15), we change variable from φ to w [33], giving where the prime now indicates derivative with respect to x and we have defined: If we temporarily ignore the functional dependence of s on w, and assume that µ I (R) is known, Eq. (C1) is linear and can be solved as The knowledge of φ + , φ − then allows us to fix the integration constant C.
The perturbative strategy then goes as follows. At O(1/R 0 ), we can approximate φ + and φ − to be the binodals, i.e. φ ± φ 1,2 and, at this level of approximation, the chemical potential is also known µ I (R) =μ . We can thus obtain w(φ) from (C3). Using this value of w(φ), we can then obtain S 0 and S 1 . Once these are known, we can solve the simultaneous equations (16) and (19) to obtain φ ± , along with µ I (R) at next order O(1/R). Using these new values of φ ± we again calculate w(φ), S 0 and S 1 , and again solve (16) and (19) to get φ ± at O(1/R 2 ). We cannot go further because (C1) is obtained discarding terms of order O 1/R 2 which would otherwise make the equation for w non-autonomous. However, we already have all the results we need. This perturbative solution was implemented numerically to give the results in Fig. 3.
Appendix D: Effect of K1 = 0 As stated in the main text, setting K 1 = 0 does not affect qualitatively any of the conclusions of our work. All our results can be explicitly generalized to this case, but we restrict ourselves here to a few statements about how K 1 enters the calculations.
First of all, the chemical potential in spherical symmetry now reads where we recall that K(φ) = K + 2K 1 φ and we have defined ν = λ − K 1 − ζ 2 . Let us first consider the calculation of the binodals φ 1,2 . This goes similarly to the one presented in Section II except that ψ takes the form while g K1 is the solution to ∂g K1 /∂ψ K1 = ∂f /∂φ that we do not report here explicitly. For K 1 = 0, φ 1,2 can thus be obtained solving (8), (10) with ψ and g replaced by ψ K1 and g K1 . Similarly, when considering a droplet of radius R, the values of φ ± (R) are also affected by K 1 = 0. Again, the calculations are similar to those in Section II and III but with ψ and g replaced by ψ K1 and g K1 .
Generalising the argument leading to Ostwald ripening presented in Section III to the case of K 1 = 0, we find that (27)-(29) still hold after replacing σ with where S 2 , S 3 are parameters that depends on the shape of the interface and are defined by (φ(y)) φ 2 (y)dy .
Conditions on λ, ζ, K 1 such that σ K1 < 0 can be obtained from (D4), leading again to reversed Ostwald ripening. We leave the detailed study of the phase diagram in terms of λ, ζ, and K 1 to future studies.

Appendix E: Evolution of two droplets
In this Appendix, we adapt the argument of Section III to explicitly deal with the initial condition with two droplets as in Fig. 4(b), and predict the full time evolution of their radii R 1,2 (t). We believe these results to be interesting by themselves, as we have found nowhere in the literature -not even the equilibrium literature -an adaptation of the argument leading to Ostwald ripening able to track the evolution of two droplets. Moreover, it will turn out clearly that the quantity used in Section III has to be interpreted as the supersaturation among distant droplets, and it is thus negative when σ < 0. This is a subtle but important point as it implies, as already stressed in Section III, that < 0 when σ < 0.
First we consider droplet 1 whose radius is R 1 (t). We assume it to be stationary, and thus its density inside to be φ + (R 1 ), while the density just outside has to be φ − (R 1 ). We furthermore assume that, at infinity, droplet 1 "sees" the outer density of droplet 2, i.e. φ − (R 2 ). Following the derivation in Section III, but without using the approximation φ ± φ 1,2 , the time evolution of R 1 (t) is then (for dimension d = 3): where ∆φ(R 1 ) = (φ + (R 1 ) − φ − (R 1 )); we recognise f (φ − (R 1 )) to be the chemical potential just outside the interface of droplet 1 and f (φ − (R 2 )) to be the chemical potential which droplet 1 "sees" at infinity. Note that in d = 2, there will be a logarithmic correction but this does not change the phenomenology. Similarly, we assume droplet 2 to be stationary, and thus to be formed by the densities φ ± (R 2 ), but to "see" the outer density of droplet 1 at infinity. The time evolution of R 2 (t) is then: It is now important to recall that we were actually able to compute φ ± (R) in Section II, see Fig. 3. Using these values, we can solve (E1,E2) simultaneously to obtain R 1 (t) and R 2 (t). As stressed in the main text, these predictions where compared to direct numerical simulation of AMB+ in Fig. 4(b), finding good qualitative agreement. Quantitative agreement is not expected, since the calculation leading to (E1,E2) assumes an effective medium picture, whereby the presence of the second droplet breaks the spherical symmetry in the neighborhood of the first droplet. Moreover, it was assumed that the two droplets are separated by a distance R 1,2 , a condition which is not really met in simulations due to computational limitations. Finally, for a precise quantitative comparison, one shall take care of the fact that simulations are performed in d = 2 with periodic boundary conditions.

Contribution due to quorum sensing
We first consider the quorum sensing part and give the explicit expressions of f qs , K qs , K qs 1 and λ qs in terms of microscopic parameters.
We first expand v[ρ] at lowest order in gradients and retain only terms allowed by rotational symmetry: v[ρ] = v(ρ)+β(ρ)|∇ρ| 2 +γ(ρ)∇ 2 ρ. Moreover, we Taylor expand v(ρ), β(ρ), γ(ρ) with the notation where higher order terms are irrelevant at the order considered. From (39) and (G1), we obtain Expanding in powers of ρ leads to (41,42), where the local part of the free energy is f qs = ρ(log ρ − 1) + 1 2 ρ dϕ log(T + τ v 2 (ϕ)) (G3) and the other parameters are set by 2. Contribution due to the two-body potential U We finally compute the effect of U to get the form of µ r in (43) and derive (44)(45). We first expand the convolution U ρ in gradients of the density where B, C are given in (45,46). It is then easy to see that (43) holds, with ζ cg given in (44) and where we have defined (G12)