Stability phase diagram of active Brownian particles

Pin Nie,1,2 Joyjit Chattoraj,1,3 Antonio Piscitelli,1,4 Patrick Doyle,2,5 Ran Ni,6,* and Massimo Pica Ciamarra 1,4,† 1School of Physical and Mathematical Science, Nanyang Technological University, Singapore 2Singapore-MIT Alliance for Research and Technology, Singapore 3Institute of High Performance Computing, Agency for Science, Technology, and Research, Singapore 4CNR-SPIN, Dipartimento di Scienze Fisiche, Università di Napoli Federico II, I-80126 Napoli, Italy 5Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA 6School of Chemical and Biomedical Engineering, Nanyang Technological University, Singapore


I. INTRODUCTION
Many biological and synthetic systems of self-propelled particles exhibit a transition from a homogeneous state to one in which a gas-and a liquid-like phase coexist [1][2][3]. While diverse physical processes might be responsible for the observed transition, the bare presence of motility is enough to induce it [4]. Indeed, motility-induced phase separation (MIPS) occur in systems of particles whose interactions are purely repulsive and do not promote the alignments of the self-propelling directions. The prototypical simulation model is the active Brownian particle (ABP) model, which consists of spherical self-propelled particles interacting via excluded volume forces, and subject to thermal noise [1,[5][6][7][8][9].
In active systems, two particles colliding head-to-head, or nearly so, severely slow down their motion, reducing the local pressure. This pressure drop may seed a positive-feedback mechanism leading to the formation of a dense cluster of active particles, and hence to phase separation [10,11]. A similar scenario occurs in granular systems, where the pressure drop is due to the dissipative nature of the interparticle collisions [12,13]. A homogeneous system of active particles is not always unstable toward phase separation, as there are physical processes that promote the homogeneous phase, opposing the above instability mechanism. The balance between the mechanisms promoting phase separation and those promoting the homogeneous phase sets the limit of stability of the homogeneous phase in the motility-density plane.
The rotational diffusion plays a role because, before a collision seeds the growth of a cluster, the two colliding particles may change their self-propelling direction and swim away [4,7,[14][15][16]. This process gives rise to a flux of particles from the dense to the less dense phase promoting the homogeneous phase. By balancing this flux and the reverse flux of particles migrating toward the denser phase, which is controlled by the activity, Redner et al. [7,16] predicted a low-density coexistence line of ABPs in good agreement with numerical results but did not predict the location of the critical point and the upper coexistence line. The other possibility is that the translational rather than the rotational noise promotes the homogeneous phase. This scenario is suggested by a continuum equation for the evolution of the coarse-grained density and polarization fields [1,8,[17][18][19], which is formally related to a thermodynamic approach aiming to map active Brownian particles into an equilibrium system [11]. This scenario predicts a U-shaped spinodal line in the activity-density plane resulting from a diffusive instability. The prediction correctly reproduces the divergence of the lower spinodal line at a finite density, as well as the existence of a critical point. However, this approach underestimates [4] the minimum value of the activity at the critical point by a factor 10. The limitations [4] of theoretical approaches based on the rotational rather than on the translational diffusivity in predicting the phase diagram of ABPs suggests that additional processes promoting phase separation may exist.
In this paper, we demonstrate a physical process promoting the homogeneous phase driven by the motility of the particles. Motility, therefore, promotes and opposes phase separation at the same time. We formalize this and the other mechanisms promoting and opposing phase separation in a collisional framework and predict the spinodal line of ABPs. Our prediction favorably compares to both two-and threedimensional numerical simulations, for different values of the control parameters. Furthermore, we demonstrate that friction tunes the features of the motility-induced phase diagram, as it controls the instability mechanisms we have uncovered.

II. KINETIC MODEL
We develop a kinetic model to predict the spinodal line of ABP particles. In this model, the dynamics is described by the following overdamped equation of motion: Here sD r and D t = D r σ 2 /3 are the rotational and the translational diffusion coefficients, γ r = γ σ 2 3 , η is a Gaussian white noise variable with η = 0 and η(t )η(t ) = δ(t − t ), F a is the magnitude of the active force acting on the particle, and F i = f i j are the forces arising from the interparticle interactions. In the absence of interaction and noise, particles move with velocity v a = F a /γ , and do not rotate. The control parameters are the volume fraction φ and the Peclet number Pe ≡ v a D r σ = v a σ 3D t , with σ the average particle diameter. For Brownian spheres, s = 1 in Eq. (2). We develop our theoretical model for arbitrary values of s, to allow for a stringent numerical test of our theoretical predictions.
We theoretically determine the spinodal as the limit of stability of a homogeneous system toward the growth of density fluctuations. Within the spinodal region, the system is unstable as the diffusivity is D < 0, so that the flux of particles induced by a concentration gradient, J = −D∇ρ, enhances the gradient in a positive feedback mechanism. The limit of stability of the homogeneous phase is thus D = 0, or equivalently J = j g − j s = 0, where j g = ρτ −1 g and j s = ρτ −1 s are the particle fluxes promoting and suppressing density fluctuations, respectively.
The flux of particles promoting the growth of density fluctuations results from the interparticle collisions inducing a sensible drop in the local pressure. In a homogeneous system, only long-lasting interparticle collisions resulting from headto-head collisions induce such a drop, as we demonstrate in the Appendix B. Most of the collisions, therefore, do not promote density fluctuations but rather slow down the particles, endowing them with an effective velocity, v e . Collisions promoting phase separation have, therefore, a frequency τ −1 g ∝ φv e . Previous results have demonstrated that in the homogeneous phase, the effective particle velocity is v e = v a (1 − φ φ * ). This density dependence has been related to the pair-correlation function anisotropy [17], and rationalized in terms of the collision rate [19]. The estimation of the effective velocity allows that of the typical inverse agglomeration time, where a is a constant of order 1, and hence allows the estimation of j g . In the limit of stiff particles φ * correspond to the close packing volume fraction. The inverse agglomeration time vanishes for φ → 0, due to the absence of nearby particles, as well as for φ → φ * . In this limit particles are stuck, and particles self-propelling in opposite directions are not able to meet and promote a density fluctuation. We remark that our estimation concerns the agglomeration time in a homogeneous system at volume fraction φ. This time differs from that needed, in a phaseseparated state, by a gas particle to join a cluster [4,7,[14][15][16]. In particular, in the gas phase far from the critical point the volume fraction dependence of v e is negligible.
The flux of particles promoting phase separation is contrasted by fluxes promoting the homogeneous phase. One of these fluxes is driven by the rotational diffusivity of the particles, as two colliding particles might resolve their collision rotating their self-propelling direction. This physical process is the same allowing, in a phase-separated state, particles on the rim of an active cluster to escape from it [4,7,14,15]. The inverse timescale of this rotational detaching mechanism is with b a constant of order 1. In principle, particles may also resolve their collision by diffusing away, giving rise to a flux promoting phase separation driven by the translational diffusivity. However, the phase diagram of ABPs appears [1,18] insensitive to D t , for reasons we rationalize later on. We therefore do not consider any stabilizing flux associated with D t . The balance of τ g and τ rd allows us to predict a spinodal line, Pe ∝ φ −1 (1 − φ/φ * ) −1 . This prediction captures the numerically observed U shape of the spinodal line, and the divergence of Pe in the φ → φ * limit. However, according to this prediction, the spinodal line also diverges in the φ → 0 limit; hence, regardless of the volume fraction, a homogeneous system should become unstable and phase separate as Pe increase. Conversely, previous results indicate that only a system with volume fraction larger than a threshold φ m becomes unstable. There is, therefore, an additional stability mechanism in ABPs, which should be relevant at high Pe and small φ.
We identify this mechanism considering that two particles may resolve their collision without any change in the orientation of their self-propelling direction, but rather by sliding off each other [16,20], effectively rotating around their center of mass. This mechanism is illustrated in the top row of Fig. 1. This physical process is promoted by the activity, which sets the scale of the particle velocity. In a crowded environment, this process is hindered by the density, which slows down particle motion. We therefore assume particles to slide off each other with a velocity proportional to the effective active velocity. Thus, the inverse timescale associated with this sliding-detaching mechanism is where c is a constant of order 1. By balancing j g = ρτ −1 g and j s = ρτ −1 with φ m = c a and A = bφ * a . Given that a, b, and c are of order 1, so are the φ m and A, in both 2D and 3D. The critical point is at The prediction of a vertical spinodal line in the absence of rotational motion, s = 0, agrees with previous investigations [1,18].

III. FIXING φ * AND φ m
Our theoretical prediction of Eq. (6) has three free parameters. We have, however, independently estimated both φ * , the jamming volume fraction, as well as φ m , the spinodal line in the absence of rotational noise, for the numerical model we consider in the following. We find φ * 0.879 (0.645) in 2D (3D), and φ m 0.25 (0.345) in 2D (3D), as we detail below. Hence, we are left with a theoretical prediction with a single free parameter, the scaling amplitude A.

A. Jamming volume fraction, φ *
The jamming volume fraction φ * is a protocol-dependent quantity, which activity pushes toward its maximum value. To determine φ * , we cyclically compress and decompress the system across the expected jamming transition, in the absence of motility and noise. The volume fraction varies in steps of 10 −3 , and the energy of the system is minimized after every change of volume fraction, via the conjugate gradient protocol. Figure 2 illustrates subsequent compression curves, in both two and three dimensions. The pressure converges after a few cycles to an asymptotic curve, which grows linearly for φ > φ * . From these results, we estimate φ φ * 0.879 in two spatial dimensions, and φ * 0.648 in three dimensions.

B. Phase separation in the absence of rotational noise
In the absence of rotational noise, s = 0, our theoretical model predicts the spinodal line to be Peclet independent, φ = φ m . Previous numerical results have investigated this limit, confirming this theoretical prediction [1,18]. The investigation of the stability phase diagram for s = 0 thus allows estimating φ m .
We report our numerical results in Fig. 3, for both twoand three-dimensional systems. We recovered the Peclet independence of the spinodal line. Deviations from the theoretical predictions occur at small Peclet number, as in this limit thermal diffusivity start being relevant.
From this investigation, we estimate φ m 0.25 in two spatial dimensions, and φ m 0.345 in three spatial dimensions.

IV. NUMERICAL VALIDATION IN TWO AND THREE SPATIAL DIMENSIONS
We investigate the phase diagram of ABPs, whose dynamics is governed by Eqs. (1) and (2). We work in the hard-sphere limit, using the interparticle interaction and the parameter detailed in Appendix A, and consider systems with up to N = 32 000 particles in 2D, and up to N = 64 000 in 3D. We determine the coexistence line evaluating the position of the peaks of the local density distribution, as discussed in Appendix C 1. Investigating the dynamics of phase separation [7,19], and specifically the time dependence of both the fraction of the volume in the low-density phase and the characteristic size of the density fluctuations, we identify the state points undergoing spinodal decomposition. See Appendix C 2 for details on the phase separation dynamics at different state points. Figure 4(a) summarizes our results for standard (s = 1) ABPs in two dimensions (2D). In the figure, stars identify the coexistence line, squares points phase separating via spinodal decomposition, and circles points where either nucleation or no phase separation occurs. This phase diagram is consistent with those previously reported in the literature as concerns the critical value of the Peclet number, the typical values of the volume fractions, the shape of the coexistence line, and the location of the spinodal region. Our theoretical prediction of Eq. (6), represented as a full black line, correctly delimits the spinodal region. Analogous results for the three-dimensional (3D) case are in Fig. 4 While these results support our model, they do not clarify whether our approach or the continuum one better captures the physics of ABPs. Indeed, a U-shaped spinodal line, which at low density diverges at a finite density, has also been predicted within a continuum description. In this approach the coarsegrained density is found to evolve according to a diffusion equation with effective diffusivity [1,8,17,18] where v cg is the coarse-grained velocity along the polarization direction, so that D = 0. Assuming [1,8,17,18] the coarsegrained velocity to behave as the effective single-particle , and interpreting D as density-independent particle diffusivity, this approach predicts a spinodal line Pe cont we illustrate this in Fig. 4(a). This parameter-free prediction largely underestimates the critical Peclet number, as previously noticed [4]. However, treating φ * [18] (or φ m [8,9]) as a free parameter, and allowing for a scale factor possibly associated with the density dependence of the particle diffusivity, the prediction of the continuum model becomes comparable to ours, for s = 1.
The two theoretical predictions, however, differ as concerns the dependence of the spinodal line on the rotational diffusivity, sD r in our formalism. Indeed, the spinodal line scales linearly with s, according to our theoretical prediction of Eq. (6), while it scales as s 1/2 according to the continuum model, Eq. (9). This consideration makes compelling the investigation of the s dependence of the phase diagram that we have performed in 2D. We compare the s = 1 phase diagram with those obtained for s = 1/2 and s = 10, in the (Pe/s, φ) plane, in Fig. 5. In the figure, open symbols correspond to s = 1, full ones to s = 1, and the full black line is as in Fig. 4(a). The phase diagrams in the (Pe/s, φ) plane are almost indistinguishable. This finding indicates that the spinodal line scales linearly with s, and strongly supports our theoretical model. Furthermore, we notice that the continuum approach predicts a vertical phase boundary in the absence of translational noise [18], regardless of the rotational noise. Conversely, a vertical phase boundary occurs with no rotational noise [1,18], as predicted by Eq. (7).

V. ROTATIONAL VERSUS TRANSLATIONAL DIFFUSIVITY
Our model, which neglects the role of translational diffusivity, successfully rationalizes the motility-induced phase diagram. We rationalize why the translational diffusivity is irrelevant and the limit of validity of this result, by comparing the diffusion coefficient of the passive suspension to the activity-induced effective diffusion coefficients. For the diffusivity of the passive suspension, we find (2D) D p (φ) = D t (1 − φ/φ d ) with φ d φ * , in the volume fraction range we have considered. This result is consistent with the expectation for the low-density behavior of Brownian particles [21]. We associate two diffusion coefficients to the active suspension, FIG. 6. Peclet number dependence of rescaled diffusion coefficients associated with the motion parallel and perpendicular to the self-propelling direction of each particle. by describing particle motion in the homogeneous phase as resulting from a sequence of steps alternatively taken from distributions corresponding to two different stochastic processes, describing the motion in between collisions and during a collision.
The stochastic process describing motion in between collisions is that of a persistent random walk, with persistence time 1/sD r . The corresponding diffusivity D is evaluated, following previous works [1,19], considering the steps to where t c and t mf are the mean duration of a collision and the mean time between collisions. At low density In the above equation, we have taken into account the contribution of the thermal diffusivity, which is divided by a factor d accounting for the fact that D is effectively a one-dimensional diffusivity. The stochastic process describing the motion resulting from the steps performed during the collisions is that of a simple random walk, with step size ∝ σ and step frequency 1/(t c + t mf ) ∼ 1/t mf ∝ v a φ/σ , so that We numerically validate these theoretical predictions by decomposing the instantaneous velocity of particle i in the components parallel and perpendicular to its self-propelling The time integration of these velocities defines a normal and a tangential displacement, dt, from which we estimate the diffusion coefficients, D ,⊥ = lim t→∞ r 2 ,⊥ (t ) /2t. Figure 6 shows that these numerical estimates compare well with the theoretical predictions. The theoretical predictions work well at Pe 1. Importantly both diffusion coefficients, and in particular D ⊥ which describes a physical process promoting the homogeneous phase, grow with the Peclet number, and are much larger than the diffusivity of the passive suspension. This result explains why the diffusivity of the passive suspension does not influence the motility-induced phase diagram.
However, our theoretical prediction fails at small Pe, where the collisional description of the dynamics is no longer appropriate. Hence, the thermal diffusivity may influence phase separation if the critical point is at Pe < 1, which may occur at very small values of s. In this limit, our theoretical prediction breaks down. Indeed, for s = 0 our model predicts a vertical phase boundary, and hence a motility-induced phase separation, also in the limit of vanishing motility. With no motility, however, no MIPS occurs and the system behaves as a thermal one [22]. Understanding how the U-shaped spinodal line evolves into a vertical phase boundary as the rotational diffusivity decreases is an interesting avenue of research we leave for the future.

VI. FRICTIONAL CONTROL OF THE PHASE DIAGRAM
The sliding-detaching mechanism occurs when two colliding particles resolve their collision without their selfpropelling directions rotating. Interparticle interactions which induce a torque on the particles, therefore, suppress this mechanism, as illustrated in Fig. 1, bottom row. These interactions may result from the shape of the particles, if these are elongated [23,24], and from lubrication forces in wet systems [25]. In dry systems, frictional forces also induce torques and hence suppress the sliding-detaching mechanism in dry active matter. Recent results have demonstrated that frictional forces are also present in colloidal hard-sphere suspensions at high enough shear rates, where they give rise to the discontinuous shear thickening phenomenology [26][27][28][29]. This frictional dependence is rationalized assuming that the hydrodynamic layer, which would give rise to diverging normal forces, breaks down at a characteristic length set by the particle surface asperities. Hence, frictional forces may play a role also in experiments of wet colloidal scale active particles, at high enough Peclet number. The tunability of the frictional interaction in colloidal systems [28], therefore, may allow us to control the motility-induced phase separation of these systems.
We investigate the influence of static friction on the motility-induced phase separation in three-dimensional numerical simulations, resorting to the frictional Mindlin model, as described in the methods section. We find Coulomb's friction coefficient μ to influence the lower spinodal line, which is critically affected by the sliding-detaching mechanism, not the upper spinodal line. At each value of the Peclet number, a homogeneous system becomes unstable toward phase separation at a volume fraction φ s (Pe, μ) which exponentially decreases with μ, approaching a limiting value, as in Fig. 7(a). Consistently, the spinodal region widens on increasing the friction coefficient and the Peclet number, as in Fig. 7(b).
The combined effect of friction and Peclet number is rationalized considering that the frictional forces, whose strength scales as μv a ∝ μPe, can be disrupted by thermal ones, which have a constant magnitude, through an activated process. This activated dynamics naturally explains the exponential decay of φ s (Pe, μ). The activation probability decreases as the Peclet number increases, so that the spinodal φ s (Pe, μ) approaches φ = 0, at all finite values of the friction coefficient, making the frictionless limit a singular one. These results demonstrate that friction modulates the phase diagram of ABPs, possibly guiding future experiments. Importantly, this modulation occurs as friction suppresses the sliding-detaching mechanism, indirectly confirming the relevance of this mechanism.

VII. CONCLUSIONS
While not representing any experimental system faithfully, active Brownian particles emerged as a prototypical model exhibiting a motility induced phase separation, and are the standard benchmark for statistical physics of active matter theories. In ABPs, motility promotes phase separation, as clarified by the MIPS [4] or, equivalently, by a mechanistic approach [11]: due to the presence of motility, colliding particles are much slower than the others, so that collisions induce a local pressure drop which seeds phase separation. Previous theoretical approaches argued that stochastic processes related to the translational [1,8,[17][18][19] or rotational [7,14,16] diffusivity conversely promote the homogeneous phase.
Here we have demonstrated that besides these previously identified mechanisms, there exists a physical process induced by motility that promotes the homogeneous phase. This process dominates over the others at high motility. In this limit, the balance of two motility-driven processes fixes the spinodal line, which thus becomes motility independent. In the opposite limit of small motility, the rotational diffusivity becomes relevant, and the phase diagram becomes motility dependent. We have formalized these mechanisms in a kinetic approach and predicted the spinodal line of ABPs up to a scaling amplitude of order 1. The estimation of this constant remains an open problem.
We have explicitly checked that in standard ABPs the translational diffusivity is negligible with respect to diffusivities induced by the collisions. This explain why the translational diffusivity does not affect the phase diagram. This scenario, however, certainly changes in the limit of small rotational diffusivity, where thermal and activity-induced effects may compete.
The sliding-detaching mechanism we have uncovered involves the coordinated motion of colliding particles. Coordinated motion is, by definition, suppressed in active Ornstein-Uhlenbeck [30] and Monte Carlo models [31]. Accordingly, for these models we expect that the spinodal line does diverges in the φ → 0, as observed. Friction also suppresses the sliding-detaching mechanism, by inducing the rotation of the self-propelling directions of the particles. Our investigation of the effect of friction indicates that this could be used to control the features of the motility-induced phase diagram, also in light of recent experimental results [28]. The frictional dependence is qualitatively rationalized considering that frictional forces scale as μPe. The coexistence of this force scale and of stochastic forces then makes the dynamics of frictional systems an activated one. The quantitative rationalization of the effect of friction on the phase diagram remains, however, an open problem. ACKNOWLEDGMENTS P.N., J.C., and M.P.C. acknowledge support from the Singapore Ministry of Education through the Academic Research Fund (Tier 2) MOE2017-T2-1-066 (S) and from the National Research Foundation Singapore, and are grateful to the National Supercomputing Centre (NSCC) of Singapore for providing computational resources.

APPENDIX A: NUMERICAL DETAILS
We perform numerical simulations of ABPs in two and three spatial dimensions. We use an interparticle interaction model borrowed from the granular community, to model frictional particles; the frictionless model is recovered setting to zero Coulomb's friction coefficient, μ.
The interparticle interaction force has a normal and a tangential component, F i j = f n i j + f t i j . The normal interaction is a purely repulsive Harmonic interaction, f n i j = k n (σ i j − r i j ) (σ i j − r i j )r i j , (x) is the Heaviside function, σ i j = (1/2)(σ i + σ j ), r i j = r i − r j , and r i is the position of particle i. The tangential force is f t i j = k t ξ i j , where ξ i j is the shear displacement, defined as the integral of the relative velocity of the interacting particle at the contact point over the duration of the contact, and k t = 2 7 k n . In addition, the magnitude of tangential force is bounded according to Coulomb's condition: We work in the hard-sphere limit considering stiff particles, the maximum relative deformation of a particle being 10 −4 for the range of parameters we have considered. This makes our results insensible to k n , but our numerical investigation more computationally costly than previous ones.
For the frictionless case, μ = 0, and the equation of motion is as in Eqs. (1) and (2). In the presence of friction, a torque T i γ r arising from the frictional interparticle interaction is added to Eq. (2).
Simulations [32] are done integrating the equation of motion via the overdamped Langevin algorithm, with integration time step 2 × 10 −8 D −1 r .

APPENDIX B: AGGLOMERATION TIMESCALE AND EFFECTIVE PARTICLE VELOCITY
Our kinetic model requires the estimation of the agglomeration timescale τ g , which is the average time a particle waits before experiencing a collision promoting the formation of a cluster. This timescale depends on the typical particle velocity. We have assumed this to be an effective density-dependent velocity, rather than the bare particle velocity v 0 . Here we provide data supporting this assumption.
If all collisions promote the formation of a cluster, then the timescale of interest is the mean free time which, in the ABP context, depends on the density and the bare particle velocity v 0 . In ABPs, however, only the long-lasting collisions are able to sensibly slow down particle motion, hence reducing the local pressure and possibly seeding the formation of a cluster.
To evaluate how may collisions could potentially lead to the formation of a cluster, we have investigated the dynamics of two-particle collisions, as a function of the impact parameter b and of the relative angle of the self-propelling directions of the two particles, θ . The inset of Fig. 8(a) defines these quantities. In these simulations, there is no rotational noise.
Collisions able to sensibly slow down particle motion are those with large t coll and small d. Hence, we consider the ratio t coll /d as a proxy of how likely a collision acts as a cluster seed. Panel (c) shows that this clustering ability is strongly FIG. 11. Evolution of the fraction of the overall volume with a density smaller than the low-density coexistence density, for an N = 32 000 particle system, in two dimensions. In panel (a) φ = 0.35, while in panel (b), φ = 0.4. In both panels, Pe = 20. At φ = 0.35 (a) the low-density fraction only increases after a period of transition. This indicates that phase separation occurs via nucleation. At φ = 0.4 (b) the fraction grows quickly at early times, indicating that the phase separation proceeds via spinodal decomposition. Afterward, the system coarsens. The growth of the length scale of the density fluctuations is compatible with the expected λ ∼ t 1/4 law, as illustrated in the inset. The bottom panels illustrate maps of the coarse-grained density distribution, which consistently suggests that phase separation occurs via different processes at the different volume fractions. Note the different timescales. Since the clustering ability sharply peaks around some characteristic values of the impact parameters, only a tiny fraction of all collisions induce the formation of a cluster. Because of this, in between two collisions promoting the formation of a cluster, a particle experiences many other collisions. These collisions slow down the particle motion endowing the particles with an effective velocity. We argue that this effective velocity, rather than the bare one, sets the agglomeration timescale.

APPENDIX C: DETERMINATION OF THE SPINODAL REGION
The spinodal line is a mean-field concept, and in finite systems, the separation between nucleation and spinodal decomposition is not sharp. Nevertheless, on increasing the system size, the crossover between the two different phase separation mechanisms allows for a meaningful operative identification of the state points where phase separation occurs via spinodal decomposition.
To identify the state points within the spinodal region, we have first identified those that phase separate after relatively short simulations, investigating the distribution of the coarsegrained density, as described next. Then, for a relevant subset of those points, we have investigated the dynamics of phase separation, to distinguish between nucleation and spinodal decomposition, as described in Appendix C 2.

Coarse-grained density
We determine whether a system is homogeneous or phase separated investigating the probability distribution of the coarse-grained density, ρ cg (r), or equivalently of the coarsegrained volume fraction, φ cg (r) = ρ cg (r) v , with v the average particle volume. Following Ref. [33], we define ρ cg (r) by convoluting the number density i δ(r − r i ) with f (r) = Z exp[−1/(1 − r 2 /w 2 )], with w = 3.5σ and Z a normalization factor. Figure 9 illustrates that in simulations starting from a homogeneous configuration, ρ cg (r) evolves until it converges to a steady-state distribution, which in the figure is bimodal. In the range of parameters we have considered, convergence occurs within 100D −1 r . We consider a point in the (φ, Pe) plane to be phase separated if the probability distribution of the local density is bimodal, and if doubling the system size yields consistent results. In two dimensions, we consider systems with N up to 32 000, in three dimensions with N up to 64 000. Figure 10 illustrates example distributions.
We evaluate the coexisting densities by locally approximating the distribution via Gaussian functions. The results of these fits are illustrated in the figure. The coexisting volume fractions are in Fig. 1 of the main text.
Note that while the two panels of Fig. 10 refer to volume fractions which are close to the critical ones, the height of the two peaks is sensibly different. This difference occurs as the coexistence curve is extraordinarily flat and asymmetrical close to the critical point, in particular in three dimensions.

Dynamics of phase separation
We investigate the dynamics of phase separation to rationalize whether a system undergoes separation via spinodal nucleation rather than via nucleation. First, we consider the time dependence of the percentage of the total volume with local volume fraction smaller than that of the coexisting low-density phase, α(t ) = V [φ cg φ coex (t )]/V tot . If phase separation proceeds via spinodal decomposition, then α(t ) quickly varies at short times, the homogeneous system being unstable. Conversely, if phase separation occurs via nucleation, then α(t ) only starts varying after a period of transition, corresponding to the nucleation time. Figure 11 illustrates the result of this investigation, in two dimensions. Panels (a) and (b) refer to different state points that have the same Pe and differ in volume fraction by φ = 0.05. At φ = 0.35, phase separation is seen to occur via nucleation, while at φ = 0.4, it occurs via spinodal decomposition. The associated snapshots of the local density field confirm this interpretation.
Then, we investigate the time evolution of the typical length scale λ of the density fluctuations. We define λ as the distance at which the correlation function of the coarsegrained density first becomes zero. The inset of panel (b) shows that λ grows as a power law with time at φ = 0.4 and Pe = 20. The growth exponent is compatible with 1/4, the expected exponent for the coarsening exponent in twodimensional systems with conserved order parameter. This result further supports our interpretation, namely that at φ = 0.4 and Pe = 20 the system phase separates via spinodal decomposition.
We provide more examples of systems phase separating via spinodal decomposition, for state points close to our identified spinodal line, in both two and three dimensions, in Fig. 12.