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

Minimal models of self-propelled particles with short-range volume exclusion interactions have been shown to exhibit 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 relation in MIPS, which provides a functional relationship between the size of a condensed drop and its surrounding vapor concentration. As a result, the coarsening dynamics of MIPS at vanishing supersaturation follows the classic Lifshitz-Slyozov scaling law at the late stage.


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][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][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][9][10][11][12][13][14][15][16][17][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. 20 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. 18 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 are where i is an integral index enumerating the particles in the system, v i cos y i x + sin y i ŷ with the angle y 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, Z is the drag coefficient and f a is the constant active force that drives the particles in the system. In particular, u f a /Z 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: 21 UðrÞ ¼ 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.

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 y C p. This results from the fact that only left-going particles remain at the wall. Just outside the wall, the distribution is highly peaked at y just below p/2 and just above 3p/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.
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 W (y). Note that since we are dealing with a 2D system, the unit of pressure is To further analyse w W (y), one can perform dimensional analysis to conclude that where F(y) is a function dependent only on y, and r N is the particle concentration far from the wall. To obtain the exact functional form of F(y), one needs to solve a set of two coupled differential equations under mixed boundary conditions, 22 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, F(y) can be readily estimated numerically (left column of Fig. 1), which allows one to obtain the following equation: FðyÞ cos ydy ¼ f a 2 r 1 2ZD r : 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.

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 C1/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 C 17 in Fig. 1(f)). This represents the accumulation of particles with an orientation highly centred at around y C p, 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 r 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 2 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).

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. 27 In the case discussed here, since the pronounced minimum of hv x i hcos yi is easy to locate (indicated by the red broken line) and its location also marks the onset of the increase in Q yy Àhcos(2y)/2i (the yy component of the nematic order parameter Q), 28 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. 18 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.

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. (c) The yy component of the nematic order parameter Q vs. x, where Q yy Àhcos(2y)/2i 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. 27 The working definition proposed here is that the interface is set to be at the pronounced minimum of |hv x i|. View Article Online 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 1801 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 2 r v /(2ZD r ) (eqn (6)); this pressure is balanced by the pressure exerted by the condensed phase: ar c f a /p + 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), ar c f a /p C 30, P r C 230 and r v f a 2 /(2ZD r ) C 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 , 30 and (ii) the concentration ratio between the condensed phase and the dilute phase (r c /r 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. 11 The result also supports the notion that the Péclet number (Pe), usually defined as Pe p f a /(aZD 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 2 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 0 inside the interfacial region, i.e., where hv x (x 0 )i is varying (Fig. 2). Since hv x (x 0 )i o 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 0 + Dx). For the particle to remain still, the sum of these forces must be countered by the repulsive forces coming from the left. Therefore, 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) C aP p (x), eqn (10) leads to ' a 2 dP p ðx 0 Þ dx : 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 r(x): For equilibrium fluids, this is of course the virial expansion, where c 0 = 0 and c 1 = k B T. 27 Since our system is fundamentally not in equilibrium (no translational Brownian motion, i.e., k B T = 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 = 0 = c 1 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 2 r 2 . Note that c 2 4 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. 32 Here, to order O(r 2 ), eqn (12) then leads to Remarkably, the simulation result shown in Fig. 2(b) indeed seems to vindicate eqn (14).

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.

Gibbs-Thomson relationship: dimensional analysis
In equilibrium phase separation, the Gibbs-Thomson (GT) relationship dictates that the concentration f R right outside a droplet (of the condensed phase) of radius R is where f 0 is the supersaturation concentration, i.e., the threshold concentration beyond which phase separation occurs, and n ¼ 2gv k B T is the capillary length with g 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. 33 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/Z, 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 R *, and the density inside the drop (the condensed phase) by r c , then by dimensional analysis we have where F is some unknown scaling function dependent on its two dimensionless arguments. If we now assume that F is regular with respect to the second argument in the sense that a Taylor series expansion exists (around u/D r R = 0), then the ratio above can be re-expressed as where H and K are now just dimensionless functions of the first argument (u/D r a), i.e., R-independent. In terms of r N *, eqn (17) can be re-written as whereñ Ku HD r , 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.

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 (C0.91), ã should be greater than the particle's Fig. 3 (a) A schematic of the interfacial condition at a curved interface of curvature R À1 . 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 2s 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 y (blue arrow) is located at the position (r, j). The angle c is equal to the difference between the orientation y and the azimuthal coordinate j.
(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.  Fig. 2(a) suggests that the concentration is around 0.72 at the interface, which indicates that ã C 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 2s range depicted. Based on the diagram shown in Fig. 3, a simple trigonometric exercise leads to As expected, a lower curvature leads to a smaller escape range (smaller s R ).
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, j, y), 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 c y À j, 34 and study instead the reduced distribution function z R ðr; cÞ Ð 2p 0 p R ðr; j; j þ cÞr cos jdj. On the drop's periphery, the corresponding reduced distribution function is denoted by w R (c). Since the periphery is assumed to be infinitely thin, w R is only a function of c and is therefore dimensionless. In addition, I have assumed that drops of all sizes have the same interior concentration r c , w R (c) and are thus related to the r c as follows: Ð dcw R ðcÞ ' 2ãRr c . 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 c i is between s R and s R . The simulation results are shown in (Fig. 4). Away from the interface, it is observed that the concentration rapidly reaches a stationary value r R Ã 1 r Ð dcz R ðr; cÞ for, say, r 4 R + 20 ( Fig. 4(a)). As expected from previous discussion, the vapor concentration r 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 R * decays to r N * (the vapour concentration as R -N) like R À1 , which, as we have seen, is consistent with the Gibbs-Thomson relationship in equilibrium systems.

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: 18 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 hv x i (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).
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 a denotes the rate of particles coming at the interface, and thus contributing to the growth of h, while a 0 (1 + b 0 k(y)) denotes the rate of the particles escaping, with the effect of the local interface curvature (k(y) = q 2 h/qy 2 ) 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, a has to be identical to a 0 . From now on, I will focus exclusively on the hydrodynamic Fig. 4 (a) The variation in the vapour concentration r R ðrÞ 1 r Ð dcz R ðr; cÞ 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 R * a short distance away from the interface. (b) r R */r N * vs. R. The curve decays to 1 like R À1 , 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 À1 . Note that the constant r N * is estimated numerically from the Gibbs-Thomson relationship. See A.2 for simulation details.
limits (large temporal and spatial scales). So let us coarse grain h by defining a new coarse-grained height function h(y): where a { l { L y and l is large enough that Gaussian as a result of the central limit theorem. The EOM of h(y) is then where g i are now Gaussian noises such that hg i (y,t)i = 0 (26) In the long-wavelength limit, the fluctuating term abkg 2 B q 2 h/qy 2 -0 and so the only relevant fluctuations come from the Gaussian fluctuations a(g 1 À g 2 ). Therefore, in the hydrodynamic limits, eqn (25) is exactly the Edwards-Wilkinson model, 35 with the effective surface tension given by ab/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. 18 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. 18 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 F EW (Á) (Fig. 6). In a 2D system where the interface is a line, a = 1/2 and b = 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 -N) B L 1/2 recently observed numerically. 18 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 (abk/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.

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.

Flat interface
Utilizing the sharp interface model for MIPS in 3D, the pressure balance equation in eqn (8) becomes 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. where on the R. H. S. the active force contribution corresponds to the swim pressure of active particles in 3D in the vapour phase, 23 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 r c /r v C 1, we arrive at which is very similar to the expression in 2D (eqn (9)).

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

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.

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 hR(t)i scales as 15,33,37 t 1/3 . In addition, the droplet size distribution should approach asymptotically the universal size distribution obtained by Lifshitz and Slyozov. 33 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. 40 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.

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+Dt i o 0, then it is reset to zero, and if x t+Dt i 4 L x , then x t+Dt i = L x and y t+Dt i is an angle chosen at random. The parameters of the simulations are: Dt = 10 À5 , a = 1, Z = 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(p/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 R then the particle becomes part of the condensed drop periphery, and the orientation follows the update: while the position remains the same until c t |y t À j t | o s R (Fig. 3(b)), in which case the particle leaves the condensed phase. If r t+Dt i 4 R + L r , then r t+Dt i = R + L r and y t+Dt i is an angle chosen at random.
For distinct annulus geometry, the density r R (r) is normalised by r c , which is the density of particles at the inner circular wall. This is based on the assumption that r c is the same irrespective of drop sizes.
The parameters of the simulations are: Dt = 10 À3 , Z = 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: where g t i are either 0 or 1 chosen with equal probability, and b = 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 /2.