Prediction of Ordered and Disordered States in Colloidal Dispersions BY IAN SNOOK* AND WILLIAM VAN MEGEN Department of Applied Physics, Royal Melbourne Institute of Technology, Melbourne, Victoria, Australia Received 28th May, 1975 Excess pressures and radial distribution functions of a monodisperse colloid are calculated from a statistical mechanical model using the Monte Carlo method to evaluate the configurational integrals. Particular regard is paid to the effects of variation of volume fraction and background electrolyte concentration on the structure of the system. Numerous experimental studies have been reported 1-4 on the structure of colloidal dispersions of spherical particles. These studies indicate the existence of regions of order (solid-like regular arrangement of particles) and regions of disorder (liquid-like, showing no long ranged order).In view of this experimental activity we have attempted a description of this behaviour in terms of the forces between the colloidal particles. Recently we have developed a model for colloidal dispersions based on classical statistical mechanics which explains the main features of the compression of, for example, an aqueous suspension of polystyrene spheres. 9 To further demonstrate the versatility and usefulness of the application of models based on classical statistical mechanics, this paper presents the results of calculations on spherical colloidal particles dispersed in aqueous electrolyte. Systematic variation of the background electrolyte concentration and volume fraction, with all other conditions held constant, predicts some interesting structural effects similar to recent experimental observation.An outline is given of the method and interparticle potentials used and the results are presented and discussed. THEORY As outlined in earlier work 59 we calculate the equilibrium properties of an electrostatically stabilized colloidal dispersion using classical statistical mechanics, the Monte Carlo method being used to evaluate numerically the multidimensional integrals in the canonical ensemble. The Monte Carlo method is adequately de- scribed elsewhere 7-9 and the details and problems associated with its application to colloidal dispersions are discussed in ref.(6). In particular we evaluate the internal energy E, the pressure P and the radial distribution function g(rI2), given by ; E 31 -=-NkT 2+NkT-<w PV 1 -= l--(YN}NkT 3NkT 216 I. SNOOK AND W. VAN MEGEN and where p = N/V is the number density and 1exp(-ON/kT)dr, . . . dr,N! d2)(r,,r,) = -(N-2)! eXp(-@N/kT) dr, . . . dv, (4) is the two-particle distribution function. Vis the volume, T the absolute temperature, k Boltzmann's constant, N the number of particles and DN and YNare the total potential energy and the virial term, respectively. Furthermore, the weighted aver- ages in eqn (1) and (2) are given by and where Uijis the pair potential and rijthe interparticle distance. The interparticle potential Uijused consists of the usual lo double layer repulsion URplus the van der Waals attraction UA,which for spherical particles of radius a are given by; uR(rij)= 271caI,biIn(1 +e-'''j) (7) and In the above expressions, z = k-a, uij = rij/a-2, I,bo is the surface potential, E the dielectric constant of the background, A the Hamaker constant, and l/k-the usual characteristic thickness of the double layer, given by where No is Avagadro's number, IZ the electrolyte concentration, e the electronic charge and z the valence of the ions in the background medium (assuming a sym-metrical electrolyte).Since, in this work, we are only attempting to produce the main features of the experimental results, the use of such simple pair potentials is fully justified.Further-more, their use substantially reduces computer tinie requirements. CALCULATIONS AND RESULTS Unfortunately, in the available experimental papers insufficient data is given pertaining to the exact Eature of the dispersions used. Hence, in our calculations we have chosen the following parameters which appear to be representative : a = 5.95 x loA7m, t+bo = 0.06 V, A = 2.5 x J, z = 1, T = 300 K and E = 80 E~ (E~is the permittivity of free space). Equilibrium properties have been calculated for the above system with background electrolyte concentrations of 1, 5 x lo-', lo-', 5 x and mol m-3. 21 8 PREDICTION OF COLLOID STABILITY rlo FIG. 1.-The interparticle potential U/kTas a function of the particle separation in units of (J = 2a for n = mol m-3 (-), n = 5 x mol m-3 (---), n = 10-1 mol m-3 (-.-.-), n = 5 x 10-1 mol m-3 (-..-..-) and y1 = 1 mol m-3 (..--..).5 4 0 A Q I'I-* 4 bQ-2 A 0 -.-IA I ElI -.-I El h -IC 70 60 50 40 30 20 10 0 4 FIG.2.-The logarithm of the reduced pressure, P* = PV/NkT,as a function of the volume fraction mol m-3, (-t-) y1 = 10-1 mol m-3, (m) n = 5 x 10-1n = mol m-3, (A)n = 5 x4, for ; (0) (:#: andm-3mol ) n = 1 mol m-3. I. SNOOK AND W. VAN MEGEN Fig. 1 shows the total potential energy as given by eqn (7) and (8) for the five electrolyte concentrations. Note the increase in the repulsive range of the potential as the electrolyte concentration is reduced. I 1.2 1.4 1.6 rlu Frc.3.-The radial distribution function g(r) for n = 1 molm-3 with ---, 4 = 5 o/o; --I-, ,$ = 15 %; ....., 4= 25y--.--.-, ,$ = 35 -..-..-, 4= 50 %. UY 0, r.1. FIG. 4.-The radial distribution function g(r) for /I = lo-' mol nr3 with ----, ,$ = 57;;---,4= 15 ;<; ....., 4 = 25 %; -.-.-, 4 = 35 x. In most cases the initial configuration was taken as the usual face centred cubic structure 6* at the number density pof the system NN p=-=-v L3' PREDICTION OF COLLOID STABILITY where L is the length of the cell side, this being related to the usual percentage volume fraction 4 by As in our earlier work, the initial few hundred thousand configurations were rejected to allow the system to equilibrate. Then, at least 700 000 further configura- tions were used to obtain the reported averages.rlu FIG.5.-The radial distribution function g(r),as given by eqn (3) and (4) as a function of the particle separation, for the lowest electrolyte concentration of n= mol m-3 with 4 = 5 % (-) and 4 = 15 %(----). 12-8-Nl 4-04 0 10 20 30 40 50 60 70 4 FIG.6.-The number of nearest neighbours N,, as given by eqn (12), plotted against 4with (x) nn= 1 mol m-3, (0)= 5 x lo-’ rnol m--3,(+) n= lo-‘ mol m-3, (A) n= 5 x lo-* rnol m-3 and n= mol m-3. (The dashed lines are to facilitate reading the results and are not meant (0) to indicate a continuous variation of Nl with d). I. SNOOK AND W. VAN MEGEN 221 Fig. 2 shows the resulting reduced pressure P* = PV/NkT as a function of 4 for the five values of the electrolyte concentration.In each case P*increases quite drama- tically with volume fraction up to a 4 of around 70 %. Beyond this point the system is no longer stable as is evidenced by a sharp drop in the calculated pressure. In fact, these latter pressures are negative due to the infinitely deep potential well [see eqn (7) and (S)] at very small surface to surface separations of the particles. Note that in the results for the lowest electrolyte concentrations the rate of pressure increase with 4 is reduced, from moderate to large values of 4. This effect can be rationalized by noting that the pressure is directly related to the pair virial rdU/dr and that for the lower electrolyte concentrations this virial is smaller.Fig. 3 to 5 show some typical radial distribution functions [given by eqn (3) and (4)] for ranges of 4 covering the transition from disorder to order. Note, for example in fig. 3, g(r) indicates an ordered “ solid-like ” structure for 4 = 50 %, whereas at lower volume fractions no such longer ranged structure is evident. For 4 > 50 %, of course, ordered systems were also obtained with the first peak becoming sharper (and higher) with increasing 4. 60-. *. 0. *. 0. 40 0 00 00 0 DO 002olI 0 00 00 ‘Ol 0-t -2 -I 0 log (4 FIG.7.-A summary of which systems were found to be ordered (0)and disordered (0). For all systems (4 = 5,15,25,35,50,60,70 % and n = 1,sx lO-l, 5 x mol m-3) the number of nearest neighbours Nl was determined, as shown in fig.6. This quantity is determined from the usual symmetrized expression l1 N1 = 8nprmaxr”g(r)dr 0 where r,,, is the position of the first peak in g(r). Note the shift to lower 4 in the transition from order to disorder as the electrolyte concentration is reduced. This is readily explained by the fact that a decrease in electrolyte concentration increases the double layer thickness, hence the effective particle diameter and effective volume fraction also increase. For a perfect hexagonal close packed system of spheres one expects N1 to be 12 as we obtained for higher volume fractions. However, for q5 just beyond the order/ PREDICTION OF COLLOID STABILITY disorder transition, calculated values for N, are slightly less indicating the structure obtained is not quite a perfect close packed solid.Furthermore, since the first peak in the radial distribution functions are not perfectly symmetrical about r,,,, the actual values of N, are approximate. Finally in fig. 7 results of ordered and disordered systems are summarized. The features of this graph seem to agree remarkably well with the experimental results of Hachisu et aZ.,4 except, of course, the region of coexistence of the ordered and dis- ordered phases is not obtained in our work. The reasons being that in our computer experiment the system is not subjected to a gravitational field. Furthermore, the small system of 32 particles as we have treated it, although quite adequate to evaluate bulk properties of systems of particles with short ranged interparticle potentials, is l2quite inadequate for a realistic treatment of two phases in eq~ilibrium.~* DISCUSSION The Monte Carlo method as it stands can, in principle, be used to study phase transitions (e.g.disorder to order in a colloidal dispersion) and even two phases in equilibri~im.~ However in practice, a very large number of particles is needed because of the density gradients and the long range correlations existing at the phase transition. Thus phase transitions are studied approximately in several ways ;for example by keeping the system in one phase beyond its region of stability (e.g. the single occupancy model of Ree et a1.)’ or by placing an external potential on the system.12 These calculations still require large numbers of particles and many more hours of computer time than the normal Monte Carlo method.In this study we merely find the region in which a phase transition occurs, the most stable phase being determined by starting the system initially in either phase and finding the phase in which the system finally settles (after around a million con- figurations). This approach is similar to the early studies of systems of particles with hard sphere square-well and Lennard-Jones 12-6 potential^,^^ and like these studies, our calculations yield only the region of the phase transition to within a range of 5 to 10 % in 4. However, we expect this is quite sufficient for an initial study just as it was in the case of molecular systems,13 particularly as there is little precise experi- mental data on the phase transition region anyway.It may be noted that the region of the phase transition occurs at lower and lower values of 4 as the background electrolyte concentration is decreased. This packing fraction is based on the size of the particle, not on the size of the particle with its attached double layer. The phase of the system was determined from the radial distribution function. For instance, from fig. 4 one can readily observe that the dispersion with an electrolyte concentration of 10-1 mol m-3 is disordered when 4 5 25 % and ordered when 4 2 35 %. This can be checked by calculating the number of nearest neighbours (NJ. Approximate calculations show N1 to be close to 12 for ordered systems and much less than this value for disordered systems.In all cases the initial configuration was taken as the face centred cubic configura- tion. However, attempts were also made to obtain the ordered structures of the n = mol m-3 and 4 = 15, 25 % systems from a random (disordered or liquid- like) initial configuration as generated for the corresponding $ but higher (n = 1 mol m-3) electrolyte concentration. Although ordered structures were obtained in these runs, the exact details of those runs starting with the face centred cubic structure could not be reproduced within a realistic number of configurations (about one million). This situation reminds one of the random hard sphere packing experiments of the last decade in which it was found impossible to obtain packing fractions or I. SNOOK AND W.VAN MEGEN volume fractions in excess of about 65 %.14 Thus it appears impossible to generate a perfect close packed solid randomly, or the probability of doing so is negligible. 9* l3In practical terms, one generates a glassy A very eloquent explanation of this effect has been given by Wood who states that at very high densities (and we are here discussing dispersions of high effective density) trajectories in phase space are located within a pocket around the point corresponding to the initial configuration and the probability of a trajectory to another pocket, around a point corresponding to the close packed hexagonal configuration, is extremely small.At all but the lowest electrolyte concentration, the onset of coagulation was indicated at volume fractions beyond 70 %. Again, the exact point of this phase transition is difficult to locate for reasons mentioned above. Furthermore, the large negative pressures obtained for these partially coagulated states are indeed question- able. This is due to the unrealistic form of the inter-particle potential ;approaching minus infinity as the particle surface to surface separation approaches zero. In reality this narrow infinite well is truncated by the Born repulsion of at least a mono-layer of ions absorbed on the surface of the particles. Since in our model we have not accounted for the structure of this inner region of the double layer, it would be unrealistic to accept any actual values generated by the Monte Carlo calculation for the partially coagulated state.For the lowest electrolyte concentration (n = mol m-3) coagulation is evident beyond 4 = 60 %. The long range of the repulsive potential for this case results in systems (at least at moderate to large 4) with a high free energy. This enhances the coagulation of the primary particles. The authors thank the Royal Melbourne Institute of Technology Computer Centre for the use of their facilities, Dr. R. 0. Watts for his valuable assistance with the Monte Carlo program and Dr. Andrew Homola for his fruitful comments. G. W. Brady and C. C. Gravatt, Jr., J. Chem. Phys., 1971, 55, 5095. L. Barclay, A.Harrington and R. H. Ottewill, Kolloid-2.2. Polyrnere, 1972, 250, 655. P. A. Hiltner and I. M. Krieger, J. Chem. Pliys., 1969, 73, 2386. S. Hachisu, Y.Kobayashi and A. Kose, J. Colloid Interface Sci., 1973, 42, 342. I. Snook and W. van Megen, Chem. Phys. Letters, 1975,33, 156. W. van Megen and I. Snook, J. Colloid Interface Sci., in press. 'W. W. Wood, Physics of Simple Liquids, ed. H. N. V. Temperley, J. S. Rowlinson and G. S. Rushbrooke (North-Holland, Amsterdam, 1968). R. 0. Watts, Rev. Pure Appl. Chem., 1971, 21, 167. F. H. Ree, Physical Chemistry, An Adzanced Treatise, ed. D. Henderson (Academic Press, New York, 1971), vol. VIIIA, chap. 3. lo H. R. Kruyt, Colloid Science (Elsevier, Amsterdam, 1952), vol. 1. l1 C. J. Pings, Physics of Simple Liquids, ed. H. N. V. Teniperley, J. S. Rowlinson and G. S. Rushbrooke (North Holland, Amsterdam, 1968). l2 J. K. Lee, J. A. Barker and G. M. Pound, J. Chem. Phys., 1974, 60, 1976. l3 W. B. Streett, H. J. Raveche and R. D. Mountain, J. Chem. Phys., 1974, 61, 1960. l4 J. D. Bernal and J. Mason, Nature, 1960, 118, 910. (PAPER 5/991)