Interface stability, interface fluctuations, and the Gibbs–Thomson relationship in motility-induced phase separations

Abstract

Minimal models of self-propelled particles with short-range volume exclusion interactions have been shown to exhibit the signatures of phase separation. Here I show that the observed interfacial stability and fluctuations in motility-induced phase separations (MIPS) can be explained by modeling the microscopic dynamics of the active particles in the interfacial region. In addition, I demonstrate the validity of the Gibbs–Thomson relationship in MIPS, which provides a functional relationship between the size of a condensed drop and its surrounding vapor concentration. As a result, the late-stage coarsening dynamics of MIPS at vanishing supersaturation follows the classic Lifshitz–Slyozov scaling law.

---

1 Introduction

Phase separation is a ubiquitous phenomenon in nature and is manifested by the partitioning of a system into compartments with distinct properties, such as the different particle densities in the two co-existing phases in the case of liquid–vapour phase separation. Phase separation under equilibrium dynamics is a well-investigated physical phenomenon.^1,2 Recently, signatures of phase separation have been reported in non-equilibrium systems consisting of active particles.^3–19 It is therefore a natural question to ask to what degree we can extend our knowledge of equilibrium phase separation to the phase separation phenomenon observed in active systems. In the case of minimal models of active particles with simple volume exclusion interactions, the phenomenon of motility-induced phase separations (MIPS) has received considerable interest.^3,5,8–18 In particular, the idea of an effective surface tension in motility-induced phase separations (MIPS) has been advocated.^15,18 At the gas–liquid interface in thermal equilibrium, surface tension results from the pulling of molecules at the interface due to their attractive interactions.²⁰ In a system of active particles with purely repulsive interactions, it is unclear how such “pulling” can occur as the particles can only push. To probe what happens at the interface, I have studied the microscopic dynamics of the active particles in the interfacial region using a combination of simulation and analytical methods. Specifically, using mean-field type arguments, I will demonstrate how pressure balance is achieved between the condensed phase and the dilute (vapor) phase, and how the Gibbs–Thomson relationship arises in a system where a circular condensed drop co-exists with the vapour phase. Furthermore, by incorporating the stochastic nature of particle dynamics, I will explain the scaling between the interfacial width and the system size recently observed in MIPS.¹⁸

1.1 Motility-induced phase separation

I will first focus on a minimal model system that exhibits MIPS in two dimensions (2D)—a collection of self-propelled particles with excluded area interactions that undergo rotational fluctuations. Specifically, the dynamical equations arewhere i is an integral index enumerating the particles in the system, v_i ≡ cos θ_ix̂ + sin θ_iŷ with the angle θ_i (with respect to the x-axis) being the orientation of the i-particle, g_i(t) is a noise term with Gaussian probability distribution with zero mean and unit variance, D_r sets the magnitude of the rotational fluctuations, U(·) corresponds to the potential function for short-ranged area exclusion interactions, η is the drag coefficient and f_a is the constant active force that drives the particles in the system. In particular, u ≡ f_a/η is the constant speed of a particle when it is not within the area exclusion zone of another particle. Previous numerical work has indicated that phase separation in this minimal system occurs as u increases, but the actual form of U is unimportant.^8,11,14 For instance, U could be in the form of a Weeks–Chandler–Andersen potential:²¹This will be the particular form of the potential function utilized in this work. Also, the time and length units will be set as having a = 1 and D_r = 3. Note that I will focus exclusively on the non-equilibrium dynamics of the system and therefore translational Brownian motion is ignored.

2 Flat interface

2.1 Point particles

To understand the microscopic dynamics at the interface, it is instructive to first look at a system of point particles, i.e., the interaction potential U is zero. Even this simple system distinguishes itself from the equilibrium system in that aggregation will spontaneously occur in the proximity of a force-absorbing but frictionless wall (left column of Fig. 1). In other words, the particles are free to slide and rotate at the wall, but they cannot penetrate the wall. Further complex patterns are revealed when one looks at the particles' orientation distribution as well as the position distribution (Fig. 1(e)). At the wall, most of the particles are left-going, as indicated by the high concentration of orientation at around θ ≃ π. This results from the fact that only left-going particles remain at the wall. Just outside the wall, the distribution is highly peaked at θ just below π/2 and just above 3π/2, which reflects the particles' orientation after they move away from the wall. The orientation anisotropy decays as a particle moves away from the wall.

Fig. 1 Steady-state configurations of active particles confined by a force absorbing wall on the left: point particles (left column) and repulsive particles (right column). (a and b) A snapshot of the system at the end of the simulations with the orientations depicted by the blue arrows. The red circles depicted in (b) are of diameter a = 1. The wall at x = 0 is perfectly force absorbing (see A.1 for simulation details). (c and d) The histograms show the horizontal distributions of the particles. (e and f) The colourmaps show the deviation from the mean in the particles' orientations at different horizontal positions. The colour scale corresponds to the measure: where i is the row index and h_i(x) is the frequency. The simulation parameters are: f_a = 100, η = 1, a = 1, D_r = 3, and A = 25/6.

In this system, the pressure acting on the wall can be expressed as

where the orientation distribution function of the particles at the wall per unit length is denoted by χ_W(θ). Note that since we are dealing with a 2D system, the unit of pressure is [force]/[length].

To further analyse χ_W(θ), one can perform dimensional analysis to conclude that

where [scr F, script letter F](θ) is a function dependent only on θ, and ρ_∞ is the particle concentration far from the wall. To obtain the exact functional form of [scr F, script letter F](θ), one needs to solve a set of two coupled differential equations under mixed boundary conditions,²² whose solution consists of a series of Mathieu functions. Unfortunately, the expansion coefficients in the series are not analytically tractable and so an analytical expression is lacking. However, [scr F, script letter F](θ) can be readily estimated numerically (left column of Fig. 1), which allows one to obtain the following equation:The second expression is equivalent to the swim pressure^23–26 of a system of active particles in 2D, which I have obtained here numerically.

2.2 Repulsive particles

Remarkably, much of what we have seen in the point particle case remains true when we add mutually repulsive interactions to the particles. When a force absorbing but frictionless wall constitutes the left boundary of a semi-infinite system, phase separation occurs where the condensed phase is located close to the wall (Fig. 1(b)). Inside the condensed phase, the orientation is isotropic (Fig. 1(f)). The reason behind the isotropy is that the impeded motility of the particles make them stay for a duration much longer than the orientation decoherence time ≃1/D_r. However, note that the particles' locations are not frozen in time as shown by the black particles in Fig. 1(b), which were the first column of particles next to the wall at the beginning of the simulation. As a particle moves further to the right, it first encounters a layer of left-going particles (shown by the bright red patch centred at x ≃ 17 in Fig. 1(f)). This represents the accumulation of particles with an orientation highly centred at around θ ≃ π, analogous to the accumulation of active point particles at the wall, except that the particles are now spread over a range of x positions due to volume exclusion interactions. Further rightwards, we encounter a pattern of two escape trajectories away from the interfacial region indicated by the two red-yellow branches emerging from the red patch. As a particle moves further away to the right, the orientation becomes isotropic again. From this discussion, it is clear that in the bulk of the condensed and vapour phases, the corresponding orientation distributions are both isotropic, while in the interfacial region separating them there is a high level of orientation anisotropy.

Let us now calculate the force exerted on the wall by the active force of these particles. Since in the condensed phase, the orientation is isotropic, the pressure felt by the wall due to these active forces is

where ρ_c is the concentration in the bulk condensed phase. Due to the orientation isotropy, the expression here scales like f_a instead of f_a² in the point particle system (eqn (6)). Besides the active force contribution, the wall will also feel additional forces arising from the repulsive interactions, which, we will see, constitute an important contribution in achieving a pressure balance in the interfacial region in MIPS (Section 2.4).

2.3 Locating the interface

The concentration variation across the two phases shown in Fig. 2(a) is similar to a typical equilibrium phase separation. What distinguishes MIPS is the high orientational anisotropy between the phases (Fig. 2(b) and (c)). As for equilibrium fluids, the location of a sharp interface between the two phases can be defined somewhat arbitrarily.²⁷ In the case discussed here, since the pronounced minimum of 〈v_x〉 ≡ 〈cos θ〉 is easy to locate (indicated by the red broken line) and its location also marks the onset of the increase in Q_yy ≡ −〈cos(2θ)/2〉 (the yy component of the nematic order parameter Q),²⁸ which signifies the escape of particles from the condensed phase, it is a convenient choice for the interface location. In other words, this convention implies that right outside the interface of the condensed phase, there is a layer of particles travelling preferentially along the interface, as indicated by the peak in Q_yy (Fig. 2(c)). The active forces of these escaped particles are potentially the cause of the emergence of a negative surface tension according to its mechanical definition.¹⁸ The definition of the interface location has of course no physical significance, but this does provide a working definition useful for the sharp interface model discussed below.

Fig. 2 (a) Particle concentration as a function of x. The same figure as in Fig. 1(d). (b) The mean horizontal component of the particles' orientations 〈v_x〉 vs. x (+ symbols, left y-axis), and ρ∂_xρ vs. x (○ symbols, right axis). (c) The yy component of the nematic order parameter Q vs. x, where Q_yy ≡ −〈cos(2θ)/2〉 so that a high Q_yy signifies that the orientations of the particles are preferentially pointing up or down. Similar to equilibrium fluids, defining the location of a sharp interface is somewhat arbitrary.²⁷ The working definition proposed here is that the interface is set to be at the pronounced minimum of |〈v_x〉|.

2.4 Interface stability

2.4.1 Pressure balance at a sharp interface. In this section we will see how pressure balance can be achieved in MIPS. Note that the discussion in this section amounts to a simplified exposition of that in ref. 23 and 29. Its presentation here is for the self-containedness of this paper and will help in understanding the approximations used in the later sections.

I will start by discussing a sharp interface model. In this drastically simplified model, let us imagine that the phase separated system is partitioned by a sharp interface where in the vapour phase, the concentration is low enough for the system to behave like a system of active point particles, and in the condensed phase, the orientation distribution is isotropic. One can imagine such a system by first rotating Fig. 1(e) by 180° and then collating it to Fig. 1(f) on the left. As calculated before, the pressure exerted by the vapour phase on the sharp interface is f_a²ρ_v/(2ηD_r) (eqn (6)); this pressure is balanced by the pressure exerted by the condensed phase: aρ_cf_a/π + P_r, where the first term comes from the active force (eqn (7)) and the second denotes the pressure arising from the repulsive force due to the area exclusion interactions. In other words, pressure balance is achieved if

For the simulation parameters used in Fig. 1 (with the units set by a = 1 and D_r = 3), aρ_cf_a/π ≃ 30, P_r ≃ 230 and ρ_vf_a²/(2ηD_r) ≃ 330. We thus see that in the simulated system, the L. H. S. and the R. H. S. of eqn (8) are of the same order of magnitude, indicating that the pressure balance conditions are qualitatively satisfied. In addition, we see that much of the active force coming from the vapour phase is used to compress the condensed phase via the pressure P_r. Note that although eqn (8) provides a pressure balance condition for the system, it does not mean that any system satisfying this condition is stable as it may still be unstable against fluctuations. This is not dissimilar to equilibrium fluids where pressure balance together with chemical potential balance is needed to achieve phase stability.

The pressure balance condition in eqn (8) already allows for crude estimation of the minimal active force required for MIPS. I assume for simplicity that at the onset, (i) the condensed phase is not very compressed and so we can ignore P_r,³⁰ and (ii) the concentration ratio between the condensed phase and the dilute phase (ρ_c/ρ_v) is of order 1, which leads us to the minimal force requirement below for MIPS:

The above condition comes from the pressure balance condition at the interface alone. Interestingly, eqn (9) reproduces the same scaling as obtained by Redner et al. via a different approximation.¹¹ The result also supports the notion that the Péclet number (Pe), usually defined as Pe ∝ f_a/(aηD_r) in this context, is a key control parameter in MIPS.^11,31

From this sharp interface model, we can now see why the condensed phase with a high density of active and repulsive particles can remain stable against a backdrop of dilute concentration of active particles – the active particles in the vapour phase impact an active pressure that scales as f_a² directed towards the normal of the interface, while the countering active pressure from the condensed phase scales as f_a. The quadratic dependence on the active force is due to the fact that only particles pushing against the interface remain at the interface, while particles oriented away from the interface leave. The escaping particles thus open up space for other particles that serve to push against the interface. From this perspective, the low concentration in the vapour phase is paramount for the stability of MIPS, otherwise the particles oriented away from the interface may be unable to leave effectively.

2.4.2 Force balance in an interface of finite width. Let us now go beyond the previous sharp interface picture and see what happens within the interfacial region from the view point of particle dynamics. Ignoring fluctuations, the stability of the interface means that if a particle happens to be lying at the interface it will, on average, remain put. In our minimal model, a particle can only move due to two reasons: (i) its own active force driving it to move in the direction dictated by its orientation, and (ii) a repulsive force that pushes it away from its neighbours if it is less than a unit distance away from them. While the second force is common in both active and passive (equilibrium) systems, the active force is unique to non-equilibrium systems. Consider now a particle located at x₀ inside the interfacial region, i.e., where 〈v_x(x₀)〉 is varying (Fig. 2). Since 〈v_x(x₀)〉 < 0, the active force will on average drive this particle to the left. Moreover, there are repulsive forces coming from the neighbouring particles on the right hand side f_r(x₀ + Δx). For the particle to remain still, the sum of these forces must be countered by the repulsive forces coming from the left. Therefore,

f_r(x₀ − a/2) = f_a〈v_x(x₀)〉 + f_r(x₀ + a/2). (10)

Since the repulsive forces come from the repulsive potential function U, let us replace the repulsive force by the pressure P_p(x) (subscript p for passive) resulting from the corresponding system with the same particle configuration and interaction potentials, but with the active force omitted. Since f_r(x) ≃ aP_p(x), eqn (10) leads to

f_a〈v_x(x₀)〉 ≃ a[P_p(x₀ + a/2) − P_p(x₀ − a/2)] (11)

In principle, P_p(x) depends on the exact configuration of particles in the system, but if one adopts the simplifying assumption that the passive pressure depends solely on the particle concentration, one can then increase P_p(x) with respect to concentration ρ(x):

P_p(x) = c₀ + c₁ρ(x) + c₂ρ(x)² + [scr O, script letter O](ρ³). (13)

For equilibrium fluids, this is of course the virial expansion, where c₀ = 0 and c₁ = k_BT.²⁷ Since our system is fundamentally not in equilibrium (no translational Brownian motion, i.e., k_BT = 0), there is no guarantee that the same would apply here. But let us assume that such an expansion is possible in our system, and that P_p develops purely due to the repulsive interactions between particles. We thus expect that c₀ = 0 = c₁ because as the concentration goes to zero, there would not be any pairwise interactions. Therefore, the first non-trivial term in the expansion is c₂ρ². Note that c₂ > 0 since P_p arises purely due to the repulsive interactions. One could also incorporate active pressure into the analysis as, for example, done by Winkler et al.³²

Here, to order [scr O, script letter O](ρ²), eqn (12) then leads to

Remarkably, the simulation result shown in Fig. 2(b) indeed seems to vindicate eqn (14).

3 Circular interface

We have seen in the previous section how pressure balance is achieved at a flat interface. However, previous 2D simulation studies have shown that similar to equilibrium phase separation, if the condensed phase in MIPS does not span the system size, the condensed phase is circular. Here, we will see how the curvature of the interface affects the particle dynamics at the interface, and its consequence in terms of the coarsening dynamics. We will first study the emergence of the Gibbs–Thomson relationship using dimensional analysis.

3.1 Gibbs–Thomson relationship: dimensional analysis

In equilibrium phase separation, the Gibbs–Thomson (GT) relationship dictates that the concentration ϕ_R right outside a droplet (of the condensed phase) of radius R iswhere ϕ₀ is the supersaturation concentration, i.e., the threshold concentration beyond which phase separation occurs, and is the capillary length with γ being the surface tension and v being the volume of the molecule. Since the concentration in the vapour phase outside a big drop is lower than that outside a small drop, a diffusive flux is set up which transfers the material from the small droplet to the big droplet. This is the Ostwald ripening mechanism that dominates the phase separation kinetics at the late stage for systems with a small supersaturation.³³

I will now discuss why the GT relationship would arise naturally in our active system. In the minimal MIPS system considered, the only parameters in the dynamical equations are the free roaming speed u = f/η, the rotational diffusion coefficient D_r, and the length scale of the short range area exclusion interaction a. Denoting now the vapour density far from a drop of radius R by ρ_R*, and the density inside the drop (the condensed phase) by ρ_c, then by dimensional analysis we have

where [scr F, script letter F] is some unknown scaling function dependent on its two dimensionless arguments. If we now assume that [scr F, script letter F] is regular with respect to the second argument in the sense that a Taylor series expansion exists (around u/D_rR = 0), then the ratio above can be re-expressed aswhere H and K are now just dimensionless functions of the first argument (u/D_ra), i.e., R-independent. In terms of ρ_∞*, eqn (17) can be re-written aswhere , which may be termed as the effective capillary length. In the large R limit, eqn (18) becomes exactly the GT relationship in eqn (15). This analysis provides an intuitive reason why one would naturally expect the GT relationship to emerge as the drop radius increases in MIPS.

3.2 Gibbs–Thomson relationship: numerics

I will now test eqn (18) by simulating a coarse-grained model of MIPS. Let us first consider what happens to an active particle in the condensed phase at a curved periphery (Fig. 3(a)). For such a particle, I assume that it occupies a zone of radius ã (shown in red) sandwiched by two zones occupied by two neighbouring particles (light blue). Note that since the concentration at the interface may not reach the level of optimal packing concentration (≃0.91), ã should be greater than the particle's diameter a. Indeed, Fig. 2(a) suggests that the concentration is around 0.72 at the interface, which indicates that ã ≃ 1.2. To incorporate the effects of the neighbouring particles on the pink particle, I assume that as a result of the caging effect, the particle can only move out of the droplet if its orientation is within the 2σ range depicted. Based on the diagram shown in Fig. 3, a simple trigonometric exercise leads toAs expected, a lower curvature leads to a smaller escape range (smaller σ_R).

Fig. 3 (a) A schematic of the interfacial condition at a curved interface of curvature R⁻¹. The particle is assumed to occupy a zone of diameter ã and the particle can leave the droplet if its orientation is within the escape range of 2σ_R indicated. (b) A schematic showing a droplet (blue) in the condensed phase of the radius R located at the origin co-existing with the dilute medium (vapour phase). An active particle (red circle) with orientation θ (blue arrow) is located at the position (r, φ). The angle ψ is equal to the difference between the orientation θ and the azimuthal coordinate φ. (c) A snapshot of a simulated system with 1000 active point particles (red dots with orientations indicated by blue arrows) in an annular system with inner radius R = 100 and outer radius R + L_r = 150.

To analyse how the variation in the escape orientation range affects the phase separated system at the steady-state, I have considered here a system with one condensed drop of radius R co-existing with the vapour phase (Fig. 3(b)). I have denoted the particle distribution function in the vapour phase by p_R(r, φ, θ), where the first two arguments correspond to the particles' locations and the last argument corresponds to the particles' orientations. Due to rotational symmetry, one can eliminate one angular argument by introducing the variable ψ ≡ θ − φ,³⁴ and study instead the reduced distribution function . On the drop's periphery, the corresponding reduced distribution function is denoted by χ_R(ψ). Since the periphery is assumed to be infinitely thin, χ_R is only a function of ψ and is therefore dimensionless. In addition, I have assumed that drops of all sizes have the same interior concentration ρ_c, χ_R(ψ) and are thus related to the ρ_c as follows: .

To study the distribution functions, I assume again that the vapour concentration is low enough that pairwise repulsive interactions can be ignored, and simulate the dynamics of active point particles, i.e., non-interacting active particles, in an annular geometry of inner radius R and outer radius R + L_r (Fig. 3(b)). As in the linear case, the particles' orientations are randomised if they reach the outer circular boundary, while if the i-th particle reaches the inner boundary, its positions will remain fixed until its orientation is within the escape range, i.e., until ψ_i is between σ_R and σ_R. The simulation results are shown in (Fig. 4). Away from the interface, it is observed that the concentration rapidly reaches a stationary value for, say, r > R + 20 (Fig. 4(a)). As expected from previous discussion, the vapor concentration ρ_R* decreases with R since a flatter interface leads to a narrower escape range, which leads to a smaller outflux of particles from the condensed phase. Fig. 4(b) shows that ρ_R* decays to ρ_∞* (the vapour concentration as R → ∞) like R⁻¹, which, as we have seen, is consistent with the Gibbs–Thomson relationship in equilibrium systems.

Fig. 4 (a) The variation in the vapour concentration away from the interface of droplets with sizes R = 20, 60, 100. The concentration decays rapidly from the interface and reaches a stationary value ρ_R* a short distance away from the interface. (b) ρ_R*/ρ_∞* vs. R. The curve decays to 1 like R⁻¹, which is consistent with the Gibbs–Thomson relationship (eqn (18)). Blue crosses represent the simulation results and the red curve depicts the function 1 + 1.84R⁻¹. Note that the constant ρ_∞* is estimated numerically from the Gibbs–Thomson relationship. See A.2 for simulation details.

4 Fluctuating interface

I have so far ignored fluctuations in the interfacial profile. In reality, the interface of course fluctuates, which is already discernible from the spatially constrained system shown in Fig. 1(b). In particular, previous simulation results point to the scaling law:¹⁸

w_L² ∼ L, (20)

where w_L is the steady state interfacial width:with h̄ being the average position of the interface. Here, the symbol h̃(y) denotes the location of the interface, i.e., the location of the peak of 〈v_x〉 (Fig. 2). To understand the scaling observed, let us consider the effects of the interface curvature on the particle exchange dynamics. Although the previous section focuses only on a circular interface, i.e., the curvature is positive, one can easily extend eqn (19) to allow for a concave interface as well (Fig. 5). The physical motivation behind the formula is the same, a particle at a highly convex portion of the interface will have a wider escape orientation range (red particle in Fig. 5) than a particle at a highly concave interface (green particle).

Fig. 5 Schematic depicting a wavy interface where the condensed phase is depicted in blue. The location of the interface (purple) is given by the function h̃(y). Due to the caging effect from neighbouring particles, the red particle at the interface will have a higher chance of escaping compared to the green particle because the escape orientation range (grey area) is bigger.

Since the fluctuations ultimately come from the fluctuating dynamics of particle exchange at the interface, we need to model the steady-state dynamics of h̃ stochastically. The simplest equation of motion (EOM) for the interface that incorporates both the effects of stochasticity and curvature-modified outflux is

where α denotes the rate of particles coming at the interface, and thus contributing to the growth of h̃, while α′(1 + β′κ(y)) denotes the rate of the particles escaping, with the effect of the local interface curvature (κ(y) = ∂²h̃/∂y²) taken into account. The noise terms are b_i(x) which are Markovian, spatially independent and are either 0 or 1 with equal probability. Since the interface does not move in the steady state by assumption, α has to be identical to α′. From now on, I will focus exclusively on the hydrodynamic limits (large temporal and spatial scales). So let us coarse grain h̃ by defining a new coarse-grained height function h(y):where a ≪ [small script l] ≪ L_y and [small script l] is large enough that become Gaussian as a result of the central limit theorem. The EOM of h(y) is thenwhere g_i are now Gaussian noises such that

〈g_i(y,t)〉 = 0 (26)

In the long-wavelength limit, the fluctuating term αβκg₂ ∼ ∂²h/∂y² → 0 and so the only relevant fluctuations come from the Gaussian fluctuations α(g₁ − g₂). Therefore, in the hydrodynamic limits, eqn (25) is exactly the Edwards–Wilkinson model,³⁵ with the effective surface tension given by αβ/2. Note that the effective surface tension here is always positive, which is a requirement for having a stable interface. As such this effective surface tension is distinct from the mechanical definition of surface tension, which has been shown to be negative in MIPS.¹⁸ As mentioned in Section 2.2, the negative tension around the interface is likely to be caused by the particles escaping from the condensed phase that are now travelling close to being parallel to the interface. Consistent with our definition of the location of the interface (Fig. 2), these escaped particles are not considered to be part of the condensed phase and are thus ignored in our discussion of the interfacial fluctuations. In other words, the effective surface tension derived here is distinct from the mechanical definition of the surface tension discussed by Bialké et al.¹⁸

Given that our stochastic model is equivalent to the Edwards–Wilkinson model, the temporal and steady-state dynamics of interfacial width are known to follow the scaling form:

for some scaling functions [scr F, script letter F]_EW(·) (Fig. 6). In a 2D system where the interface is a line, α = 1/2 and β = 1/4.^35,36 As shown in Fig. 6, these expectations are confirmed with direct simulation of a discretised version of the original EOM in eqn (22). This model thus provides an analytical argument supporting the steady state scaling w_L(t → ∞) ∼ L^1/2 recently observed numerically.¹⁸

Fig. 6 Interface fluctuations as measured by the interface width . The curve collapse of systems with difference linear dimension L upon rescaling is predicted using the EW model. Inset plot: The interface width at the final time w_L(t_f) shows a linear dependence with , where L is the system size. Blue crosses represent the simulation results and the red line is a guide for the eyes. See Appendix A.3 for simulation details.

To summarise this section, I have incorporated the caging effect as discussed in Section 3 into the modelling of the stochastic dynamics of particle exchanges at the interface. The model equation is then shown to be equivalent to the Edward–Wilkinson model in the hydrodynamic limits. In particular, the emergence of the effective surface tension term (αβκ/2) from the particle dynamics at the interface also explains why the interface is flat when both phases span the system, and circular when one phase does not span the system.

5 MIPS in 3D

I have so far analysed the interfacial properties in MIPS in 2D using a combination of analytical and numerical methods. Here I will extrapolate the obtained results to MIPS in 3D.

5.1 Flat interface

Utilizing the sharp interface model for MIPS in 3D, the pressure balance equation in eqn (8) becomeswhere on the R. H. S. the active force contribution corresponds to the swim pressure of active particles in 3D in the vapour phase,²³ and on the L. H. S., the first term comes from the active contribution to the force assuming again that the orientation is isotropic in the condensed phase. With regards to the minimal active force required for MIPS, using the approximations that f_r is negligible and that ρ_c/ρ_v ≃ 1, we arrive atwhich is very similar to the expression in 2D (eqn (9)).

5.2 Spherical interface

The dimensional analysis presented in Section 3.1 also includes spherical drops in 3D. Therefore, if the assumption that the escape range decreases with the curvature of the drop, then the Gibbs–Thomson relationship should emerge in the large drop limit (eqn (18)). In particular, we again expect the MIPS coarsening kinetics to be equivalent to the equilibrium scheme at low supersaturation.^15,33,37

5.3 Interface fluctuations

For MIPS in 3D, the interface is two dimensional and so there are two principal curvatures. If we adopt the natural assumption that the escape range now depends on the mean curvature, then the theoretical analysis presented in Section 4 applies in 3D in a straightforward manner. As a result, keeping only the linear terms will again lead to the Edwards–Wilkinson model in the hydrodynamic limits.

6 Discussion and outlook

In this paper I have investigated the microscopic dynamics of active particles in the interfacial regions in MIPS using a combination of simulations and analytical arguments, and demonstrated (i) how interface stability is achieved, (ii) why the GT relationship emerges in MIPS, and (iii) how interface fluctuations scale with the system size. Therefore, I have shown that all the observed “surface tension” related phenomena found in MIPS result from the microscopic dynamics of the active particles. More specifically, I have demonstrated that pressure balance in MIPS is achieved because of the orientation anisotropy in the region, which leads to a high active force directed towards the condensed phase. By incorporating the caging effects of neighbouring particles in the peripheral of a condensed drop, I have shown how the Gibbs–Thomson relationship emerges naturally in MIPS, which dictates that the larger the condensed drop, the smaller the vapour concentration outside the drop. If the supersaturation level is small, the GT relationship leads to diffusive transfer of active particles from small drops to larger drops. As a result, the late-stage coarsening kinetics in an active phase-separating system should follow the temporal scaling as in equilibrium phase separation: i.e., the average droplet size in the system 〈R(t)〉 scales as^15,33,37t^1/3. In addition, the droplet size distribution should approach asymptotically the universal size distribution obtained by Lifshitz and Slyozov.³³ Lastly, motivated by the same caging effects, I have proposed a stochastic model that describes the interfacial fluctuations in MIPS. In this model, the probability of particles leaving the interface is assumed to be dependent on the interface curvature. An analytical argument was then provided to show that the proposed model belongs to the same universality class as the Edwards–Wilkinson model.

There are a number of future research directions that are of interest. For instance, phase separation may play a role in cytoplasmic re-organisation during asymmetric cell division.^38,39 How the activity in the cytoplasm due to the many motor proteins contributes to such re-organisation via phase separation awaits more attention. Moreover, the fact that active phase separation occurs naturally raises the question of the existence of the critical point as in the equilibrium case. The critical transition in incompressible active fluids has recently been shown to give rise to a novel universality class.⁴⁰ Furthermore, if a critical transition exists in MIPS, will the critical exponents be identical to those in the equilibrium case which belong to the Ising universality class? This question awaits further investigation.

Appendix A Simulation details

A.1 Particle dynamics simulation

For the simulation results reported in Section 2 (Fig. 1 and 2), I numerically integrate the Langevin equations for each particle in the bulk of the system of size L_x × L_y by using the following updated equations:where g^t_i are Gaussian distributed random variables with zero mean and unit variance, and U is given by the Weeks–Chandler–Andersen potential shown in eqn (3). The periodic boundary condition is enforced in the y, direction; while in the x direction, if x^t+Δt_i < 0, then it is reset to zero, and if x^t+Δt_i > L_x, then x^t+Δt_i = L_x and θ^t+Δt_i is an angle chosen at random.

The parameters of the simulations are: Δt = 10⁻⁵, a = 1, η = 1, f_a = 100, D_r = 3, and A = 25/6 for repulsive particles and A = 0 for point particles. The system has width L_x = 50 and height L_y = 10 sin(π/3), with 300 particles initialized in the configuration of a hexagonal lattice (with spacing 1) next to the left boundary, and random orientations. Two hundred million time steps are evolved to equilibrate the system and then data are collected in the subsequent two hundred million time steps.

A.2 Point particles in an annular geometry

For the simulation results reported in Section 3 (Fig. 3), the system is now annular with inner radius R and outer radius R + L_r. The same updates as in eqn (31) and (32) are included for the point particles in the bulk of the system. Concerning the boundary conditions, if r^t ≤ R then the particle becomes part of the condensed drop periphery, and the orientation follows the update:while the position remains the same until ψ^t ≡ |θ^t − φ^t| < σ_R (Fig. 3(b)), in which case the particle leaves the condensed phase. If r^t+Δt_i > R + L_r, then r^t+Δt_i = R + L_r and θ^t+Δt_i is an angle chosen at random.

For distinct annulus geometry, the density ρ_R(r) is normalised by ρ_c, which is the density of particles at the inner circular wall. This is based on the assumption that ρ_c is the same irrespective of drop sizes.

The parameters of the simulations are: Δt = 10⁻³, η = 1, f_a = 100, D_r = 3, and A = 0. The system has 1000 particles initialized with random orientations and positions. Twenty million time steps are evolved to equilibrate the system and then data are collected in the subsequent twenty million time steps. Simulations are performed for R = 20, 40, 60, 80 and 100, while L_r is always 50.

A.3 Fluctuating interfaces

To simulate interface fluctuations according to eqn (22), I discretise the interface (a line) into L sites with height values h_i where i = 1,…,L. Periodic boundary conditions are enforced. The updates are performed as follows:

h^t+Δt_i = h^t_i + [1 + β(h^t_i−1 + h^t_i+1 − 2h^t_i)]g̃^t_i − [1 − β(h^t_i−1 + h^t_i+1 − 2h^t_i)]g̃^t_i, (34)

where g̃^t_i are either 0 or 1 chosen with equal probability, and β = 0.1. The system sizes simulated are L = 160, 200, 240, 280, 320 and 360. The total number of time steps simulated for each system size is 5L²/2.

Acknowledgements

I thank Christoph Weber (Max Planck Institute for the Physics of Complex Systems), Fernando Peruani (Université Nice Sophia Antipolis) and Richard Sear (University of Surrey) for stimulating discussions and for their comments on the manuscript.

References

1. C. Domb, J. L. Lebowitz and M. S. Green, Phase Transitions and Critical Phenomena, Academic Press, 1983.
2. A. J. Bray, Adv. Phys., 2002, 51, 481–587.
3. J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2008, 100, 218103.
4. M. E. Cates, D. Marenduzzo, I. Pagonabarraga and J. Tailleur, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 11715–11720.
5. F. Peruani, T. Klauss, A. Deutsch and A. Voss-Boehme, Phys. Rev. Lett., 2011, 106, 128101.
6. J. P. D. Clewett, K. Roeller, R. M. Bowley, S. Herminghaus and M. R. Swift, Phys. Rev. Lett., 2012, 109, 228002.
7. F. D. C. Farrell, M. C. Marchetti, D. Marenduzzo and J. Tailleur, Phys. Rev. Lett., 2012, 108, 248101.
8. Y. Fily and M. C. Marchetti, Phys. Rev. Lett., 2012, 108, 235702.
9. J. Schwarz-Linek, C. Valeriani, A. Cacciuto, M. E. Cates, D. Marenduzzo, A. N. Morozov and W. C. K. Poon, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 4052–4057.
10. M. E. Cates and J. Tailleur, EPL, 2013, 101, 20010.
11. G. S. Redner, M. F. Hagan and A. Baskaran, Phys. Rev. Lett., 2013, 110, 55701.
12. J. Bialké, H. Löwen and T. Speck, EPL, 2013, 103, 30008.
13. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
14. J. Stenhammar, A. Tiribocchi, R. J. Allen, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2013, 111, 145702.
15. R. Wittkowski, A. Tiribocchi, J. Stenhammar, R. J. Allen, D. Marenduzzo and M. E. Cates, Nat. Commun., 2014, 5, 4351.
16. A. Y. Grosberg and J. F. Joanny, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 32118.
17. J. Stenhammar, R. Wittkowski, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2015, 114, 18301.
18. J. Bialké, J. T. Siebert, H. Löwen and T. Speck, Phys. Rev. Lett., 2015, 115, 98301.
19. R. P. Sear, 2015, arXiv:1503.06963.
20. A. Marchand, J. H. Weijs, J. H. Snoeijer and B. Andreotti, Am. J. Phys., 2011, 79, 999–1008.
21. J. D. Weeks, D. Chandler and H. C. Andersen, J. Chem. Phys., 1971, 54, 5237–5247.
22. See eqn (15), (17) and (18) in C. F. Lee, New J. Phys., 2013, 15, 055007.
23. S. C. Takatori, W. Yan and J. F. Brady, Phys. Rev. Lett., 2014, 113, 28103.
24. X. Yang, M. L. Manning and M. C. Marchetti, Soft Matter, 2014, 10, 6477–6484.
25. S. A. Mallory, A. Šarić, C. Valeriani and A. Cacciuto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 89, 52303.
26. A. P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri, M. Kardar and J. Tailleur, Nat. Phys., 2015, 11, 673–678.
27. J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, Academic Press, 2006.
28. P. G. de Gennes and J. Prost, The Physics of Liquid Crystals International Series of Monographs on Physics, Oxford University Press, USA, 1995.
29. A. P. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri, M. E. Cates and J. Tailleur, Phys. Rev. Lett., 2015, 114, 198301.
30. An anonymous referee has indicated that numerical simulation shows that P_r is about a third of the other terms at the onset of MIPS. However, eqn (9) remains valid as P_r does not depend on the Péclet number at the onset, as shown in the preceding paper.
31. J. Stenhammar, D. Marenduzzo, R. J. Allen and M. E. Cates, Soft Matter, 2014, 10, 1489–1499.
32. R. G. Winkler, A. Wysocki and G. Gompper, Soft Matter, 2015, 11, 6680–6691.
33. I. Lifshitz and V. Slyozov, J. Phys. Chem. Solids, 1961, 19, 35–50.
34. A. Pototsky and H. Stark, EPL, 2012, 98, 50004.
35. S. F. Edwards and D. R. Wilkinson, Proc. R. Soc. London, Ser. A, 1982, 381, 17–31.
36. A. L. Barabási and H. E. Stanley, Fractal Concepts in Surface Growth, Cambridge University Press, 1995.
37. J. Yao, K. R. Elder, H. Guo and M. Grant, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 8173–8176.
38. C. Brangwynne, C. Eckmann, D. Courson, A. Rybarska, C. Hoege, J. Gharakhani, F. Jülicher and A. Hyman, Science, 2009, 324, 1729–1732.
39. C. F. Lee, C. P. Brangwynne, J. Gharakhani, A. A. Hyman and F. Julicher, Phys. Rev. Lett., 2013, 111, 088101.
40. L. Chen, J. Toner and C. F. Lee, New J. Phys., 2015, 17, 042002.

---