Active Brownian Equation of State: Metastability and Phase Coexistence

As a result of the competition between self-propulsion and excluded volume interactions, purely repulsive self-propelled spherical particles undergo a motility-induced phase separation (MIPS). We carry out a systematic computational study, considering several interaction potentials, systems confined by hard walls or with periodic boundary conditions, and different initial protocols. This approach allows us to identify that, despite its non-equilibrium nature, the equations of state of Active Brownian Particles (ABP) across MIPS verify the characteristic properties of first-order liquid-gas phase transitions, meaning, equality of pressure of the coexisting phases once a nucleation barrier has been overcome and, in the opposite case, hysteresis around the transition as long as the system remains in the metastable region. Our results show that the equations of state of ABPs quantitatively account for their phase behaviour, providing a firm basis to describe MIPS as an equilibrium-like phase transition.

Active matter made of self-propelled particles can be found in a wide variety of contexts, both in living and synthetic systems, ranging from cell populations to bacteria suspensions, animal groups or colloidal artificial swimmers [1]. The fundamental difference between active and 'passive' matter made of thermally agitated constituents, is that the microscopic dissipative dynamics of the former breaks detailed balance and, as such, evolves out-ofequilibrium. The intrinsic out-of-equilibrium nature of active matter manifests strikingly in the presence of different kinds of interactions. Self-propelled spherical particles accumulate in regions of space where their velocity decreases as a consequence of collisions (between particles or with external obstacles). Simple models that capture in a minimal way the competition between self-propulsion and steric effects, like the so-called Active Brownian Particles (ABP) model, have provided much insight into the generic behavior of such systems. For instance, at high enough densities and activities, a purely Motility-Induced Phase Separation (MIPS) generically takes place, leading to the coexistence of an active low density gas with a high density drop in the absence of attractive forces [2][3][4][5][6][7]. This out-of-equilibrium transition is reminiscent of equilibrium liquid-gas de-mixing. It is thus tempting to extend the thermodynamic description of first order phase transitions in terms of, for instance, equations of state, to ABPs. However, this poses several fundamental difficulties since no thermodynamic variable is, in principle, well defined in this context. Much effort has been recently devoted to this question from a statistical mechanics perspective: the notions of effective temperature [8][9][10][11][12] and chemical potential [13] have been introduced, and special attention has been paid to the notion of pressure in active fluids [14][15][16][17][18][19][20]. In particular, it has been shown that an equation of state exists for isotropic ABP.
In spite of the above-mentioned efforts, the equality of pressure at coexistence implied by the very existence of an equation of state has not been observed so far (neither in experiments nor simulations) and the possibility to construct from the equations of state an effective thermodynamic description of MIPS is still a matter of debate. In order to move forward in our fundamental understanding of active systems, it is crucial to provide a consistent interpretation of pressure measurements that allows to confront them with the theory. In this Letter we clarify this issue and provide a full description of the phase behavior of ABPs in terms of its equations of state.
The equations of state of self-propelled Janus colloids have been experimentally measured in the absence of phase coexistence [17]. Several numerical studies have measured the pressure in systems undergoing MIPS, but the way to interpret the results raises some serious conceptual problems [14,18,19]. For instance, an abrupt pressure drop at the vicinity of MIPS has been reported recently and it has been argued that it constitutes a distinctive feature of this transition, hindering the analogy with equilibrium phase separation [18,19] and in apparent contradiction with theoretical predictions [16,21]. A pressure loop generically appears in phase separating finite systems in equilibrium, an effect that can be suppressed by the Maxwell construction to extract the thermodynamic behavior. However, this construction is violated for ABP [16,21], so there is no direct way to understand the phase diagram of the system from the simulated equations of state.
We show here that the equations of state of ABP across MIPS are consistent with the first-order phase transition scenario and allow to characterize its phase behaviour, providing a consistent thermodynamic interpretation of the numerical data. We bring out the existence of a metastability region: an hysteresis around the coexistence pressure is found as the system is quenched to the coexistence region from 'below' (from a homogeneous state) or 'above' (from a phase separated state) MIPS.
We analyze both open and confined systems to show that all the conflicting results found in previous simulation studies are due to the presence of a large nucleation barrier that can be easily bypassed by including a nucleation core (a wall). Then one can generate a motility induced phase separated state at densities well below those previously reported. To show these results, we perform constant-density simulations (N V T ensemble) of Active Brownian Hard-Disks and ABP interacting with a WCA potential.
To be specific, we consider N disks of diameter σ in a 2d volume V = L x × L y , with φ = πσ 2 N/4V the packing fraction. The model is defined by the following equations of motion for each particle at position r i = (x i , y i ) and with orientation n i (t) = (cos θ i , sin θ i ): where ξ i and ν i are zero-mean unit-variance Gaussian noises, and v 0 is a constant self-propulsion velocity. The force F i acting on particle i comes from inter-particle interactions, 2 j =i f ij , and external potentials, F ext i = −∇ i w(r i ). The units of length and time are given by σ and τ = D −1 θ , respectively, and fix D θ = 3D 0 /σ 2 . Together with the packing fraction, we identify two dimensionless parameters that control the phase behaviour of our system: the Péclet number, Pe = v 0 /σD θ , which quantifies the strength of selfpropulsion, and the effective particle stiffness Γ = µ σv0 , which naturally arises from the potential energy scale (in units of k B T = µ/D 0 = 1). This parameter quantifies to what extent particles become effectively softer as their activity increases, an effect surprisingly disregarded in the literature.
In order to focus on the essential features determining the phase behavior of the system, we compare a suspension of Active Brownian Disks with a Weeks-Chandler-Andersen potential (AB-WCAD): u(r) = 4 σ r 12 − σ r 6 + 1 4 with a upper cut-off at r = 2 1/6 σ and such that Γ > 0.03; with one composed by infinitely hard active particles, i.e. u(r) = 0 if r > σ and u(r) = ∞ otherwise [22]. The hard sphere fluid is the most studied model in liquid theory: it was the first fluid to be simulated with molecular dynamics [23] and several analytical developments can be carried to predict its thermodynamic properties [24][25][26]. In both cases we carry out BD simulations to study their evolution in the absence of any external forcing. It is not straightforward to deal with the singular nature of hard core interactions in Brownian dynamics (BD) simulations, commonly used in the context of ABP. We therefore use a variant of Evendriven Brownian Dynamics (ED-BD) for hard spheres [27] that accounts for self-propulsion [28] [39].
The phase diagrams shown in Fig. 1 can be readily constructed from the coexisting densities. ABP with different stiffness exhibit the same qualitative phase behaviour: above a critical Péclet number they undergo MIPS if its average density is large enough. Increasing the potential stiffness one observes both a decrease of the coexisting packing fractions and an increase in the structure, as depicted in Fig. 3 (d). This feature implies that the structure of the high-density phase is more compact for stiffer potentials. In the limit of hard disks, the high-density branch of the binodal quickly saturates at close packing, indicating that the competition between self-propulsion and excluded volume generates a nearly perfect hexagonal crystal at high enough activities, as shown in Fig. 3 (d) (e.g. for Pe= 60, the average density of the condensate is φ high = 0.98 φ cp , in coexistence with a low density phase at φ low = 0.089).
We identify the pressure, P , by an extension of the virial theorem that accounts for active forces [18,29]. After projecting eq. (1) on r i and averaging over the noise, we get P = ρk B T + P int + P a , where ρk B T is the thermal ideal gas pressure, is the standard virial expression of the collisional pressure, and is the active contribution to the pressure or 'swim' pressure, which can be mechanically interpreted as a body force [14,15,30]. This yields an ideal gas law The total pressure P defined this way is a state function for spherical ABP [16] and its definition can be extended to unconfined systems with PBC where P is equal to the internal bulk pressure [18].
Following usual practice, we study the phase behaviour of ABP by letting the system evolve from a homogenous initial state at some density φ towards its stationary state corresponding to a given value of Pe. This procedure corresponds to a quench from Pe = 0, below MIPS, to Pe > 0. However, in order to gain insight in the intrinsic nature of the pressure, we compare the values obtained for a periodic system and one confined in the x-direction by two walls at x = 0 and x = L x , and periodic boundary conditions (PBC) in the y-direction [40].
The P (φ; Pe) equations of state computed numerically using eq. (2-3) in confined and PBC geometries, following high-Pe quenches at fixed φ, are shown in Fig. 2 (a-c). If the system does not phase separate, for Pe<Pe c , the total pressure of the AB-HD fluid with PBC coincides with the pressure computed in confined conditions (see Fig. 2 (a)). At higher activities, for PBCs the pressure drops abruptly in the vicinity of MIPS [18,19], while in the confined system the pressure displays a monotonic increase with density and remains roughly constant in the coexistence region. As shown in Fig. 2 (a-b), the pressure jump becomes more pronounced as the activity and the stiffness of the particles increases (note the similarity between AB-HD with Pe = 30 and AB-WCAD with = 250 and Pe = 100). The formation of a macroscopic drop, due to MIPS, coincides with the location of the pressure drop, φ n (see [22]). Particles are dramatically slowed down as they aggregate, making v φ , and therefore the active pressure, to drop (see eq. (3)). The low-density branch of the binodal, φ low , is well below the density φ n at which the dense phase appears.
Above Pe c , and at low densities, below the onset of MIPS, the equations of state collapse into a single curve which is very well fitted by a first order virial expansion, P = P 0 (1 + B 1 φ) [18]. For a system with confining walls, B 1 decreases with Pe (see [22]) in consistence with the viewpoint that self-propulsion induces an effective particle attraction [2,12,17]. Surprisingly, for PBC, B 1 is independent of Pe, indicating that the activity-induced attraction saturates above Pe c . As we show below, this is due to the fact that in the absence of walls the system evolves along a metastable branch as the density is increased and remains in its pure low-density phase above φ low , until it reaches the MIPS threshold, φ n , when a dense cluster nucleates.
Interestingly, the equations of state in confined and PBC geometries converge at large densities to a value close to the expected coexistence pressure. In equilibrium, finite-size systems display a pressure loop, which, in contrast to a van der Waals loop, is thermodynamically stable and due to the formation of an interface between the two phases [31,32]. Accordingly, P (φ) shows a peak when the dense phase develops. In simulations, this discontinuous change is typically smoothed out by interface fluctuations. The sharp decrease in P for ABPs implies that interfacial fluctuations are considerably suppressed at high Pe, suggesting a large interface tension (AB-HDs might generate an interface between a dilute gas and a closed packed crystal). Hence, we can interpret the pressure drop for PBC systems as a finite size effect that should vanish in the thermodynamic limit. For finite systems the equation of state is affected both because of an interfacial contribution to the pressure and a shift in the location of nucleation [41]. These two effects are iconsistent with the data shown in Fig. 2 (c). The active pressure in systems of size N 4000 converges above φ n , meaning that there is no interfacial finite-size contribution to the pressure after phase separation. It is less clear though whether the pressure jump is reduced when increasing N [19]. We push thermodynamic ideas even further, and claim that the system sizes used so far are not large enough to observe how the nucleation approaches the binodal while the pressure drop vanishes. The idea behind this claim is that MIPS involves a large critical nucleus: our systems turn out to be too small to phase separate close to the binodal and thus remain metastable in a sub-region of the coexistence region [42]. The convergence of φ n towards the binodal as the system is made larger is due to the fact that the probability to spontaneously form a cluster of size m > m c , (where m c is the critical nucleus size) increases with system size. However, if N m c is not guaranteed, then we cannot expect to observe how P (φ) approaches its infinite-size behavior. In the following we show this important claim.
In order to probe metastability, we initialize the open system deeply in the coexistence region. As we expand the system nucleation is avoided since an interface between a dense and low density phase is already present from the very beginning of the simulation. We performe such a low-φ quench at fixed Pe by letting AB-HDs evolve from a steady-state at φ = 0.50, Pe = 30 to lower packing fractions. The equation of state we obtain is depicted in blue in Fig. 2 (d). We find that: (i) the system remains phase separated at much lower densities than the nucleation point found by high-Pe quenches (see Fig. 3 (a)); (ii) the pressure remains roughly constant down to φ ≈ φ low . We thus found hysteresis around MIPS, a typical signature of first-order phase transitions.
Accordingly to classical nucleation theory (CNT), in the absence of a preferential site (homogeneous nucleation), phase separation can only be triggered by a rare event: the spontaneous formation of a critical nucleus of size larger than m c ∝ γ/∆G homo (where γ is the tension interface and ∆G homo their free energy difference).
Since our system is intrinsically out-of-equilibrium, CNT cannot be directly applied. However, borrowing ideas from equilibrium systems, Redner et al. have developed a theory analogous to CNT to describe the kinetics of phase separation in ABP [33] which provides a good theoretical description of several simulation results [34]. Within this framework, the critical nucleus m c ∝ φ cp / ln 2 (Pe φ low ) where φ low is predicted by the theory. To make a comparison with simulations, an expression φ(A, φ high , φ low ) for the average packing fraction in terms of the theoretical predictions is established, where A is a structural parameter completely determined by the cluster size distribution P m (see eq. (11) in [33]). The nucleation barrier is thus controlled by two main ingredients: the location of the binodals and the structure of the clusters in the metastable region. We computed P m for AB-HD and AB-WCAD with different stiffness and found a roughly identical distribution (see Fig. 3 (c)), meaning that the nucleation barrier is mainly controlled by the location of the binodals. For stiffer potentials the structural difference between the coexisting phases is more severe. In the hard-disk limit a nearly perfect crystal coexists with a very dilute gas in the high activity regime, pushing the binodals to its extreme values, φ low → 0 and φ high → φ cp . Therefore, the critical nucleus is expected to be very large at high Pe. As shown in Fig. 3 (d), the crystalline order of the dense phase is suppressed by softening the particles. In turn, φ low is larger for softer potentials such that m c is made smaller than in the hard limit case, thus qualitatively explaining the reduction of the pressure drop for softer disks (see Fig. 2 (b)). To be specific, Redner et al. predict m c ≈ 5000 at Pe= 30 and φ = 0.30 for AB-WCAD with = 1 in our units (see Fig. 2 in [33]). This value lies quite close to the binodal in their case while for AB-HD it falls deep into the coexistence region since the lowdensity branch of the binodal is shifted to much lower packing fractions in the hard disk limit. We did not intend to make a direct quantitative comparison between the theory and our simulations, but, at this level, we are able to insure that for hard ABPs, MIPS features a large nucleation barrier that discourages attempts to observe how the system escapes from metastability with BD simulations. To give further support to this idea we turn now our attention into the system in presence of hard walls.
Confinement facilitates nucleation because of wetting. As shown in Fig. 3 (b), self-propelled particles accumulate at walls which thus act as natural nucleation seeds [12,20,35]. The free energy associated with (heterogeneous) nucleation in a confined system is ∆G het = F (α) ∆G homo , where F (α) = (2 + cos α)(1 − cos α) 2 /4 is a geometric function of the contact angle α between the wall and the dense phase. The adsorption of ABPs into layers give α = 0 (pure wetting) and therefore, by extension of CNT, a vanishing nucleation barrier ∆G het = 0. This equilibrium-like description is consistent with the absence of a pressure drop in the presence of walls and confirming our overall interpretation of the equations of state across MIPS (see Fig. 2 (b)).
We have carefully examined the pressure of ABPs using different potentials, topologies and preparation protocols in order to show that the equations of state are fully consistent with the classical (equlibirum) first-order phase transition scenario. The equations of state of ABP do not exhibit any fundamental difference to an equilibrium system showing phase coexistence, besides: (i) the absence of a Maxwell construction on P ; (ii) the extreme structural difference between the two coexisting phases, giving rise to a large nucleation barrier. Over-all, our work results on a systematic way to interpret the equations of state of ABP using equilibrium-like concepts. We quantitatively confirm the debated analogy between MIPS and equilibrium phase separation, and, as such our work should represent an important step towards the construction of an effective thermodynamic description of active systems.