Thermodynamics predicts a stable microdroplet phase in polymer-gel mixtures undergoing elastic phase separation

We study the thermodynamics of binary mixtures with the volume fraction of the minority component less than the amount required to form a flat interface and show that the surface tension dominated equilibrium phase of the mixture forms a single macroscopic droplet. Elastic interactions in gel-polymer mixtures stabilize a phase with multiple droplets. Using a mean-field free energy we compute the droplet size as a function of the interfacial tension, Flory parameter, and elastic moduli of the gel. Our results illustrate the role of elastic interactions in dictating the phase behavior of biopolymers undergoing liquid-liquid phase separation.


I. INTRODUCTION
Membraneless compartmentalisation in cells that are driven by phase-separation processes due to changes in temperature or pH, and maintained by a non-vanishing interfacial tension, is one of the most exciting recent biological discoveries [1][2][3][4]. These membraneless compartments composed bio-molecular condensates have been implicated in important biological processes such as transcriptional regulation [5], chromosome organisation [6] and in several human pathologies e.g. Huntington's, ALS etc. [7]. Self-assembly processes that lead to organelle formation however need to be tightly regulated such that the phase separated droplets do not grow without bound and remain small compared to the cell size. Understanding the regulatory processes that controls droplet size in cellular environments is therefore a crucial interdisciplinary question. The two candidate mechanisms proposed for arresting droplet growth are (i) incorporation of active forces that break detailed balance [8,9], and (ii) non-equilibrium reaction mechanisms which couple to the local density field [4]. Although biological systems are inherently out of equilibrium, an estimation of diffusion constant of bio-molecules indicates that non-equilibrium effects are negligible on length-scales beyond microns and timescales beyond microseconds. Hence, the framework of equilibrium thermodynamics can be readily applied to analyse biological phase separation in cells [10].
For synthetic polymer mixtures, in the absence of active processes, droplet growth is limited by the elastic interactions of the background matrix that alters the thermodynamics of phase separation [11][12][13]. Recent experiments on mixtures of liquid PDMS and fluorinated oil in a matrix of cross-linked PDMS show the dependence of the droplet size on the nucleation temperature and the network stiffness [14,15]. Despite theoretical attempts [16,17] a complete understanding of elasticity mediated arrested droplet growth is still lacking.
In this paper, we develop a consistent thermodynamic formalism to compute the equilibrium radius of the droplet of the minority phase in (a) binary polymer, and (b) a polymer-gel mixture, using mean-field theories utilising the Flory-Huggins [29] and the Flory-Rehner [30] free energies, respectively. A parallel tangent construction for droplets, used to obtain the densities of coexisting phases is presented. This procedure is a generalisation of the common tangent construction for flat interfaces and in the thermodynamic limit allows us to compute the equilibrium radius of a single droplet. For phase separation processes in mixtures with a gel component, elastic interactions limit droplet growth stabilising a phase with multiple droplets, in the correct parameter regime [31].

II. MODEL
Consider a binary mixture of a gel and a solvent, close to but below the gelation temperature, where the gel-formation and the phase-separation are competing processes. The physical system considered is different from the recent experimental [14,15] and theoretical studies [16,17]. The experiments have been performed on mixtures of liquid PDMS (uncrosslinked PDMS polymers) and flourinated oil in a matrix of cross-linked PDMS, thus it is a ternary system. The two previous theoretical attempts [16,17] however approach this by describing the thermodynamics of a binary mixture (oil and uncrosslinked PDMS) in the background of the elastic matrix (crosslinked PDMS), where the volumefraction of the matrix does not enter the calculation. The matrix only provides an elastic background in which the phase separation of the binary mixture occur. On the other hand, the elastic matrix is considered in reference [14,15], but the translational entropy of the gel has been explicitly put to zero. However, this is a contentious issue, as we discuss later in the manuscript, and it leads to unstable solutions for a binary mixture of a gel and a solvent. References [14,15] does not encounter this issue as they do not perform the parallel tangent construction, which is a condition that arises from the minimisation of the free-energy, and they bypass this by assuming that the dispersed microdroplets of the solvents can be described as an ideal gas.
The thermodynamic formalism to understand phase separation is as follows: an unstable mixture of composition φ 0 splits into two coexisting phases in a slab-like geometry respecting volume and mass conservation, with the equilibrium configuration being a minimum of the free energy ( Fig. 1(a)). The volume fraction of the two coexisting phases are φ in , and φ out respectively, with V d denoting the volume occupied by phase with density φ in , and F b (φ) is the Helmholtz free-energy per unit volume (in units of k B T a 3 ). The solvent fraction is f = V d V , and the free-energy density F(φ) of the planar configuration ( Fig. 1(a)) is given by, where F s = 2γV −1/3 , corresponds to the surface energy with γ being the surface tension, V the volume of the system considered, and λ a Lagrange multiplier that enforces the mass conservation constraint. A calculation of the equilibrium thermodynamics proceeds via minimising the free energy in Eq. (1) w.r.t to the independent quantities φ in , φ out , f , and λ. The constrained minimization of the free-energy function in Eq. (1) w.r.t. φ in , φ out and f leads to the common tangent construction where µ(φ), and Π(φ) refers to the exchange chemical potential and the osmotic pressure of the phases respectively. Eq. (2) ensures chemical, and mechanical equilibrium (see SI). Thermal equilibrium is ensured as calculations are carried out in a constant temperature ensemble. We obtain coexistence volume fractions φ in and φ out from Eq. (2). The solvent fraction f is obtained by minimising the functional w.r.t λ, i.e ∂F(φ in , φ out , f, λ)/∂λ = 0, which yields, f = φ0−φout φin−φout . For a planar interface, the surface energy term does not explicitly depend on the solvent volume fraction f . Consequently, the minimisation conditions lead to four uncoupled equations (SI) and a knowledge of the coexistence volume fractions φ in and φ out is enough to determine f . As evident from Eq. (1), the effect of the surface energy term vanishes in the thermodynamic limit, i.e., as volume V → ∞. In contrast, a spherical droplet geometry introduces a non-trivial coupling among the minimisation conditions and a knowledge of the volume, V , of the system is required to obtain the equilibrium configuration.
FIG. 1: A (i) common tangent (solid pink line) and a (ii) parallel tangent (dashed green line) construction for planar interfaces (a) and droplets (b) (with volume V d ) for a binary polymer mixture. A Flory-Huggins functional with χ = 1.2χc, NA = 100, NB = 200 is used. Coexistence volume fractions inside φin and outside φout droplet approach the values obtained for a flat interface φα, and φ β in the thermodynamic limit V → ∞.

A Spherical Droplet
Spherical droplets of the minority phase arise in finite systems when the volume fraction is less than a critical value [32,33]. The thermodynamics in such situations differ from the common tangent construction and leads to the classical Gibbs-Thomson relations [4]. Fig. 1(b) shows an unstable system of volume fraction φ 0 , that phase separates into a background matrix of volume fraction φ out and a single droplet of radius R of volume fraction φ in in a finite box of volume V . Assuming an ansatz of a phase separated mixture comprising of N spherical droplets of identical radius R, (referred to as the microdroplet phase henceforth), the solvent fraction is given by f = N ( 4 3 πR 3 /V ). The free energy of the micro-droplet phase is therefore where F s (f ) = N V 4πR 2 γ, accounts for the interfacial energy between the droplet and the background phase. By imposing the mass conservation constraint and expressing the surface energy in terms of the solvent fraction f , the free energy per unit volume is given by The surface energy of the droplet depends on the solvent fraction f on account of the its spherical shape. The equilibrium conditions therefore lead to four coupled equations, involving the yet unknown system volume V . The chemical and mechanical equilibrium conditions for the micro-droplet phase involving the coexisting densities translates to, µ(φ in ) = µ(φ out ) and Π(φ in ) = Π(φ out ) + 2γ( 4πN 3f V ) 1/3 , where the extra term in the pressure equation accounts for the Laplace pressure acting across the interface. We carry out a minimisa-tion procedure akin to the planar interface to obtain the solvent volume fraction f , and the coexistence volume fractions inside and outside the droplet, φ in and φ out respectively for a given box volume V . In the absence of elastic interactions the equilibrium phase corresponds to a single droplet of the minority phase, i.e. N = 1 in Eq. (3). The radius of the drop is determined in terms of the coexistence densities and is given by where ν = 3(φ0−φout) 4π(φin−φout) 1/3 , and L = V 1/3 is the length of the cubic box. We apply the framework to compute the radius of the minority phase droplet of a binary polymer mixture described by a Flory-Huggins free energy in the thermodynamic limit i.e. V → ∞, performing our calculation for different box volumes V . The surface tension γ for the micro-droplet phase is taken to be the same as that of a planar interface. The thermodynamics of binary polymer mixtures is well described by the Flory-Huggins free-energy where N A , and N B are the lengths of A and B polymers respectively, and χ is the mixing parameter. For χ > χ c , where χ c is the value of the mixing parameter at criticality, the mixture is unstable and spontaneously phase separates into low and high volume fraction phases determined by the minimisation conditions. We consider an unstable polymer mixture with N A = 100, N B = 200, and an initial composition φ 0 = 0.35, and χ = 1.2χ c . Fig. 1(a) shows the common tangent construction which yields the coexistence volume fractions φ α = 0.235 and φ β = 0.885 for a flat interface. If the amount of material is not enough, the minority phase forms a droplet whose coexistence volume fractions outside φ out and inside φ in are determined by the parallel tangent construction ( Fig. 1(b)) as a function of the box volume V . A combination of the parameters φ 0 , φ α and φ β determines that the fraction of the solvent-rich phase f ≈ 0.17. The equilibrium phase is a single drop. To obtain the coexistence volume fractions and the droplet radius in the thermodynamic limit, we perform parallel tangent constructions for cubic boxes of lengths L = 160, . . . 10 4 using Eq. (4). Fig. 2 shows a finite size scaling analysis of the droplet radius R in units of the box-size L, (R/L) as a function of 1/L. The thermodynamic limit 1/L → 0 corresponds to the y-intercept R/L ≈ 0.34 for the Flory parameters listed above. The numerical derivative of R/L w.r.t. L approaches zero in this limit (Fig. 2 inset). The solvent fraction f , is also a function of the systems size (f ∼ (R/L) 3 ). The coexistence densities calculated from Eq. (4) are functions of L and can be quantified in terms of their deviation from the coexistence volume fractions for a planar interface, i.e.,φ in = (φ in − φ β )/φ β and φ out = (φ out − φ α )/φ α . As shown in Fig. 2 φ out → φ α , and φ in → φ β in the thermodynamic limit.

A microdroplet phase
The Helmholtz free-energy per unit volume of the micro-droplet configuration of a gel-solvent mixture, with N droplets (see Fig. 4 inset), is given by where F b (φ) is the Flory-Huggins free energy given by . We consider a situation where the strand length of the gel, N B , is considered to be finite in these calculations (N B = 25 and N A = 1). The reason for this, and not letting N B → ∞, is based on stability arguments and is discussed in the SI and we set χ(T ) = 1.38χ c in our calculations. The surface-energy per unit volume in Eq. (5),  5) can be expressed as a function of the solvent fraction, f , (see SI) and is given by To incorporate the effects of the finite stretch-ability of the gel, we adopt the Gent model [34,35]. The elastic free energy density has the form, with λ's corresponding to the strains in the radial, azimuthal, and polar directions, J m ∼ 10 6 is the stretching limit of the network, and G is the shear modulus. The shear modulus is related to the microscopic parameters via the relation, Surface energy Fs, increases, the elastic energy F el , decreases, whereas, the bulk free energy F b (φ)is almost independent of N as shown in (b). Panel (c) shows the dependence of the droplet radius R and the number density of the droplets n and the shear modulus of the gel G. The inset of panel (a) shows the free-energy as a function of the number of droplets for a system with no elastic interaction. Thus a single macrodroplet is the stable phase.
, where n dry and R 0 are the average cross-link density and the mesh size of the dry gel respectively [36] (see SI). Due to the volume-preserving nature of the deformation, λ r = 1/λ 2 and λ φ = λ θ = λ and its magnitude is bounded, i.e., 0 < J/J m < 1 [34]. The energy minimisation conditions w.r.t the independent variables as outlined earlier, leads to a modified equilibrium conditions: These conditions lead to a set of coupled equations that we solve numerically to yield the four unknown variables, φ in ,φ out ,f , and λ, associated with each droplet number, N . A geometrical interpretation of these equations lead to the construction of parallel tangents.
We substitute the equilibrium values of the coexistence volume fractions and solvent fraction into the original free-energy expression in Eq. (5), to obtain a free energỹ F (N ), as a function of the number of droplets N . The minimisation ofF (N ) w.r.t N yields N m , the optimal number of droplets of the micro-droplet phase. Fig. 3(a) shows the free-energyF (N ) = F g (N )−F g (1) (Eq. (5)) as a function of the number of droplets, once the coexistence volume fractions have been obtained for a cubic box of side L = 200 and the surface tension γ = 1.67 × 10 −3 (in units of k B T /a 2 ). It is evident that this is a convex function, with a well defined minimum occurs around N m ≈ 23. The inset shows the contrasting behaviour ofF (N ) for a binary polymer mixture. In the absence of elastic interactions, surface tension dominates the thermodynamics and a phase with a single droplet is the equilibrium state corresponding to the free energy minimum. The convex nature of the free energỹ F (N ) arises from a balance between the surface, elastic, and bulk free energies of the micro-droplet phase. As the number of droplets N increases, the surface energy monotonically increases on account of the increase of the total interfacial area. In contrast, the elastic energy monotonically decreases as a function of N , since an increase in the number of droplets translates to smaller sized drops and less deformation of the gel matrix. The elastic free energy has a lower bound corresponding to a minimum droplet of size R/L ∼ a, length of a monomer. The combined effect of these two contributions to the free energy therefore stabilizes the micro-droplet phase. The bulk free energy is nearly independent of N . Fig. 3(b) shows the variation of the different components of the total free energy as a function of the number of droplets N , while Fig. 3(c) shows the variation of number density n = N m /V , and droplet radius R as a function of the shear modulus G. The shear modulus G is tuned by varying the mesh size, R 0 , of the gel. We compute the number density by minimizingF (N ) w.r.t. N and determine the drop radius using R(N m ) = ν/N 1/3 m L for a given shear modulus G. As shown, the radii of the droplets decrease (and hence the number density n increases commensurately) as the gel becomes stiffer.
The convex nature ofF (N ) as a function of N is independent of the system size L as shown in Fig. 4(a). Fig. 4(b) shows the dependence ofF (N ) as a function of the surface tension, γ = 0.0025, . . . 0.004, while keeping the shear modulus of the gel-solvent mixture fixed at G = 1.9 × 10 −4 (in units of k B T /a 3 ). The free energy minimum shifts to smaller values of N m with increasing surface tension as shown in Fig. 4(b). The inset of Fig. 4(b) shows that for γ < γ c ≈ 4.0 × 10 −3 , a microdroplet phase is the equilibrium configuration, with N m , monotonically increasing with decreasing γ. Fig. 4(c) shows that the equilibrium number density of droplets n = N m /V , and the droplet radius R(N m ) have reached a thermodynamic limit and are independent of the system size L. Fig. 4(d) shows the phase boundary demarcating regions of a stable macro-droplet and dispersed micro-droplet phases in the γ − G plane. The mean field phase-boundary (symbols) qualitatively agrees with the scaling results [31] (red dashed line) for softer gels while significant deviations are observed for stiffer ones. The mean-field phase boundary (symbols) is now a function of the gel-strand length N B , a variable that is associated with the network heterogeneity of the system. Such quenched disorder dramatically modifies the equilibrium thermodynamics of gel networks.  Fig. 4(d), shows the contour-plot of the dimensionless ratio between the surface energy and the elastic energy, h/α, has been shown in the γ-G plane, where α is equal to 2.5 (see SI for a discussion on this). Also shown is the phase boundary from the mean field theory calculations (inverted triangles, the inverted triangle and the dashed line are sim- ilar to that presented in Fig. 4(d)). Simple scaling arguments would suggest that the phase boundary would occur at h/α equal to unity (see the dashed line in Figure  5 (a)) and we observe that for small values of the shear modulus, G, this is indeed the case. However, as the value of G increases deviations between the mean-field phase boundary (inverted triangles) and the h/α equal to unity increase. In order to facilitate comparison with present and future experiments, we have studied how the equilibrium number of droplets evolve as a function of a tuning parameters (shear modulus or surface tension in this case) as one crosses the phase boundary along the principal directions in the phase plane. Panel (b) shows the transition from a dispersed micro-droplet to a single macro-droplet as one crosses the phase boundary while keeping G fixed and increasing γ. For γ < γ c , the dependence of the number of droplets on the surface tension follows the linear relationship, N m = 31.8(1 − γ γc ). Similarly, panel (c) shows the transition from a single macro-droplet to a dispersed micro-droplet state when one keeps γ constant and increases G and here the dependence of the number of droplets on the elastic modulus again follows a linear dependence N m = 39.2( G Gc − 1).
The linear dependence of the number of droplets on the elastic modulus of the matrix is a result of the mean-field theory calculations (and not an assumptions as in [16]) and has been observed in the experiments [14]. In summary, we consider phase separation in an elastic medium, where the background matrix influences the equilibrium thermodynamics of the mixture. Previous studies consider the background matrix as an inert phase [16,17]. For composition regimes where the solvent is a minority phase and there is a dearth of material to form a flat interface, solvent-rich droplets coexist with the majority phase. We demonstrate, via a mean-field theory that the dispersed micro-droplet phase is indeed a thermodynamic minimum for a binary gel-solvent mixture. A competition between surface tension and network elasticity stabilizes this phase. When the surfacetension exceeds a critical value, a single macroscopic droplet is the stable thermodynamic phase. Though the Flory-Huggins functional has been used to describe polymer mixtures, our results are generic and valid for any bistable potential [37].

III. DISCUSSION
A macroscopic gel would possess intrinsic heterogeneities in the mesh size resulting in different values of N B in different that leads to a micro-droplet phase with different sized droplets [38]. Thus our mean-field theory needs to be extended to incorporate a distribution of mesh sizes, i.e. P (N B ). Assuming that the disorder correlation length is ψ, our mean-field theory is applicable for length scales ≤ ψ from which the coexistence densities φ in and φ out can be obtained. The coexistence densities are functions of N B , a parameter in the FH free-energy. Thus a distribution of mesh sizes, P (N B ), leads to a distribution of coexistence volume fractions (akin to "mosaic state" in spin-glass models [39]) within the sample, which can be computed using the formalism presented here. The differing coexistence volume fractions in different parts of the sample corresponding to different local values of N B would result in additional surface energy cost between domains that has not been considered in the present calculation. However, this would have effect on the thermodynamics of the mixture gel-solvent mixture [13]. Upon investigating the slope of the common-tangent for the bulk free-energy, F b (φ), for different values of N B , we infer that at a constant temperature (and hence constant χ(T )) the effect of increasing N B leads to the lowering of the slope of the common tangent. Thus, a heterogeneous mesh-size would result in random slopes of the common tangent to F b (φ). The situation is analogous to the behaviour of random-field Ising models, where the relative depth of the bistable free-energy is set by the value of the field h(x) [40].
The effect of network disorder and its relation to the thermodynamics of random field Ising models would be studied in a future work. Elastically mediated phase transitions admit a third thermodynamic phase, where the gel network partially wets and intrudes the solvent rich droplets [31]. A variational calculation allowing for polydisperse droplets and their associated wetting behaviour is currently underway and will be reported elsewhere.
We place our work in context of previous work in this exciting area. The importance of elastic interactions in modifying the equilibrium state of phase separating system was first discussed in context of a ternary system with the elastic network and a polymer interacting with a solvent [14,15]. The stability of a droplet phase is argued along the lines of classical nucleation theory, using the Gibbs free energy (Eq. 1 of [15]) to relate the work done by an expanding drop against the pressure exerted by the bounding polymer network. The droplet is identified as a dilute solvent and the ideal gas equation is used to determine the chemical potential difference ∆µ = k B T ln φ φsat . Based on this formalism (and the Eqs. (1)-(9) of the SI) the authors argue that when φ < φ sat , the mixture is stable, independent of elasticity. When φ > φ cond , the mixture is unstable. While in the interim region φ sat < φ < φ cond a microdroplet phase is stabilised. This experimental situation is closely modelled by Kothari et al. [17] who focus on the kinetics of a three-component system written in terms of the volume fractions of liquids A (uncrosslinked part of the background gel) φ A , part of liquid B φ B that resides within the gel, and φ D , the part of liquid B which exists in droplet form. Our model bears resemblance with the model free energy proposed by Wei et al. [16], though differing significantly in detail. Perhaps the work that is most relevant to the present study is the beautiful scaling theory backed by simulation data by Ronceray et al. [31]. We believe that our work is the first calculation against which these results can be compared. In fact, the schematic phase diagram (Fig 2 of [31]) can be derived from the thermodynamic treatment presented in the present manuscript. In addition, deviations from the scaling theory can also be captured within our model. We hope that our work will prompt careful experimental and theoretical studies in this area. Lastly, we note that our thermodynamic formalism does not capture the exciting non-equilibrium effects [38]. A time-dependent Ginzburg-Landau formalism based on the free energy form explored in this article that incorporates network inhomogeneity, and adhesion of droplets to gel matrices will be explored in a future study. We hope that our theoretical work will instigate experimental work on binary gel-polymer mixtures towards a complete understanding of this fascinating problem.

Author Contributions
B.M., and B.C. designed the research. B.M. and S.B. contributed equally to this work. B.C. obtained funding for the research. All authors contributed to the paper.

Conflicts of interest
The authors declare that no competing interests exist.

Supplementary Information Thermodynamics of droplets undergoing liquid-liquid phase separation
The planar interface: Consider the geometry of the plane interface in the inset (a) of Figure 1 of the main manuscript. The Helmholtz free-energy per unit volume of this system can be written in the following form, where F b (φ) is the free-energy per unit volume of the bulk, f is the faction of the solvent phase, γ is the surface tension, V is the box volume, and φ 0 is the initial composition and the Lagrange multiplier λ ensures mass conservation. Since we treat the total volume V as a parameter there are four unknowns, φ in , φ out , f and λ in the above free-energy. Minimisation w.r.t these four unknowns leads to the following equations, These three equations can be solved to yield the three unknown variables φ in , φ out , and λ. It should be noted that upon rearranging the three above equations, one arrives at the familiar common-tangent conditions : µ(φ in ) = µ(φ out ) and Π(φ in ) = Π(φ out ), where µ(φ) = F b (φ) is the chemical potential and the osmotic pressure is similarly given by Π(φ) = φF b (φ) − F b (φ). These two conditions ensure chemical and mechanical equilibrium, respectively.Thermal equilibrium is ensured as the Helmholtz free-energy is defined in a constant temperature ensemble. Once we know these, the solvent fraction can be found out from the fourth equation ∂F b (φ)/∂λ = 0, which yields, f = φ0−φout φin−φout . Note that these four equations are decoupled as first three equations do not involve the solvent fraction f . The situation is different for a spherical interface and that introduces a non-trivial coupling which we discuss in detail. Once these unknowns are determined, we are free to take the thermodynamic limit, which ensures that the effect of the interface term vanishes as V → ∞. The equilibrium configuration is characterised by two coexisting phases with a planar interface as shown in the inset of Fig. (1) of the main manuscript. The interfacial tension between the coexisting phases has the form, /∂φ| φ1 is the freeenergy after subtracting the common tangent, and k(φ) is the energetic cost associated with spatial variations of order parameter φ [41]. k(φ) has dimensions of ∼ a 2 , where a is the microscopic Kuhn length. The spherical interface: For a system with spherical interface (see inset of (b) of Figure 1 of the main manuscript), the Helmholtz free-energy per unit volume of the droplet phase has the following form : Minimising with respect to the four unknowns result in the following equations, The first two equations imply the equality of chemical potentials : µ(φ in ) = µ(φ out ) and upon substituting the value of λ from the third equation into the first and second one arrives at the second condition : Π(φ in ) = Π(φ out ) + 2γ( 4πN 3f V ) 1/3 . Now by substituting f from the last equation one ends up in the two equilibrium conditions expressed in terms of the coexistence volume fractions, φ in and φ out , and they are solved numerically to yield the coexistence densities for a given box volume V and surface tension γ.

The Choice of the Bulk Free-Energy in Presence of Elastic Interactions
The choice of the exact functional form of the bulk free-energy, F b (φ), is dictated by the form of the equation arising from the condition of equilibrium of the osmotic pressure ( derived below) in the solvent rich and the solvent depleted phases : We will describe how we arrived at Eq. (13) in the text below, however first let us relate the bulk freeenergy to the osmotic pressure, Π b (φ). The osmotic pressure is related to the free-energy via the expression, It is clearly evident that in a F b (φ) with a single minima the tangent at φin lies above the tangent at φout thus leading to Π b (φin) < Π b (φout) (unstable solutions). Stable solutions with Π b (φin) > Π b (φout) can only be found for functions with two minima (lower panel).
. This equation can be rearranged to the form : y(φ) = φF b (φ 0 ) − Π b (φ 0 ). Upon substituting φ = 0, one obtains the intercept to the vertical axis occurs at (0, −Π b (φ 0 )). Note the negative sign as it has an important role to play in the subsequent discussion. The equation which we are solving to determine the equilibrium volume fractions is Eq. (13). In addition to this, the equality of the exchange chemical potentials implies that the tangent lines at the coexistence volume fractions are parallel. This osmotic pressure equation implies that the equilibrium coexistence volume fractions should be such that Π b (φ in ) > Π b (φ out ). This is evident as the second term on the right hand side of the above equation is a positive quantity as it is equal to 2γ/R and similarly we have verified that the last term, , it implies that tangent line at φ in must lie below the tangent line at φ out . In order to demonstrate this refer to Figure (6), where a pair of parallel tangents are drawn at two coexistence volume fractions φ in and φ out .
For the form of the bulk free-energy and the parameter values used in this Figure (6) and also in our calculations presented in this manuscript, the tangents intercepts the vertical axis at negative values, which this implies that the sign of both Π b (φ in ) and Π b (φ out ) are positive. In the upper panel we have a free-energy where entropy term associated the gel, 1 set to zero by putting N B explicitly equal to ∞. Note that since F b (φ) has a single minimum, in the upper panel of Figure (6), the parallel tangents have been constructed at a stable φ out (F b (φ) > 0) and an unstable φ in (F b (φ) < 0). Here the tangent at φ in is located above the tangent at φ out and thus would leading to an unstable solution where Π b (φ in ) < Π b (φ out ). This thus implies that no stable solutions can be obtained for the phaseseparated configurations and thus only the mixed state with uniform order parameter φ 0 would be stable. On the other hand, if one has a high but finite N B (equal to 25) is shown in panel (b), the tangent at φ in lies below the tangent at φ out and thus it leads to stable configurations where Π b (φ in ) > Π b (φ out ). Note that for bulk free-energies with finite N b , at both φ in and φ out , F b (φ) > 0. This thus proves that stable solutions for the phase-separated configurations (as admitted by Eq. (13) of the main manuscript) can only appear if the bulk freeenergy admits two minima, which, in turn, can only occur if the entropy associated with the gel is small, but finite, which is brought about by a high but finite N B . The way out of this situation is to consider forms of F b (φ) where the translational entropy of the gel is finite but small. Physically, when one considers a gelling mixture, which has been thermally quenched to prepare the gel, the thermal disorder would result in parts of the sample which would have strong gelation resulting in high values of N B , which is the polymerisation index of the gel strands. Coexisting with these regions, would be those where the local value of N B is smaller. Thus we consider a form of the bulk free-energy, F b (φ), where the translational entropy of the gel has been divided by a high value of N B . We have performed our calculations for various values of N B and we find that the basic result, which is the stabilisation of the dispersed micro-droplet phase arising via a competition between surface and elastic energies remain valid for computations performed for all values of N B . Thus the free-energy, describing the bulk gel-solvent mixture, that has been considered is given by, Upon analysing the above form of F b (φ) we find that depending on values of χ(T ) and N B , this function can admit both two minima (where stable solutions, Π bulk (φ in ) > Π bulk (φ out ), are possible to find and would thus stabilise the dispersed micro-droplet phase) and one minima and one maxima free-energy landscapes (the mixed phase would be the only stable phase for these parameter values). The left panel of Figure ( closer to unity moves even closer when N B is increased. This makes computations at a finite precision difficult and thus to avoid this we have performed computations at values of N B ranging between 12 to 50 and present results only for N B equal to 25. Upon increasing the temperature (decreasing χ(T )), one loses the unstable region, where F b (φ) < 0 and the free-energy thus becomes one with a single minima. Thus we have performed our computations at χ(T ) = 1.38χ c (see the point marked by a dot in the right panel of Figure (7)), where spontaneous phase separation is possible.
The elastic free-energy : The total elastic free energy associated with accommodating a single solvent droplet of radius, R, inside the gel of mesh size, R 0 is given by [16], To incorporate the effects of the finite stretchability of the gel, we adopt Gent model [34,35], where the elastic free energy density is of the form, W (λ) = − GJm 2 ln 1 − J Jm , where J = λ 2 r + λ 2 θ + λ 2 φ − 3, J m ∼ 10 6 is an upper-limit of the stretching and G is the shear modulus of the network gel. The shear modulus is related to the microscopic parameters via the relation, , where n dry is the cross-link density of the dry-gel and R 0 is the mesh size of the dry gel [36]. The mesh size R 0 is given by where N m is the number of monomers along the backbone of the dry gel, between two cross-links (N m = 16 in our subsequent calculations). The parameter b is the linear dimension of an effective polymeric monomer and following Tanaka [36] we take b ∼ 5a, where a is the Kuhn length or the smallest length-scale associated with the polymer. Due to the volume-preserving nature of the deformation, one has λ r = 1/λ 2 θ and λ φ = λ θ = λ [34] and the deformation obeys the following bound : 0 < J/J m < 1.
The upper limit of the above integral signifies the droplet-gel interface and the lower limit, of unity, signifies a region far away from the centre of the droplet where the stress fields have decayed and the gel there is completely unstressed. As the droplet radius R is related to the solvent fraction, f , via the relation f = (N/V ) 4 3 πR 3 , the elastic free-energy per unit volume is thus expressed in terms of the solvent volume fraction, f . In those situations when the upper limit of the integral, R/R 0 , is less than unity, it means that the droplets do not deform the elastic network and thus F el (f ) = 0. Thus the final form for the elastic energy per unit volume is, where the normalisation by the volume of the gel inside the box is evident.
The primary set of parameters on which the thermodynamic phase of the system depends are : the surface energy γ, the shear modulus of the gel, G, which again depends on the mesh-size of the gel, R 0 . The presence or absence of the dispersed micro-droplet phase depends on the relative weights of the elastic and the surface-energy interactions [31]. In the limit where one has a single macroscopic solvent droplet of radius R → ∞ inside the gel the associated elastic energy per unit volume can be written in the form, (17) In the limit of large droplet radius, R, the elastic energy per unit volume can be cast in the form, f el (R) = αG, where α ∼ 2.5 is a dimensionless constant. In the limit where there are micro-droplets of solvent dispersed inside the gel, the surface tension becomes important. The surface energy per unit volume of the droplets is given by f surf (R 0 ) = 3γ R0 . The ratio of the surface and the elastic energies per unit volume is given by, where, the dimensionless elasto-capillary number is given h = 3γ R0G . If h < α, then the thermodynamic stable state is that of dispersed micro-droplets in the gel, while if h > α, the stable phase is one with a single macroscopic droplet. In the subsequent calculations we choose a value of γ such that h ∼ α and we go on to see whether we indeed observe a dispersed droplet phase in a more detailed microscopic mean-field theory calculations. The value of γ is set to 1/600, unless in the set calculations where the characteristics of the droplet phase is investigated by varying γ while keeping G fixed.
Thus the total free-energy of the dispersed droplet phase in the background gel-matrix has the following form, when every term is expressed as a function of the solvent volume fraction, f , Upon minimising the above free-energy w.r.t the four unknowns one has the following equations, Again, by substituting the expression for λ into the first two equations and by identifying that f = φ0−φout φin−φout , the above equations can be recast into two equilibrium conditions which impose chemical and mechanical equilibrium, respectively, The two equations, Eq. (22), have been solved numerically, via the parallel tangent construction, for the two coexistence densities provided one inputs the values of the surface tension and the box volume and the composition φ 0 . The numerical solution of these two equations proceeds in the following manner : we ensure that the chosen pair of points φ in and φ out , have local tangents those are parallel (µ(φ in ) = µ(φ out )). Since the value of the composition, φ 0 , is an input, the value of the solvent fraction, f , is readily computed via f = (φ0−φout) (φin−φout) . As a result, the difference in osmotic pressure, Π b (φ in ) − Π b (φ out ), becomes a function of the solvent fraction, f , and the surface tension and elastic terms are intrinsic functions of f . One then plots the three terms in the appearing in the mechanical equilibrium condition as a function of the solvent volume fraction, f , and search for the intersection of these curves (see Fig. (8)). The chosen value of φ 0 is 0.45, which lies between the two local minima of the free-energy, F b (φ) and the shape of Π b (φ in ) − Π b (φ out ) can be understood in the following way : low value of f implies φ out is close to φ 0 and thus the slopes of the parallel tangents are maximum. This implies that the vertical distance between the tangents is large and this translates to a large value of Π b (φ in ) − Π b (φ out ) in the plateau region for low f . At high values of f , the tangents are closer to the common tangent of F b (φ) and this implies smaller vertical separation between them and thus translates to a small value of Π b (φ in ) − Π b (φ out ). The shape of this function, Π b (φ in ) − Π b (φ out ), depends on the temperature : at high temeperatures (low χ(T )), the height of the plateau is smaller and the knee region is more pronounced. At lower temperature (high χ(T )), the plateau occurs at a higher value and it drops almost vertically at large f . The sum of the surface and elastic contributions, Π surf (N )+Π el (N ), has two contributions, with the surface contribution, 2γ R(N ) , dominating at low f and the elastic contribution dominating at larger f . This combination also depends significantly on the number of droplets, N , and as a result value of f at the intersection and consequently φ in and φ out becomes a function of N . As a result, the procedure of determining f , φ in and φ out is then repeated for all number of droplets.
After each of these four unknown variables, φ in ,φ out ,f , and λ have been determined from the computation associated with each droplet number, N , the equilibrium values are substituted back into the original free-energy expression (see Eq. (19)), which we set out to minimise, we get a new free-energy,F (N ), which is a function of the number of droplets N . We explore the properties of F (N ) to investigate its shape and whether it admits a single or a multiple droplet minimum.