Fluid–fluid coexistence in an athermal colloid–polymer mixture: thermodynamic perturbation theory and continuum molecular-dynamics simulation

Using both theory and continuum simulation, we examine a system comprising a mixture of polymer chains formed from 100 hard-sphere (HS) segments and HS colloids with a diameter which is 20 times that of the polymer segments. According to Wertheim's first-order thermodynamic perturbation theory (TPT1) this athermal system is expected to phase separate into a colloid-rich and a polymer-rich phase. Using a previously developed continuous pseudo-HS potential [J. F. Jover, A. J. Haslam, A. Galindo, G. Jackson, and E. A. Muller, J. Chem. Phys. 137, 144505 (2012)], we simulate the system at a phase point indicated by the theory to be well within the two-phase binodal region. Molecular-dynamics simulations are performed from starting configurations corresponding to completely phase-separated and completely pre-mixed colloids and polymers. Clear evidence is seen of the stabilisation of two coexisting fluid phases in both cases. An analysis of the interfacial tension of the phase-separated regions is made; ultra-low tensions are observed in line with previous values determined with square-gradient theory and experiment for colloid–polymer systems. Further simulations are carried out to examine the nature of these coexisting phases, taking as input the densities and compositions calculated using TPT1 (and corresponding to the peaks in the probability distribution of the density profiles obtained in the simulations). The polymer chains are seen to be fully penetrable by other polymers. By contrast, from the point of view of the colloids, the polymers behave (on average) as almost-impenetrable spheres. It is demonstrated that, while the average interaction between the polymer molecules in the polymer-rich phase is (as expected) soft-repulsive in nature, the corresponding interaction in the colloid-rich phase is of an entirely different form, characterised by a region of effective intermolecular attraction.


Introduction
Colloid-polymer mixtures are present in many substances commonly encountered in everyday life, such as foods and paints, as well as being pervasive in biological systems. They have been found to exhibit a variety of scientifically interesting or practical behaviour [1], and their importance is reflected in the large body of work that has been devoted to their study. Experimental investigations of the phase behaviour of colloid-polymer mixtures have been performed on a variety of systems ranging from comparatively simple latex-polystyrene and silica-polydimethylsiloxane colloidal dispersions to biological systems containing proteins or DNA [2][3][4][5][6][7][8][9][10][11][12][13]. Phase separation has been described by many authors [3,4,8,[14][15][16][17][18][19][20][21][22][23][24][25][26][27]; see the excellent book of Lekkerkerker and Tuinier [1] for a comprehensive historical review. Fluid-solid coexistence is observed when the radius of gyration of the polymer is small compared to that of the colloid; here one can draw an analogy between the behaviour of the colloids and that of a system of pure hard spheres (HSs). On the other hand, by increasing the ratio of the size of the polymers to that of the colloids one can induce a departure from HS-like phase behaviour; when the polymers are sufficiently large, fluid-phase separation between low-density and high-density colloidal phases is observed; at higher densities solid colloidal phases can also be found.
An important concept in understanding the phase behaviour of colloid-polymer mixtures is that of the depletion interaction. Asakura and Oosawa were the first to recognise this idea in two seminal papers [28,29] published over 50 years ago, and the concept was later reprised by Vrij and co-workers [14,30]. Each colloid is imagined to be surrounded by a so-called depletion layer. When two colloids lie in close proximity the depletion layers overlap and the volume available to the polymer chains increases; as a consequence, the free energy of polymer chains is minimised when the colloids are close together. The overall behaviour can be described in terms of an effective attraction between the colloids, which can give rise to aggregation or phase separation into polymer-rich and colloid-rich fluid phases, often referred to as colloidal vapour and colloidal liquid (respectively). Thorough explanations of the depletion interaction are provided in, for example, references [1,[31][32][33][34][35]; depletion forces have been the focus of a number of studies, including [36][37][38][39][40][41][42][43]. In particular, Binder et al. [35] have provided a detailed appraisal of the Asakura-Oosawa (AO) model, and its scope. Accordingly, we provide here only a brief summary.
The AO model is defined through the interactions between the species in the colloid-polymer mixture. Both colloid and polymer are modelled as spherical; the colloids are hard, while the polymers are hard with respect to the colloids but are treated simply as soft, interpenetrating spheres with respect to each other. This model, or the related AO approximation, in which the polymers are treated as ideal chains, has formed the basis of much of the subsequent theoretical work on colloid-polymer systems. An important example in the context of our current paper is the work of Gast and co-workers [16,18], who used thermodynamic perturbation theory [44,45] to predict the phase behaviour of an AO model colloid-polymer mixture; they were able to account for the change from fluid-crystal to fluid-fluid coexistence with increasing polymer size, observing the transition at a size ratio q = R g /R c 1/3, where R g is the radius of gyration of the polymer and R c is the radius of the colloid. Gast et al. [16] discussed at some length the limitations to their approach brought about by the neglect of three-and higher-body interactions. Later Lekkerkerker et al. [46] introduced the free-volume approach for the determination of the phase behaviour, which captures some of the effects of many-body interactions, as well as avoiding the assumption of Gast et al. that the polymer concentration is the same in the coexisting phases thereby incorporating partitioning of the polymer between phases.
To date, the phase behaviour of colloid-polymer mixtures has been addressed in a number of theoretical studies, in addition to those already mentioned; see, for example, references [34,43,[47][48][49][50][51][52][53][54][55][56][57]. Among these, we highlight the studies [43,[53][54][55]58] based on Wertheim's first-order thermodynamic perturbation theory (TPT1) [59,60]; these deal with the polymer at the monomer segment level [61], in contrast to studies based on the AO model or the AO approximation. While it is gratifying that theoretical studies can be used to explain the phase separation seen experimentally, to assess the quality of the underlying models one still requires accurate verification of theories by direct molecular simulation.
Simulations of prototypical colloid-polymer systems have been conducted at different levels of sophistication, e.g., using one-component [48,62] and two-component [37,38,[63][64][65][66][67][68][69][70][71][72] descriptions; see reference [35] for a review of simulation studies in relation to the AO model. In the approach pioneered by the Hansen group in Cambridge [38,[73][74][75][76][77], the 'polymer as soft colloids' model is adopted; the simulation of sufficiently large systems is facilitated by a coarse-grained approach in which the polymer-polymer interactions are incorporated using an effective potential that acts between the centres of mass of the chains. Vink and Horbach [78] have also introduced 'cluster moves' that allow tractable simulation of the AO and related models using grand canonical Monte Carlo [70,[79][80][81]. An alternative approach to simulate large systems has been lattice Monte Carlo, using an explicit (monomer) description of the polymers ( [36,82]). However off-lattice simulations of a system with an explicit model are rare. Starr et al. [83] simulated a system consisting of Lennard-Jones (LJ) chain polymers and nanoparticles, represented by LJ spheres bonded together in an icosahedron, however these authors did not observe phase separation; rather, they reported that the nanoparticle clustering was akin to equilibrium polymerisation. Lo Verso et al. [84] have simulated the interesting and related system of polymer-brush-coated spherical nanoparticles, but the nature of their model precludes phase separation of the type that is of interest in our current work. More recently, Meng et al. [72] simulated systems of LJ-chain polymers and nanoparticles conceptualised as multiple LJ centres grouped together to form larger colloidal spheres; the potentials of interaction between the colloidal nanoparticles, and between nanoparticles and polymer LJ segments were obtained via integration of all the LJ sites in the nanoparticles. Although these authors examined the phase transitions in the system, these were of solid-fluid nature (rather than fluid-fluid), except when monomeric LJ particles were used as solvent in place of the polymers. Also of particular relevance to our current study is the series of papers by Panagiotopoulos and co-workers [71,[85][86][87]. Mahynski et al. [85] performed continuum simulations to complement their more extensive study [71] in which the polymer chains were constrained to lie on a lattice. However, in these two contributions the authors were primarily interested in the so-called protein limit, where the polymer radius of gyration exceeds the colloid radius. By contrast, in our current work, we consider a system corresponding to the colloid limit, which is characterised by a size ratio q = R g /R c < 1. Subsequently, Mahynski et al. have carried out continuum simulation studies of colloid-polymer mixtures wherein both colloids and polymer-chain segments were first [86] represented using repulsive LJ potentials of Weeks, Chandler and Anderson (WCA) [88] type and, in their most recent study [87], using more-general attractive potentials. Although the size ratios considered in these latter studies correspond to the colloid limit, the authors were concerned with the role played by the polymers in the crystallisation of the colloids. Our primary concern is a continuum simulation study of a fluid-fluid phase-separated mixture. We focus our analysis on a colloid-polymer system in which the colloid diameter is ∼20 times that of the polymer-segment diameter. The computational expense of simulating such a system in the protein limit (rather than the colloid limit) is prohibitively high, due to the considerable length of the polymer chains that would be required.
Although there is a widely held view that the physics of the colloid limit is now well understood [56,[89][90][91], aspects in which this understanding is incomplete have been highlighted in some recent publications. For example, Shvets and Semenov [92] have pointed out that the mechanism of depletion stabilisation is, as yet, not fully understood; more recently, Binder et al. [35] have suggested that the success of various generalised free-volume theories, such as that of reference [93], in accounting for the experimental phase boundaries should not be taken as a validation of the underlying models for the effective colloid-colloid interactions. Most striking of all, in a quite fascinating study of a mixture of 3-methacryloxypropyl trimethoxysilane (TPM) colloids with polyethylene oxide depletant, characterised by q ∼ 0.4, Feng et al. [94] have observed hitherto unsuspected re-entrant phase behaviour, demonstrating that the conventional picture of the colloid-polymer depletion systems misses an important feature; weak temperaturedependent adsorption of the polymer on the surface of the colloid can give rise to a temperature dependence of the depletion interaction. Although an investigation of such behaviour is beyond the scope of our current investigation, it serves to illustrate that there yet remains much to learn concerning the physics of the colloid limit. From our own perspective, and more pertinent to our current work, as a final observation in this context we note that in theoretical approaches based on the AO model, it is usual to assume that the effective interaction between colloids and polymers is the same throughout the system, even though the environment in which the polymers find themselves is substantially different in the colloidal vapour and liquid phases. Continuum molecular simulations provide a vehicle for testing the validity of this assumption, as well as providing a means for testing our understanding of the physics of colloid-polymer mixtures in the colloid limit in general.
In our current contribution, we use the TPT1 approach of Wertheim [53,59,60] to identify a phase point well inside the fluid-fluid coexistence region for mixtures of model hard-core colloid and polymer systems. The same mixture is then studied using continuum molecular-dynamics (MD) simulation.
The remainder of this paper is laid out as follows: in Section 2 we describe the models used in our work, provide a short description of the Wertheim TPT1 approach, and set out our simulation technique. In Section 3 we present in brief our theoretical calculations of the fluid-phase diagram of our colloid-polymer mixture, and describe in detail the results of our simulations of the mixture. Finally, we present our conclusions in Section 4.

Models and theory of the colloid-polymer system
Very accurate versions of the statistical associating fluid theory (SAFT), which incorporate Wertheim's TPT1 treatment of chains [59], have been developed recently for molecules formed from segments interacting through the Mie potential (SAFT-VR and SAFT-γ Mie) [95,96]. The analytical SAFT-γ Mie equation of state can be employed to develop coarse-grained force fields for use in large-scale molecular simulation of a variety of systems [97][98][99][100][101] including surfactants [102,103], long chain molecules and polymeric systems [104]. In future work we plan to examine colloidal suspensions represented with more-realistic coarse-grained interactions based on the SAFT-γ Mie approach. However, in this first contribution, we consider an athermal colloid-polymer mixture, using models based on the HS potential. Any observed effective attractive interactions can then confidently be ascribed to the depletion phenomenon, without a need to decouple any attractions incorporated within the potential model. The focus here is on the fundamental nature of the phase separation observed in a well-characterised model system, to highlight the physical features that lead to the underlying behaviour.

Athermal mixture of hard spheres and tangent hard-sphere chains
The underlying model in this work is that of HSs (representing the colloids) and flexible chains of tangent HS segments (representing the polymers); the relative diameters of the spherical moieties are chosen such that the diameter of the colloid is much larger than that of the spherical monomeric segments of the polymer chain. The focus of our study is on a colloid-polymer system characterised by a colloid-tomonomer HS diameter ratio σ c /σ m = 20, and a polymer chain length m p = 100. An important descriptor of colloidpolymer systems is the ratio, q, between σ p , the effective 'diameter' of the polymer, and σ c , the diameter of the colloid; q = σ p /σ c = R g /R c , where R g is the radius of gyration of the polymer and R c = 10σ m is the radius of the colloid. Based on the analysis of reference [105], one expects the radius of gyration of 100-segment tangent HS polymer chains to lie in the interval (5.8σ m R g 7.3σ m ), depending on the density, which would correspond to 0.58 q 0.73. Thereby, our chosen colloid-polymer system corresponds to the so-called colloid limit, which is characterised by q < 1. It is also important to note, however, that the size ratio for our system lies above the threshold value (q ∼ 0.5) below which traditional thermodynamic perturbation theory is inapplicable as suggested by Pelissetto and Hansen [90]. The phase behaviour of the colloid-polymer system is examined both using the Wertheim TPT1 relations and direct molecular simulation. In our simulations the (discontinuous) HS interaction is modelled using a (continuous) pseudo-hard-sphere (PHS) potential [106]. This approach has the advantage that it allows straightforward and efficient simulation using 'off-the-shelf ' MD software, and the deployment of parallel computation with no loss of accuracy. The PHS potential employed is a cut-and-shifted Mie potential in the spirit of the WCA treatment [88], and is expressed as a function of interparticle separation, r, as where is an energy parameter representing the depth of the potential well. Unlike the true HS potential, this interaction is not athermal and one is required to set a temperature. It has been previously shown [106] that an excellent representation of (athermal) HS systems is provided at a reduced temperature of T * = k B T/ = 1.5, where k B is the Boltzmann constant and T the absolute temperature. At this condition, the size parameter, σ , effectively matches the HS diameter, i.e., no density correction is required. Furthermore, to facilitate a direct comparison with the Wertheim TPT1 description [53,59,60], no additional bonding, bending, or torsion intramolecular interactions are considered when modelling the tangent HS chains and the distance between connected segments in the polymer is kept fixed at the value of σ . It will be convenient to define a scale of reduced (dimensionless) pressure. The appropriate reduced pressure for our model system is defined in terms of the energy scale as P * = P σ 3 c / . One should note that the well depth enters the definition of the reduced pressure because of the need to make a connection with the PHS potential model employed in the molecular simulations; its relationship to the natural choice of reduced pressure of the athermal system within the TPT1 description P * TPT1 = P σ 3 c /(k B T ) follows simply from P * TPT1 = P * /T * where T * = k B T/ = 1.5 throughout.

Wertheim first-order perturbation theory (TPT1)
We provide here only a brief description of the important features underlying the Wertheim TPT1 approach; a much more detailed exposition of the use of the algebraic equation of state in relation to the colloid-polymer system can be found in reference [53]. The potential of interaction between the spherical moieties is the HS potential: where r represents the distance between the centres of the spheres, and σ is the contact distance or hard-core diameter.
For the interaction between a colloid and a monomeric segment of the polymer chain, the contact distance is obtained simply from the arithmetic mean, (σ c + σ m )/2. Following the TPT1 approach [53,59,60], the Helmholtz free energy, A, of a binary mixture of N c colloids and N p polymers (comprising N m = N p m p monomer segments) is written as a sum of contributions: A ideal is the ideal free energy of a mixture of N = N c + N p colloid and polymer particles; A HS is the residual free energy of a binary mixture of N c HSs of diameter σ c and N m HSs of diameter σ m ; A chain is the contribution to the free energy of bonding together the N m HS segments into N p chains of length m p . The HS free-energy contribution A HS is obtained as a function of reduced densities using the expression of Boublík [107] and Mansoori et al. [108]. This expression has been found to provide for an accurate description of binary mixtures of HSs of significantly different size (corresponding to diameter ratios of up to 20:1 [109], the ratio adopted in the current work) over a wide range of compositions. Detailed expressions of each of the individual terms in Equation (3) can be found in reference [53]. The thermodynamic properties of the system are obtained using standard relations, and the fluid-phase equilibria are determined by ensuring that the pressure, P = −(∂A/∂V) N, T , and chemical potential of each component, μ i = (∂A/∂N i ) V ,T ,N j =i for i = c, p, are equal in the coexisting phases; the temperature plays a trivial role for athermal systems of this type. A limitation of the TPT1 approach is that information about the structure of the chain is lost, because the m-body distribution function for the tangent polymer chain is approximated by that of the reference HS monomeric mixture: g(1. . .m p ) = g(1, 2). . .g(m p − 1, m p ) = g(σ p ) m p −1 [53]. As a consequence, differences in the conformation of the chain cannot be treated at this level of approximation [110]. Further, there is no explicit description of the end-to-end distance or the radius of gyration of the chain. The second virial coefficient obtained with the TPT1 treatment scales linearly with the chain length m p [111,112] instead of the expected m 3ν p , where ν is the Flory exponent; [55] moreover, the use of TPT1 leads to ν = 1/2 [55], the ideal-chain value, whereas for self-avoiding chains the Flory exponent takes the value ν ∼ 3/5 in the dilute regime. As a result the swollen or dilute regime is not described accurately with the TPT1 approach [110]. Indeed, Mahynski et al. [71] have found that coexistence densities obtained using TPT1 are systematically lower than those obtained in their fine-grained lattice simulations; they attributed this inconsistency to the implicit incorrect scaling of TPT1 in the dilute regime. Within our current study, however, this issue will be of no serious consequence, since the theory is used simply to locate a state point well inside a dense twophase (liquid-liquid) region of the colloid-polymer mixture phase diagram.

Results and analysis of model colloid-polymer system
Although the use of the PHS WCA(50,49) potential represents an efficient route to perform molecular simulations of systems with the marked size asymmetry and high connectivity inherent in colloid-polymer systems, the large number of interaction sites involved makes such simulations computationally demanding. Consequently, even if one  [53,59,60]. The red dashed curve represents the polymer-rich phase; the green continuous curve represents the colloid-rich phase. The filled circle represents the selected state point at which the initially demixed and the initially mixed simulations are performed. The triangle and square represent the polymer-rich and colloid-rich coexisting phases, respectively, as used for the single-phase simulations (described in Section 3.2.2); the asterisk denotes the predicted critical point of the mixture.
uses efficient parallelised MD software such as DL_POLY [113], the currently available hardware still restricts the size/number of simulations that can be performed. By contrast, calculations with the algebraic TPT1 expressions require a very small fraction of the computer time. Accordingly, to identify a suitable phase point for simulation, an initial examination of the phase space is made using TPT1.

Theoretical examination of the phase space and identification of phase point for simulation
The fluid-phase equilibria of our colloid-polymer system characterised by σ c /σ m = 20 and m p = 100 are determined using TPT1. Representations of the phase behaviour in terms of the reduced pressure vs. mole fraction and reduced pressure vs. packing fraction planes are shown in Figure 1(a) and 1(b), respectively. The polymer mole fraction is defined as x p = N p /N, and the overall packing fraction (a reduced, or dimensionless, density) of the system as where V is the system volume. The packing fractions of the individual species, η c and η p , are obtained analogously, including only the first or second parenthesised term in Equation (4), respectively. It would be usual to define a scale of reduced pressure for such a system as ; in order to avoid confusion we reiterate here that in our current work we adopt instead a scale related to the energy scale of our PHS potential, P * = P σ 3 c / , where in our case P * TPT1 = P * /T * and T * = k B T/ = 1.5. The choice of thermodynamic state (pressure and overall composition) in which to carry out the direct molecular simulations is made keeping in mind the following constraints: (1) the pressure of the system should be significantly higher than the lower critical point of the mixture, corresponding to a point well inside the binodal boundary; (2) the pressure of the system should be sufficiently low so that the fraction of colloids in the polymer-rich phase is not so small as to provide unacceptable statistics; (3) the estimated volume of each coexisting phase should be of the same order of magnitude (for a fixed overall composition); (4) the liquid phase should be thermodynamically stable (relative to the solid).
It is well known that TPT1 (in common with other analytical theories of fluids) provides a classical mean-field description in the vicinity of the critical point. Theoretical calculations usually deviate from experiment or simulation in this region; in this case the calculated critical solution pressure is expected to be lower than that obtained from simulations (or experiments) of equivalent systems. Constraint (1) thereby prevents the accidental selection of a phase point in the single-phase region. Based on constraints (1-3), an appropriate state point for simulation corresponds to a reduced pressure of P * = 10, a polymer mole fraction of x p = 0.67, and an overall packing fraction of η total = 0.28 (at a reduced temperature of T * = 1.5). The sole use of TPT1 takes no account of the possibility that a liquid phase may be metastable relative to the solid phase, since solid phases are not explicitly considered in the development of the liquid-state perturbation theory. By comparison with the phase diagrams provided in figure 12 of reference [31] our chosen state point, although close to the triple point, is expected to correspond to a fluid-fluid and not fluid-solid region, thereby we expect constraint (4) to be fulfilled as well. Accordingly, we select this state point chosen for our simulations, as indicated in Figure 1 by a filled circle.

Molecular-dynamics simulation
The dimension of the system is set to 5 × 5 × 10σ 3 c (L x × L y × L z ) in units of the colloid diameter. The number of molecules is then obtained from the packing fraction of the chosen state point; for the selected dimensions this corresponds to N p = 256 polymers (totalling N m = 25, 600 spherical segments) and N c = 128 colloids. All simulations are performed at a reduced pressure of P * = 10 and reduced temperature of T * = 1.5, and are carried out using standard periodic boundary conditions, with the Verlet leap-frog algorithm [114] as the time-integration method. The reduced time step is t * = tσ −1 √ /ω, where ω = 1 is the mass of a PHS segment. For the polymer chains, the bond distance is kept fixed using the SHAKE algorithm [115] with a tolerance of 10 −5 σ m . Full details of the simulation technique are provided in the appendix to reference [106].
Different types of simulations of the colloid-polymer system with σ c /σ m = 20 and m p = 100 are carried out: an initially fully mixed system; an initially demixed system; and two separate simulations of the bulk phases corresponding to the compositions and densities of the coexisting fluid phases as determined with the TPT1 approach. We describe first the simulations of the phase separation, and then discuss the single-phase simulations.

Study of the phase separation
Two simulations of the colloid-polymer mixture are performed independently, one starting from an initial configuration in which pure colloid and polymer phases are brought into contact, and the other starting with a fully mixed system (cf. Figure2(a) and 2(b), respectively). The aim is to reach a stable phase-separated system starting from opposite extremes, thus establishing proper equilibration. Both systems are first simulated under isobaric-isothermal (NPT) conditions to equilibrate the initial configurations, establishing the appropriate pressure, and then in the canonical (NVT) ensemble; except where stated explicitly, results given correspond to canonical simulations.
Snapshots of the final configurations with both approaches are presented in Figure 3. Visually, colloid-rich and polymer-rich regions are clearly identifiable, although the system does not appear to be separated into two welldefined bulk regions, with a clearly demarcated interface. This does not necessarily indicate that phase separation has not occurred. For example, it may be that one is unable to stabilise an interface for a system of this size. Estimates have been made of the interfacial width of AO colloidpolymer mixtures. Using density-functional theory, Brader et al. [116] estimated the width to be ∼σ c close to the triple point while Bryk [55] suggested the width to be a few colloidal diameters wide, and increasing with increasing chain length of the polymer. Vink et al. [79] obtained interfacial widths ∼7-9σ c using grand canonical Monte Carlo simulation (see figure 3 of reference [79]). Since our simulations are periodic, two interfaces are required in order for the system to stabilise into well-defined colloid-rich and polymerrich regions. For interfacial widths of ∼5σ c then even the relatively large system size chosen in our current work will be insufficient to stabilise two phase-separated regions with two corresponding interfaces because the longest side of the simulation cell, 10σ c , would be no greater than the width of two adjacent interfaces. Moreover, due to the low interfacial tension of colloid-polymer systems, the interfaces are characterised by large local fluctuations not only in width [79] but also in composition [55,79,116]. Accordingly, even when phase separation has occurred, a clear profile of bulk regions separated by an interface is unlikely to be seen for all but the largest systems.
To establish whether or not phase separation has occurred, careful analysis of the polymer-rich and colloid-rich regions apparent in Figure 3 is required. One cannot simply sample properties along the long (z) axis of the simulation cell, however, since the molecular distributions (polymer and colloid) in the x and y directions within a given sampling volume are no longer uniform. A different approach is required to account properly for these differences in the distribution of the local density of each species. To this end, an algorithm is implemented that is based on the procedure [117] used to sample phase coexistence in systems where the spacial domains of the phases are not fully segregated [118]. In essence, the algorithm consists of sampling the local compositions by randomly analysing small sub-regions of given configurations; a histogram of the local compositions exhibits bimodality if phase separation is present.
We incorporate two important modifications of the procedure to account for the highly asymmetric nature of the colloid-polymer system. In the original scheme [117] onedimensional profiles are generated by sampling cuboidal sub-systems, randomly distributed in the simulation cell.  This is not suitable for use in the study of systems with such a large asymmetry in particle/segment size; the colloid spheres are so large that the diameter of a colloid may be of similar order of magnitude to the width of the sampling volume. A consequence of this is that a sampling volume could contain portions of many colloidal particles, yet be counted as devoid of colloid should none of the colloid centres lie within the sampling volume. Accordingly, using Equation (4) directly (but with V representing the sampling volume) to obtain the local packing fraction leads to unsatisfactory results for the local colloid packing fraction, which is thereby returned as a very noisy function. Here, we adapt the algorithm to account for the portions of the colloids (and only those portions) that lie within the sampling volume, irrespective of where the centres of the colloids lie.
For computational simplicity, spherical sampling volume elements are employed, thus the portion of a colloid to assign to a volume may be calculated as the intersection of two spheres, which is a simple mathematical object referred to as a spherical lens; the volume of a spherical lens depends only of the radii of the spheres and the distance between their centres. The specific procedure employed to calculate the volume of a spherical lens is presented in Appendix A.
A radius of r samp = 12.5σ m is chosen for the sampling sphere, based on the best compromise between a large sampling sphere, which more efficiently samples the colloidrich regions, and a smaller sampling sphere, which more effectively samples the polymer-rich regions. (A discussion of the effect of the sampling-sphere radius on the resulting density distributions is provided in Appendix A.) The analysis is conducted for the last 5 × 10 5 time steps of each simulation; the resulting probability distributions of the density, P(η total ), are displayed in Figure 4. It can be seen that both runs (initially mixed and initially demixed) lead to very similar average density distributions, featuring a sharp peak corresponding to the packing fraction, η pol total , of the polymer-rich phase and a broader peak corresponding to the packing fraction, η col total , of the colloid-rich phase. This is significant as it indicates that the simulations have  The dashed curve represents the initially mixed simulation, while the continuous curve represents the simulation in which polymers and colloids were initially demixed (see Figure 2). The average is accumulated over the last 5 × 10 5 time steps in each case, using 10 5 samples per time step, with a sampling-sphere radius of r samp = 12.5σ m . The vertical green dashed and red dot-dashed lines represent the densities of the two phases coexisting predicted with the Wertheim TPT1 approach.
converged to configurations with essentially the same density distribution even though, as can be discerned from the snapshots of the initial configurations in Figure 2, the initial density distributions of the two simulations were entirely different. The average distributions obtained from both the initially demixed state and the initially mixed state indicate a total packing fraction of η pol total ∼ 0.029 for the polymer-rich phase. Although the peaks corresponding to η col total are very broad the maxima are almost coincident at a total packing fraction of η col total ∼ 0.43 for the colloid-rich phase. The location of these in relation to the pressure-density phase diagram (Figure 1(b)) is illustrated in Figure 5 (the 'error bars' reflect the broad nature of the two peaks, estimated as half their widths at half height). It is interesting to note that the packing fraction of η col total ∼ 0.43 is significantly lower than the value of 0.49 where the pure HS system is known to exhibit a transition from a fluid to a solid phase [119].
It is striking that the maxima of the average density distributions (i.e., the most-probable densities) obtained from the very different starting configurations coincide closely with the densities of each of the phases (polymerrich and colloid-rich) calculated using the Wertheim TPT1 approach; note that the positions of these maxima are found pol total (filled red square), and colloidrich, η col total (filled green square), phases for the colloid-polymer system with a size ratio of σ c /σ m = 20 and polymer length m p = 100 at the reduced pressure P * = 10. The coexisting packing fractions determined from the two peaks in Figure 4 are seen to coincide closely with the predictions of the Wertheim TPT1 approach (dashed and continuous curves); the 'error bars' reflect estimates of the widths of the two peaks (see the text for details). As in Figure 1, the dashed curve represents the polymer-rich phase, the continuous curve represents the colloid-rich phase and the asterisk denotes the critical point of the mixture.
to be essentially independent of the sampling-sphere radius employed (see Appendix A). This close correspondence between the theoretical and simulated values is not in line with the findings of Mahynski et al. [71], who obtained a systematic disagreement between coexistence densities for a colloid-polymer system calculated with TPT1 and those obtained in their simulations on a fine-grained lattice. It is important to note, however, that in our current simulations, the polymers are not in the dilute regime (in contrast to those of Mahynski et al.), since the overall packing fraction of even the lower density polymer-rich phase, η pol total ∼ 0.029, is in excess of 0.001, the value at which the crossover from dilute to melt behaviour for HS chains is observed [105]. The Flory exponent for swollen HS chains is found to be ν dilute ∼ 0.60, whereas the exponent for HS chains at η = 0.029 is expected to be ν ∼ 0.53 [105]; this is much closer to the ideal-chain exponent of ν TPT1 = 0.5 implicit to the TPT1 approach [55]. Moreover, as will be discussed in Section 3.2.2, the polymers behave as though the system is at an effective packing fraction in excess of 0.029, whereby their behaviour will be characterised by an exponent still closer to the ideal-chain value. It may also be significant that Mahynski et al. focused their study on the so-called protein limit, i.e., the regime in which q = R g /R c > 1 whereas the mixture examined in our current study lies in the regime corresponding to the colloid limit where q < 1.
The close agreement between the TPT1 predictions and direct simulation obtained in our current work not only provides mutual corroboration of the two approaches, but also lends strong support to the suggestion that phase separation has been achieved in both the initially demixed and initially mixed simulations. Moreover, one can thus be confident that the systems have reached equilibrium, allowing for the fact that (as previously discussed) the size of the system may be too small to allow two well-defined bulk regions separated by clearly demarcated interfacial regions to develop. This confirms that, regardless of the starting conditions (fully mixed or completely demixed), a mixture of our model colloids and polymers at P * = 10 with x p = 0.65 and η total = 0.28 is unstable and will phase-separate into a colloid-rich and a polymer-rich phase.
To verify that the phase separation observed is indeed a fluid-fluid separation, the mean squared displacement (MSD) of the colloids is assessed over the final stages of both simulations. Through the last 5.5 × 10 5 time steps (i.e., ∼8% of the overall length of the simulation), the MSD is found to increase linearly with time. Moreover, the colloids have moved (on average) more than 20 times their diameter during the course of the simulation. One can therefore conclude that the colloids are behaving as a fluid in each phase, and that neither a solid nor a glassy state has been formed. This conclusion will be confirmed by the additional analysis provided in the next section.
Before moving on to discuss our examination of the individual phases, it is interesting to analyse our two-phase system in terms of the interfacial tension that is implicit from the information accrued thus far, and to make comparisons with existing work. Adopting the approach of Binder [120][121][122] the interfacial tension γ of our two-phase colloid-polymer system can be estimated from the height P of the peaks in the probability distributions of the density (measured relative to the 'interfacial' states), which is related to the surface free energy F surf and interfacial tension through where A represents the area of the interface. From Figure 4 we estimate P ∼ 0.01, and we take A ∼ 25σ 2 c , corresponding to the smallest cross section of our simulation box. We thus obtain an estimate of the dimensionless interfacial tension as Moncho-Jordá et al. [52] used square-gradient theory to compute the interfacial tension for systems of colloids with interacting and non-interacting polymers in a good solvent. In their figure 4, these authors provide values of the interfacial tension as a function of the difference η in the packing fraction of the two phases, and the polymer/colloid size ratio q. From this figure, one can read off γ R 2 c /(k B T ) ∼ 0.05 (where R c is the colloid radius) for η ∼ 0.41 and q = 0.6 (corresponding to the system in our current work; see Section 3.2.2 for our estimate of q ∼ 0.6), which yields an interfacial tension of γ * ∼ 0.2, in remarkably good agreement with our estimate of γ * ∼ 0.18.
To make a comparison with an experimental system, one first needs to define a length scale. Aarts et al. [123] studied the free fluid-fluid interface in a phase-separated colloid-polymer dispersion using laser scanning confocal microscopy. The size of the colloids was measured as σ c = 142 nm, and the size ratio q = 0.6, the same as in our current study. Using σ c = 142 nm to define a length scale for our simulations, and assigning an ambient temperature (T ∼ 300 K), our simulated interfacial tension corresponds to γ ∼ 40 nN m −1 , which is in order-of-magnitude agreement with the ultra-low interfacial tensions reported in reference [123]. The agreement we obtain with this study (and the theoretical work of Moncho-Jordá et al. [52]) not only further supports the conclusion that our colloid-polymer mixture is indeed phase-separated, but also adds further insight into the absence of a clearly defined interface in our simulations. Aarts et al. [123] were able to directly observe thermal capillary waves in their study; from their figure 2, the amplitudes and wavelengths of the capillary waves can be seen to be of the order of a few microns, corresponding to a few tens of colloid diameters. Recalling that the longest dimension of our simulation cell is 10σ c , one can thereby again conclude that a very much larger simulation cell would be required in order to stabilise a well-defined interface between the two phases.

Studies of the separate polymer-rich and colloid-rich phases
To properly understand the colloid-polymer system in terms of its structural and thermodynamic behaviour, it is useful to examine each of the coexisting phases separately. Based on the evidence provided by the probability distribution of the density (cf. Figure 4(b)) and quantified in Figure  5, we henceforth assume that the two phases exactly match the theoretical predictions from TPT1 [59] and, accordingly, the density and composition of each phase are taken to be those indicated by TPT1 (see Figure 1(a) and 1(b)).
(Although this assumption results in very slightly higher densities in each phase than those indicated by the peaks of the distributions in Figure 4, one thereby avoids possible complications arising from the statistical uncertainty in the compositions of the two phases that stems from the small Table 1. Details of the simulations of the separate polymer-rich and colloid-rich phases as predicted with the TPT1 approach for the colloid-polymer system with a size ratio of σ c /σ m = 20 and polymer length m p = 100 at the reduced pressure of P * = 10; L is the length of the (cubic) simulation cell. number of colloids present in the simulation.) Using this information, separate cubic simulation boxes are constructed, each representing one of the two phases; details are given in Table 1. Simulations of each phase are performed in the canonical (NVT) ensemble.
To examine the dimensions of the polymers, we calculate their average radius of gyration, R g , in each of the two phases. R g is defined as the root-mean-square distance of the monomers from the centre of mass (CoM) of the chain: The sum is over the entire chain length of the polymer m p (in the current work m p = 100 segments); R CoM represents the position of the centre of mass of the chain (assuming all segments have the same dimensionless mass), and is given by and r i similarly represents the position of the ith segment. The probability distributions obtained for the radii of gyration are presented for the polymer-rich and the colloidrich phases in Figure 6. Both distributions exhibit a peak at R g /σ m ∼ 6, so that the size ratio for our system is q = R g /R c ∼ 0.6. One can discern from Figure 6 that the difference between the most-probable R g values of the polymer in the two phases is less than 6%. It is important to note that, although this difference is small, this does not mean that the structure of the polymers in the two phases will be similar. Indeed, the similarity in the dimensions of the polymers may itself seem a surprising result, considering the difference in the packing fractions of the polymer-rich (η pol total = 0.0427) and the colloid-rich phase (η col total = 0.4581). Even though the overall packing fraction of the colloid phase is only ∼3% lower than the value at which pure HSs undergo a liquid-solid transition, this suggests that the polymers are able to occupy the interstitial regions, which are relatively large empty regions when compared to the volume occupied by the monomers.
An order-of-magnitude calculation of the relative volumes involved is helpful. In the close-packed limit, in which the colloids form a tightly packed fcc lattice, the packing fraction is η fcc = 0.74; in other words 26% of the space is free. In an fcc unit cell with four colloids, given that σ c = 20σ m , there is an empty volume equivalent to that of 496 monomeric segments of the polymer. In the colloidrich phase the ratio of polymers to colloids is 1/8 (see Table 1), so the number of polymers that would correspond to four colloids is ∼0.5, comprising ∼50 monomers. In other words, even in the case of close-packed colloidal particles, there would be space for ∼50 monomers in a unit cell; of course, the actual density of the colloids is substantially lower than that of a close-packed system, with even more volume available to the polymers. Using Equation (4) but replacing V with the empty volume that has just been discussed, a 'local' polymer packing fraction of η p, local = 0.013 would be obtained. However, the polymers and colloids share the same space and the colloids are not stationary, so such a calculation is rather crude. In reality, the polymers behave as though the system is at an effective packing fraction which will be higher than η p, local = 0.013 but much less than η col total = 0.4581. It has been shown that for (pure) HS polymer chains of degree 50 ≤ 2000 the crossover from the swollen (dilute) to the melt regime takes place between η ∼ 0.001 and η ∼ 0.17 [105]. If the effective packing fraction of the polymers in the colloid-polymer mixture lies within this interval, then one way to establish its value is to perform simulations of pure polymers at varying packing fractions in this range, comparing the probability distributions of the dimensions obtained from these simulations with those of the mixture. Accordingly, also displayed in Figure 6 are distributions obtained from the simulation of the corresponding pure-polymer systems at packing fractions η p = 0.05 and η p = 0.075. Gratifyingly, the former lies close to the data corresponding to the simulation of the polymer-rich phase of the colloid-polymer mixture, for which the packing fraction η pol total = 0.0427, but centred on a slightly smaller value of R g , consistent with the slightly higher packing fraction of the pure-polymer simulation. The distribution obtained from the simulation at η p = 0.075 coincides quite closely with the data corresponding to the simulation of the colloidrich phase of the mixture, indicating that the polymers in this phase behave as though the system is at an effective packing fraction close to, but a little below this value.
Frenkel [124] has pointed out that hard-core spherical particles with an attractive range spanning less than 20% of their diameter do not exhibit a stable vapour-liquid transition, since the vapour-liquid critical point drops below the triple point. This concurs with the 'rule of thumb' for the depletion interaction between colloids in colloid-polymer mixtures, prescribed by Likos [32], that when the range of attraction is more than 20% of the range of the repulsive cores, the liquid phase is stable and the phase diagram features two fluid phases and one solid phase; such behaviour has been seen experimentally, e.g., by Calderon et al. [3]. For polymers in the swollen (dilute) regime, one would expect that the thickness of the depletion layer around the colloids is close to R g and that the range of the depletion interaction for our system is ∼2R g = 12σ m = 0.6σ c . [1,125] This would correspond to a range of the attractive depletion well which is in excess of the threshold value (∼σ c /5), below which the colloidal liquid phase is expected to become metastable relative to the solid [124]. However, since the effective polymer packing fraction corresponds to the crossover (or semi-dilute) region between swollen chains and the melt, the range of the depletion interaction is determined not by R g but, rather, by a correlation length ξ , often referred to as the blob size of the polymer [33,[125][126][127]. Using information pertaining to the dimensions of HS polymer chains (see table III of reference [105]) together with a simple scaling argument [105], the blob radius ξ of the 100segment HS polymer chains at η p = 0.075 can be estimated to be of the order of 3σ m to 4σ m , whereby the range of the interaction would be approximately 6σ m to 8σ m σ c /3 (see Appendix B for a description of the scaling argument). This value still exceeds the theoretical threshold below which the liquid phase is expected to be metastable, and corresponds to a depletion thickness to colloid radius ratio of q s ∼ 0.3-0.4. The phase diagram for a system with this value of q s may be expected to correspond roughly to that provided in figure 1.1 of reference [34], which features 'gas-liquid', 'gas-solid' and 'fluid-solid' phase boundaries. It is noticeable that the solid phase is characterised by η c 0.54, corresponding to the well-known coexistence packing fraction for pure HS, and significantly in excess of the total packing fraction (polymers + colloids) of the denser of the two phases we find in our simulations. One can also discern that the packing fractions of the liquid phase along the gas-liquid binodal lie between η c ∼ 0.27 (at the critical point) and η c ∼ 0.46 (at the triple point), spanning our value of η col total = 0.4581 for the colloid-rich phase. Taken together, these observations provide strong corroboration of our earlier conclusion that the phase transition we observe is fluid-fluid, and not fluid-solid in nature.
It is known that HS chains in simulations of pure polymers at densities corresponding to those of our current work (η p = 0.05 and η p = 0.075) are not, on average, extended so that polymers behave as though spherical [105,128]. The close coincidence between the distributions relating to the pure polymers and those relating to the corresponding simulations of the colloid-polymer mixture suggests that, on average, the polymers in the mixture are also spherical and not extended, even in the colloid-rich phase. This implies that although the polymer chains are mutually interpenetrable, effectively they behave as spherical entities that are almost impenetrable by the colloids -rather in the spirit of the AO model. This is a helpful conclusion as it implies that one may usefully probe the structure of the fluid phases using molecular distribution functions of the CoM.
The radial distribution function, g(r), is defined as the ratio of the local number density of particles or molecules from a central particle to the bulk number density, ρ. It may be evaluated as [129] g(r) = n(r) where n(r) is the number of particles in a shell bounded by the distances r and r + δr from the given particle, and the angular brackets represent the ensemble average over all particles. The internal structure of the two phases is presented in Figure 7 in the form of the colloid-colloid, g cc (r), colloid-polymer, g cp (r), and polymer-polymer, g pp (r), radial distribution functions; the colloid-rich phase is represented in (a), and the polymer-rich phase in (b). As one might expect, the colloid-colloid distribution function g cc (r) in the colloid-rich phase (the green curve in Figure 7 number of colloids in the polymer-rich phase the statistics relating to the colloids in this phase are poor (to keep the simulation tractable, the simulation contains only eight colloids). Consequently, g cc (r) is too noisy to draw any firm conclusions about the behaviour of the colloids in this phase. We consider next the colloid-polymer distribution function g cp (r). In the colloid-rich phase (cf. Figure 7(a)), g cp (r) has its first peak almost exactly at σ cp = (σ c + σ p )/2, with subsequent weaker peaks at 2σ cp and 3σ cp , just as one would expect for a liquid-like radial distribution function for mixtures of spheres of diameter σ p and σ c . This corroborates the earlier conclusion that the polymers behave as spherical entities from the perspective of the colloids. Moreover, the location of the first peak at σ cp is a further indication of the relative impenetrability of the polymers to the colloids. The widths of the peaks provide an indication of the relative softness or hardness of the colloid-polymer interaction; although broader than those of g cc (r), it is interesting to note that the peaks are nevertheless quite sharp. In the polymer-rich phase (cf. Figure 7(b)) g cp (r) is rather noisy. Even so, it is consistent with a low-density vapour-like radial distribution function for mixtures of spheres of diameter σ p and σ c , once again supporting the earlier interpretation that the colloids effectively 'see' the polymers as spherical particles.
The different forms of the polymer-polymer distribution function g pp (r) in the two phases are of particular interest. As expected, g pp (r) is essentially liquid-like in the colloid-rich phase. However one can see clear evidence of interpenetration of the polymer spheres. This is indicated firstly by the trend of g pp (r) as r approaches zero; for impenetrable spheres, g(r) would be zero for r < σ p , whereas g pp (r) remains non-zero within the polymeric core. The next indicator is that the first major peak does not appear exactly at σ p (as would be the case for fully impenetrable spheres) but, rather, at a slightly smaller value of r; the displacement of the first peak alone reflects that the interpenetration of the soft polymer spheres occurs (by its nature) at small inter-polymer separations. All of the peaks in g pp (r) are broader (and noisier) than the corresponding peaks in g cc (r) or g cp (r); this is consistent with a softer intermolecular interaction between the polymers than that between polymers and colloids, and that between colloids. In the polymer-rich phase, g pp (r) has a form somewhat characteristic of a 'vapour' state. Once again, there is evidence of interpenetration of the spherical polymeric cores since, as in the colloid-rich phase, g pp (r) remains non-zero at values of r substantially lower than the 'contact distance", σ p = 2R g .
The radial distribution function can be related to the potential of mean force, ω, via the expression  Figure 7 for the colloid-polymer system with a size ratio of σ c /σ m = 20 and polymer length m p = 100 at the reduced pressure P * = 10. The potentials relating to the colloid-rich phase are given in (a), while those relating to the polymer-rich phase are given in (b). The continuous green curves represent the colloid-colloid ω cc (r), the black curves represent the colloid-polymer ω cp (r), and the red curves represent the polymer-polymer ω pp (r). Note that the statistics are too poor to extract ω cc (r) for the polymer-rich phase. ω(r) = −k B T ln g(r). Accordingly, one can extract information relating to the average interactions between the particles from the appropriate radial distribution functions; indeed, in the low-density limit ω tends to the effective pair potential between particles. Since our simulations are not carried out at low densities, a rigorous approach such as that followed by Bolhuis et al. [74] should, strictly speaking, be used to extract the effective pair potentials from g(r) in our simulations. Bolhuis et al. solved the Ornstein-Zernicke [130] equation with the hypernetted-chain closure [131] using a self-consistent iterative procedure to obtain effective CoM-CoM pair potentials for their lattice-model polymers. In our current work, we seek only to examine the general features of the interactions, so an examination of the potential of mean force is sufficient.
In Figure 8 we display the potentials of mean force relating to the colloid-colloid, colloid-polymer, and polymerpolymer interactions extracted from the corresponding g(r) plots displayed in Figure 7. The colloid-colloid ω cc (r) within the colloid-rich region is characterised by a deep trough at σ c : this indicates the expected (depletion-induced) effective attractive interaction between the colloids. The statistics are too poor to extract a meaningful ω cc (r) effective interaction for the polymer-rich phase. Although less pronounced, there is also a clear well in the colloidpolymer ω cp (r) in the colloid-rich phase, indicating an ef-fective attraction between the polymers and the colloids. Perhaps most surprising is that there is also a deep trough in the polymer-polymer ω pp (r) in the colloid-rich phase, similarly indicating that an effective attraction exists between the polymer molecules in this phase. It appears from these observations that the depletion interaction is not restricted to the colloids; the presence of polymer molecules not only results in a depletion interaction between the colloids but, in concert, the presence of the colloids appears to induce a depletion interaction between the polymers.
By contrast, in Figure 8(b), there is no clear evidence of an effective attractive interaction between the molecules, although one can perhaps discern a small well in ω cp (r) at r just greater than σ cp . The differences between the potentials of mean force displayed in (a) and (b) clearly indicate that both the colloid-polymer and polymer-polymer interactions in the colloid-rich phase differ significantly from their respective counterparts in the polymer-rich phase. This is a particularly interesting conclusion in relation to the effective polymer-polymer interaction, considering the apparent similarity in the dimensions of the polymers in the two phases (cf. Figure 6). Such a difference in the effective interactions between the polymer chains in the two phases is not generally taken into account in theoretical treatments of colloid-polymer systems. It is interesting to compare our ω pp (r) with the effective CoM intermolecular potentials obtained by Bolhuis et al. [74] in their lattice simulations of self-avoiding-walk polymers. To facilitate comparison with reference [74], and to highlight the difference between the polymer-polymer interactions in the two phases, ω pp (r) for each phase is plotted again in Figure 9, though here normalised relative to the radius of gyration of the polymer in the corresponding phase. Also given in the figure are two of the effective CoM-CoM potentials obtained by Bolhuis et al. [74] from their lattice simulations. One corresponds to the potential between isolated chains of 100 segments, and another to chains of 500 segments in the semi-dilute regime at a density corresponding approximately to the effective packing fraction of polymers in the colloid-rich phase of our current simulations. For consistency with these simulations, one might expect our results to lie approximately between these two curves. In the polymer-rich phase ω pp (r) ∼ 2k B T at r ∼ 0, dropping smoothly to zero at r ∼ 2 to 3 R g ; this is indeed consistent with the effective CoM-CoM potentials of reference [74]; interestingly, along much of its length our curve almost superimposes with that corresponding to the isolated 100-segment chains. However, the behaviour of ω pp (r) in the colloid-rich phase is in stark contrast to the curves obtained by Bolhuis et al. [74].
One can further probe the polymer behaviour by looking at the site-site (colloid-colloid, colloid-monomer, and monomer-monomer) radial distribution functions. The site-site radial distribution functions obtained from our MD simulations of the colloid-rich and polymer-rich phases are presented in Figure 10; note that since the colloids are modelled as single spheres, the colloid-colloid functions represented here are identical to the molecular g cc (r) of Figure  7. In addition, for comparison, distribution functions relating to simpler HS systems are also included in the figures: g HS (r) of a pure (one-component) system of diameter σ c at the same packing fraction as the mixture; binary HS mixture of large and small spheres corresponding to the colloidmonomer g cm,HS mix (r) and colloid-colloid g cc,HS mix (r) cases. The HS mixture mimics the colloid-polymer mixture except that the monomer HS segments comprising the polymer are not connected (i.e., the mixture comprises simply large HS of diameter σ c and small HS of diameter σ m ); it thus has (by construction) the same density and (segment) composition as the phase of the actual mixture.
It is evident from Figure 10(a) that the presence of the monomers, whether unconnected or joined into polymer chains, has almost no influence on the structure of the colloids in the colloid-rich phase, since g cc (r) for this phase coincides almost exactly with both g cc,HS mix (r) and g HS (r). Although this is consistent with the earlier conclusion that the polymers occupy interstitial regions between the colloids, it is still a surprising result, since it appears to indicate that any depletion interaction experienced by the colloids is, at most, weak. It is interesting, in particular, to draw attention to the especially close match of g cc (r) and g cc,HS mix (r). Hitherto, we have given no consideration to the so-called 'depletion stabilisation' phenomenon: under certain conditions, the depletion layer can give rise to an effective repulsive interaction between the colloids. Shvets and Semenov [92] have suggested that the nomenclature 'depletion stabilisation' is misleading, preferring the term 'long-range free polymer induced' (FPI) repulsion. These authors have shown that while one expects only the classical polymer-mediated depletion attractive interactions when the polymers are in the dilute (swollen) regime, an FPI repulsion may be expected for colloid surface-surface separations in excess of 10ξ , when ξ 0.5R g . Recalling our earlier estimate of ξ ∼ 3σ m to 4σ m and that R g 6σ m , one might therefore expect to see some evidence of an FPI repulsion in our simulations at colloid separations in excess of 30σ m . On close inspection of g cc (r) one does indeed see that the first trough in the distribution, at ∼30σ m , is indeed shifted marginally to longer r in comparison to g HS (r); this is reflected in the potential of mean force (cf. Figure 8), with a peak at r slightly in excess of 30σ m , and one could g (r ) +1 colloid-monomer (colloid-polymer mix) g (r ) +1 colloid-sphere (binary HS mix) g (r ) pure HS system at same density g (r ) colloid-colloid (colloid-polymer mix) g (r ) colloid-colloid (binary HS mix) (a) distribution function g (r ) pure HS system at same density g (r ) +1 colloid-sphere (binary HS mix) g (r ) +1 colloid-monomer (colloid-polymer mix) g (r ) colloid-colloid (colloid-polymer mix) g (r ) colloid-colloid (binary HS mix) (b) Figure 10. Site-site pair radial distribution functions in (a) the colloid-rich and (b) the polymer-rich phases for the colloid-polymer system with a size ratio of σ c /σ m = 20 and polymer length m p = 100 at the reduced pressure P * = 10, compared with those for equivalent HS systems at the same density. Broken curves represent binary HS mixtures of the same density and (segment) composition as the phase (i.e., comprising the same numbers of HS colloids and (unconnected) HS polymer monomer segments; see text for details). The continuous curves represent the colloid-polymer systems; the green curves indicate colloid-colloid g cc (r) (reproduced from Figure 7) and the violet curves indicate colloid-monomer g cm (r). For clarity, the curves representing g cm (r) are shifted up by unity. The filled grey circles indicate g HS (r) of a pure (one-component) system of HSs of diameter σ c at the same packing fraction, η total , of the phase.
easily be misled into interpreting these observations as evidence of an FPI repulsion interaction. However, the almost exact coincidence of g cc (r) and g cc,HS mix (r) demonstrates that this cannot be the case; the origin of the FPI repulsion is a chain-end effect [92] whereas, from this comparison, one can quite clearly see that there is no impact of any sort of chain connectivity on the structure of the colloidal fluid.
Although the structural influence on the colloids is negligible, that on the monomers is not. This is shown by the disagreement at short range between g cm (r) and g cm,HS mix (r), indicating that the connectivity of the segments constrains the manner in which they order in the vicinity of the colloids; evidently the connected segments do not fill so closely the space between the colloids as unconnected HS of the same diameter. A similar short-range effect is manifest relating to these two functions in the polymer-rich phase, illustrated in (b). Unfortunately, as already discussed, due to the necessarily small size of the system there are insufficient colloids present in the polymer-rich phase to provide good statistics, so that the distribution functions in (b) are rather noisy. There is a pronounced difference between g cc (r) and the colloid-colloid g cc,HS mix (r) of the binary HS mixture, while neither coincides closely with g HS (r) from the pure-HS simulation. Accordingly, it is tempting to conclude not only that the presence of the monomers influences the structure of the colloids in this phase, but also that the connectivity of the monomers into chains is important. This would be consistent with the earlier conclusion that, effectively, the polymer behaves as spheres that are almost impenetrable by the colloids. However, one cannot rule out the possibility that the differences are due to no more than noise.
As the final component of this analysis, we examine the site-site distribution functions relating to the monomers comprising the polymer chains; monomer-monomer distribution functions for the colloid-rich and polymer-rich phases are presented in Figure 11. The features in these curves are dominated by the connectivity of monomers belonging to the same chain. This gives rise to a δ function at r = σ m , which is manifest in the figures as an extremely sharp peak (the height and width of which are determined by the bin size used in the accumulation of the statistics).
A second peak appears at r = 2σ m , and there is a further discontinuity in the gradient at r = 3σ m , after which the function decays smoothly. These features are clearly discernible in Figure 11(a), in which g mm, internal (r) is displayed.
Only monomers belonging to the same chain are counted in the accumulation of g mm, internal (r); accordingly the function decays to zero (rather than unity), at r ∼ σ p (σ pol p ∼ 14.5σ m for the polymer-rich phase; σ col p ∼ 13.8σ m for the colloidrich phase). One can discern a small but clear difference in the structure of the chains in the two phases at short r, as is highlighted by the inset in the figure; g mm, internal (r) is systematically larger in the colloid-rich phase. However, the internal structure of the chain for r 10σ m appears to be almost identical regardless of the phase.
In Figure 11(b) the overall g mm (r) is also presented; this function takes into account all of the monomer segments, regardless of whether or not they belong to the same chain. The features of these curves are still dominated by chain connectivity; they exhibit qualitatively similar characteristics to those in (a). However, one can discern from (b) that g mm (r) in the polymer-rich phase remains above unity well in excess of r = σ p (see the inset graph in (b)); this indicates local 'clustering' of polymer molecules in the polymer-rich phase, as one might intuitively expect given that the polymer chains interpenetrate. In addition, the plots in (b) reveal further information concerning the interpenetration of the polymer molecules. The difference seen between the curves in the polymer-rich and colloid-rich phases is partly due to the difference in the densities of polymers in the two phases.
However, the considerable separation of the two curves at short r (i.e., r < R g ) demonstrates that there is considerably more interpenetration of the chains in the polymer-rich phase than in the colloid-rich phase. This is to be expected, since the density of polymers in the polymer-rich phase is much higher than that in the colloid-rich phase.
That there is, nevertheless, interpenetration of chains in the colloid-rich phase is confirmed by comparing g mm (r) in (b) with its counterpart g mm, internal (r) in (a); for convenience, the tail of the latter is reproduced in the inset as the black dashed curve. The difference between the two curves at r σ p is due to monomers from other chains penetrating the chain to which the reference monomer segment belongs. Interpenetration of chains also explains why g mm (r) reaches its asymptotic value of unity well before σ p .
Altogether, the information in Figure 11 reveals that the internal structure of the polymer chains in the two phases is quite similar (thereby providing some enlightenment into the success of simple theoretical approaches, such as that of Asakura and Oosawa [29]) but that there are important differences which are manifest in the different degrees of interpenetration of the polymers in the two phases; the presence of the colloid molecules in the colloidrich phase influences the penetrability of the polymers by each other.  Figure 11. Monomer-monomer radial distribution functions of the polymers in both phases for the colloid-polymer system with a size ratio of σ c /σ m = 20 and polymer length m p = 100 at the reduced pressure P * = 10: (a) internal monomer distribution g mm, internal (r) for segments on the same chain; and (b) overall monomer distribution g mm (r). The green curves represent results for the colloid-rich phase while the red curves correspond to those of the polymer-rich phase. The dashed black curve in the inset of (b) represents the internal monomer distribution for the colloid-rich phase (i.e., the green curve from (a)); comparison of this curve with the overall distribution reveals the extent of interpenetration of the polymer chains in this phase (see the text for a discussion). The separation of the green and red curves in (b) indicates that the polymer chains in the polymer-rich phase interpenetrate far more than those in the colloid-rich phase.

Discussion and conclusions
Using MD simulation we have successfully demonstrated fluid-phase separation in a model colloid-polymer system comprising HS polymer chains (formed from 100 tangentially bonded HS segments) and HS colloids with a diameter 20 times that of the polymer segments. The system is found to exhibit coexistence between a colloidrich phase and a polymer-rich phase. Under the conditions of our simulations, the radius of gyration of the polymers is ∼6 monomer diameters, so that the size ratio q = R g /R c ∼ 0.6, whereby our system lies at the so-called colloid limit (q < 1). The densities of the two phases are in accordance with predictions made using Wertheim's TPT1; this agreement between the simulation and TPT1 provides mutual corroboration of each approach. A useful outcome of the theoretical predictions is an identification of the compositions of the two coexisting phases, particularly that of the low-density polymer-rich phase. Although this information is also available from the simulations, the statistical uncertainty associated with the simulated value of the composition of the polymer-rich phase is inevitably high due to the small number of colloids present.
While TPT1 provides no insight at the molecular level, the molecular simulations allow one to probe the structure of the fluid mixture and, thereby, develop insight into the effective interactions between the species in the mixture. Our detailed examination of the two coexisting phases has been especially revealing. Various indicators, such as the molecular (CoM) colloid-polymer distribution function, g cp (r), suggest that as expected, from the point of view of the colloids, the polymers behave (on average) as soft, impenetrable spheres of radius σ cp ∼ (σ p + σ c )/2. The polymers are interpenetrable by each other, as can clearly be seen from both the polymer-polymer and monomer-monomer distribution functions, g pp (r) and g mm (r). The latter further reveal that the interpenetration is greater in the polymer-rich phase; this is intuitively sensible since the density of polymers in this phase is higher. Comparing the colloid-colloid distribution function (g cc (r)) with the analogous function (g HS (r)) accumulated in a simulation of pure HSs (colloids) surprisingly reveals that the structure of the fluid with respect to the colloids is virtually unaffected by the presence of the polymers, indicating that any depletion interaction between the colloids due to the presence of the polymers is rather weak. However, the presence of the colloids does affect the structure of the polymer molecules. Based on an analysis of the polymer dimensions, the polymers in the colloid-rich phase behave as though the polymers are at an effective packing fraction of η p ∼ 0.075, although the actual polymer density in this phase is much lower. This effective density is higher than the density of the polymer-rich phase. The interpenetration of polymer molecules would be expected to increase with density in systems of pure polymer whereas, as already noted, the interpenetration of polymers in the colloid-rich phase is demonstrated by analysis of g mm (r) to be lower than that in the polymer-rich phase. The increase in effective density, allied to a decrease in interpenetration, demonstrates that the structure of the polymers is affected by the colloids. Denton and co-workers have performed simulation studies of crowding of polymer coils in mixtures with nanoparticles in the colloid [132] and the protein [91] limits. These authors observed that polymer coils, modelled as effective penetrable ellipsoids, both contracted and (in the protein limit) became more spherical when crowded by impenetrable, spherical nanoparticles. The analysis in our current study precludes any assessment of possible change in the relative sphericity of the polymer chains brought about by the presence of the colloids; however, the high effective density of the polymers we observe in the colloid-rich phase is consistent with the compaction of the polymer coils observed by Denton and co-workers.
Radial distribution functions are closely related to potentials of mean force and, thereby, provide information relating to the effective interactions between molecules in the low-density limit. The interactions between polymers in the polymer-rich phase are seen to be consistent with the effective potentials obtained in previous work [74]; however, the analogous interactions in the colloid-rich phase turn out to be entirely different, featuring a region of effective attraction between the polymer chains that is absent in the polymer-rich phase. This could, perhaps, be evidence of a depletion interaction between the polymers; such an interpretation is strengthened by an analysis of the interaction between polymers in a pure-polymer simulation at the effective packing fraction η p ∼ 0.075, where no evidence of attraction is seen. Based on our other observations, the difference between the polymer-polymer interaction in the two (colloid-rich and polymer-rich) phases may be due to the manner in which the colloids affect the interpenetration of the polymer molecules.
It is widely regarded that the physics of the colloid limit is now well understood [56,[89][90][91] (although various authors have recently highlighted some shortcomings in the current understanding [35,92,94]). However, in theoretical treatments of colloid-polymer systems, it is usual to assume that the effective interaction between polymers is the same at the different thermodynamic states, even when going beyond a simple AO treatment. It is clear from our study that this assumption is not necessarily valid, although the importance of this is not clear since a good representation of colloid-polymer systems have been obtained without taking account of the observed difference between the polymer-polymer interactions in colloid-rich and polymer-rich environments. For example, in simulations using their coarse-grained potential [74] to represent the polymer-polymer interaction, Bolhuis et al. [77] obtained excellent agreement with the experimental volumefraction (φ c vs. φ p ) phase diagram of a colloid-polymer system with a size ratio of σ p /σ c = 0.67. Although their potential is (weakly) density dependent, it was obtained from earlier bulk simulations of self-avoiding-walk polymers with no colloids present and, thereby, cannot take into account any influence of the colloids on the effective polymer-polymer interaction such as that uncovered in our current work. Accordingly, the incorporation of a more representative polymer-polymer interaction into a theoretical approach, allowing for the difference between colloid-rich and polymer-rich phases, would be very interesting.  Figure A1. The effect of the sampling-sphere radius, r samp , on the probability distribution of the total packing fraction, P(η total ), obtained from 50,000 samples taken from (a) the initially demixed and (b) the initially mixed states of the colloid-polymer system with a size ratio of σ c /σ m = 20 and polymer length m p = 100 at the reduced pressure of P * = 10. The sharp peak at low density, corresponding to the polymer-rich phase, is well defined for all values of r samp . The shape of the broader peak at higher density, corresponding to the colloid-rich phase, is sensitive to the size of the sampling sphere, becoming sharper at higher values of r samp , although the position of the maximum is largely unaffected.  Figure A2. The sampling-sphere radius, r samp , does not affect the comparison of the initially demixed (continuous curves) and initially mixed (dashed curves) simulations, in terms of the calculated probability distributions of the total packing fraction, P(η total ), which is here reaffirmed to be essentially the same for the two systems, irrespective of the value of r samp . The conclusion that both simulations have converged to configurations with the same distribution of density is thereby shown not to rely on the chosen value of r samp .