Tuning the phase diagram of colloid-polymer mixtures via Yukawa interactions

Theory that predicts the phase behavior of interacting Yukawa spheres in a solution containing nonadsorbing polymer is presented. Our approach accounts for multiple overlap of depletion zones. It is found that additional Yukawa interactions beyond hard core interactions strongly affect the location and presence of coexistence regions and phase states. The theoretical phase diagrams are compared with Monte Carlo simulations. The agreement between the two approaches supports the validity of the theoretical approximations made and conﬁrms that, by choosing the parameters of the interaction potentials, tuning of the binodals is possible. The critical end point characterizes the phase diagram topology. It is demonstrated how an additional Yukawa interaction shifts this point with respect to the hard sphere case. Provided a certain depletant-to-colloid size ratio for which a stable colloidal gas-liquid phase coexistence takes place for hard spheres, added direct interactions turn this into a metastable gas-liquid equilibrium. The opposite case, the induction of a stable gas-liquid coexistence where only ﬂuid-solid was present for hard spheres, is also reported.


I. INTRODUCTION
A colloidal dispersion can be regarded as a collection of big atoms [1,2]. This colloid-atom analogy was used by Einstein in his theory for Brownian motion of particles suspended in a solvent [3]. The later experimental verification by Perrin [4] of the colloidal barometric height distribution following Boltzmann's law constituted the starting point of systematic studies on the dynamics and equilibrium properties of colloidal systems [5]. The effective interaction between colloidal particles mediates the physical properties of a colloidal dispersion, the hard sphere (HS) model being the simplest and widely used approach [6][7][8].
Perrin's colloid-atom analogy also holds for the phase behavior. In fact, nanoparticle dispersions exhibit the same phase states as those found in atomic systems (i.e., gas, liquid, and solid). Colloidal interactions beyond the pure hard-core interaction influence the topology of the phase diagram [9]. Attraction between the colloidal particles can induce instability of the fluid phase, leading to a dilute fluid (colloidal gas, G) coexisting with a denser fluid (colloidal liquid, L). For a range of attraction greater than 1/3 of the particle diameter, G-L coexistence can be stable [10][11][12][13]. The phase behavior of a colloidal dispersion can be tuned, for instance, by adding nonadsorbing polymers. This leads to an effective attraction between the colloidal particles, as shown by Asakura and Oosawa [14]. The strength and range of attraction between nanoparticles can be controlled via the depletant * r.tuinier@tue.nl concentration and the depletant-to-colloid size ratio (q D ). The depletion interaction is used to rationalize both fundamental as well as practical issues [15][16][17][18][19].
Two different regimes for depletion interaction can be distinguished. For q D 0.15 a collection of spheres mediated by depletants is pair-wise additive [20]. Otherwise, multiple overlap of depletion zones needs to be accounted for. The thermodynamical consequences of assuming pair-wise addition in the latter case have been analyzed and discussed in detail [20,21]. Describing the mixture by an effective pairpotential approach has been widely applied [21][22][23][24][25][26][27]. Contrary to the usually complex or computationally based approaches followed, the method reported here is simple, yet accurate, and constitutes a systematic tool for predicting phase diagrams of colloid-polymer mixtures. Moreover, the theory accounts for multiple overlap of depletion zones.
Much attention has been paid to the phase behavior of HS plus depletants [9,20,[28][29][30][31] when multiple overlap of depletion zones occurs (q D 0.15). However, limited theories are available for dispersed colloidal particles with realistic interactions and added free polymers (acting as depletants). Available approaches focused on short-range repulsive colloids [16,[32][33][34][35]. In real systems, more involved (soft) interactions are often present. For example, in paint and food dispersions, particles are often charged and/or partially (de)stabilized by adsorbing polymers [36].
We show in this work how additional direct interactions between the colloidal particles dramatically modify the thermodynamic equilibrium of the colloid-polymer mixture. The direct interactions between the spherical colloidal particles are described via hard-core Yukawa potentials [37][38][39][40] (HCYs). The nonadsorbing polymers are treated as penetrable hard spheres [41] (PHSs). In the results section, phase diagrams of HCY spheres in a sea of PHSs are presented and compared with Monte Carlo (MC) simulations. It is shown how additional direct interactions dramatically modify the topology of the phase diagrams, as well as the location of the phase transition lines within the theoretical framework of free volume theory (FVT) [10].

II. THEORY
The present section summarizes the tools required for the calculation of the thermodynamic properties of colloid-polymer mixtures, which are used to calculate phase diagrams. The main modification with respect to the theory for HS plus depletants [10] is presented to some extent, while further details are shown in the Appendixes. Furthermore, the pair potentials used as input for the MC simulation are also introduced.

A. Pair potentials and second virial coefficient
The interaction potential between colloidal spheres is described in this research via a hard-core plus Yukawa (HCY) [40,42] attraction or repulsion. The Yukawa interaction enables us to mimic a wide range of interactions, since both the range and strength can be tuned: it could represent, for instance, a screened double layer repulsion or a Van der Waals attraction between the colloids [32,[43][44][45]. Since the focus is on the effect of the direct interaction (the additional Yukawa interaction) on the phase behavior of colloid-polymer mixtures, the depletants are treated in the simplest way. Hence, the depletants are described as penetrable hard spheres (PHSs): they do not interact with each other, but have a hard-core repulsive interaction with the colloidal spheres. The PHS approximation can be used to describe depletion effects of dilute polymer solutions. The method developed can, however, be extended towards polymers in θ and good solvent conditions at elevated concentrations [34,[46][47][48][49][50][51].
Due to radial symmetry of all interaction potentials, it is convenient to work with the dimensionless distance between the centers of two colloidal spheres x = r/σ , with r the center to center distance and σ the HS diameter (σ = 2R, with R the HS radius). The depletion thickness is δ, which is identified here with the PHS radius. Furthermore, the relative range of the Yukawa interaction is characterized by q Y = 2/κσ = 1/κR, with screening length κ. When the Yukawa interaction mimics a double layer repulsion, its screening length is related to the salt concentration and dielectric properties of the solvent. The depletant to colloid size ratio q D = δ/R defines the relative range of the depletion interaction.
The HCY pair potential between the colloidal spheres is written in terms of x, q Y , and , and normalized by the factor β = 1/(k B T ) (with k B the Boltzmann's constant and T the temperature): In order to calculate the phase equilibria it is convenient to consider the colloid-polymer mixture to be in equilibrium with an external reservoir that does not contain colloidal particles. In this reservoir, as explained in the next subsection, the PHS concentration (φ R d ) fixes the fugacity of the depletants. In this way the PHSs lead to an effective depletion attraction, expressed in Eq. (2) in terms of q D and φ R d [41]: The strength of the Yukawa potential is set by the (dimensionless) contact potential , defined such that > 0 implies a HCY attraction and < 0 a HCY repulsion. For = 0 the HCY reduces to the HS interaction. Similarly, the strength of the depletion potential depends on the relative polymer concentration in the external reservoir, φ R d . When φ R d = 0 the depletion pair potential reduces also to the HS one. Some insight into the total interaction potential can be obtained by taking the sum of the HCY and depletion contributions (superposition approximation): . The quantity U T provides insight into the governing forces between the colloidal spheres that determine the dispersion stability. Figure 1 presents the depletion, HCY, and total potentials for two different sets of parameters with the same |β |. Note how the (local) attractive or repulsive nature of U T sensitively depends on the key parameters {q D ,φ R d ,q Y , }. It should be also clarified that this pair potential does not account for the multibody nature of the depletion interaction [21,28,52,53].
Depending on the particular values for the additional Yukawa potential considered, the effective excluded volume between colloidal particles increases (repulsions) or decreases (attractions) with respect to the hard sphere case. Such effect is accounted for in the second osmotic virial coefficient (B 2 ). Its general expression in spherical coordinates can be simplified for potentials that are hard core within the colloidal sphere's volume using the dimensionless distance x: where the colloidal sphere volume v c has been introduced (v c = πσ 3 /6). This leads to higher (repulsions) or lower (attractions) B 2 with respect to the hard sphere case (B HS 2 = 4v c ). As for the pair potentials, the excluded volume depends on the total set {q D ,φ R d ,q Y , } (as can be seen in Fig. 2). It is further convenient to define the effective excluded volume per colloidal particle volume as B * 2 = B 2 /v c . For attractive interactions, the Vliegenthart-Lekkerkerker (VL) criterion [12] states that at the colloidal gas-liquid critical point B 2 has a unique value: However, it has been recently reported that the VL criterion does not hold when considering the gas-liquid (G-L) critical point (CP) for spherical particles that attract one another via the depletion interaction [53]. This effect is re-evaluated in here in view of the results presented for the phase behavior. Furthermore, B 2 is also studied for repulsive HCY spheres in a sea of polymeric depletants.

B. Basic elements of free volume theory (FVT)
Free volume theory (FVT) is used here for the calculation of the phase diagrams of mixtures containing colloidal particles and depletants in a background solvent. Within this theory, a system (S) is considered at given temperature and held in osmotic equilibrium with a reservoir (R). System and reservoir are separated by a hypothetical membrane permeable to the depletants and the solvent but impermeable to the colloidal particles. The system of interest contains colloidal particles, depletants, and solvent, while the reservoir only contains solvent and depletants. This establishes the semigrand potential (N c ,V ,T ,μ d ) that depends on the temperature (T ), the depletant chemical potential (μ d ) and the number of colloidal particles (N c ) in a volume V . The grand potential can be split into a canonical contribution plus a contribution arising from the presence of depletants, which is related to φ R d : where F (N c ,V ,T ) is the Helmholtz free energy of the suspension of interacting colloidal particles in S in the absence of depletants [so, formally, As in R only depletants are present, φ R d fixes the depletant activity (hence the state). The chemical potential of the depletants in the reservoir and in the system are equal in equilibrium. The colloid-depletant mixture (S) is characterized by the colloid packing fraction (η) and the depletant concentration in the system (φ S d ). The state of S depends on the particular interaction between the colloids, the colloidal packing fraction and the depletant concentration.
As established in original FVT, (i) the depletants are ideal (no interactions between them), and (ii) the available free volume in S for the depletants does not depend on φ S d (the depletant concentration in the system).
This allows rewriting the semigrand potential in dimensionless units: where α is the free volume fraction for the depletants in S, R is the osmotic pressure of depletants in R. The tilde over the magnitude indicates dimensionless units. The chemical potential of the colloidal particles ( μ c ) and the total osmotic pressure of colloidal particles plus depletants in the system ( T ) follow (in dimensionless units) as Note here that μ c relates to the chemical potential of different colloidal phases in S, as colloidal particles do not trespass the fictional membrane. Since the depletants behave ideally, the osmotic pressure in R ( R ) follows van't Hoff's law: The depletant concentration in the system (φ S d ) is defined in terms of the free volume fraction (α) and the concentration in the reservoir as This free volume fraction and the HCY contributions to FVT are quantified in Appendixes A and B respectively.

C. Calculation binodals and critical end point
Provided that all contributions to the grand-canonical potential are known (see Appendixes A and B), calculation of phase diagrams of HCY interacting colloids plus PHSs using FVT is possible. A phase diagram may contain various characteristic quantities defined below. The limiting value of the colloidal gas-liquid (G-L) coexistence is characterized by the critical point (CP), which satisfies where μ and T are the chemical potential and the osmotic pressure of the colloidal suspension. The triple point (TP) at which colloidal gas and liquid phases coexist with a solid phase satisfies the conditions Phase coexistence binodals are defined by the family of points that satisfy with {i,j } the colloidal solid (S), liquid (L), gas (G), or fluid (F) phase. For sufficiently short-range attractions, G-L coexistence is metastable [10,11]. The critical end point (CEP) marks the critical point that coincides with the triple point [13,50], defined by

III. RESULTS AND DISCUSSION
This section focuses on the phase diagrams obtained for mixtures of HCY spheres and PHSs. First, the binodals for colloidal HS in a sea of PHSs are recalled, as a point of reference for later phase diagrams. Subsequently, the effect of added Yukawa interactions on phase diagrams is discussed. Phase coexistence is considered for variable ranges and strengths of the Yukawa interaction and depletant to colloid size ratios. Further, the CEP that characterizes the topology of the phase diagram and its dependence on the additional Yukawa interaction are reported. The value of the second osmotic virial coefficient at the CP is discussed for both attractive and repulsive HCY spheres plus depletants. Finally, the phase diagrams predicted using FVT are compared with MC simulation results. in close agreement with computer simulations [54]. Colloidal G-L phase coexistence is obtained for φ R d values above the critical point and is characterized by a characteristic U shape. For short-range depletion attractions (q D 0.33) only colloidal fluid-solid (F-S) coexistence is observed. Then the G-L coexistence is metastable. For long-range depletion attractions (q D 0.33), colloidal G-L coexistence can take place. The transition between these two regimes is defined by the critical end point (CEP) [13,50], as discussed later in this paper. Two binodal curves are compared with simulation results as an illustration of the accuracy of FVT [28]. These phase diagrams constitute the reference for the ones in the next sections, where added Yukawa interactions are accounted for.

B. Yukawa contributions to FVT
Before presenting phase diagrams of mixtures of HCY interacting spheres in a sea of PHSs, the influence of the added Yukawa interaction on FVT is discussed. As shown in Fig. 4, added Yukawa attractions lower the osmotic pressure of the colloidal dispersions while repulsions increase with respect to the pure HS interaction. This can be explained by the fact that an additional weak repulsion increases the effective excluded volume of the colloidal particle [32]. An inverse effect is expected for weak attractions. This has a significant influence on the stability of G-L and F-S coexistence regions, inducing different topologies with respect to the HS of the phase space at given set of interaction (q D ,q Y , ).
In Fig. 5 the free volume fraction α for depletants in a suspension of HCY spheres is plotted for different colloidpolymer size ratios (q D ). Results are shown for HS and for both attractive and repulsive HCY spheres. The difference in α between HS and HCY is small, but increases with q D and with increasing strength (| |) of the Yukawa interactions. Further effects are quantified in Appendix B.

C. Complete phase diagrams and critical points
In this section, phase diagrams are presented for colloidal spheres with HCY interaction plus PHS. As expected from the pair potentials (see Fig. 1), the key parameters that determine phase coexistence are the strength (β ) and range (q Y ) of the attractive or repulsive HCY interaction, the depletant concentration (φ R d ), and the relative size of the depletant (q D ). In Fig. 6 phase diagrams calculated for a long-range Yukawa interaction (q Y = 2) combined with short-range depletion attraction (q D = 0.15) are presented.  attractions are strong enough (β > 0.25). Such G-L coexistence leads to a narrow liquid window for (β = 0.5), and for higher contact potentials the whole phase space gets unstable (demixing). While HCY repulsions hardly affect the F-S phase in {0.49 η 0.55}, additional HCY attractions broaden this F-S coexistence region. Although obvious (but often ignored), it is concluded here that added direct repulsions between colloidal spheres stabilize (widen) the stability regions of colloid-polymer mixtures while direct attractions decrease the stable regions.
Phase diagrams obtained as a result of long-range depletion (q D = 1) and short-range HCY interactions (q Y = 0.15) are presented in Fig. 7. Independently of the strength and nature of the HCY interaction, the characteristic U shape of G-L coexistence is always present for the combination {q D ,q Y } = {1,0.15}. There is a remarkable increase of the width between the two vertical curves corresponding to the F-S coexisting phase at low φ R d when direct attractions are incorporated (i.e., the HS F-S coexistence broadens due to added direct attractions). Even though the Yukawa attraction is short-ranged, the binodals are shifted tremendously. The main trends (an increased stability region for more repulsive interactions and decreased stability for attractions) of the shifts in the phase lines are similar to those in Fig. 6. There is always G-L coexistence at η < 0.4 for these sets of {q Y ,q D , } values; see the inset of Fig. 7.
These two limiting situations indicate how additional Yukawa interactions modify the colloid plus depletant phase diagram with respect to the HS case. Based on the results presented above, the following scenarios are expected: colloidal G-L coexistence is observed at relatively low |β | when it is present in the hard sphere case (q D 0.33). On the other hand, for q D 0.33 only F-S coexisting phases are expected.
Next, different combinations of Yukawa and depletion interactions are explored by systematically varying the strength of the Yukawa interaction for two different sets of ranges of the depletion and Yukawa interaction (Fig. 8). The most striking observation is that the direct Yukawa interactions change the topology of the phase diagrams with respect to the HS case, inducing a stable G-L coexistence when only F-S is present for HS [ Fig. 8(b)] and vice versa [ Fig. 8(a)]. This topological change with respect to the HS case depends on the particular HCY interaction. In the bottom panels of Fig. 8, phase diagrams are presented using the depletant concentration in the colloid-polymer mixture (φ S d ). As can be appreciated, the topology in this phase space is significantly more sensible when direct attractions are added as compared to the case that direct repulsions between the colloidal particles are added. This is in line with the experimental observation that added repulsions between colloidal particles stabilize the suspension.
The following general trends can be withdrawn from computed phase diagrams plotted in Figs. 6-8. Additional Yukawa attractions between colloidal particles lower the depletant concentration required for phase coexistence. This applies independently of the nature of the coexistence phase (G-L or F-S) and the ranges of the interactions considered. Equivalently, added repulsions increase the required depletant concentration for phase coexistence. The latter is in qualitative agreement with computer simulations and experimental results for charged nanoparticles in polymeric solutions [16,32,34,35]. Not surprising, the difference in the phase diagram with respect to the HS cases is more dramatic for larger values of |β | and q Y . Sufficiently strong, short-ranged Yukawa attractions turn the G-L coexistence region eventually metastable, inducing a thermodynamically stable F-S coexistence even at low colloid volume fractions. This shows that it is possible to modify the phase diagram topology of colloid-polymer mixtures by tuning the Yukawa pair potential. Such tuning of the direct colloid interactions in experimental systems can be achieved for example by modifying the double layer repulsion [55,56] on brush repulsion [57,58] or in near-critical solvent mixtures [59].
The width of the colloidal G-L region also varies with the set {q D ,q Y ,β }. The onset of the G-L region is marked by the CP. It is noted that the G-L critical point can be metastable as it might lay inside the more stable F-S coexistence phase. At a fixed contact potential ( ), its position varies depending on the range of the Yukawa interaction. In Fig. 9 the position of the CP is plotted at a given |β | for a collection of q Y values. As expected, increasing the range of repulsive interactions increases the depletant concentration needed in order to obtain a stable G-L coexisting phase. When extra attractions between the colloidal particles are considered, φ R d (CP) is lowered with respect to the HS case, as expected. This produces either a G-L or a F-S phase coexistence, depending on the phases present in the = 0 case and the strength of the contact potential. For strong and long-ranged enough attractions, the coexistence phase extends practically over the whole phase space. All curves collapse at φ R d = 0, providing the q Y value at which a suspension of purely attractive-HCY spheres gets critical at = 0.5 ({η,q Y } = {0.15,1.66}). Note the inflexion point at β = 0 (HS case), which is particularly visible for small q D .

D. Critical endpoint
The main indicator of the topology of the phase diagram is the critical end point (CEP). It has been demonstrated that the CEP is correlated with an effective range of attraction close to 1/3 of the particle diameter [11]. The CEP has been determined both for HS in a sea of PHSs [50] and attractive HCY fluids [13]. For a suspension of attractive HCY spheres [13], the CEP is q In Fig. 10 binodals are presented for the set {q Y ,q D } = {0.26,0.33}, the corresponding CEP of purely attractive HCY fluids, and HS in a sea of PHSs. The difference for HCY spheres in a sea of PHSs with respect to the HS case is small and follows the trends previously described. All binodals in Fig. 10 exhibit a shape of the fluid branch that is critical (or nearly critical), and a narrow G-L coexistence arises for repulsive Yukawa interactions above |β | = 0.5: due to combined short-ranged indirect attractions and short-ranged repulsions, G-L coexistence takes place, while it was metastable for HS plus PHSs.
The calculated CEP values are summarized in Fig. 11, and further details are given in Appendix C.
For short-range Yukawa attractions (q Y < 0.26), q cep D increases with |β | with respect to the HS case. Due to a short-range direct HCY attraction that destabilizes colloidal G-L coexistence, the range of the depletion interaction needs to be increased in order to still achieve a stable G-L coexistence. In contrast, if the range of the HCY attractions is long enough (q Y > 0. 26

E. Second virial coefficient at the critical point
It is worthwhile to analyze the value of the second virial coefficient at the critical point, given the tools for systematically calculating the CP from FVT. In Fig. 12 the normalized second osmotic virial coefficient (B * 2 = B 2 /v c ) is plotted as a function of q D for various q Y values at a fixed contact potential (|β | = 0.5). Especially for sufficiently long-ranged repulsive potentials B * 2 exceeds the HS (B * ,o 2 = 4) value for small q D . As expected from previous research [53] the Vliegenthart-Lekkerkerker (VL) criterion [B * 2 (CP) VL = −6] does not hold for a mixture of colloidal spheres and PHSs even when the added direct attractions are sufficiently long ranged. This is due to the indirect nature of the depletion interaction: the VL criterion is based upon direct attractions. As expected, indirect attractions lead to different physical behavior. As a consequence, when q D is high enough, B * 2 at the critical points can be smaller than specified by the VL criterion even for additional strong, long-ranged direct attractions (as can be seen for q Y = 1.5 in Fig. 12). Not surprisingly, low values of q Y exhibit trends closer to the HS case than high ones.   (NVT). The details of the simulations performed and the calculation of the free energies are presented in Appendix D.
Phase diagrams obtained using the free energies obtained from MC simulations for HCY spheres plus the considered depletion potential are presented as the data points in Fig. 13. The trends observed match the ones of FVT for HCY spheres plus PHSs (solid curves). The F-S coexistence arising from the HS case is broader at high colloidal packing fraction than predicted by FVT at sufficiently high {q D ,φ R d }. This (most likely) reflects to the lack of accounting for multioverlap of the depletion zones in the MC simulations, which is accounted for within FVT. In all cases studied, the MC-predicted phase diagrams quantitatively match with the FVT results when  following the perturbation approach. Moreover, the results for HS + PHSs differ from the original calculations [20], where no G-L coexistence was found for q D = 0.4. Experimentally, G-L coexistence has been reported for q D > 0.3 [50,60]. Future improvements of the simulation method are possible, mainly accounting for the multibody nature of the depletion interaction [28], without presuming the FMSA free energy (i.e., performing two λ integrations for the free energy) or using binary mixtures of hard-spheres and penetrable hard spheres in a canonical or grand-canonical ensemble.

IV. CONCLUSIONS
The quantitative match between the Monte Carlo-generated phase diagrams and those arising from free volume theory implies that the tuning knobs that determine the phase diagram of hard-core Yukawa (HCY) spheres in a sea of penetrable hard spheres (PHSs) have been identified correctly within our theoretical framework. The depletant concentration and relative depletant size plus the range and strength of the direct interactions allow controlling the stable phase regions of colloid-polymer mixtures. Deviation from the hard sphere (HS) case takes place as the direct HCY interactions become stronger, inducing different phase coexistence regions than those present for a pure suspension of HS in a sea of polymeric depletants. Here we have identified how the critical end point (CEP) is shifted due to combined interactions, leading to values far from the pure HCY or depletion interactions. The CEP determines the types of coexisting phases present given a particular set interaction potentials.
In summary, the results shown in this paper pave the way towards a better understanding of complex colloidpolymer mixtures, providing the necessary tools required when considering soft interparticle interactions.

ACKNOWLEDGMENTS
We thank M. S. Feenstra for her efforts in preliminary work. We acknowledge the Netherlands Organization for Scientific Research (NWO TA Grant No. 731.015.025) for financial support. In additional, we acknowledge DSM and SymoChem for support. Professor M. Dijkstra is thanked for fruitful discussions in the preparation of the MC simulations. We are indebted to Dr. H. Friedrich for the suggestions in the post-processing of these simulations. Dr. M. Vis is thanked for useful discussions, and Professor G. de With is thanked for a critical reading of an earlier version of the manuscript.

APPENDIX A: SCALED PARTICLE THEORY
Following Widom's insertion method [61], the free volume fraction for depletants in S can be related to W , the work required for bringing a depletant into the colloidal suspension: where V free is the ensemble-averaged free volume for depletants in S. Using the approximations made in original FVT, V free is not affected by the presence of the depletants and thus V free can replaced by V free 0 , the undistorted free volume. The work of insertion (W ) is calculated by scaling the size of the depletant inserted from 0 to its final size by a tuning parameter λ (δ ↔ λδ). The actual work of insertion (λ = 1) is calculated by connecting the λ → 0 and λ → ∞ limits via expansion of W in a series of λ. This is the core idea of scaled particle theory (SPT) [62][63][64][65].
In the limit λ 1, depletion zones do not overlap and the volume of depletants is negligible so α only depends on geometrical factors, hence not on the interaction potential between the colloidal particles. This results in the following expressions for α and W : In the opposite limit (λ 1), W can be approximated as the work required for creating a cavity of volume (4π/3)(λδ) 3 against the osmotic pressure of colloidal spheres S c : Within original FVT, the osmotic pressure of the colloidal suspension results from pure HS contributions. Here additional Yukawa interactions are accounted for. As a consequence, the total osmotic pressure of the colloidal particles is given as where 0 c and Y c are the HS and additional Yukawa contributions.
The connection between the two limiting cases given in Eqs. (A2) and (A3) provides Taking λ = 1, multiplying with β, and using q D = δ/R gives As = βv c , and collecting the HS terms: Inserting of Eq. (A1) into Eq. (A4) finally yields: where the function Q(η; q D ) arises from the hard spheres plus depletants contribution [9,10]: The additional Yukawa interactions are accounted for in the last term of Eq.
The quantity Y (η; q Y , ) is expressed in Eq. (B6). SPT provides an essential part of the information required for calculating the phase behavior. As follows from Eq. (6), the free energy of the pure colloidal suspension ( F ) is a key contribution to the grand potential ( ), and it is specified below.

APPENDIX B: FMSA CONTRIBUTION TO FREE VOLUME THEORY
The free energy ( F c ) of the colloidal particles interacting via hard-core Yukawa is described as consisting of a hard core plus an additional Yukawa contribution [13,40]: The pure HS contributions to the free energy ( F HS ) are well known [66,67]. For a fluid of HS, an accurate expression up to η ≈ 0.5 follows from the Carnahan-Starling (CS) [67] equations of state: while for a fcc crystalline solid phase the Lennard-Jones-Devonshire (LJD) [66] equation of state reads where η cp = π/(3 √ 2) ≈ 0.74 is the close packing fraction of HS. The value 2.1306 derives from Monte Carlo simulations of the pure HS system [68].
Tang et al. [37] derived an expression for the free energy of a collection HCY spheres via a first-order mean spherical approximation (FMSA). The HCY potential allows an analytical solution of the Orstein-Zernike integral upon using the mean spherical closure approximation. This leads to analytical expressions for the radial distribution function and the direct correlation function up to first order in inverse temperature. These expressions result in a Yukawa contribution to F beyond F HS ( F = F k HS + F Y , with k denoting the fluid or solid state). Tang's approach can be extended to Multi-Yukawa potentials, and has been successfully applied to study the interactions between charged colloidal particles [69] and Lennard-Jones pair interactions [40]. This Yukawa contribution to the free energy can be written in a van der Waals form [40,50]: in which the functions γ 1 (η; q Y ) and γ 2 (η; q Y ) can be expressed in terms of the auxiliary functions L Y (η; q Y ) and Q Y (η; q Y ): and Using Eqs. (B1) and (B4), the Yukawa contribution to the HS osmotic pressure can be written as FIG. 14. Relative percentage difference between the free volume fraction for PHSs in a HCY sphere suspension compared to the HS case using STP for q Y = 1.0 and |β | = 0.5 for a collection of q D . Dotted curves hold for attractions, and solid ones for repulsions.
For small values of |β | the free energy is symmetric with respect to the HS case ( = 0). For large values of {|β |,q Y } the HCY free energy becomes asymmetric due to the |β | 2 term present in Eq. (B4). The consequences of this symmetry break are discussed below. From the free energy of the pure colloidal suspension, its osmotic pressure and chemical potential can be calculated using standard thermodynamic relations: In Fig. 14 the absolute percentage difference % |α HCY − α HS | is plotted for the same values as in Fig. 5. When comparing the effects of repulsive to attractive Yukawa interactions there is no full symmetry around the abscissa due to the nonsymmetric Yukawa free energy defined by the FMSA [37] [see Eq. (B4)]. Due to the lowering in osmotic pressure caused by additional attractions (see Fig. 4), the work of insertion lowers for bringing a depletant into the suspension and the free volume available for the depletants increases with respect to the HS case. Likewise, for added repulsions the free volume for the depletants is decreased with respect to the HS case. becomes small: in this case the G-L coexistence phase extends over practically the entire phase space (see Fig. 7 for η below the F-S coexistence region). For very high {q Y ,β } values the numerical solution for the CP provides a negative value of φ R d , which induces the large unstable coexistence region observed for instance in Figs. 6 and 7 (dashed curves, attraction).
Second, the effect of HCY repulsions is examined. The intricate combinations of depletion-mediated attraction and direct Yukawa repulsion give rise to an interesting behavior. For small values of |β |, the following scenario is observed. When q Y < 0.26, the q cep D value is lowered with respect to the HS (β = 0) case with decreasing β . For q Y > 0.26, the q cep D value is slightly higher than for the HS scenario.  Upon exceeding a certain value for β (dependent on the particular q Y considered), q cep D is lower again with respect to the HS case for all q Y values. Further insight about this counterintuitive behavior is provided in Fig. 16, where the depletant concentration in the reservoir at the CEP is plotted as a function of β for the same collection of q Y as in Fig. 15. The trends followed for φ R,cep D are straightforward. For β < 0, the required depletant concentration for achieving G-L coexistence increases with q Y and |β |. Opposite to this effect, added attractions β > 0 make G-L coexistence more favorable and hence take place at lower depletant concentrations. Not surprisingly, the topology of the phase diagrams presented (Figs. [3][4][5][6][7][8][9][10]  Further insight of the CEP is withdrawn from the nature of the total pair interaction that allows a stable G-L colloidal coexistence. Figure 17 presents the total potential at the CEP for a given repulsive contact potential and varying range of repulsion [normalized by the contact potential, as U T (x)/U (x = 1)]. The depletant concentration and the depletant to colloid size ratio increase as the range of the interaction increases (inset of Fig. 17). This produces a lowering in the total (attractive) contact potential as repulsions become stronger. However, the local repulsive potential increases while the attractions close to colloidal particle's surface becomes stronger. This indicates a competition between long-range direct repulsion interaction and the indirect depletion attraction which makes possible colloidal G-L coexistence for repulsive colloidal particles.
In view of the computed values for the CEPs, it can be concluded that systems with very strong repulsions (both in range and contact potential) exhibit a G-L coexisting phase if the depletant concentration is high enough and the depletant to colloid size ratio is above approximately 0.33. For added Yukawa attractions, colloidal G-L coexistence becomes metastable at lower colloid-polymer size ratio with increasing direct attractions (both in range and contact potential) when the range of the Yukawa attraction is above the one of the pure attractive Yukawa fluid.

APPENDIX D: DETAILS ON THE MONTE CARLO PHASE DIAGRAM CALCULATION
In a simulation box at constant volume and temperature, a collection of 256 HS interacting via a depletion potential (U D ) plus a Yukawa potential (U HCY ) was employed (NVT ensemble). The pair potentials as specified in Eqs. (1) and (2) were used as input for the MC routine. The simulation approach is based on the method of Dijkstra and co-workers [20,32,70].
A collection of {η,φ R d } state points in the phase diagram is considered for each set of {q D ,q Y ,β } (hence, the total interaction U T is set at each state point). The free energy of the ensemble at each {η,φ R d } is estimated using the λ-integration method for the depletion interaction, with a ten-point Gauss-Legendre integration (for further details of the integration method, see [20,32,70] and references therein). The MC simulations ran for 32 000 lattice sweeps at each λ value, enough for the system to equilibrate in (practically) all state points studied. The considered free energy of the system was averaged over the last 3200 MC cycles at each state point. The collected energy can be associated with the (dimensionless) free energy of the system: HCY contri.
where the angular brackets ( . . . ) indicate ensemble averages. Even though already proven [13], we checked the feasibility of the FMSA. In Fig. 18, the MC-averaged free energies in the absence of depletants (φ R d = 0) are compared with the analytical expressions. The match between the MC-averaged free energies and the analytical FMSA expressions invites us, indeed, to perform the λ integration only over the depletion pair potential.
On this basis, the depletion pair-wise additive contribution to the free energy ( F Dep ) can be estimated by subtracting the analytical FMSA free energy (which, for HS, reduces to CS or LJD equations of state) to the total free energy of the system: FIG. 20. 3D representation of the phase diagram in Fig. 3. This allows treating this depletion contribution as a perturbation from the FMSA free energy, as a simple third-order polynomial fit in η for each φ R d (provided {q D ,q Y ,β }). The chemical potential and osmotic pressure are obtained applying standard thermodynamic to the total free energy obtained as a sum of the FMSA free energy and the fitted depletion free energy: It is noted that this method does not account for multiple overlap of depletion zones [20,28], which becomes especially relevant for q D > 0.4 [21]. Using the same approach done with the analytical FVT phase diagrams but in a canonical ensemble, colloidal G-L and F-S coexistence points are calculated from the simulation data (see Sec. III F and Fig. 19).

APPENDIX E: 3D PLOTS OF THE PHASE DIAGRAMS
In Figs. 20-22 of this short Appendix, an alternative visualization of some of the phase diagrams is presented.