Active Brownian Equation of State: Metastability and Phase Coexistence

TL;DR

This study demonstrates that the phase transition in active Brownian particles undergoing motility-induced phase separation exhibits equilibrium-like properties, including pressure equality and hysteresis, despite being a non-equilibrium system.

Contribution

It provides the first systematic computational evidence that ABP phase behavior mimics equilibrium first-order transitions, across various potentials and boundary conditions.

Findings

Pressure equality at coexistence confirms first-order transition.

Hysteresis observed in metastable regions.

Equations of state quantitatively describe phase behavior.

Abstract

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.

Equations

$${\dot{\mathbf{r}}}_i=v_0{\mathbf{n}}_i+\mu {\mathbf{F}}_i+\sqrt{2D_0}{\mathbf{\xi }}_i,{\dot{\theta }}_i=\sqrt{2D_{\theta }}{\nu }_i$$

$$P_{int}=\frac{1}{4V}\sum _i\sum _j⟨{\mathbf{f}}_{ij}\cdot ({\mathbf{r}}_i-{\mathbf{r}}_j)⟩$$

$$P_a=\frac{v_0}{2\mu V}\sum _i⟨{\mathbf{n}}_i\cdot {\mathbf{r}}_i⟩=\frac{\rho v_0{v}_{\varphi }}{2\mu D_{\theta }},v_{\varphi }=\frac{1}{N}\sum _i⟨{\mathbf{n}}_i\cdot {\dot{\mathbf{r}}}_i⟩$$

Full text

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 Marchetti et al. (2013). 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-of-equilibrium. 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 Tailleur and Cates (2008); Cates and Tailleur (2015); Fily and Marchetti (2012); Redner et al. (2013); Stenhammar et al. (2014); Levis and Berthier (2014). 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 Loi et al. (2008); Palacci et al. (2010); Szamel (2014); Levis and Berthier (2015); Marconi and Maggi (2015) and chemical potential Dijkstra et al. (2016) have been introduced, and special attention has been paid to the notion of pressure in active fluids Takatori et al. (2014); Yang et al. (2014); Solon et al. (2015); Ginot et al. (2015); Winkler et al. (2015); Patch et al. (2017); Speck and Jack (2016). 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 Ginot et al. (2015). Several numerical studies have measured the pressure in systems undergoing MIPS, but the way to interpret the results raises some serious conceptual problems Takatori et al. (2014); Winkler et al. (2015); Patch et al. (2017). 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 Winkler et al. (2015); Patch et al. (2017) and in apparent contradiction with theoretical predictions Solon et al. (2015, 2016). 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 Solon et al. (2015, 2016), 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 ($NVT$ ensemble) of Active Brownian Hard-Disks and ABP interacting with a WCA potential.

To be specific, we consider $N$ disks of diameter $\sigma$ in a 2$d$ volume $V=L_{x}\times L_{y}$, with $\phi=\pi\sigma^{2}N/4V$ the packing fraction. The model is defined by the following equations of motion for each particle at position $\bm{r}_{i}=(x_{i},y_{i})$ and with orientation $\bm{n}_{i}(t)=(\cos\theta_{i},\sin\theta_{i})$:

[TABLE]

where $\bm{\xi}_{i}$ and $\nu_{i}$ are zero-mean unit-variance Gaussian noises, and $v_{0}$ is a constant self-propulsion velocity. The force $\bm{F}_{i}$ acting on particle $i$ comes from inter-particle interactions, $\bm{F}^{int}_{i}=-\frac{1}{2}\nabla_{i}\sum_{j\neq i}u(r_{ij})=\frac{1}{2}\sum_{j\neq i}\bm{f}_{ij}$, and external potentials, $\bm{F}^{ext}_{i}=-\nabla_{i}w(r_{i})$. The units of length and time are given by $\sigma$ and $\tau=D_{\theta}^{-1}$, respectively, and fix $D_{\theta}=3D_{0}/\sigma^{2}$. Together with the packing fraction, we identify two dimensionless parameters that control the phase behaviour of our system: the Péclet number, $\mbox{Pe}={v_{0}}/{\sigma D_{\theta}}$, which quantifies the strength of self-propulsion, and the effective particle stiffness $\Gamma=\frac{\epsilon\mu}{\sigma v_{0}}$, which naturally arises from the potential energy scale $\epsilon$ (in units of $k_{B}T=\mu/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\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}+\frac{1}{4}\right]$ with a upper cut-off at $r=2^{1/6}\sigma$ and such that $\Gamma>0.03$; with one composed by infinitely hard active particles, i.e. $u(r)=0$ if $r>\sigma$ and $u(r)=\infty$ otherwise SM . The hard sphere fluid is the most studied model in liquid theory: it was the first fluid to be simulated with molecular dynamics Alder and Wainwright (1957) and several analytical developments can be carried to predict its thermodynamic properties Hansen and McDonald (1990); Baxter (1968); Carnahan and Starling (1969). 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 Even-driven Brownian Dynamics (ED-BD) for hard spheres Scala et al. (2007) that accounts for self-propulsion Ni et al. (2013)111The constraint that two particles cannot overlap can be implemented by simply discarding displacements that do not obey this condition - this corresponds to the standard Monte Carlo (MC) procedure used in Levis and Berthier (2014) - or by performing elastic collisions when two particles are at contact, as we do here. The MC dynamics with persistence used in Levis and Berthier (2014) gives rise to clusters with self-limiting size, while using ED-BD we find clusters that coarsen by MIPS, as generically observed in previous BD simulations of ABP. .

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 $\phi_{\text{high}}=0.98\,\phi_{\text{cp}}$, in coexistence with a low density phase at $\phi_{\text{low}}=0.089$).

We identify the pressure, $P$, by an extension of the virial theorem that accounts for active forces Winkler et al. (2015); Falasco et al. (2016). After projecting eq. (1) on $\bm{r}_{i}$ and averaging over the noise, we get $P=\rho k_{B}T+P_{int}+P_{a}$, where $\rho k_{B}T$ is the thermal ideal gas pressure,

[TABLE]

is the standard virial expression of the collisional pressure, and

[TABLE]

is the active contribution to the pressure or ’swim’ pressure, which can be mechanically interpreted as a body force Takatori et al. (2014); Takatori and Brady (2015); Yang et al. (2014). This yields an ideal gas law $P_{0}=\rho k_{B}T[1+{v_{0}^{2}}/{(2D_{0}D_{\theta}})]$. The total pressure $P$ defined this way is a state function for spherical ABP Solon et al. (2015) and its definition can be extended to unconfined systems with PBC where $P$ is equal to the internal bulk pressure Winkler et al. (2015).

Following usual practice, we study the phase behaviour of ABP by letting the system evolve from a homogenous initial state at some density $\phi$ towards its stationary state corresponding to a given value of Pe. This procedure corresponds to a quench from $\text{Pe}=0$, below MIPS, to $\text{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 222We use, for convenience, an external confining potential $w(x)=\left[\left(\frac{1}{x}\right)^{12}+\left(\frac{1}{L_{x}-x}\right)^{12}\right]$. Then the wall-particle force acts along the $x$-axis: $\bm{F}^{ext}(\bm{r}_{i})=\bm{F}^{ext}(x_{i})=-\partial_{x_{i}}w(x_{i})\bm{u}_{x}$, where $\bm{u}_{x}=(1,0)$ is the base vector..

The $P(\phi;\,$Pe$)$ equations of state computed numerically using eq. (2-3) in confined and PBC geometries, following high-Pe quenches at fixed $\phi$, are shown in Fig. 2 (a-c). If the system does not phase separate, for Pe$<$Pec, 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 Winkler et al. (2015); Patch et al. (2017), 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 $\text{Pe}=30$ and AB-WCAD with $\epsilon=250$ and $\text{Pe}=100$). The formation of a macroscopic drop, due to MIPS, coincides with the location of the pressure drop, $\phi_{\text{n}}$ (see SM ). Particles are dramatically slowed down as they aggregate, making $v_{\phi}$, and therefore the active pressure, to drop (see eq. (3)). The low-density branch of the binodal, $\phi_{\text{low}}$, is well below the density $\phi_{\text{n}}$ at which the dense phase appears.

Above $\text{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}\phi)$ Winkler et al. (2015). For a system with confining walls, $B_{1}$ decreases with Pe (see SM ) in consistence with the viewpoint that self-propulsion induces an effective particle attraction Tailleur and Cates (2008); Ginot et al. (2015); Marconi and Maggi (2015). Surprisingly, for PBC, $B_{1}$ is independent of Pe, indicating that the activity-induced attraction saturates above Pec. 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 $\phi_{\text{low}}$, until it reaches the MIPS threshold, $\phi_{\text{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 Mayer and Wood (1965); Schrader et al. (2009). Accordingly, $P(\phi)$ 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 333For a first-order phase transition in equilibrium, the interface free energy density vanishes as $\Delta f\sim 1/\sqrt{N}$ and the nucleation point is shifted by $\phi_{\text{n}}\sim N^{-1/3}$ Mayer and Wood (1965).. These two effects are iconsistent with the data shown in Fig. 2 (c). The active pressure in systems of size $N\gtrsim 4000$ converges above $\phi_{\text{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$ Patch et al. (2017). 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 444Small systems should rather be thought in terms of clusters as discussed in McGinty (1973); Lee et al. (1973); Reguera et al. (2003) for equilibrium systems. The convergence of $\phi_{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\gg m_{c}$ is not guaranteed, then we cannot expect to observe how $P(\phi)$ 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-$\phi$ quench at fixed Pe by letting AB-HDs evolve from a steady-state at $\phi=0.50$, $\text{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 $\phi\approx\phi_{\text{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}\propto\gamma/\Delta G_{\text{homo}}$ (where $\gamma$ is the tension interface and $\Delta G_{\text{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 Redner et al. (2016) which provides a good theoretical description of several simulation results Richard et al. (2016). Within this framework, the critical nucleus $m_{c}\propto\phi_{\text{cp}}/\ln^{2}(\text{Pe}\,\phi_{\text{low}})$ where $\phi_{\text{low}}$ is predicted by the theory. To make a comparison with simulations, an expression $\phi(A,\phi_{\text{high}},\phi_{\text{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 $\mathcal{P}_{m}$ (see eq. (11) in Redner et al. (2016)). 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 $\mathcal{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, $\phi_{\text{low}}\to 0$ and $\phi_{\text{high}}\to\phi_{\text{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, $\phi_{\text{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}\approx 5000$ at Pe$=30$ and $\phi=0.30$ for AB-WCAD with $\epsilon=1$ in our units (see Fig. 2 in Redner et al. (2016)). This value lies quite close to the binodal in their case while for AB-HD it falls deep into the coexistence region since the low-density 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 Fily et al. (2014); Marconi and Maggi (2015); Speck and Jack (2016). The free energy associated with (heterogeneous) nucleation in a confined system is $\Delta G_{\text{het}}=F(\alpha)\,\Delta G_{\text{homo}}$, where $\ F(\alpha)=(2+\cos\alpha)(1-\cos\alpha)^{2}/4$ is a geometric function of the contact angle $\alpha$ between the wall and the dense phase. The adsorption of ABPs into layers give $\alpha=0$ (pure wetting) and therefore, by extension of CNT, a vanishing nucleation barrier $\Delta G_{\text{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. Overall, 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.

Acknowledgments

DL warmly thanks Ran Ni, Gabriel Redner, Alexandre Solon, Thomas Speck, Caleb Wagner and Thomas Voigtmann for useful exchanges, and in particular Julien Tailleur for his comments and suggestions. DL acknowledges funding from a Marie Curie IE-Fellowship (G.A. no 657517). I.P. acknowledges MINECO and DURSI for financial support under projects FIS2015-67837- P and 2014SGR-922, respectively.