Phase behavior and surface tension of soft active Brownian particles

Abstract

We study quasi two-dimensional, monodisperse systems of active Brownian particles (ABPs) for a range of activities, stiffnesses, and densities. We develop a microscopic, analytical method for predicting the dense phase structure formed after motility-induced phase separation (MIPS) has occurred, including the dense cluster's area fraction, interparticle pressure, and radius. Our predictions are in good agreement with our Brownian dynamics simulations. We, then, derive a continuum model to investigate the relationship between the predicted interparticle pressure, the swim pressure, and the macroscopic pressure in the momentum equation. We find that formulating the point-wise macroscopic pressure as the interparticle pressure and modeling the particle activity through a spatially variant body force – as opposed to a volume-averaged swim pressure – results in consistent predictions of pressure in both the continuum model and the microscopic theory. This formulation of pressure also results in nearly zero surface tension for the phase separated domains, irrespective of activity, stiffness, and area fraction. Furthermore, using Brownian dynamics simulations and our continuum model, we showed that both the interface width and surface tension, are intrinsic characteristics of the system. On the other hand, if we were to exclude the body force induced by activity, we find that the resulting surface tension values are linearly dependent on the size of the simulation, in contrast to the statistical mechanical definition of surface tension.

---

1 Introduction

Active-matter systems consist of “active” components (e.g. self-propelled nanorods, molecular motors) that locally consume energy to move, exert forces or perform chemical reactions, thus being inherently out of equilibrium. Properties of active matter such as adaptation, responsiveness, and self-healing may enable the development of novel materials and technologies.^1–12 However, to develop these next-generation technologies, a deeper theoretical understanding and description of active systems is needed. In the past two decades both mechanical and thermodynamic approaches, predicated on our understanding of equilibrium matter, have provided great insight towards an understanding of active systems, which are out of equilibrium and violate detailed balance.^13–18 However, in its current state, active matter has no complete theory, no “real gas model” which can predict emergent behavior based on system parameters.

Simulations are a useful platform for testing active-matter theories by allowing the calculation of properties inaccessible or difficult to obtain experimentally, as well as the investigation of a broad parameter space. Here, we focus on the active Brownian particle (ABP), a model subject to the overdamped Langevin equations of motion in which a particle propels itself at an intrinsic speed while rotating randomly in time.^19–21 One of the most surprising and interesting behaviors observed with the ABP model is motility-induced phase separation (MIPS), where the system undergoes a first-order phase transition into dense and dilute (gas-like) phases induced by the activity of the particles in absence of an attractive potential.^20,21

Though it is mainly activity and density that have been shown to induce MIPS,²² there are other parameters that we expect would influence the resulting structure after a phase transition. The degree of particle softness has been shown to influence macroscopic properties of colloidal suspensions in equilibrium systems,²³ with investigations both via theory²⁴ and experiment.^25–27 Specifically, the influence of particle softness has been shown to alter the flow of the liquid phase,²⁸ the conditions for glass formation (as well as its aging process²⁹) and requires the reconsideration of the relevant driving forces (e.g. the source of entropy) that determine phase behavior.²³ Moreover, experimentalists have demonstrated a great degree of control in synthesizing colloids of a specific softness, e.g. by functionalizing colloids with polymers of different lengths and densities.^30,31 Thus, a question arises, how does the rich behavior accessible by varying softness in Brownian colloidal systems transfer to active matter?

So far, no experimental studies have systematically investigated the effect of particle softness in active matter. Levis et al. computationally examined the effect of particle softness and obtained a phase diagram relating activity and softness for four distinct repulsive strengths.³² They found that making particles softer made the dense phase denser, and increased the threshold activity at which phase separation occurred. Additionally, various types of isotropic potentials have been examined: the Yukawa potential for soft particles³³ or different strengths of the WCA potential^17,34 as well as anisotropic interactions e.g. Janus interactions for Janus particles.³⁵

Most studies focus on different parameters that control the onset of MIPS. Here we focus, instead, on the structure of the dense phase and its interface with the gas phase after MIPS has occurred. The dense phase exhibits two spatial regions with distinct characteristics: a bulk dense phase and a dense–dilute interface.³⁶ The bulk dense phase has constant density, whereas, the dense–dilute interface exhibits a monotonically decreasing density from the dense to the dilute phase density,³⁶ resembling that of typical equilibrium liquid–gas interfaces.^37–39 The stability of the dense phase is dictated by the balance of incoming and outgoing flux of particles from the cluster's surface. Incoming particles from the dilute phase are initially oriented towards the dense phase until rotational diffusion causes the particle's body axis to no longer be perpendicular to the cluster's surface.²⁰ Alone, this would result in a rough interface lacking orientational alignment.⁴⁰ However, particles which bump into a rough interface will gradually move into convex regions of the surface, smoothing the interface and promoting local alignment.^41–43 This gives rise to a dense–dilute interface with a high degree of polarization of the body axes towards the cluster's center of mass, resulting in aligned body forces at the interface. To determine how these aligned body forces play a role in the mechanical stability of the steady state, we must first understand the momentum equation and its components.

In Brownian suspensions and molecular liquids the stress due to interparticle interactions is computed using virial formulae, which involves a volumetric integral of interparticle force moment.^44,45 This definition of stress recovers the Cauchy stress in continuum mechanics, σ, defined as a second-rank tensor that relates the traction vector, F̂ (force per unit area) on a surface with normal vector n̂ as F̂ = σ·n̂. Similarly the trace of the stress tensor, defined as pressure, computed from determining the force per unit area of the surface and evaluating the volumetric integral yields the same results.

Though the equivalence of interpreting physical processes from both a mechanical (microscopic) standpoint and a statistical mechanical (continuum-level) perspective applies for equilibrium systems, there has been ongoing debate about the appropriate microscopic formulation of stress in active suspensions, that is also consistent with a continuum definition. Brady and coworkers used the virial formulation to compute the average pressure within a domain containing ABPs and showed that the change in the direction of swimmers due to interactions with the neighboring ABPs reduces the effective diffusivity of the swimmers and, thus, reduces the entropic stress. They referred to this activity-induced modification to pressure as swim pressure, and used this quantity to predict the onset of MIPS in ABPs.⁴⁶ Consequent studies have shown that the pressure defined as the force per unit area on the boundaries of the computational domain is dependent on the detailed interactions of the particles with the boundary,^47,48 leading to an argument that pressure is not a state variable in active systems.

Surface tension, γ, similar to stress, is a surface quantity and is defined as the energy required for creating a unit area of the interface.⁴⁹ Kirkwood and Buff⁵⁰ showed that, similar to stress, the surface tension in molecular liquids can be formulated as integrals of interparticle forces over both phases and the interface. This formulation is consistent with the continuum definition for equilibrium systems.⁵¹ Other studies have found that using a pressure formulation that contains swim pressure and deploying Kirkwood and Buff formulation of surface tension results in extremely negative surface tension.^36,52–55

More recent studies^56–58 have argued that these inconsistencies can be resolved if the swim pressure is not included in the stress calculations and instead the effect of particle activity in ABPs is modeled through a body force, due to the net alignment of ABPs at the interface, in the continuum limit. Particularly, Omar et al. showed that ignoring the swim pressure term leads to negligible surface tension in the dense–dilute interface of phase separated ABPs.

In this paper, we analytically and computationally investigated the effect of softness for monodisperse active Brownian particles across a range of activities (Pe = 3v_pτ_rσ where v_p is the intrinsic particle velocity, τ_r is the rotational frequency, and σ is the particle diameter) and system area fractions, using the ABP model, see Table 1 for parameters details. We build upon the work of Levis et al.³² by deriving analytical formulae that predict the resulting steady state structure of soft ABP systems. Focusing on the dense phase after MIPS has occurred, we describe two analytical approaches, a microscopic and a continuum one, built from few assumptions (average interaction between particles, hexagonal-close packing structure). We derived analytical expressions for the lattice spacing and the area fraction of both the bulk dense and dilute phase. Then, in concert with kinetic theory, we obtained formulae for the cluster radius and the interparticle pressure. We found great agreement between analytical predictions and simulation results. To relate the microscopic pressure to the macroscopic pressure in the momentum equation, we explored the effect of particle softness, activity, and area fraction on surface tension. Consistent with the finding of ref. 57 and 58, we found that the swim pressure should not be included in the definition of the point-wise stress, and that the particle activity leads to a body force in the momentum equation near the interface. With these modifications, we derived a continuum approach for calculating the pressure arising from the aligned body forces at the interface (which approximately equals the interparticle pressure of the bulk dense phase), the interface width (which we find to be intrinsic to the system irrespective of varying particle softness, activity, or area fraction), and the surface tension. One important implication of our results is that across a range of parameters (softness, activity, system area fraction or size), the surface tension was found to be nearly zero and, therefore, play a negligible role in mechanically sustaining the steady state.

Table 1 Definitions of important parameters in our analytical derivations in Section S2 (ESI) and the values they take within our simulations in Section S3 (ESI)

Simulation parameters	Definition	Value
Particle diameter	σ	1.0
System size	N	10^5
System area fraction		[0.45,0.65]
Interparticle interaction strength	ε	[10^−4,10^0]
Rotational frequency	τ r = Dr^−1
Ratio of active to thermal forces (Péclet number)		[0,500]
Ratio of active to pair potential forces		[10^0, 10^6]

---

The structure of the paper is as follows. In Section S2 (ESI†) we present our microscopic theoretical framework. In Section S3 (ESI†), we outline the simulation model and details for the systems studied here. We describe our results in Section S4 (ESI†) showing comparisons between our analytical predictions and simulation results. In addition, we write down a continuum-theory approach and compare results with microscopic theory and simulations. Finally, we end with conclusions and outlook in Section S5 (ESI†).

2 Theory

Consider a colloidal particle with a stiff repulsive core that is functionalized with a weakly repulsive polymer brush. In the extremely rarified case (a dilute gas), this colloid does not interact with neighbors and has a resting diameter, σ, when only thermal forces are present. However, upon increasing the system density, the functionalized colloid deforms due to an increasing number of interparticle interactions and has an effective diameter (less than σ) defined by the distance to its nearest neighbor, ‖r_i‖. With this kind of experimental system in mind we develop our analytical model of soft ABPs.

Now, consider a system of ABPs which has undergone MIPS and is at steady-state: there is a dense phase and a dilute (gas-like) phase (Fig. 1a). In what follows, we will be focusing on the dense phase and will be calculating the lattice spacing, area fraction, cluster radius, interparticle pressure, and the surface tension based on a small number of assumptions, discussed first.

Fig. 1 (a) Schematic representation of the cluster with three color-coded regions based upon the distance from the cluster's center of mass (‖r‖): the bulk phase (green), the interface (yellow), and the dilute phase (red) with the cluster radius and interface width labeled (r_c and h respectively). The local area fraction of the bulk of the dense phase (for 0 ≤ ‖r‖ ≤ r_c − h) is constant (ϕ_d). In the interface (r_c − h ≤ ‖r‖ ≤ r_c), the local area fraction (ϕ_i) decreases from the bulk phase area fraction, ϕ_d, until it reaches the area fraction of the dilute phase outside the cluster, ϕ_g. (b) Use of the pair force in computing the total force on a reference particle as a vector sum (eqn (3)) over neighboring particles in a HCP structure, where p̂_i indicates the ith particle's orientation, θ_i represents the angle between p̂_i and the separation unit vector, r̂, and a represents the average interparticle separation distance in the bulk of the dense phase, the lattice spacing.

In the dense phase, the velocity of particles is negligible compared to the velocity of particles in the dilute phase. Thus, we assume that the particle hydrodynamic drag forces (which are proportional to the particle velocity) are also negligible. In accordance with the previous and our own findings from simulations, we also assume that the particles in the dense phase are arranged in a hexagonally close-packed (HCP) lattice²⁰ (see Fig. 1b and Fig. S1, ESI†), despite any local particle flows occurring throughout the dense phase.⁹⁰ The lattice spacing (a) between neighboring particles will be determined through a balance of the neighbors' active forces (which compress a particle) and the repulsive forces (which resist particle overlap). First, we distinguish two regions within the dense phase: the bulk, which includes the majority of particles in the dense phase, and whose constituent particles' body axes exhibit no orientational alignment on average, and the interface between the dense and dilute phases, where particles possess a high degree of orientational alignment towards the cluster's center of mass. Note that the existence of a bulk dense phase and a dense–dilute interface is supported by previous works^59–61 as well as our own simulation results presented in Section S4 (ESI†) and Fig. 5. As such, the compression of particles within the bulk results from the aligned particles at the edge of the dense phase, i.e. the interface pushing inward towards the center of the cluster. As a result, the bulk particles' effective diameter is smaller than the resting diameter. The effective diameter, which is equal to the interparticle separation of immediate neighbors, has little variability within the bulk dense phase. Therefore, we assume that each particle within the bulk dense phase has a constant distance apart from its nearest neighbors, equal to the lattice spacing, a.

Based on these assumptions, our first aim is to analytically compute the different structural and mechanical parameters of the dense phase, including its lattice spacing, cluster radius and interparticle pressure, at a variety of activities and stiffnesses (repulsive strengths). Our particles interact through the Weeks–Chandler–Andersen (WCA) potential

which provides repulsion at distances up to slightly greater than the resting particle diameter and is zero beyond that distance . Here σ defines the resting particle diameter when only thermal forces are present, ε determines the interaction strength, and r_i,j = ‖r_j − r_i‖ is the center-to-center separation between two particles i and j. The interparticle force applied by particle j on particle i is the gradient of this potential, F^WCA(r) = −∇_rU, and is given by:where r_i,j = r_j − r_i and r̂ = r_i,j/‖r_i,j‖ is the relative separation unit vector. Note that particle stiffness is modulated via the interaction strength, ε, where larger ε corresponds to stronger repulsive forces and, in turn, stiffer particles.

Let us begin by computing the active force exerted by an isolated pair within the dense phase; see Fig. 1b. The average force applied by particle 2 on particle 1 can be generally expressed as

where P(p̂₁,p̂₂) is the probability density function of observing particles 1 and 2 at orientations p̂₁ and p̂₂, respectively; see Fig. 1b. We know that in the bulk dense phase the orientations of the particles are independent of each other and are uniformly distributed: P(p̂₁,p̂₂) = P(p̂₁)P(p̂₂) = (1/2π)². The pair force is nonzero only when the relative motion of the pair causes overlap. Therefore, since the large active force (F^a) dominates over translational Brownian fluctuations, the pair force, F^pair, is only a function of activity and interparticle forces. In the absence of any orientational anisotropy, the pair force acts only along the line of centers of the particles, r̂, and is equal to the projection of the pair relative velocity in the r̂ direction:

Substituting P(p̂₁,p̂₂) = (1/2π)² and eqn (4) into eqn (3), we simplify to find the average pair force experienced by an isolated pair of particles in a dilute system

This calculation, however, ignores the effect of surrounding “bath” particles on 〈F^pair〉, which may dominate the pair interactions at large area fractions. Thus, instead of using the prefactor 4/π², we assume the general form 〈F^pair(a)〉 = βF^ar̂.

Since the dense phase particles are assumed to be static, the pair active force (〈F^pair〉) and the repulsive interparticle force (F^WCA) must be equal, giving a force balance equation that enables the determination of the lattice spacing in the dense phase (a):

where is the ratio of active to interparticle forces. We, then, proceed to use simulation results to compute β. Fig. 2 shows the simulation values of a vs. F* using eqn (6) for a wide range of Pe and ε. The data collapse into a single curve. The dashed lines shows the fit from eqn (6) for β = 2.0, which is in excellent agreement with simulation results (discussed in detail in the ESI,† see Section S1). To simply things further, we neglect the second term on the RHS of eqn (6) and compute σ/a as

Fig. 2 Lattice spacing (a) of the HCP phase for variably soft ABPs at distinct, constant, potential well depth (ε, evenly spaced on log-scale, see legend) for both simulation (points of varying shape and color, see legend) and with the best fit using β = 2.0 (discussed in detail in the ESI,† see Section S1) with eqn (6) (dashed lines) and eqn (7) (dotted lines). Increasing activity (Pe ∝ F*) or softness (decreasing ε) corresponds to a shorter lattice spacing. Varying system area fraction (different hatching, see legend) has negligible effect on the lattice spacing at constant softness and activity. Furthermore, all data collapses to a single curve in log–log scale. We find the simplified eqn (7) with β = 2.0 reliably agrees with that derived analytically in eqn (6) and measured via simulation.

We find the simplified eqn (7), plotted as a dotted line using β = 2.0 in Fig. 2, is in almost equally strong agreement with our simulation data as eqn (6). Note that β = 2 is approximately 5 times larger than the computed value, when the effect of bath particles is neglected (4/π²), indicating the dominant role of bath particles in determining the effective pair interactions; this is to be expected in this range of area fractions (discussed in detail in the ESI,† see Section S1).

Knowing the lattice spacing a enables us to determine the area fraction of both the dense and dilute phases, which allow for the calculation of three important quantities: (i) a binodal of the dense and dilute phase area fractions, (ii) the number of particles in the dense phase and (iii) the radius of the cluster at steady-state. The dense phase area fraction can be calculated from:

where is the area fraction of disks in a HCP lattice. Note that because our particles are soft they can compress so that the lattice spacing (a) is smaller than the particle diameter (Fig. 2) and the dense phase area fraction is greater than close packing, ϕ_cp. In order to obtain the area fraction of the dilute phase we follow Redner et al.²⁰ and define rates of adsorption on to (k_in) and desorption from (k_out) the dense-phase cluster:where n_g is the number density of the dilute phase, ν_p is the swim velocity of a single particle, D_r is the rotational diffusion rate, and a is the lattice spacing of the dense phase, as predicted by eqn (6). Redner et al.²⁰ used the fitting parameter κ = 4.5 explaining it stems from the observation that particle desorption occurs in avalanche-like events. At steady state, the rate of adsorption is equal to the rate of desorption, k_in = k_out, which gives

Multiplying eqn (11) by the area of a particle (A_p = πσ²), we can present this quantity in terms of typical input parameters, namely the area fraction and activity,

where is the ratio of active to thermal forces.

So far, we have calculated the area fraction of both dense (eqn (8)) and dilute (eqn (12)) phases. We can also compute the number of particles in the dense phase in terms of ϕ, ϕ_g, and ϕ_d given the system size (N), simulation box area (A), and the system area fraction (ϕ = NA_p/A), which are all known inputs for our simulations. Further simplification using eqn (8) and (12) enables us to calculate N_d (Lever rule) based upon our physical input parameters of ϕ, Pe, and ε, (discussed in detail in the ESI,† see Section S2),

The area of the dense phase cluster, A_d can be expressed as a function of the effective particle diameter (a), the number of dense-phase particles (N_d), and the packing fraction of the HCP lattice (ϕ_cp):

Next, we compute the cluster radius (r_c) as a function of activity (Pe), softness (via a), system area fraction (ϕ), and resting particle diameter (σ):

We now seek to compute the interparticle pressure within the dense phase using the virial formulation,

where r_i is the position vector between the centrally-tagged reference particle and its ith neighbor, is the number density of the dense phase, [capital Pi, Greek, circumflex]^P_d is the interparticle pressure on a single particle in the dense phase, the superscript P denotes interparticle interactions. Using eqn (14) and F^WCA(a) = 2F^a simplifies eqn (16) to

The form of pressure in eqn (17) is not immediately intuitive. The difficulty arises since we are studying the system after MIPS, where the particles are forming a crystalline phase. Our simulations and theory show that the pressure (and other variables) of the crystalline phase is independent of the initial area fraction for ϕ ≥ 0.45, even though ϕ is a determinant of the onset of MIPS. Below this area fraction, we did not observe the transition to the crystalline phase. Thus, the dimensionless pressure in our system is only a function of activity (F^a) and interparticle interactions (ε, σ). To make eqn (17) more intuitive, we now explicitly present the pressure in its dimensionless form. To do so, we substitute eqn (7) in for a in eqn (17). The form of eqn (17) is unlike the typical equations of state, where the pressure is expressed in terms of particles' area fraction (volume fraction in 3D), ϕ, and Pe.^{46,47,58,67,72,91} Note that we are interested in the structure and mechanics of the system after MIPS has occurred, where the particles in the dense phase are arranged in a hexagonal crystal configuration. As shown in Fig. 2, the hexagonal lattice spacing, and thus the pressure, of the dense phase is independent of the initial value of the area fraction for ϕ ≥ 0.45, even though the onset of MIPS depends on ϕ. Below this area fraction, we did not observe MIPS and the formation of the crystalline phase. Instead, as shown in eqn (7) the lattice spacing is entirely determined by F* and σ. Substituting eqn (7) into eqn (17) gives an expression for pressure in terms of the dimensionless quantity F* and Π₀ = F^a/σ wthout any explicit dependency on ϕ:

Recall that we are interested in the limit of Pe ≫ 1, beyond the MIPS critical point. Considering this limit in eqn (12) and (13) gives ϕ_g → 0 and N_d → N i.e. all particles will be adsorbed to the dense phase. Similarly, taking Pe ≫ 1 and writing eqn (14) in terms of r_c, we find that the radius of the dense phase scales linearly with the lattice spacing and the dimensions of the simulation box :

To sum up, we have given analytical expressions for the macroscopic mechanical variables, including interparticle pressure in the dense phase, as well as microstructural variables, including the lattice spacing, a, and the radius of the dense phase, r_c, in terms of activity Pe, softness ε, resting particle diameter σ and area fraction ϕ. Our only assumptions were that the dense phase forms an HCP lattice, that the activity is high and dominates over Brownian motion, and the only interparticle forces we consider are from immediate neighbors. To test the validity of our analytical calculation, next we compare our predictions against the results from Brownian Dynamics simulations.

3 Simulation methods

We simulate N = 10⁵ spherical active Brownian particles (ABPs), of diameter σ, confined to a two-dimensional simulation box and subject to periodic boundary conditions. Each particle's activity is modulated by varying the active force (F^a = ξv_pp̂ with drag ξ from an implicit solvent§), which is applied along a body axis, or swim direction, defined through unit vector p̂ = (cos φ,sin φ), where φ is the angle between the body axis and the positive x-axis. The active force is varied via the intrinsic particle velocity (v_p). A particle's motility, or activity, is quantified by the (dimensionless) Péclet number (Pe = 3v_pτ_r/σ). We ensure our results are not influenced by finite-size effects from the periodic boundary conditions by testing different system sizes, see ESI,† Fig. S5 and S9.

Particles translate and rotate in accordance with Brownian dynamics,

where, r_i provides the ith particle's position in two dimensions, Λ_i, and Γ_i represent zero-mean unit variance Gaussian noise, , , where α and β denote Cartesian coordinates, and is the repulsive interparticle force from the WCA potential (eqn (1)). Drag (ξ = 3πησ) and translational/rotational noise (D_t = k_BTξ, D_r = 3D_tσ²) are set according to system temperature (T) and particle size (resting particle diameter σ) for a given fluid with dynamic viscosity (η). In addition to sudden orientation changes from collisions, the body-axis reorients randomly over time according to a characteristic timescale τ_r = D_r⁻¹ = 1/3. Note that our implementation and results apply also to systems where random translational and rotational motion do not stem from the system temperature. That is to say, the emergent phenomena we observed result from the relationship between a particle's velocity and rotational frequency (τ_r = D_r⁻¹), which is encoded in the persistence length (l_p = v_pτ_r). Thus, the rotational frequency need not be explicitly thermal in origin (and could be set to any reorientation rate that reflects a particular system).

Interaction between particles is described through a WCA potential (eqn (1)). We implement several values of the repulsive well depth (ε) to study the effect of softness on the structure and mechanics of the dense phase. Recall that at a fixed interparticle distance, a larger value of ε produces a greater repulsive force i.e. hard particles have a larger value of the repulsive well depth than soft particles (ε_hard > ε_soft). For a constant repulsive strength (ε), increasing the activity (and therefore the active force) results in greater particle overlap (Fig. 3). Despite using a constant value of softness for particles in a given simulation, there is a distribution of effective particle diameters resulting from a different degree of particle compression. Depending on the environment of any given particle, the degree of compression will vary according to its neighbors' orientations and the resulting compressive forces acting on it. An experimental analogue to this implementation of the interparticle potential could be a colloid functionalized with a polymer brush where distinct repulsive strengths can be viewed as a brush of different length and density.²³

Fig. 3 Force at varying interparticle separation distance (r_i,j for the WCA potential (eqn (1), black) for constant repulsive well-depth (ε = 0.1), showing sensitivity to activity. For a head-on collision of active particles (Pe_red < Pe_blue < Pe_green) points indicate the maximum particle deformation while lines indicate range of interparticle distance available in other types of collisions. Colored spheres demonstrate the maximum deformation at each activity and wire mesh overlay shows the extent of particle overlap (for particle diameter σ = 1.0).

We used the molecular dynamics package HOOMD-Blue^62–64 to simulate N = 10⁵ monodisperse active particles for a simulation time interval of τ = 300τ_r ensuring that steady state had been reached. We varied: system area fraction (ϕ = 0.45, 0.55, 0.65), particle activity (Pe = 50–500 in steps of 50), and the potential well depth resulting in different softness (ε = 1, 10⁻¹, 10⁻², 10⁻³, 10⁻⁴k_BT). We focus on systems that phase separate via MIPS into dense and dilute phases (see Fig. 4). To overcome kinetic limitations of cluster formation, we instantiate small circular clusters (discussed in detail in the ESI,† see Section S3). We note that the steady-state composition of a cluster is independent of its initial seeded size (see ESI,† Fig. S1 and S2). As in the analytical approach, the steady-state dense phase is comprised of a bulk domain in the interior and an interface that separates the bulk dense phase from the gas phase.

Fig. 4 (a) Corresponding local area fraction of the dense phase from simulation data is computed as a local bin and shows sound agreement with analytical approach. The dense phase becomes more densely packed (and, in turn, the gas phase becomes more dilute) via increasing particle softness or activity (at constant softness). (b) Analytically computed cluster radius from eqn (15) at system area fraction of ϕ = 0.65 with simulations at various softness (color). Strong agreement between the simulation values and the analytical results are seen for all other system area fractions tested (ϕ = 0.45 and ϕ = 0.55).

4 Results

4.1 Properties of the dense phase: gas, bulk, interface

The dense phase cluster is highly dynamic, i.e. it frequently changes size (see ESI,† Fig. S3–S6) and shape (see ESI,† Fig. S7–S10), with some parts breaking up into smaller clusters and merging back, similar to ref. 20 and 22. As the activity is increased the fluctuations of the interface are decreased leading to a more stable shape (see ESI,† Fig. S7). In the analysis that follows, we are concerned with the dynamics of the largest length-scale/wavelengths (i.e. cluster radius) and sufficiently large activities that lead to crystallization of the dense phase. Specifically, we are not concerned with interface fluctuations at shorter wavelengths (the amplitude of these fluctuations is less than 1% of the cluster radius, see ESI,† Fig. S7–S10).

How do the properties of the dense phase, such as cluster size and pressure, change when the particles become softer? Levis et al. computationally calculated the phase diagram for ABPs at various particle softnesses and found a softness-dependent binodal.³² They showed that soft particles undergo MIPS at a smaller critical cluster size than hard particles; however, due to a lower nucleation barrier, softer clusters could more easily destabilize and break apart, necessitating larger activities or area fractions for sustained phase separation.³² Here, we explore a broader range of repulsive strengths (ε), activities (Pe) and system area fractions (ϕ) and compare both with our analytical predictions and the observations of Levis et al.³²

Our simulation results, in agreement with Levis et al.,³² show that softer particles pack more densely and therefore shift the dense phase of the binodal to higher densities at constant activity, (Fig. 2). At fixed softness, increasing particle activity also makes the dense phase denser and reduces the lattice spacing (Fig. 2). Increasing the system area fraction (ϕ) at values greater than the critical area fraction has negligible influence on the area fraction of the dense phase, as we see nearly perfect overlap of points with constant activity (Pe) and softness (ε) (Fig. 4a). Note that the theoretical predictions of ϕ_d from eqn (8) are in good agreement with the simulation results, as we would expect given that the lattice spacing, a, was computed accurately in our theoretical model (see Fig. 2).

The cluster radius at each softness changes little with activity, see Fig. 4b, and remains roughly unchanged with system concentration at high activities (Pe > 150), see ESI,† Fig. S11 for r_cvs. Pe of ϕ = 0.45 and 0.55. The theoretical predictions of cluster radius given by eqn (15) are in good agreement with simulation results (Fig. 4b).

Our analysis of simulations has thus far treated the cluster as a single entity. But as mentioned earlier, it is useful to distinguish two regions within the dense phase: a bulk and an interface. We define the interface as the region, where particles are orientationally aligned (pointing towards the center of mass of the cluster), and the bulk as everywhere else in the dense phase (where there is no orientational alignment), see Fig. 5(a)–(c). The interface width shows a weak dependence on softness (Fig. 5(e)) but is mostly found to be constant over all activities (see ESI,† Fig. S12), area fractions (see ESI,† Fig. S12), and simulation box sizes (see ESI,† Fig. S13) signifying that this measured interface width is an intrinsic quantity of the system. We will approximate the regions as a function of the distance from the center of mass (‖r‖): bulk ‖r‖/r_c ≈ [0,0.8] and interface ‖r‖/r_c ≈ [0.8,1.0].

Fig. 5 Orientational alignment towards the cluster's center of mass (a–c), local area fraction (d–f), and interparticle pressure (g–i) calculated using the virial formulation of pressure (eqn (16)). Data is radially binned and measured over twenty 18° conical surfaces per time step. As evident in the x-axis, the radial location (‖r‖) is normalized by the measured cluster radius (r_c) of each conical surface, and all data is averaged over time at steady-state (for at least 50τ_r) and all conical surfaces (twenty conical surfaces per time step). Similarly, the measured local area fraction and pressure are normalized by dividing through by their analytically predicted values (eqn (8) and (17)), respectively. Data shows the effect of both system area fraction (a, d and g, see legend), softness (b, e and h, see legend), and activity (cf. i, see legend). For each column, the parameter being varied is in the above legend while the other two system parameters are held constant (Pe = 350, ε = 1.0, or ϕ = 0.65).

The bulk maintains a constant average local area fraction, [small phi, Greek, macron]_local, approximately equal to that predicted by our theory, ϕ_theory in eqn (8), see Fig. 5d–f. Therefore, the area fraction is nearly constant for the majority of the dense phase, supporting the assumption of a constant lattice spacing for our analytical approach. In addition, the bulk phase exhibits no orientational alignment, , of the body forces, p̂, towards the cluster's center of mass, r, where r = ‖r‖, see Fig. 5a–c. Utilizing the virial formulation of pressure, we can calculate the interparticle pressure (eqn (16)) where we see a non-spatially varying pressure experienced throughout the bulk (Fig. 5g–i), signifying an equal degree of compression among bulk particles. Similarly, as softness (Fig. 5) or activity (Fig. 5) increases, so too does the interparticle pressure within the bulk, enabling a greater degree of particle compression (a trend that is captured, through a, in our analytical formulation of pressure as well, eqn (17)).

However, still within the dense phase cluster, we find a thin surface layer where the local area fraction begins to decrease from the bulk phase area fraction, ϕ_theory, until reaching that of the dilute phase, ϕ_g, see Fig. 5d–f. The decreasing area fraction at the interface results in a drop in the interparticle pressure (as particles are now at distance greater than a from one another, Fig. 5g–i). This monotonically decreasing area fraction between two phases resembles the density profiles of typical equilibrium liquid–gas interfaces,^37–39 thus we will henceforth label this surface layer as the dense–dilute interface. The body axes of particles becomes aligned within the interface (pointing towards the interior of the dense phase, Fig. 5a–c). As the body-axis simply dictates the direction of a particle's active force, we find that the interface exhibits an inwardly aligned active force density (Fig. 5a–c), compressing the bulk particles and giving rise to the opposing, outward interparticle pressure from the bulk of the dense phase. We claim this region of both sharply decreasing area fraction and large inwards orientational alignment is a thin dense–dilute interface layer of width h (see ESI,† Fig. S12), motivating us to create a mathematical definition to accurately identify the start (r_c − h) and end (r_c) of the interface (discussed in detail in the ESI,† see Section S4).

4.2 Surface tension and momentum transport within dense–dilute interface

We showed in Section S2 (ESI†) that the pressure in the dense phase is . Also, we know that the pressure in the dilute gas phase is negligible compared to the pressure in the dense phase. The transition from dense to dilute phase properties – including pressure, surface tension, and particle alignment – occur through a thin, dense–dilute interface. Force balance dictates that the jump in force per unit area (traction) across the interface must balance against the force induced by the interface itself. Assuming that the interfacial forces are entirely due to surface tension, the equation describing this force balance reduces to

ΔF̂_I = 2γκ_mn̂ + (I − n̂n̂)·∇γ, (22)

where ΔF̂_I = F̂_d − F̂_g is the jump in force per unit area across the interface, n̂ is the normal unit vector of the surface pointing outwards that is separating the dense and dilute phases, γ is the surface tension and κ_m = 1/r_c is the mean curvature of the interface. The first and second terms on the right hand side represents the force jumps along the normal and tangential directions of the surface, respectively. Note that the tangential component becomes negligible, compared to the normal direction, in our system. The ABP model predicts that phase separated domains coarsen with time, ultimately leading to a single cluster that scales with the dimension of the simulation box, , as derived in eqn (19). On first examination, one may think of using Young–Laplace equation, which is the form eqn (22) takes in stationary drops, to determine the surface tension of the active drop as γ_active = (Π_d − Π_g)r_c/2, (discussed in detail in the ESI,† see Section S5). Here, we assume that the gas pressure is negligible and so γ_active ≈ Π_dr_c/2. Substituting eqn (17)–(19) to compute a positive surface tension at high activities, Pe ≫ 1:
The computed surface tension through this formulation is independent of, ε, ϕ and Pe. Within this formulation, the predicted active surface tension (γ_active) is linearly increasing with the dimensions of the simulation box () without limit. However, the surface tension should be an intrinsic property of the system and, thus, must be independent of system size.

A closer examination of this formulation reveals why it cannot be used to measure surface tension. First, note that in this treatment the surface forces arise from the the net inward orientation of active particles normal to the cluster boundary (see Fig. 6), leading to a pressure jump across the surface. This is true whether the interface is flat or curved. In contrast, in mechanical and thermodynamic formulations of the surface tension for gas–liquid interfaces of passive systems,⁴⁹ we have the pressure coexistance condition i.e. the pressures in the gas and liquid phase are equal. The surface tension arises due to tangential interfacial forces along the boundary that resist the increase in the surface area.⁵⁰ In other words, the net alignment of particles at the interface of ABP clusters are not acting to minimize surface area; instead they arise in response to normal stress gradients between two phases, including pressure gradients.

Fig. 6 Steady-state ABP system with Pe = 500, ϕ = 0.55, and ε = 1.0 corresponding to a simulation frame at τ = 186τ_r. The bulk (green), interface (yellow), and dilute (red) phases are labeled according to the average interface width for the system, h ≈ 25. Particles are binned and the average orientation per bin is plotted as the arrows. It is evident that there is essentially zero alignment in the bulk while the interface is highly aligned towards the interior of the cluster. The orientation of the gas is highly random as particles are freely moving with minimal interactions. Our observations for the general alignment trends are in agreement with Fig. 5a–c.

To resolve these inconsistencies, we separate this alignment term from the calculations of surface tension and explicitly include it in the momentum equation (eqn (22)) as a body force, n(r)F^ap̂(r), while still maintaining the term modeling force jump due to surface tension, ΔF̂_I. Besides aligning at the gas-dense interface, particles can become locally aligned within the bulk dense phase, forming particle flows, vortices, and oriented grain-like domains.^88–90 However, in our systems, local alignment is short-lived and short-ranged. Given that the particles within the dense phase move as a rigid body with no relative motion with respect to the fluid, we can neglect the hydrodynamic forces and stresses. In this limit the momentum equation in a liquid–gas interface reduces to

n(r)F^ap̂(r) + f_I − ∇Π = 0, (23)

where Π is the macroscopically consistent (true) pressure of the system. Note that we have rewritten the interfacial force (ΔF̂_I) as a body force in the momentum equation: f_I = ΔF̂_Iδ(x − Γ(r)), where δ(x − Γ) Γ(r) is the Dirac delta function ensuring the surface tension term is localized to the interface region, defined by Γ(r).¶

Assuming that the thickness of the interface is much smaller than the radius of the dense phase, h/r_c ≪ 1, and that surface tension is spatially constant, the momentum equation in the radial direction across the interface simplifies to

where x = 0 and x = h specify the boundaries of the interface residing at the end of the bulk phase and at the cluster radius, respectively, 0 < x_I < h is the approximate position of the interface; and is the projection of the active force in the radial direction (see Fig. 5a–c), signifying the net orientational alignment towards the cluster's center of mass.

Multiplying both sides of eqn (24) by dx and integrating across the interface gives an expression for computing surface tension:

Fig. 7 shows the computed value of surface tension from eqn (25), utilizing simulation data for α(x) (Fig. 5a–c), n(x) (Fig. 5d–f), Π_d (total interparticle pressure calculated by eqn (17) in Fig. 8), and h (see ESI,† Fig. S12). Note that the surface tension is made dimensionless through dividing by γ_active = r_cΠ_d/2. As it can be seen, the computed surface tension fluctuates around zero without any apparent dependency on softness, activity, and area fraction (discussed in detail in the ESI,† see Section S4). This finding is in line with those of Omar et al. who also found the surface tension to be nearly zero.⁵⁸ In addition, we find the surface tension is approximately independent of simulation box size (see ESI,† Fig. S13).

Fig. 7 The non-dimensional surface tension, (2γ_true)/(Π_dr_c), calculated viaeqn (25) using values for α(x) (Fig. 5a–c), n(x) (Fig. 5d–f), r_c (Fig. 4b), and Π_d (hollow circles from Fig. 8) measured from simulation. At all activities, γ_true remains approximately constant near zero with a slight bias in the positive direction (discussed in detail in the ESI,† Section S6). The inset shows the normalized surface tensions averaged over softness (ε) and area fraction (ϕ) at each activity with error bars corresponding to a single standard deviation. In the inset, all surface tension measurements (colored) are fitted (dashed line) such that we do not bias low activity where fewer systems undergo MIPS. The line of best fit is found to be approximately constant near zero while being encompassed in the standard deviation at most activities.

Fig. 8 Interparticle pressure computed from the analytical pair-force approach (dashed lines, eqn (17)) at distinct particle softness (color) and averaged over the steady state lasting for τ ≥ 50τ_r. Simulation data calculated via the microscopic approach (eqn (16), hollow circles) at ϕ = 0.65 demonstrates good fit for stiff interparticle potential. The quality of fit decays with decreasing stiffness. Simulation data calculated via the continuum approach (eqn (26), plus markers) at ϕ = 0.65 show good agreement with both eqn (17) and (16), demonstrating the possibility to accurately calculate pressure using either a microscopic or a continuum-based approach. In addition, increasing particle activity and softness correspond to smaller clusters (see Fig. 4b) with a higher interparticle pressure.

Now that we have established that interfacial forces are negligible compared to gradients of stress normal to the boundary, we can substitute γ_true ≈ 0 into eqn (25) and compute the pressure in the dense phase by evaluating the following integral:

where the values of α(x) (Fig. 5a–c) and n(x) (Fig. 5d–f) are provided by simulations.

Having discussed both the virial formulation for calculating the interparticle pressure within the bulk dense phase and the continuum formulation for calculating the pressure arising from the aligned body forces at the interface, we proceed by calculating the total pressure experienced by each region in addition to the resulting pressure equivalence. We start by measuring the total interparticle pressure experienced by each particle in the bulk dense phase (Hollow circles in Fig. 8) from its nearest neighbors using the virial formulation of pressure (eqn (16)). The total interparticle within the bulk dense phase agrees excellently with our analytical predictions (dashed line in Fig. 8), which are linearly increasing with activity and have a slope that increases with softer particles, signifying the greater degree of compression for softer particles in the bulk dense phase. Secondly, using data from Fig. 5a–f and eqn (26), we obtain the total pressure from the aligned body-forces in the interface (Plus markers in Fig. 8), which we found to be approximately equal to the interparticle pressure of the bulk dense phase (Hollow circles in Fig. 8) at every activity and softness. As the pressure of the gas phase is negligible, this finding satisfies a steady-state force balance: the aligned active body-forces at the cluster interface are offset by compressing the particles in the cluster interior to an equilibrium separation, providing an outward interparticle pressure which balances this directed active body-force.

The agreement of interparticle pressure with the true pressure of the system in the macroscopic scale strongly supports the argument given by Omar et al. that the swim pressure, introduced in earlier studies by Takatori and coauthors,⁶⁷ should not be included in point-wise definition of the true stress in the continuum scale and the true stress can be computed using the same processes as in passive systems. Of course, particle activity does change the stress indirectly through generating a body force due to the net alignment of particles and density gradients across the interface.

Next, we ask what determines the thickness of the dense–dilute interface. Our simulations at the same ϕ, ε and Pe at different simulation box sizes show that, unlike the cluster radius, the interface thickness, h, remains unchanged. Previous studies show that the thickness of the boundary layer that forms by accumulation of ABPs near the walls scales inversely with Péclet number.⁵⁶ In contrast, the interface thickness in our case is independent of Pe, as shown in Fig. 9.

Fig. 9 The ratio (h/h_cont) of the interface width (h, see ESI,† Fig. S12) calculated via the method described in Section S4.1 (ESI†) and the interface width (h_cont, see ESI,† Fig. S16) calculated through the continuum method (eqn (27)). At all activities, h is less than h_cont by at most ≈15%. The inset shows the width of the interface measured via simulation (h). When considering the dimensionless interface thickness (h/a. See ESI,† Fig. S12), where a decreases with both activity and softness (see Fig. 2), the interface width consists of more particles for both more active systems at constant particle stiffness (ε) and softer particle systems at constant activity (Pe). The system area fraction (ϕ) has a negligible influence on the interface width, similar to its role in the surface tension.

What, then, determines the interface thickness? Moving radially outwards from the center of mass of the dense phase, the start of the dense–dilute interface is marked by a decrease in the pressure (Fig. 5g–i) and density (Fig. 5d–f), and an increase in the alignment, α (Fig. 5a–c); whereas, the end of the interface is marked by the pressure dropping to nearly zero and α undergoing a sharp decrease from its maximum to match the dilute pressure. Following these observations, we rewrite eqn (25), using the following change of variables:

[small phi, Greek, tilde] = ϕ/ϕ_d, [small alpha, Greek, tilde] = α/α_max, x̃ = x/h,

where α_max ≈ 0.45 is the maximum value of α(x) from simulations (see ESI,† Fig. S15). Applying these change of variables, dropping the surface tension contribution and integrating eqn (25) across the interface gives an expression for the interface thickness in terms of the pressure in the dense phase:where the term, h_cont (see ESI,† Fig. S16), denotes calculation of the interface thickness using the momentum equation in continuum length-scale, and
Given that the integral term I only contains scaled variables, [small phi, Greek, tilde] < 1, [small alpha, Greek, tilde] < 1, we expect it to be independent of Pe, ε and ϕ. The numerical evaluation of this integral using simulation results (see ESI,† Fig. S17) confirms this assumption where I ≈ 3.0 for all Pe, ε, and ϕ. Similar to the surface tension (γ_true), the interface width (h_cont) is found to be approximately independent of the simulation box size (see ESI,† Fig. S18).

Fig. 9 shows the ratio of the interface thickness measured from simulation, h as detailed in Section S3 (ESI†), to the calculated value of interface thickness from eqn (27), h_cont, vs. Pe for different values of ε and ϕ. As can be seen, the ratio remains close to 0.9 for all values of Pe, ϕ and ε. The close agreement between the continuum calculations and simulation results is yet another observation in agreement with formulating the pointwise true pressure in the continuum scale as the pressure that arises from interparticle forces and negligible surface forces.

5 Conclusions

A lot of studies have focused solely on the process of phase separation but not the resulting steady-state dense phase. Therefore, in this paper, we characterize and predict the properties of the dense phase itself, such as the area fraction, lattice spacing, and size. To do so, we developed a simple, microscopic analytical approach which relies on (1) the approximation of an HCP dense phase and (2) that each particle interacts with each of its neighbors with an average pair force. The microscopic, analytical approach demonstrates reasonable accuracy in reproducing the trends in simulated data for area fraction of the dense and dilute phases and size of the dense phase. These results generalize to ABP systems at any particle softness, activity, simulation box size, or area fraction. Though we utilized the WCA potential to determine interparticle interactions, we fully expect our construction to apply to other short-range repulsive potentials. An experimental validation of these results is certainly viable. We expect that similar synthetic principles to induce phoresis in hard-sphere colloids, e.g. decomposition of oxygen in a hydrogen peroxide solution at the particle surface,⁶⁸ can be extended to soft particles, such as polymer functionalized colloids. Alternatively, motility can be induced via polymer chains as is evidenced within cells⁶⁹ and which has caused the theoretical examination of phoretic polymer chains.⁷⁰

Though much progress has been made in understanding how this nonequilibrium phase separation gives rise to a dynamic steady-state, one looming question remains: how does the presence of activity influence stress/force generation in the continuum scale and how are these stresses/forces linked to the collective behavior of the system? Many studies have tried to explain this nonequilibrium phenomenon from a thermodynamic perspective via a mechanical equation of state; however, these theories give disagreeing results for important physical properties, such as surface tension. The main difference between these studies is whether activity gives rise to a stress that acts as either a spatially uniform state variable (the swim pressure^42,67,71,72) or a spatially varying body force density.^57,58 In the former group of studies the stress is defined as a volume-averaged quantity within the container such that there are no spatial gradients, and it is shown, through theory and simulation, that the activity induces an extra term referred to as swim pressure.⁶⁷ Upon utilizing the swim pressure to describe the system, one can accurately predict many emergent, macroscopic properties, such as determining the onset of MIPS^32,73 and explaining how active pressure being non-monotonic with activity and area fraction gives rise to a phase transition.^32,67,73

However, problems arise when we seek answers to localized phenomena. Omar et al. recently showed that including the swim pressure in the description of total pressure results in extremely negative values of surface tension,^36,52–55 in contrast to passive systems. Extension of the statistical mechanics derived for passive systems to its active counterparts is reliant on the system being homogeneous with no concentration or alignment gradients. However, our active systems demonstrate a monotonically decreasing area fraction at the highly aligned dense–dilute interface, which gives rise to this negative surface tension term when treating the volume-averaged swim pressure as the active analogue to the osmotic pressure. Though a volume-averaged treatment of pressure fails to explain surface tension, it does work on a macroscopic scale where the volume-averaged body forces cancel out, giving rise to the swim pressure in the momentum equation.

Therefore, accurate characterization of localized phenomena, like surface tension, requires a point-wise treatment of pressure as we cannot define a swim pressure at an interface where there is no volume-averaging involved. As a result, we no longer treat activity's contribution to pressure as a spatially independent variable but instead as a spatially varying body force. Activity produces aligned body forces acting on the boundary around the container or between phases, not a stress, which is reserved for that as defined in a passive system in order for our definition to always be correct. By always, we mean that though a volume-averaged approach (utilized commonly in traditional, equilibrium thermodynamics) works as an exception⁴⁷ for determining macroscopic properties in our ABP systems, a point-wise approach to pressure, specifically by treating activity as a body force to account for gradients in density and alignment at surfaces, satisfies a mechanical force balance on both the microscopic and macroscopic scale, enabling us to explain and predict emergent behavior. Omar et al. showed that if only the stress due to passive forces in ABP systems are considered in the definition of pressure while activity is treated separately as a spatially varying body-force density, the predicted surface tension, which relies on gradients of properties across the interface,⁴⁹ becomes negligible.

The results presented here confirm that of Omar et al., demonstrating that the point-wise mechanical effect of activity is to generate the gradients in concentration and alignment of particles, resulting in a body force (not stress), α(r)n(r)F^a, in the continuum level. Using the virial formulation of pressure, we have derived an analytical expression for the interparticle pressure in the bulk of the dense phase that agrees strongly with simulations. We also derived a second, continuum approach that utilizes the radial alignment, radial area fraction, and dense phase pressure from simulation to calculate the stress from aligned body forces at the interface, which approximately equals the interparticle pressure of the dense phase. As a result, the surface tension is approximately equal to zero for all softnesses, activities, area fractions, and simulation box size, demonstrating how the surface tension is an intrinsic property of the system. As such, we similarly confirmed that the interface width was an intrinsic quantity of the system, as similarly predicted by a local free-energy approach in equilibrium liquid–vapour interfaces,^74–77 with both the analytical interface width and that measured via simulation agreeing within 10% for all activities, softnesses, area fractions, and system box sizes.

While our results demonstrate the complex behavior that is accessible to monodisperse active mixtures of varying stiffness, the work presented here is only the first step. A number of interesting future directions are evident that will help us further understand the mechanism behind nonequilibrium steady states, namely those characterizing the dense–dilute interface. Although we found surface tension plays a negligible role in mechanically maintaining the steady state, it could play a role in other important physical properties, such as controlling fluctuations and particle flows at the dense phase surface.

If the interface is disturbed by an external force, there will be a local displacement of interfacial particles in the immediate vicinity that continues to travel tangentially across the interface while decaying in the process, like a wave.⁷⁸ These surface fluctuations give rise to long-range correlations in density across the interface, which are consistent with a description of the surface in terms of capillary waves that are thermally excited against surface tension or an external force,⁷⁹ enabling us to characterize the fluctuations similarly⁸⁰ and understand their role in stability, like cascading, avalanche events.^20,81 Surface tension could also mitigate long-term surface instabilities through surface flows in ABP systems as seen in other liquid–gas interfaces.⁸² In ABP systems, curvature-dependent surface tension drives sustained local tangential motion of particles on either side of the interface, suggesting a redirection of particles to heal local fluctuations and promote stability,⁵⁴ potentially maintaining the aligned body forces at interface that stabilizes the cluster. In addition, many other interfacial properties of ABP systems have been connected to that of equilibrium phases,^60,83–87 necessitating deeper study of the interface and surface tension's role in mechanical stability of non-equilibrium steady states.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. (NSF grant number DGE-1650116) and National Science Foundation, Grant No. DMR-1753148.

References

1. P. K. Ghosh, V. R. Misko, F. Marchesoni and F. Nori, Phys. Rev. Lett., 2013, 110, 1–5.
2. J. Palacci, C. Cottin-Bizonne, C. Ybert and L. Bocquet, Phys. Rev. Lett., 2010, 105, 1–4.
3. B. Van Der Meer, L. Filion and M. Dijkstra, Soft Matter, 2016, 12, 3406–3411.
4. A. Ghosh, W. Xu, N. Gupta and D. H. Gracias, Nano Today, 2020, 31, 1–34.
5. J. Patterson, M. M. Martino and J. A. Hubbell, Mater. Today, 2010, 13, 14–22.
6. M. Brambilla, E. Ferrante, M. Birattari and M. Dorigo, Swarm Intell., 2013, 7, 1–41.
7. W. Gao, X. Feng, A. Pei, Y. Gu, J. Li and J. Wang, Nanoscale, 2013, 5, 4696–4700.
8. J. Orozco, D. Vilela, G. Valdés-Ramírez, Y. Fedorak, A. Escarpa, R. Vazquez-Duhalt and J. Wang, Chem. – Eur. J., 2014, 20, 2866–2871.
9. J. Li, P. Angsantikul, W. Liu, B. Esteban-Fernández de Ávila, S. Thamphiwatana, M. Xu, E. Sandraz, X. Wang, J. Delezuk, W. Gao, L. Zhang and J. Wang, Angew. Chem., Int. Ed., 2017, 56, 2156–2161.
10. J. Li, V. V. Singh, S. Sattayasamitsathit, J. Orozco, K. Kaufmann, R. Dong, W. Gao, B. Jurado-Sanchez, Y. Fedorak and J. Wang, ACS Nano, 2014, 8, 11118–11125.
11. L. Restrepo-Pérez, L. Soler, C. Martínez-Cisneros, S. Sánchez and O. G. Schmidt, Lab Chip, 2014, 14, 2914–2917.
12. Q. Chen and B. Q. Ai, J. Chem. Phys., 2015, 143, 1–5.
13. M. E. Cates, Rep. Prog. Phys., 2012, 75, 42601.
14. M. C. Marchetti, Y. Fily, S. Henkes, A. Patch and D. Yllanes, Curr. Opin. Colloid Interface Sci., 2016, 21, 34–43.
15. C. Battle, C. P. Broedersz, N. Fakhri, V. F. Geyer, J. Howard, C. F. Schmidt and F. C. Mackintosh, Science, 2016, 352, 604–607.
16. É. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco and F. Van Wijland, Phys. Rev. Lett., 2016, 117, 38103.
17. J. Stenhammar, D. Marenduzzo, R. J. Allen and M. E. Cates, Soft Matter, 2014, 10, 1489–1499.
18. T. Kolb and D. Klotsa, Soft Matter, 2020, 16, 1967–1978.
19. L. Schimansky-Geier, M. Mieth, H. Rosé and H. Malchow, Phys. Lett. A, 1995, 207, 140–146.
20. G. S. Redner, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2013, 110, 1–5.
21. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244.
22. I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert and L. Bocquet, Phys. Rev. Lett., 2012, 108, 1–5.
23. D. Vlassopoulos and M. Cloitre, Soft Matter, 2012, 8, 4010–4013.
24. N. Gnan and E. Zaccarelli, Nat. Phys., 2019, 15, 683–688.
25. V. Nigro, R. Angelini, M. Bertoldo, F. Bruni, M. A. Ricci and B. Ruzicka, Soft Matter, 2017, 13, 5185–5193.
26. R. P. Seekell, P. S. Sarangapani, Z. Zhang and Y. Zhu, Soft Matter, 2015, 11, 5485–5491.
27. J. Mattsson, H. M. Wyss, A. Fernandez-Nieves, K. Miyazaki, Z. Hu, D. R. Reichman and D. A. Weitz, Nature, 2009, 462, 83–86.
28. A. Grand and G. Petekidis, Rheol. Acta, 2008, 47, 579–590.
29. C. Christopoulou, G. Petekidis, B. Erwin, M. Cloitre and D. Vlassopoulos, Philos. Trans. R. Soc., A, 2009, 367, 5051–5071.
30. N. A. Mahynski, S. K. Kumar and A. Z. Panagiotopoulos, Soft Matter, 2015, 11, 5146–5153.
31. D. Vlassopoulos and M. Cloitre, Tunable rheology of dense soft deformable colloids, 2014.
32. D. Levis, J. Codina and I. Pagonabarraga, Soft Matter, 2017, 13, 8113–8119.
33. Y. Fily and M. C. Marchetti, Phys. Rev. Lett., 2012, 108, 1–5.
34. G. S. Redner, C. G. Wagner, A. Baskaran and M. F. Hagan, Phys. Rev. Lett., 2016, 117, 148002.
35. M. Pu, H. Jiang and Z. Hou, Soft Matter, 2017, 13, 4112–4121.
36. J. Bialké, J. T. Siebert, H. Löwen and T. Speck, Phys. Rev. Lett., 2015, 115, 1–5.
37. J. Miyazaki, J. A. Barker and G. M. Pound, J. Chem. Phys., 1975, 64, 3364–3369.
38. J. D. Weeks, J. Chem. Phys., 1977, 67, 3106–3121.
39. G. Chapela, G. Saville, S. Thompson and J. Rowlinson, J. Chem. Soc., Faraday Trans. 2, 1977, 73, 1133.
40. A.-L. Barabási and H. E. Stanley, Fractal Concepts in Surface Growth, 1995.
41. A. Wysocki, R. G. Winkler and G. Gompper, EPL, 2014, 105, 1–6.
42. Y. Fily, A. Baskaran and M. F. Hagan, Soft Matter, 2014, 5609–5617.
43. N. Nikola, A. P. Solon, Y. Kafri, M. Kardar, J. Tailleur and R. Voituriez, Phys. Rev. Lett., 2016, 098001.
44. J. G. Kirkwood, J. Chem. Phys., 1946, 14, 180–201.
45. J. Irving and J. G. Kirkwood, J. Chem. Phys., 1950, 18, 817–829.
46. S. C. Takatori and J. F. Brady, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 91, 1–7.
47. A. P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri, M. Kardar and J. Tailleur, Nat. Phys., 2015, 11, 673–678.
48. T. Speck and R. L. Jack, Phys. Rev. E, 2016, 93, 1–10.
49. G. Navascues, Rep. Prog. Phys., 1979, 42, 1131–1147.
50. J. G. Kirkwood and F. P. Buff, J. Chem. Phys., 1949, 17, 338–343.
51. J. U. Brackbill, D. B. Kothe and C. Zemach, J. Comput. Phys., 1992, 100, 335–354.
52. U. Marini Bettolo Marconi and C. Maggi, Soft Matter, 2015, 11, 8768–8781.
53. S. Paliwal, V. Prymidis, L. Filion and M. Dijkstra, J. Chem. Phys., 2017, 147, 1–10.
54. A. Patch, D. M. Sussman, D. Yllanes and M. C. Marchetti, Soft Matter, 2018, 14, 7435–7445.
55. A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri and J. Tailleur, New J. Phys., 2018, 20, 75001.
56. W. Yan and J. F. Brady, Soft Matter, 2015, 11, 6235–6244.
57. J. M. Epstein, K. Klymko and K. K. Mandadapu, J. Chem. Phys., 2019, 150, 1–17.
58. A. K. Omar, Z.-G. Wang and J. F. Brady, Phys. Rev. E, 2020, 101, 12604.
59. S. Paliwal, J. Rodenburg, R. Van Roij and M. Dijkstra, New J. Phys., 2018, 20, 1–12.
60. A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri and J. Tailleur, New J. Phys., 2018, 20, 75001.
61. S. Hermann, P. Krinninger, D. De Las Heras and M. Schmidt, Phys. Rev. E, 2019, 100, 1–15.
62. J. A. Anderson, J. Glaser and S. C. Glotzer, Comput. Mater. Sci., 2020, 109363.
63. J. Glaser, T. D. Nguyen, J. A. Anderson, P. Lui, F. Spiga, J. A. Millan, D. C. Morse and S. C. Glotzer, Comput. Phys. Commun., 2015, 192, 97–107.
64. J. A. Anderson, C. D. Lorenz and A. Travesset, J. Comput. Phys., 2008, 227, 5342–5359.
65. C. S. Peskin, Acta Numer., 2002, 11, 479–517.
66. J. M. Andreas, E. A. Hauser and W. B. Tucker, J. Phys. Chem., 1938, 42, 1001–1019.
67. S. C. Takatori, W. Yan and J. F. Brady, Phys. Rev. Lett., 2014, 113, 1–5.
68. J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science, 2013, 339, 936–939.
69. C. P. Brangwynne, G. H. Koenderink, F. C. MacKintosh and D. A. Weitz, J. Cell Biol., 2008, 183, 583–587.
70. D. Sarkar, S. Thakur, Y. G. Tao and R. Kapral, Soft Matter, 2014, 10, 9577–9584.
71. S. C. Takatori, PhD thesis, California Institute of Technology, 2017.
72. S. A. Mallory, A. K. Omar and J. F. Brady, 2020, arXiv:2009.06092.
73. R. G. Winkler, A. Wysocki and G. Gompper, Soft Matter, 2015, 11, 6680–6691.
74. J. D. van der Waals and J. S. Rowlinson, J. Stat. Phys., 1979, 20, 197–200.
75. J. W. Cahn and J. E. Hilliard, J. Chem. Phys., 1958, 28, 258–267.
76. S. Fisk and B. Widom, J. Chem. Phys., 1969, 50, 3219–3227.
77. Y. Singh and F. F. Abraham, J. Chem. Phys., 1977, 67, 5960–5962.
78. M. Wertheim, J. Chem. Phys., 1976, 65, 2377.
79. R. Evans, Adv. Phys., 1979, 28, 143–200.
80. A. Wysocki, R. G. Winkler and G. Gompper, New J. Phys., 2016, 123030.
81. J. Stenhammar, R. Wittkowski, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2015, 114, 1–5.
82. J. B. Keller and M. J. Miksis, SIAM J. Appl. Math., 1983, 43, 268–277.
83. U. Marini Bettolo Marconi, C. Maggi and S. Melchionna, Soft Matter, 2016, 12, 5727–5738.
84. T. Speck, EPL, 2016, 114, 1–6.
85. V. Prymidis, S. Paliwal, M. Dijkstra and L. Filion, J. Chem. Phys., 2016, 145, 1–7.
86. C. F. Lee, Soft Matter, 2017, 13, 376–385.
87. E. Tjhung, C. Nardini and M. E. Cates, Phys. Rev. X, 2018, 8, 31080.
88. L. Caprini, U. M. B. Marconi, C. Maggi, M. Paoluzzi and A. Puglisi, Phys. Rev. Res., 2020, 2, 1–13.
89. L. Caprini and U. M. B. Marconi, Soft Matter, 2021, 17, 4109–4121.
90. L. Caprini, U. M. B. Marconi and A. Puglisi, Phys. Rev. Lett., 2020, 124, 78001.
91. A. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri, M. Cates and J. Tailleur, Phys. Rev. Lett., 2015, 114, 1–6.

---

Footnotes

† Electronic supplementary information (ESI) available: Details of simulations, analytical routines, and supporting figures (.png). See DOI: 10.1039/d1sm00350j

‡ These authors contributed equally to this work.

§ Note that the drag is dependent on the effective particle area which varies with a particle's effective diameter and, in turn, softness. Despite changing the kinetics, variations in drag (and, in turn, the rotational and translational diffusion coefficients indirectly) would negligibly influence our predictions as the structure of the dense phase is static and, hence, independent of hydrodynamic interactions.

¶ Note that while we have presented the interface in the continuum limit as a surface with no volume (line in 2D), in simulations these variations in properties occur over an interface with finite thickness. This is, in practice, similar to approximating Dirac delta function with a smooth and differentiable function with finite, yet small, spreading length, as done in numerical techniques such as immersed boundary method.⁶⁵ The alignment term n(r)F^ap̂(r) is analogous to the gravitational body force that appear in the formulation of the pendant drop experiment for determining surface tension.⁶⁶

---