Statistical Mechanical Approach to Phase Transitions in Colloids BY WILLIAM VAN MEGEN AND IAN SNOOK Department of Applied Physics, Royal Melbourne Institute of Technology, Melbourne, Victoria, Australia Received 6th December, 1977 The Monte Carlo method and a number of approximate theories, based on classical statistical mechanics and using DLVO pair potentials, are used to calculate the structure and the osmotic pressure of electrostatically stabilized dispersions of spherical particles. Reasonable agreement is found between the exact Monte Carlo radial distribution functions and the corresponding quantity determined from the hard sphere model. This hard sphere model provides a simple unifying theory for the structure of these colloidal solutions provided that the equivalent hard sphere diameter is chosen correctly.From the good agreement between experiment and these and previous calculations, we conclude that the DLVO pair potentials are adequate to calculate at least the structure of colloidal solutions. 1. INTRODUCTION Many works have illustrated the usefulness of and the necessity for using statistical mechanics to determine the properties of colloidal solutions [see, for example, ref. (1)-(6)]. These ideas have more recently been applied to a variety of dispersions of spherical and non-spherical particles to calculate, in particular, the conditions for the separation of a phase of regularly arranged particle^.^*"^ In all cases, such calcula- tions have reproduced at least the qualitative features of real colloidal systems.It is the aim of this discussion to examine more closely the structure of colloidal solutions and its dependence on volume fraction and temperature and electrolyte concentration of the suspending medium. Considerable emphasis has been placed on very dilute systems of small particles for which the structure is quantitatively accessible experi- mentally. 9*10 Since the application of statistical mechanics to a many-particle system presupposes the knowledge of the potential energy of the system due to the interactions among the particles," we briefly examine the assumptions on which the usual interaction poten- tials for colloidal particles (uiz. those formulated in the DLVO theory)12 are based. Statistical mechanics explicitly applied to a system with a large number of particles is exact for given interaction potentials among the particles.Hence, any discrepancy between the results of exact statistical mechanical calculations and experiments on corresponding systems can only arise from inadequacies in the formulation of inter- action potentials. To be specific we restrict the following discussion to an electrostatically stabilized dispersion of identical solid spheres in an aqueous electrolyte. In the usual formula- tion the total interaction energy, U, between two colloidal particles of separation rtJ, is regarded as the sum of the van der Waals attraction, U,, plus the double layer repulsion, V, : l2 W , j > = U*(P'ij> + U&>. (1.1)W. VAN MEGEN AND I . SNOOK 93 In the absence of external fields and assuming pairwise additivity, the total excess1 potential energy is given by @ = 2 U(rij).i < j Since the double layer repulsion dominates in most of the situations discussed below, it may be worthwhile examining the implication of the assumption expressed by eqn (1. I) and (1.2) for this contribution. The double layer repulsion is essentially coulomb repulsion, due to the net charge on the particles, screened by the intervening electrolyte. The magnitude of this screening depends on the charge and distribution of the ions in the background medium; to determine this exactly in a colloidal solution is an immensely complex problem. The assumption of pairwise additivity simplifies this problem considerably by assuming that the distribution of ions in the vicinity of an isolated pair of charged colloidal particles remains unperturbed by the approach of a third particle. This assumption may appear somewhat drastic particularly in moderately concentrated dispersions in which the double layers overlap.It must be emphasized, however, that many properties, such as the occurrence of a disorder to order (D-0) transition, are not critically dependent on the inclusion of many-particle potentials nor on the precise details of the pair potential? This statement is clearly substantiated by the excellent overall agreement between experiment and statistical mechanical calculations, based on simple phenomenological pair potentials, for both molecular systems l1 and colloidal s o l ~ t i o n s . ~ ~ * ~ ~ Calculations using simple asymp- totic forms for both the van der Waals attraction and the double layer repulsion have yielded structure factors l3 and osmotic pressures l4 in quantitative agreement with corresponding experimental values for monodisperse polystyrene latexes.Further, since the pressure is very much more sensitive to the precise form of the interaction potential, the reasonable agreement between the calculated and measured osmotic pressures in concentrated dispersions illustrates that the effect of three or more particle potentials is not large under the conditions studied so far, and that the pair potentials formulated in the DLVO theory should certainly be adequate for determination of the structure. The structure of a many particle system is quantitatively given by the radial distri- bution function, which is of further interest since it provides a vehicle for relating the pair potential to excess equilibrium properties.ll Analogously, the time-dependent distribution function, or the van Hove space-time correlation function provides the key to the dynamic ~r0perties.l~ Both distribution functions are accessible by auto- correlation of the intensity of laser light scattered from suspensions.Examination of the radial distribution function tells us whether the suspension is ordered or dis- ordered. It is this ordering or regular stacking of the particles which gives rise to the familiar iridescence observed in concentrated polystyrene latexes. In the next section we outline the essential statistical mechanics and the main ideas underlying the numerical procedure (Monte Carlo method) for evaluating the multidimensional integrals.The one drawback of this direct procedure is the enormous amount of computing time required to locate precisely a phase transition (the D-0 transition in this case). It is, therefore, useful to use a combination of a number of approximate methods, such as the cell model and perturbation theory, to locate the D-0 transition more e~pediently.~ The simplest approximate approach is the hard-sphere model which exploits the fact that the structure in dense systems, and in particular the occurrence of a D-0 transition, is essentially a packing phenomenon arising from the effective size of the particles, i.e., the actual size of the particles plus an allowance for the strongly repulsive part of the pair potential.8*16 These approxi- mate methods are also briefly reviewed in section 2.In section 3 results are presented94 PHASE TRANSITIONS I N COLLOIDS for dispersions in which the double layer thicknesses are small and large, respectively, compared with the partick radius. 2. THEORETICAL In the usual formulation of classical equilibrium statistical mechanics the radial distribution function, g(r), and the pressure, P, for an isotropic system of spherical particles are given by and in which N is the number of particles in volume V, k is the Boltzmann constant, T is the absolute temperature and N(r, r + 6r) is the number of particles within a distance r, r + 6r from a given particle.ll The pair potential U(rij) is given by eqn (1 .l) with12 UR(rlj) = 2neoeryi ln(1 + exp(-icauij), lca $ 1 (2.4) in which uij = (rij/a) - 2 is the surface to surface separation in units of the particle radius a, K is given by IC = (ST)% and all other symbols have their usual meaning.12 The angular brackets in eqn (2.1) and (2.2) represent averages over the canonical (or N, V, T ) ensemble; in the Monte Carlo (MC) method of Metropolis et a1.l’ these are explicitly evaluated by a numerical technique. For practical reasons the real system containing a very large number of particles is replaced by a periodic system.Each cell of this periodic system contains a small number (up to 108 for the present computations) of particles which are allowed to interact among themselves and also with their periodic images. This restriction does not introduce any significant error provided the range of the correlation between the particles is smaller than the length of the unit cell.Considerable difficulty arises in the vicinity of a phase transition (such as the D-0 transition) where the small number of particles is inadequate to examine two phases in equilibrium.ll It is still possible, however, to locate the volume fraction at which the disordered phase begins to order (qD) and that at which the ordered phase disorders (q0) by using an indirect approach discussed, for example, by Hoover and Reel8 to locate the melting transition in molecular solids. In this approach one constructs a reversible thermodynamic (P,q) curve connecting the ordered and disordered phases by using a single occupancy model to localize the particles and thus stabilize the ordered phase at low volume fractions.The excess free energies are determined by integrating the pressure curves of the disordered and artificially ordered systems. Subsequent equating of the Gibbs free energy and pressure for each phase locates pD and qo.W . VAN MEGEN AND I . SNOOK 95 Since the above procedure is tedious and requires comparatively large amounts of computer time it is extremely useful to develop simpler, albeit approximate, methods for evaluating the properties of many-particle systems. Unlike the MC method, these approximate techniques are only likely to be applicable under certain restricted conditions. Two very successful approximate theories are the hard-sphere perturba- tion theory (for disordered or fluid-like phases)ll and the cell model (for ordered or crystal-like phases).7* l9 Briefly, in perturbation theory7*11 the pair potential is divided into a steep, repul- sive reference part U, and a longer ranged, weaker perturbation, Up, i.e., u = u, + up.(2.7) A given thermodynamic property can then be written as a perturbation series, ex- panded about the property of the reference system, i.e., the system of particles with the pair potential U,. For example, the Helmholtz free energy A and the radial distribu- tion function, g(r), are ----+J+>+...., A - A0 A A NkT NkT NkT NkT and It is, of course, particularly useful to relate A. and go(r) to a system of hard spheres since the properties of this system are already known very accurately.For a given system (pair potential) the central problem is then to carry out suitably the division expressed in eqn (2.7) and choose an appropriate equivalent hard sphere to replace U,. This in general depends on a number of factors such as the range and softness of the pair potential, the density and the temperature.ll One method for dividing up the pair potential is given in the Barker-Henderson (BH) theory11*20 in which one writes g(r) = go@) + gl09 + * * * * * (2.9) and u, = u = o up = 0 = u r S p r > p r L P (2.10) r > p and the break point 11 is chosen so that U, is a harsh repulsive (hard-sphere like) potential and Up is a weaker slowly varying potential. The effective hard sphere diameter is given by d = /om [exp(- U,/kT) - lldr. (2.1 1) With this criterion both go and g, are usually required to obtain a good approximation for g(r). The Chandler-Weeks-Anderson (CWA) theory 11s21 presents a second method in which the pair potential is broken up as follows : and u,=u+c Up = -& = o = u (2.12) where U(p) = --e and p and d are chosen to produce a hard sphere diameter so that go(r) is as close as possible to g(r).This is achieved by making the difference between A . and Ahs zero to first order. Both the above theories give essentially equivalent results for simple liquids .I1 9 2o 2196 PHASE TRANSITIONS IN COLLOIDS In the cell model one replaces the array of particles by a regular array of potential cells which constrain the motion of each parti~1e.l~ This leads to the following expression for the canonical partition function 2; (2.13) where Eo is the static lattice energy, vf is the free volume available to each particle and all other symbols have their usual meaning. In the “ sphericallized ” or “ smeared out ” approximation, vf becomes vf = 4n r2 exp( - [ ~ ( r ) - y(O)]/kT} dr, (2.14) where C is the effective cell radius and y(r) is the particle potential energy when displaced a distance r from the cell centre.6 3. RESULTS AND DISCUSSIONS To illustrate the quantitative difference between ordered and disordered colloidal solutions, fig. 1 shows g(r) as determined by the MC method for a dispersion of spheres of radius a = 5.95 x loA7 m in a 1-1 electrolyte of concentration c = 0.1 mol m-3. For this system the pair potential is comparatively hard and short ranged [i.e., K a B 1 and we use U, given in eqn (2.4)], as opposed to systems to be discussed r l r FIG.1 .-Radial distribution function g ( r ) plotted against r in units of the particle diameter CT showing MC (-) and hard sphere (- - - -) results; (a) rp = 15%, (b) v, = 25%, (c) v, = 35%. later, for which the pair potential is soft and long ranged (rca < 1). For q of 15% and 25%, g(r) is qualitatively similar to that of a simple liquid, indicating a disordered arrangement of particles, whereas for q = 35%, g(r) indicates a long ranged crystal- like behaviour. Fig. 1 also shows the hard sphere perturbation results obtained using CWA theory (the BH results being almost identical). The hard sphere diameter is roughly at the point where U = 1 kT. Fig.2 shows Pa3/kT plotted against q, as determined by first order perturbation theory for q < q,, and the cell model for q 2 qo, for various concentrations, c, of the background electrolyte. The accuracy of these methods, by comparison with MC results, has been established in earlier work.’ Note that as c decreases so does the volume fraction pD at which the ordered phase begins to separate from the disorderedW. VAN MEGEN AND I. SNOOK 97 phase, as does the difference between qo and pD. Both these predictions are in accord with experimental observation.22 It is also interesting to note the maximum in the osmotic pressure at the D-0 transition as c is varied from 1 to mol m-3. As discussed in previous work, considerable difficulty arose in using the cell model and first order perturbation theory to locate the D-0 transition for systems in which the attractive tail in the pair potential is significant.14 In the case of c = 1 mol m-3 the depth of the secondary minimum # 1 % FIG.2.-Osmotic pressure Pa3/kT plotted against q, around the D-0 transition, for various electro- lyte concentrations; (a) c = 1 mol mW3, (b) c = 0.5 rnol mW3, (c) c = 0.1 rnol m-3, ( d ) c = 0.05 mol M-~, ( e ) c = 0.01 mol m-3. MC result indicated by *. in U is about 1 kT rendering an error in qD and qo of M 1 % in q and possibly an over- estimate of one (SI) unit in Pa3/kT. For the system under examination here, the primary maximum in the pair potential is of the order of lo3 kT and, for interparticle spacings beyond this, the effect of the van der Waals attraction is negligible compared with the double layer repulsion.At least for c < 0.5 mol m-3 the secondary minimum is less than ~ 0 . 1 kT. Thus, as c is decreased below 0.5 mol m-3 U, becomes not only longer ranged, accounting for the reduction of qD and qo, but also less steep or softer. Since the pressure is essentially the ensemble average of the gradient of the pair potential averaged over all pairs of particles in the system, the increase in the osmotic pressure at the D-0 transition can be attributed to the steepening of the pair potential with increasing electrolyte concentration. However, for c > 0.5 mol m-3 the appearance of a comparatively small secondary minimum (1 kT at c = 1 mol m-3) in the pair potential causes a sudden drop in the pressure at the phase transition.That this reduction in the osmotic pressure stems from the secondary minimum is easily confirmed by comparing the pressures, at the estimated qD (or qo), determined by the hard sphere method,8 which allows only for the positive part of the potential, and first order perturbation theory, which approximately allows for the negative tail in the pair potentials; these pressures are 2.7 x and 1.4 x N m-2, res- pectively. Interestingly, a similar effect has recently been observed by Takano and ha chis^^^ in an experimental determination of qD and qo for polystyrene latex in several concentrated electrolytes. This experiment then presents further experimental verification of the secondary minimum in electrostatically stabilized colloids.Also included in fig. 2 are the results of a complete MC determination of qo and qD for c = 0.1 mol m-3 using the approach of Hoover and Ree.18 Clearly the approximate methods, although not exact, provide a reasonable and expedient estimate of yo and qD. We turn now to systems for which the pair potential is very soft and long ranged98 PHASE TRANSITIONS I N COLLOIDS [ K a < 1, and we use eqn (2.5) for U,]. This occurs in dispersions of very small particles in dilute electrolytes. These systems are particularly attractive since g(r) can be determined experimentally by means of laser light scattering. In more concentrated dispersions of larger particles such experiments are plagued by multiple scattering effects. The effect of temperature on the structure of a system of spherical particles, a = 4.36 x m, in a monovalent electrolyte, c = 2.9 x low3 mol m-3, was recently carried out by Schaefer.1° In fig.3 we show the MC determination of 3 I FIG. 3.-Radial distribution functiong(r) plotted against r for a = 4.36 x lo-* m and c = 2.9 x 10 niol m-3, showing MC (-----) and hard sphere (- - - -) results, Q = 2a. g(r) and the BH result (go + g,) for this system. It should be noted that in this case U(d) N 7 kT and d is strongly dependent on y. Fig. 4 compares the calculated structure factor, S(K), with that measured by SchaeferlO at a temperature of 46 "C. The structure factor is related to g(r) by S(K) = 1 + - - l ] r sinKrdr, where 4nn 2. K = - sin (012) 2 1 :: 0 4 08 12 16 2 0 2 A Kb FIG. 4.-Structure factor S(K) as a function of the magnitude of the scattering vector, K, in units of the reciprocal diameter, showing MC (-) and experimental (.) results.is the magnitude of the scattering vector, 8 is the scattering angle, il is the wavelength of the scattered radiation and n is the refractive index of the medium. The agreemint is quite reasonable considering the simplicity of the interaction potential used and the expected errors in such measurement^.^ Schaefer also measured the effect of tem- perature on S(K) and found that this quantity indicates virtually no structure at 70 "C and a very pronounced structure at 40 "C (with freezing occurring at 39.5 "C). In direct contrast, our MC computations between 40 and 70 "C indicate no significant variation in S(K). The calculations assume the surface charge to remain constant,W.VAN MEGEN AND I . SNOOK 99 in which case the surface potential yo increases with increasing temperature. The resulting increase in the repulsion compensates the reduction in structure due to the Boltzmann term. The observed disappearance of the structure with temperature is readily accounted for by assuming some degradation of the ion exchange resin with the consequential increase in concentration of background electrolyte. The dramatic effect of this is shown in fig. 5 where the g(r) values are reduced to dilute gas-like shapes n L U br 2 6 0 10 12 14 r / 0, FIG. 5.-Radial distribution functions for a = 4.36 x lo-' m and c = mol m-3 (-----) and c = 2.9 x mol m-3 (- - - -). when the electrolyte concentration is increased to and 2.9 x mol m-3.Indeed, Schaefer mentions that his ion exchange resin degrades irreversibly above 60 "C. We are at a loss, however, to account for Schaefer's observation of the large increase in the magnitude of the first peak in S(K) upon decreasing the temperature below 46 "C. The only explanation we can offer is that the discrepancy may be due to experimental errors or the inadequacy in the formulation of the potential energy (or both). 4. CONCLUSIONS From the preceding discussion it is evident that the structure of monodisperse systems of spherical particles is adequately explained by statistical mechanical calcula- tions based on pair potentials given in the DLVO theory. It has also been shown that g(r) is reasonably predicted by a hard sphere approach.Considerable care must be exercised in selecting the correct hard-sphere diameter particularly for very long ranged potentials. The results are clearly very sensitive to variations in the hard sphere diameter; arbitrarily setting this at a point where the pair potential has decayed to E 1 kTis by no means sufficient. It is essential that all the parameters characteriz- ing the system are accurately known. The sensitivity of the calculated results to the surface potential, yo, has been discussed in earlier worki3 whilst the drastic variation of g(r) with electrolyte concentration has been illustrated above. The above results also indicate that the structure of colloids is insensitive to temperature variations, provided that the particle surface charge and the electrolyte concentration do not vary with temperature.L. Onsager, Ann. N.Y. Acad. Sci., 1949, 51, 627. I. Snook and W. van Megen, Colloid and Interface Sci., ed. M. Kerker (Academic Press, N.Y., 1976), vol. IVY p. 1. R. P. Keavey and P. Richmond, J.C.S. Furaday 11, 1976,72, 773. S. Marcelja, D. J. Mitchell and B. W. Ninham, Chem. Phys. Letters, 1976, 43, 353. W. G. M. Agterof, J. A. J. van Zomeren and A. Vrij, Chem. Phys. Letters, 1976, 43, 363. S. L. Brenner, J. Phys. Chem., 1976, 80, 1473. W. van Megen and I. Snook, J. Colloid Inlerface Sci., 1976, 57, 40, 47.100 PHASE TRANSITIONS I N COLLOIDS W. van Megen and I. Snook, Chem. Phys. Letters, 1975, 35, 399. J. C. Brown, P. N. Pusey, J. W. Goodwin and R. H. Ottewill, J. Phys. A , 1975, 8, 664. lo D. W. Schaefer, J. Chem. Phys., 1977,66, 3980. l1 J. A. Barker and D. Henderson, Rev. Mod. Phys., 1976,48, 587. l2 E. J. W. Verwey and J. Th. G. Overbeek, Theory of stability of Lyophobic Colloids (Elsevier, l3 W. van Megen and I. Snook, J. Chem. Phys., 1977,66, 813. l4 A. Homola, I. Snook and W. van Megen, J. Colloid Interface Sci., 1977, 61, 493. l5 S. H. Chen, Physical Chemistry, An Advanced Treatise, ed. H. Eyring, D. Henderson and W. l6 P. Jena and W. R. Smith, Chem. Phys. Letters, 1973, 21, 295. l7 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, J. Chem. l8 W. G. Hoover and F. H. Ree, J . Chem. Phys., 1967, 47,4873. l9 J. A. Barker, Lattice Theories of the Liquid State (Pergamon, Oxford, 1963). 2o J. A. Barker and D. Henderson, J. Chem. Phys., 1967, 47, 2856,4714; Accounts Chem. Res., 21 H. C. Andersen, D, Chandler and J. D. Weeks, Phys. Rev. A , 1971, 4, 1597; J. Chem. Phys., 22 S. Hachisu and Y . Kobayashi, J. Collaid Interface Sci., 1973, 42, 342. 23 K. Takano and S . Hachisu, J. Phys. SOC. Japan, 1977, 42, 1775. Amsterdam, 1948). Jost (Academic Press, N.Y., 1971), vol. VIIIA. Phys., 1953,21, 1087 1971,4, 303; Phys. Rev. A , 1971,4, 806. 1972,56, 3812.