Motility-induced phase separation of soft active Brownian particles

The non-equilibrium conditions arising from the loss of microscopic detailed balance have given rise to a wide program of research, particularly in active matter,¹ where particles self-propel by transforming energy into active motion through complex mechanisms. Non-equilibrium conditions trigger a richer phenomenology of a wide range of systems, ranging from slow dynamics relaxation, such as in the glass formation,^2–4 to systems with irreversible dynamics due to energy dissipation like in active matter.^1,5–9 The principal characteristic of the systems cataloged in the last case consists in the persistence of motility of their constituent particles, commonly due to a self-propulsion.^10–18 

The main feature of active motion is the tendency of the particles to maintain the direction of motion, called persistence. Paradigmatic examples of active motion are found in biological systems, at macroscopic scales (“dry” active matter) like in bird flocks^19,20 and fish schools^21,22 but also at mesoscopic scale (“wet” active matter) in bacterial colonies^23–25 or actin filaments.^26,27 Also in these mesoscopic lengths scales, one can find artificially designed colloidal systems like autophoretic colloids^28–30 and colloid rollers,^31,32 which exhibit the same common behavior, namely, pattern formation, non-equilibrium transitions, and anomalous fluctuations.^33–35 

A salient feature of the dynamics of these systems, observed in both experiments^36,37 and numerical simulations,^38,39 corresponds to the motility-induced phase separation (MIPS) that emerges in dense enough systems when the persistence length overpasses the mean separation length between particles. MIPS is observed in different systems of active matter, therefore is independent of the pattern of active motion (active Brownian motion, run-and-tumble, active Ornstein-Uhlenbeck, etc.). Recently, dynamic and MIPS features of active matter have been extensively studied through experimental,^24,40,41 theoretical,^6,42,43 and simulation^44–46 approaches. These studies have allowed to identify general characteristics in systems composed of self-propelled particles, as for instance, the cluster formation in active systems without attractive forces nor alignment.^29,47,48 In such cases, for systems at high enough densities and high self-propulsion velocities, particles tend to accumulate in the regions where the particles velocity abruptly decreases due to crowding effects. As a consequence of this, a fraction of particles are free to move in the remaining available space, thus forming a dilute gas and, therefore, giving place to a stable state with two separated phases. The transition from an homogeneously dense state to the phase-separated one corresponds to a purely motility-induced phase separation.^{39,44,49–54}

The phase diagram for active Brownian particles (ABP) is very rich, going from MIPS to melting,^52,55 and it has been recognized that MIPS in this kind of non-equilibrium system resembles the phenomenology of the gas–liquid transition in thermal equilibrium liquids. However, in this last case, and in contrast to active systems, attractive interactions between particles are indispensable.^56–58 Moreover, striking differences between the phase behavior of passive and active fluids have been reported. For instance, the spontaneous alignment of velocities,⁵⁹ the microphase separation,⁶⁰ negative interfacial tension,⁴⁸ and the difference between the temperature of coexisting phases.⁶¹ Thus, different thermodynamic approaches have been followed in order to determine the phase separation behavior: nucleation kinetics,⁶² effective mean-free fields energies,^63–66 pressure equations of state,^50,67–69 and effective representation through isotropic pair potentials.^70,71 Nevertheless, due to the intrinsic complexity of studying phase coexistence in both passive and active systems, such approaches are commonly carried out by means of computer simulations with different numerical techniques, such as the order parameter fluctuations and the crossing of cumulants^72,73 or the direct determination of coexisting densities.^48,58 Their main purpose is the determination of the critical density for the coexistence of phases.^48,74,75

Different studies in passive fluids^58,76 have shown that the details of the phase diagram are sensitive to the inter-particle interaction potential; slight differences in the potential can lead to significant changes of the critical point and of the shape of the gas–liquid coexistence line. For active matter, Debets and Janssen have analyzed the influence of particle dynamics on the glassy behavior of active Brownian particles;⁷⁷ Martin-Roca et al.³⁹ carried out a comparative study of the effects of the standard n = 6 Weeks–Chandler–Andersen (WCA) potential⁷⁸ and the pseudo hard-sphere (PHS)^79–81 potentials, on the local and global properties of the structure of the dense phase of MIPS, particularly in the phase space defined by the particle density and the Péclet number. Interestingly, it is reported that the liquid–gas boundary depends on how the parameters involved in the Péclet number are varied. The authors conclude that when the “self-propulsion force” is strong enough, the potential softness becomes relevant in the state diagram. In the case of the standard WCA potential, the overlapping of particles in the dense region induces a shift of the MIPS phase boundary to higher density. This effect is absent for the stiffer PHS potential. Therefore, the authors suggest that different repulsive potentials might lead to differences in the phase behavior. In this regard, it is worth mentioning that active Brownian particles are mainly modeled with the WCA potential. For this active fluid, not only the MIPS but also the liquid,⁸² hexatic,⁸³ and solid phases⁸⁴ have been studied. Accordingly, the MIPS, and specifically the phase and critical behavior, would be defined by the particle softness besides the density and self-propulsion velocities.

However, in this study, only one soft active fluid was studied and the phase diagram was determined in a qualitative way employing a local structural and dynamic criterion based on the hexatic order parameter and the non-Gaussian parameter, respectively. Unfortunately, the effect of the particle softness is not identified and the critical point of the MIPS is not quantified.

Motivated by the significant interest in active matter systems, this contribution focuses on the analysis of the influence inter-particle softness interactions on the MIPS and its corresponding critical point. To achieve this, we use soft active Brownian particles (ABP) in combination with computer simulations and methodologies based on the particle distributions, commonly employed in passive systems, to determine liquid phase diagrams. Our analysis takes into account the effects of different softness described by a WCA-like potentials with $n = 5 , … , 12$⁠, elucidating the role of these on the appearance of MIPS. We found that the appearance of MIPS shifts to higher values of the self-propulsion velocity as the potential becomes softer. In addition, we also elucidate the applicability of standard methods used in the study of phase behavior of passive fluids to the analysis of the intrinsically out-of-equilibrium active matter.

By delving into these aspects, we aim to provide deeper insights into the complex interactions that contribute to the emergence of distinct phases and critical points in active matter systems. Thus, in Sec. II, we detail the soft ABP model and the simulation scheme employed. In Sec. III, we discuss the physical conditions at which coexisting phases are formed. Then, density profiles are analyzed to determine the MIPS phase diagram and corresponding critical points for each soft fluid studied. Finally, in Sec. IV, conclusions of this work are provided.

The interaction potential of this fluids is illustrated in Fig. 1. We study five different soft ABP fluids defined by the exponents $n = 5 , 6 , 7 , 9$⁠, and 12. As one can see, if n increases, the potential becomes more alike to hard-sphere (HS) potential, also displayed in the same figure with a gray line. It is worth mentioning that depending on the particle softness, particles will exhibit a soft-core interaction range slightly smaller than the hard-core diameter, σ. This implies that particles are allowed to overlap to some extent. However, the repulsive interaction extends over ranges greater than σ. These interaction patterns are schematically illustrated in the inset of Fig. 1.

The PHS potential is also displayed in Fig. 1 with a magenta line. Although this potential shows slight deviations from the HS one, it can accurately reproduce thermodynamic⁷⁹ and transport properties^80,81 of HS; moreover, it allows to avoid any correction to the particle diameter and system density.

Regarding the simulation details, all cases considered consisted of a rectangular domain of area $A = L x L y$⁠, where $L x = 2 L y$⁠, in order to observe the phase separation in systems composed by N = 16 589 circular particles. The equations of motion were integrated with the standard Ermak and McCammon algorithm without hydrodynamic interactions,⁹⁰ and we impose periodic boundary conditions in each direction. Using the usual units of length, σ, and time, $[ D o / σ 2 ] − 1$ (one-third of persistence time), and considering equilibrium rotational diffusion,⁹¹ we express all results in dimensionless units.

Assuming that the phase behavior of ABP can be described with just two parameters, namely, the area fraction $φ = π 4 ρ$⁠, where $ρ = N σ 2 / A$ represents the number density, and the dimensionless self-propulsion speed $v p = v 0 σ / D$⁠, akin to a Péclet number, we direct our analysis toward the MIPS dependent on v_p and n. It is worth to stress that for this definition we use the particle diameter σ and the short-time diffusion coefficient D as the corresponding units to define dimensionless length $r / σ$ and time $t σ 2 / D$⁠.

Notably, the dimensionless persistence length $ℓ p = v p ( σ 2 / D )$ varies with v_p. To ensure a comprehensive exploration of the phase space, we conduct ten distinct simulations for each state of interest, originating from various starting points. The density profiles presented herein are averages derived from these simulations.

To observe phase separation, characterized by the coexistence of diluted and dense phases, in systems with relatively high total density (⁠ $φ ≈ 0.5$⁠), we systematically varied v_p across a range of values for each fluid with different degrees of softness (n). Under these conditions, a spontaneous aggregation of slower particles into a dense phase is observed, surrounded by a diluted region of faster particles. In Fig. 2, we present snapshots captured at the steady state of fluids with softness levels n = 5, 7, and 12, all possessing a self-propulsion velocity of v_p = 140.

Due to periodic boundary conditions, the emergence of two interfaces is observed; however, in contrast to passive fluids, the interface exhibits a largely fluctuating interface profile. This phenomenon can be elucidated by considering the net forces acting on particles at the interface are also fluctuating, combined with the stochastic nature of self-propulsion velocities.⁵⁴ The fluctuations of the interface, whose origin is necessarily found in the non-equilibrium nature of persistent motion, have become a topic of interest in itself;⁹² here, such effects are visualized as a rough interface as is shown in Fig. 2. Undoubtedly, the creation of the dense slab is a direct consequence of self-propulsion and, therefore, of v_p; additionally, we note that its width correlates with particle softness. As illustrated in Fig. 2, if the ABP interaction becomes more similar to the hard-sphere interaction, as studied here by increasing the n-exponent in the interaction potential [Eq. (4)], the dense phase covers a broader region of the simulation box. At the same time, the shape of the interface between the dense and dilute phases becomes more complex, i.e., the fluctuations of the interface around the average one become larger. This can be understood in terms of the relationship between the distinct contributing terms of the interaction potential. A larger repulsive range corresponds to a softer potential (see Fig. 1). Then, although the microscopic effect of a self-propulsion velocity is reflected in a closer separation distances between particles, the repulsive contributions in softer fluids prevent the formation of large dense regions, smoothing the interface.

In addition to the qualitative interface behavior presented, we quantitatively assess the density of both coexisting phases by computing the density profiles $φ ( x )$ along the x-direction. These profiles are determined relative to the system's center of mass and to mitigate interface fluctuations, we average over a minimum of ten distinct long-time simulations, as mentioned previously. In Fig. 3, we present representative $φ ( x )$ results for a fluid with softness n = 7. The illustrated density profiles correspond to systems with self-propulsion velocities $v p = 80 , 100 ,$ and 140. Additionally, the same figure depicts the particles distribution within one-half of the simulation box, offering a snapshot of the mentioned fluid with v_p = 80. For each fluid under study, the qualitative behavior of $φ ( x )$ closely resembles the behavior and features shown in Fig. 3.

Observing the density profiles one can easily recognize the regions with high and low densities, we denote the area fraction of those regions as $φ +$ and $φ −$⁠, respectively. Then, starting from high v_p values, where the dense phase is wider and the correspondingly surrounding phase is diluted, it can be noticed that, as v_p decreases, density in dense phase decreases while it increases in the surrounding phase. This behavior is similar to the one observed in passive fluids when the temperature approaches the critical value. In such a case, if coexistence densities are analyzed in terms of temperature, one can establish the phase behavior of gas, liquid, and their coexistence. Thus, by establishing an analogy between an effective temperature and the self-propulsion velocity, we can perform an equivalent description for soft ABP fluids by determining the densities of coexistence phases but as a function of the self-propulsion velocity.

Once $φ ±$ has been determined by the fit, the phase diagram is built by expressing them as a function of the self-propulsion speed v_p. Results for each soft fluid (⁠ $n = 5 , 6 , 7 , 9 , 12$⁠) are illustrated in Fig. 4. There, taking as a reference the critical point of the hard-sphere active fluid, values of $φ −$⁠, corresponding to low area fractions, are depicted toward the left (green shaded region), while high area fractions values $φ +$ toward the right one (blue shaded region).

From Fig. 4, it can be noticed that for soft particles with $n ≤ 7$⁠, the phase coexistence is characterized by large area fractions $φ +$⁠, but these values are shifted toward smaller values as n is increased, i.e., the more stiff the particle interaction became, the lower is the fraction area of the dense phase, whereas the surrounding phase increases its density. This also induces that the MIPS region found at lower self-propulsion velocities if the inter particle interaction is stiffer.

Additionally, in Fig. 4, we also include the phase results (gray crosses) obtained by Siebert et al.⁵¹ for the WCA (n = 6) potential, using the Barker-Henderson correction,⁹⁴ which implies a redefinition of the particle diameter and in consequence of the fluid density. The agreement between such an approach and our PHS results is remarkable. Although, due to the non-equilibrium nature of the system induced by activity, it is arguable that finite-size effects would be important in determining the phase behavior. Siebert et al. explore this question determining the system size-dependence of order parameter, particularly analyzing density fluctuations and the crossing cumulants. Previously, Bialké et al., using a similar approach as us, demonstrated that there are no differences between systems sizes $4900 ≲ N ≲ 52 000$⁠. In fact, our results are in good qualitative agreement with this last approach; however, the corrected Barker-Henderson diameter σ_BH used by Siebert and Bialké is quite different between them. This highlights the need for consistency between different approaches and the reason why we opt to use the PHS potential without corrections to the particle diameter.

From Fig. 5, we observe that the critical self-propulsion velocity is more sensitive to the particle softness. In this regard, for softer particles (⁠ $n ≤ 7$⁠), higher propulsion velocities are needed in order to observe a MIPS. This behavior can be seen in the results for the critical point shown in Fig. 5 for fluids with $n = 5 , 6 , 7$⁠. It is also observed that for the same softer fluids, corresponding critical densities are almost the same. Additionally, for even softer fluids, e.g., $n ≤ 4$⁠, modeled with the generalized WCA potential, we do not observe the formation and coexistence of a dense and diluted regions. We attribute this to the repulsive part of the corresponding interaction potentials, whose range is around $20 %$ of the particle size. This repulsion prevents the formation of the dense phase. On the other hand, as the interaction between particles becomes more stiffer (⁠ $n > 7$⁠) MIPS take place at lower v_p values, for instance, for n = 12 the critical velocity value is almost the half than the one obtained for the n = 5 fluid.

The critical area fraction exhibits a clear non-monotonic behavior. This phenomenon can be attributed to both the softness of the particles and their self-propulsion velocity. Within the MIPS region and at sufficiently high v_p values, particles can approach each other closely, and their interaction is primarily governed by the soft-core part of the interaction potential. As the self-propulsion velocity, v_p, decreases, interactions transition to the hard-core repulsion range. A schematic representation of these interaction regimes can be found in the inset of Fig. 1.

Therefore, softer interactions allow for denser regions with higher area fractions to form, as evidenced by the progression of this phase branch within the blue shaded region of Fig. 4. Interestingly, the dense branch corresponding to n = 9 crosses the one corresponding to the n = 12 fluid. This could stem from a competition between the effects of self-propulsion and excluded volume interactions. Moreover, it might be the underlying cause of the observed non-monotonic behavior.

We want to comment that while the binodal-like lines and critical point of the active system under study were determined by using the methods of an equilibrium description, there is still an intense ongoing research to clarify other aspects, as the universality class of the MIPS.^51,96,97 However, our results indicate that corresponding states law is not meet for these non-equilibrium systems. Furthermore, there are some controversy about the surface tension between the dense and the dilute phases^48,74,92 since it must be positive in order to address other important phenomena such surface free energy dependence on curvature in liquid drops or bubbles,⁹⁸ adsorption at fluid interface,⁹⁹ and interface tension of curved interfaces.¹⁰⁰ For the time being, such questions are out of the scope of this work, however, as previously discussed, the parameters of the interaction potential exert significant influence over the phase behavior. Consequently, it is anticipated that these parameters will similarly exert a substantial impact on surface tension phenomenon.

The hard-sphere model is usually employed as a first approach to describe a wide variety of systems ranging from simple liquids to active matter. This interaction potential can be used to characterize active colloids whose propulsion mechanism can be due to different phoretic phenomena. Living active matter, on the other hand, is far to behave like hard-spheres; therefore, one can expect a phase behavior strongly dependent on softness of the particles. Here, it is demonstrated that MIPS on soft ABP fluids is observed at higher self-propulsion velocities as particles becomes softer. This also allows us to observe dense and dilute coexistence phases with higher area fractions as shown by the corresponding binodal-like lines. The critical self-propulsion velocity turns to be very sensitive to particle softness. Moreover, critical density exhibits a clear non-monotonic behavior resulting from combined effects of self-propulsion velocity and particle softness.

We verify that this non-equilibrium system seems to not follow the corresponding states or extended corresponding states law. This does not allow to classify ABP fluids in a sort of universality class, at least for the specific motility pattern studied.

The high fraction areas observed for the dense region in soft potentials indicates a more complex distribution of particles in that region. This observation compels for the formulation of a research program (theoretical and experimental) to study such dense region which might differ from a liquid phase. Thus, the softness of the interparticle interaction would induce more interestingly separated phases.

Finally, with the confidence of using equilibrium approaches to describe, in some extent, MIPS of active Brownian particles, it would be of interest to systematically determine and describe the phase behavior of active matter systems with other motility patterns. Additionally, the clear identification of the MIPS boundary, through its critical point, will allow to perform systematic studies regarding to emerging effects near to such transition, like the formation of spatial patterns, freezing transition, and glassiness in active fluids. Work in this direction is in progress.

A.T.-C. acknowledged the postdoctoral-scholarship granted by CTIC-DGAPA-UNAM (2022). This work was supported by UNAM-PAPIIT IN112623. The authors thankfully acknowledge the computer resources, technical expertise and support provided by the Laboratorio Nacional de Supercómputo del Sureste de México, CONACYT member of the network of national laboratories. Authors also thank to C.E. López-Natarén for the valuable technical support on the use of computer infrastructure of the Instituto de Física UNAM.

The authors have no conflicts to disclose.

Alexis Torres-Carbajal: Conceptualization (lead); Methodology (lead); Software (lead); Writing – original draft (lead); Writing – review & editing (equal). Francisco J. Sevilla: Conceptualization (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article.