Eur. Phys. J.

E 4, 363370 (2001)

THE EUROPEAN PHYSICAL JOURNAL E

EDP Sciences c Societ` Italiana di Fisica a Springer-Verlag 2001

Eect of colloidal charge discretization in the primitive model

R. Messinaa , C. Holm, and K. Kremer
Max-Planck-Institut fr Polymerforschung, Ackermannweg 10, 55128, Mainz, Germany u Received 30 June 2000 and Received in nal form 28 November 2000 Abstract. The eect of xed discrete colloidal charges in the primitive model is investigated for spherical macroions. Instead of considering a central bare charge, as it is traditionally done, we distribute discrete charges randomly on the sphere. We use molecular dynamics simulations to study this eect on various properties such as overcharging, counterion distribution and diusion. In the vicinity of the colloid surface the electrostatic potential may considerably dier from the one obtained with a central charge. In the strong Coulomb coupling, we showed that the colloidal charge discretization qualitatively inuences the counterion distribution and leads to a strong colloidal charge-counterion pair association. However, we found that charge inversion still persists even if strong pair association is observed. PACS. 82.70.Dd Disperse systems: Colloids 61.20.Qg Structure of associated liquids: electrolytes, molten salts, etc. 41.20.-q Applied classical electromagnetism

1 Introduction
The electrostatic interactions in charged colloidal systems play a crucial role in determining the physical properties of such materials [1,2]. The behavior of these systems is extremely complex due to the long-range Coulomb interactions. A rst simplifying assumption is to treat the solvent as a dielectric medium solely characterized by its relative permittivity r . A second widely used approximation consists in modeling the short-range ion-ion excluded-volume interaction by hard spheres. These two approximations are the basis of the so-called primitive model of electrolyte solutions. The system under consideration is an asymmetrical polyelectrolyte made up of highly charged macroions and small counterions in solution. A further simplication can be achieved by partitioning the system into subvolumes (cells), each containing one macroion together with its neutralizing counterions plus, if present, additional salt. This approximation has been called accordingly the cell model [3,4]. The cells assume the symmetry of the macroion, here spherical, and are electrostatically decoupled. In this way one has reduced a complicated many-body problem to an eective one-colloid problem. For spherical macroions the structural charge is normally modeled by a central charge, which, by Gauss theorem, is equivalent to considering a uniform surface charge density as far as the electric eld outside the sphere is concerned. Most analytical work as well as simulation approaches rely on the above assumptions. It is well known that in the strong Coulomb coupling regime ion-ion correlations become very important, and signicant deviations from mean-eld approaches are expected. One of the eects
a

which a mean-eld theory like the Poisson-Boltzmann one cannot explain is the phenomenon of overcharge, also called charge inversion. It consists of binding excess counterions to a charged particle (macroion) so that its net charge changes sign. This has recently attracted signicant attention [514]. It may give rise to a possible mechanism for strong long-range attraction between like-sign charged colloids [12,13]. The purpose of this paper is to investigate if such a phenomenon (overcharge) depends on the way the structural charge is represented. The macroion is taken to be perfectly spherical, i.e. we neglect any surface roughness [15]. We introduce discrete charges on the macroion sphere instead of a central charge, and compare the results to those obtained with a central charge. We concentrate on the following properties in the strong Coulomb coupling: overcharging, counterion distribution and surface diusion.

2 Simulation model
2.1 Macroion charge discretization The macroion charge discretization is produced by using Nm identical microions of diameter , all identical to the counterions, distributed randomly on the surface of the macroion. Then the structural charge is Q = Zm e = Zc Nm e, where Zm > 0, Zc is the counterion valency and e is the positive elementary charge. The discrete colloidal charges (DCC) are xed on the surface of the spherical macroion. In spherical coordinates the elementary surface

e-mail: messina@mpip-mainz.mpg.de

364

The European Physical Journal E

Ey

DCC

r0 r a

Ex

the dissipation-uctuation theorem Wi (t) Wj (t ) = 6m kB T ij (t t ). For the ground state simulations the uctuation force is set to zero. Excluded-volume interactions are taken into account with a pure repulsive Lennard-Jones potential given by 12 6 4 rr0 + , for rr0 < rcut , rr0 (3) ULJ (r) = 0, for r r0 rcut , where r0 = 0 for the microion-microion interaction (the microion can be a counterion or a DCC), r0 = 7 for the macroion-counterion interaction, and rcut (= 21/6 ) is the cut-o radius. This leads to a macroion-counterion distance of closest approach a = 8. Energy and length units in our simulations are dened as =kB T0 (with T0 = 298 K), and = 3.57 , respectively. A The pair electrostatic interaction between any pair ij, where i and j denote either a DCC or a counterion, reads UCoul (r) = kB T0 lB Zi Zj , r (4)

rj (r)

Fig. 1. Schematic view of the setup: the discrete colloidal charges (DCC) of diameter are in dark grey. The radial electrostatic eld components Ex and Ey are represented. For a detailed meaning of the other symbols, see text. Note that this is a a two-dimensional representation of the three-dimensional system.

is given by
2 2 dA = r0 sin dd = r0 d(cos )d ,

(1)

and, to produce a random discrete charge distribution on the surface, we generated randomly the variables cos and . Only congurations leading to an overlap of microions are rejected. Figure 1 shows a schematic view of the setup. Note that in real physical systems like sulfonated latex spheres, no large heterogeneities are expected in the charge distribution, provided that the colloid surface is relatively regular, therefore our choice is justied. Nevertheless, the experimental situation is more complicated since other phenomena such as surface chemical reactions [16], hydration, roughness [15] and many more may be present. Here, we restrict ourselves to a simple model in order to understand the eect of macroion charge discretization, and leave the other questions for future investigations. 2.2 Molecular dynamics procedure We use molecular dynamics (MD) simulations to compute the motion of the counterions coupled to a heat bath acting through a weak stochastic force W(t). The equation of motion of counterion i reads m d2 ri dri + Wi (t) , = i U m dt2 dt (2)

where lB = e2 /4 0 r kB T0 is the Bjerrum length describing the electrostatic strength. Being essentially interested in the strong Coulomb coupling regime, we choose the relative permittivity r = 16 (lB = 10), divalent counterions (Zc = 2) and divalent DCC for the remaining of this paper. The macroion and the counterions are conned in a spherical impenetrable cell of radius R. The macroion is held xed and is located at the center of the cell. The 3 colloid volume fraction fm is dened as rm /R3 , where rm = a /2 is the colloid radius. We have xed R = 40 so that fm = 6.6 103 . To avoid image charge complications, the permittivity r is supposed to be identical within the whole cell (including the macroion) as well as outside the cell.

3 Macroion electric eld

The rst step to understand the eect of colloidal charge discretization consists of estimating the electric eld, or equivalently, the electrostatic potential generated by such a sphere in the absence of counterions. A simple graphical visualization of the eld lines is here not possible, since there is no perfect symmetry. Indeed, in the present case the electric eld becomes very anisotropic and irregular close to the sphere, which is the most interesting region where correlations are expected to be large. To describe qualitatively the eect of charge discretization on the electrostatic potential, we compute for three perpendicular directions x, y, z the resulting radial potentials Vx (r), Vy (r), Vz (r) for one given DCC random distribution as a function of the distance r a from the macroion center. The radial component of the electric eld Ei (r) = r Vi (r) has the important feature of representing the attractive

where m is the counterion mass, U is the potential force having two contributions: the Coulomb interaction and the excluded-volume interaction and is the friction coecient. Friction and stochastic force are linked by

R. Messina et al.: Eect of colloidal charge discretization in the primitive model

Zm = 50

365

(a)

component towards the sphere. The normalized radial potential Vi in the i -th direction at a distance r from the colloid center is given by
Nm 2 Vi (r) = kB T0 lB Zc j=1

-50

V(r)/(k B T0 )

1 , |rj (r)|

(5)

-100

monopole Vx Vy Vz 1 2 3

-150

r/a
Zm = 90

(b)

-100

-200

monopole Vx Vy Vz 1 2 3

-300

r/a
Zm = 180

(c)
-100

-200

-300 monopole Vx Vy Vz 1 2 3

-400

-500

r/a
Fig. 2. Radial electrostatic potential as a function of the macroion center distance r produced by the xed microscopic colloidal charges disposed on the sphere. These potentials have been measured in three perpendicular directions (x, y, z) (see Fig. 1). The isotropic case corresponds to the eld obtained with a central charge (monopole). Three structural charges are considered: (a) Zm = 50, (b) Zm = 90, and (c) Zm = 180.

where rj (r) is the vector pointing from the microion j to the point where the electric potential is computed (see Fig. 1). Physically, V (r) is the electrostatic potential interaction between a counterion and all the surface microscopic colloid charges. The monopole contribution is merely given by Vmono (r) = kB T0 lB Zm Zc . In Figure 2 we r show the electric potential for three typical bare charges, each corresponding to one given random macroion charge distribution. For all cases, one notes that in the vicinity of the surface the potential becomes very dierent from the one computed with a central charge. We carefully checked that similar results were obtained for other choices of x, y, z directions (by rotating the trihedron (ex , ey , ez )). However, if we observe the electric eld suciently far away from the colloidal surface (about one macroion diameter), the eld is almost exactly the same as the one produced by a central charge, which we term isotropic for the rest of this paper. A closer look at Figure 2 reveals that by increasing the bare charge Zm the electric eld starts to become isotropic at smaller distances from the spheres surface. This last feature can be physically easily interpreted. In fact when one increases the bare charge, one also increases the absolute number of discrete charges which has the eect of approaching the uniform continuous charge density limit (corresponding to the isotropic case). To capture the discretization eect on the surface electrostatic potential, we have measured the electrostatic potential along a circle of radius a concentric to the spherical macroion (see Fig. 1). We start from the top of a given DCC microion and move along a circle in a random direction and measure the electrostatic potential V (s) as a function of the arc length separation s from the starting point. The same formula as equation (5) has been used here. The constant monopole contribution is merely given by Vmono = kB T0 lB Zm Zc . Results are reported in a Figure 3 for the same congurations as before. It clearly shows that the electrostatic potential is strongly uctuating. More specically, the higher the structural charge Zm , the larger the oscillation frequency of the potential uctuations over the surface. This feature can be explained in terms of holes. In the very vicinity of a given DCC the potential is increased (in absolute value) in average, and around a given DCC there is a hole (depletion of charges) which tends to decrease the potential (in absolute value). The average surface of this hole is increasing with decreasing bare charge Zm (i.e. decreasing density of charged sites). In the following sections we are going to study the eect of charge discretization on the counterions distribution in the strong Coulomb coupling. For all following results we used the same random charge distributions which gave the results of Figures 2 and 3.

V(r)/(k B T0 )

V(r)/(k B T0 )

366
Z m = 50

The European Physical Journal E

80

4 Ground state analysis

(a)

-100

V(s)/(k B T0 )

-120

-140

-160

-180

10

20

30

40

50

60

s/
-180 Z m = 90

(b)
-200

In this section, we focus on counterion distribution exclusively governed by energy minimization, i.e. T = 0 K. In such a case correlations are maximal, and all the counterions lie on the surface of the spherical macroion. To avoid the trapping in metastable states, we systematically heat and cool (10 cycles) the system and retain the lowest energy state obtained in this way. Furthermore, we choose as the starting conguration the one where each DCC is exactly associated with one counterion, and each of these dipoles is radially oriented (each dipole vector and the macroion center lie on the same line). Preliminary, we checked that this method reproduces well the ground state energy and structure in simple situations where a central charge with two, three, four or ve counterions is present. The structure of these systems is well known by the Gillespie rules [17]. It turns out that in these situations no rough energy landscape (even for Zm = 180 and 90 counterions) appears and therefore the MD simulation easily nds the global minimum. It is only in the case of DCC that several energy minima are observed.

V(s)/(kB T0 )

-220

4.1 Neutral case First we consider the simple salt-free case where the system [macroion + counterions] is neutral. In order to characterize the counterion layer, we compute the counterion correlation function (denoted by CCF) g(r) on the surface of the sphere, dened as 2 g(r) = s
0 10 20 30 40 50 60
i=j

-240

-260

-280

(r ri )(r rj ),

(6)

s/
-400 -420 Z m = 180

(c)

V(s)/(kB T0 )

-440 -460 -480 -500 -520

where s = Nc /4a2 is the surface counterion concentration (Nc = Zm /Zc being the number of counterions), r corresponds to the arc length on the sphere. Note that at zero temperature all equilibrium congurations are identical, thus only one is required to obtain the CCF. Similarly, one can also dene a surface macroion correlation function (MCF) for the microions on the surface of the macroion. The CCF is normalized as follows: s
0 a

2rg(r)dr = (Nc 1) .

(7)

10

20

30

40

50

60

s/
Fig. 3. Surface electrostatic potential as a function of the arc length s along a circle of radius a concentric to the the macroion for three dierent trajectories. The monopole contribution is represented by the dashed line. The same congurations as those of Figure 2(a-c) have been used.

Because of the nite size and the topology of the sphere, g(r) has a cut-o at a (= 25.1). Therefore at large distance the correlation function diers from the one obtained with an innite planar object. The CCF and MCF for two dierent structural charges Zm (50 and 180) can be inspected in Figure 4. The CCF is computed for a system with a central charge (CC) and for the discrete colloid charges (DCC) case. One remarks that both CCF dier considerably following the nature of the colloidal charge, i.e., discrete or central (see Fig. 4). For the isotropic case (central charge) a Wigner Crystal structure is observed as was already found in references [12,13, 18]. It turns out that when we have to deal with DCC the

R. Messina et al.: Eect of colloidal charge discretization in the primitive model

Zm = 50

367

(a)
MCF (DCC) CCF (CC) CCF

g(r)
2 1

(a)
0 0 10

r/
Zm = 180

20

30

(b)

MCF (DCC) CCF ( CC) CCF

4

g(r)
2

(b)
0

10

r/

20

30

Fig. 5. Ground state structures for two values, (a) Zm = 50 and (b) Zm = 180, corresponding to the two cases of Figure 4. The colloidal surface microions are in white, and the counterions in blue. Full ionic pairing association occurs.

Fig. 4. Ground state surface correlation functions for two macroion bare charges: (a) Zm = 50 and (b) Zm = 180. The two counterion correlation functions (CCF) are obtained for discrete colloidal charges (DCC) and for the central charge (CC). To get the same distance range for CCF and the colloidal surface discrete microions correlation function (MCF), the MCF curve x-axis (r/) was rescaled by a factor a/r0 (compare the setup of Fig. 1).

4.2 Overcharge We now investigate the overcharge phenomenon. The starting conguration corresponds to neutral ground states as were previously obtained. The spirit of this study is very similar to the one undertaken in reference [12]. To produce overcharge, one adds successively overcharging counterions (OC) in the vicinity of the colloidal surface. Thus the resulting system is no longer neutral. By using Wigner crystal theory [6,19], we showed that the gain in electrostatic energy (compared to the neutral state) by overcharging a single uniformly charged colloid can be written in the following way [12,13,18]:
OC En = E cor + E mon = 3 3n 2 (n 1)n n Nc + , + (kB T0 )lB Zc 2 8Nc 2a

counterion distribution is strongly dictated by colloidal charge distribution (see Fig. 4). Ground state structures are depicted in Figure 5. It clearly shows the ionic pairing, between DCC and counterions. Also, it appears natural to call such a structure a pinned conguration. However, one can expect that the structure might become less pinned if the typical intra-dipole distance (here ) and the typical mean inter-dipole distance become of the same order. This is a case which is not discussed in the present paper. It would correspond to extremely highly charged colloids that are rarely encountered in nature. Nevertheless, we checked that even for Zm = 360 the structure is still pinned, where the average inter-dipole distance is about 2.

(8)

where E cor , which is equal to the rst term of the right member, denotes the gain in energy due to ionic correlations for n OC. The functional form of this term can be derived from the WC theory [12,13,18]. The second term on the right-hand side, E mon , is the monopole repulsion, which sets in when the system is overcharged (with

368

The European Physical Journal E

Z m = 50 - DCC Z m = 50 - CC Z m = 90 - DCC Z m = 90 - CC Z m = 180 - DCC Z m = 180 - CC Z m = 360 - DCC Z m = 360 - CC

200 100 0 -100 -200 -300 0

obtained with DCC approaches the one from the continuous case as expected. Common features of overcharging between isotropic and discrete systems are briey given here. We note that the maximal (critical) acceptance of OC (4, 6 and 8 for a central charge and 2, 4 and 6 for DCC) increases with the macroionic charge Zm (50, 90 and 180, respectively). Furthermore, for a given number of OC, the gain in energy is always increasing with Zm . Also, for a given macroionic charge, the gain in energy between two successive overcharged states is decreasing with the number of OC. Note that at T = 0, the value r acts only as a prefactor. All these features are captured by equation (8).

E/(kBT0 )

12

16

20

Overcharging counterions
Fig. 6. Electrostatic energy (in units of kB T0 ) for ground state congurations of a single charged macroion as a function of the number of overcharging counterions for three dierent bare charges Zm . CC stands for the central charge case. The neutral case was chosen as the potential energy origin. The dashed lines are produced by using equation (8).

5 Finite temperature
In this part, the system is brought to room temperature T0 . We are interested in determining the counterions distribution as well as the counterion motion within the counterion layer. The radius R of cell is again xed to 40 so that the macroion volume fraction fm has the nite value 6.6103 . Under these conditions the system is still highly energy-dominated so that at equilibrium all counterions lie on the surface of the macroion (strong condensation).

n > 1). This term will, for sucient high number n of OC, stop the process of overcharging. As before Nc = Zm /Zc is the number of counterions in the neutral state, and is a constant which was determined by using the measured OC value of E1 of our simulations. The total electrostatic energy of the system as a function of the number of OC is displayed Figure 6 for four bare charges Zm (50, 90, 180 and 360). These energy curves corresponding to discrete systems were produced by averaging over ve random DCC realizations. Again, the overcharging process is aected by the charge discretization and pinning, but it is still energetically favorable. The main eects of charge discretization are: i) the reduction of gain of energy and ii) the reduction of maximal (critical) acceptance of OC. Both points can be explained in terms of ion-dipole interaction. It is exactly this attractive ion-dipole correlation which is responsible of charge inversion for colloidal systems with discrete charges. When the rst OC is present, it is normally located in between the pinning centers, and will essentially interact with its neighboring dipoles (DCC-counterion). This interaction increases with decreasing OC-dipole separation, i.e. increasing the colloid bare charge Zm . This explains why the energy gained increases with Zm (see Fig. 6). On the other hand, the repulsion between the counterions is not fully minimized since they do not adopt the ideal Wigner crystal structure that is obtained with a central charge which in turn explains i). For a higher degree of overcharge, one has to take into account a repulsive monopolar contribution identical to E mon appearing in equation (8). Again, since for DCC structures counterions are not perfectly ordered, the attractive correlational energy is smaller (in absolute value) than the one obtained with a central charge, which in turn explains ii). Note that for very high bare charge (Zm = 360) the overcharge curve

5.1 Counterions distribution Like in the ground state analysis, we characterize the counterion distribution via its surface correlation function. At non-zero temperature, correlation functions are computed by averaging i=j (r ri )(r rj ) over 1000 independent equilibrium congurations which are statistically uncorrelated. Results are depicted in Figure 7. For both bare charges Zm (50 and 180) considered the counterions distributions are aected by charge discretization. The effect of temperature is to smooth the CCF. As expected, for the central charge case, the counterion positional order is much weaker at room temperature than in the ground state case.

5.2 Surface diusion The aim of this section is to answer the following question: do the counterions only oscillate around their equilibrium (ground state) position or do they have also a global translational motion over the sphere? To answer to this question one introduces the following quantity: x2 (t, t0 ) = 1 t t0
t t0

dt [x(t ) x(t0 )]2 ,

(9)

which is referred to as the mean-square displacement (MSD), where x(t0 ) represents the position of a given counterion at time t = t0 and x(t, t0 ) is its position at later time t. The root-mean-square displacement (RMSD)

R. Messina et al.: Eect of colloidal charge discretization in the primitive model

Zm = 50 Zm = 50

369

0.7
MCF (DCC) CCF (CC) CCF

(a)
3

(a)
0.6

RMSD/()
0 10 20 30

g(r)

0.5

0.4

0.3

r/
2 Zm = 180 5

10

15

20

25

counterion
Zm = 180

(b)
1.5

MCF (DCC) CCF (CC) CCF

(b)
4

RMSD/()
0 10 20 30

g(r)

0.5 1 0

r/
Fig. 7. Surface correlation functions at room temperature. The two CCF are obtained for discrete colloidal charges (DCC) and for the central charge (CC). (a) Zm = 50, (b) Zm = 180.

20

40

60

80

100

counterion
Fig. 8. Root mean square displacement (RMSD) for each counterion; (a) Zm = 50, (b) Zm = 180.

is dened as x(t, t0 ) = x2 (t, t0 ). Like for the surface correlation function, the MSD is measured on the spherical surface (arc length) and it has a natural upper limit (a)2 . Results for the discrete case are depicted in Figure 8 for two macroion bare charges Zm (50 and 180), where each counterion RMSD is given. In both cases, the RMSD is smaller than the typical mean inter-dipole separation, which is approximatively (s )1/2 6 for Zm = 50 and (s )1/2 3 for Zm = 180. This suggests that the motion of the counterion is purely oscillatory around its DCC pinning center. Figure 8 also shows that pinning is stronger for Zm = 50 than for Zm = 180. This is in agreement with our previous statement, where we point out that the inter-dipole distance has to be comparable (or smaller) to (than) the intra-dipole distance in order to reduce pinning eect. Thus the DCC sites do produce a noticeable energy well. One can get convinced on this point, by evaluating the electrostatic binding energy of an ionic pair Epin which is
2 Epin = kB T0 lB Zc / = 40kB T0 .

However, for much higher temperatures a liquid-like behavior is to be expected. Also, of course, the strength of the pinning can be lowered by dierent parameters: larger ions, higher dielectric constant r , monovalent ions as is captured by equation (10). For the isotropic case, we have checked that counterions have a large lateral motion and can move all over the sphere. This is obvious since in this situation there are no pinning centers.

6 Concluding remarks
We have carried out MD simulations within the framework of the primitive model to elucidate the eect of colloidal charge discretization. The main result of our study is that, in the strong Coulomb coupling, the charge inversion is still eective when the macroion structural charge is carried by discrete charges distributed randomly over its surface area. We have shown that the intrinsic electrostatic potential solely due to the surface colloidal microions varies strongly from point to point on the macroion

(10)

370

The European Physical Journal E 5. V. Perel, B. Shklovskii, Physica A 274, 446 (1999). 6. B. Shklovskii, Phys. Rev. E 60, 5802 (1999). 7. E.M. Mateescu, C. Jeppesen, P. Pincus, Europhys. Lett. 46, 493 (1999). 8. J.F. Joanny, Eur. Phys. J. B 9, 117 (1999). 9. E. Gurovitch, P. Sens, Phys. Rev. Lett. 82, 339 (1999). 10. M. Lozada-Cassou, E. Gonzlez-Tovar, W. Olivares, Phys. a Rev. E 60, R17 (1999). 11. M. Deserno, C. Holm, S. May, Macromolecules 33, 199 (2000). 12. R. Messina, C. Holm, K. Kremer, Phys. Rev. Lett. 85, 872 (2000). 13. R. Messina, C. Holm, K. Kremer, Europhys. Lett. 51, 461 (2000). 14. T.T. Nguyen, A.Y. Grosberg, B.I. Shklovskii, J. Chem. Phys. 113, 1110 (2000). 15. S. Bhattacharjee, C.H. Ko, M. Elimelech, Langmuir 14, 3365 (1998). 16. O. Spalla, L. Belloni, J. Chem. Phys. 95, 7689 (1991). 17. A theory that predicts molecular geometries using the notion that valence electron pairs occupy sites around a central atom in such a way as to minimize electron-pair repulsion. See, for example, D.W. Oxtoby, H.P. Gillis, N.H. Nachtrieb, Principles of Modern Chemistry (Saunders College Publishing, 1999), Chapt. 3, p. 80. 18. R. Messina, C. Holm, K. Kremer, submitted. 19. L. Bonsall, A.A. Maradudin, Phys. Rev. B 15, 1959 (1977).

sphere. When counterions are present, it leads to a pinned structure where every counterion is associated with one colloidal charge site. Furthermore, we observed a pure phonon-like behavior (counterion vibration with small lateral motion) found at room temperature. Future work will address the problem of valency asymmetry, that is when the colloidal charges are represented by monovalent counterions and the counterions are multivalent. This is a non-trivial situation since ionic pairing may be frustrated. Also, the case of the weak Coulomb coupling regime should be addressed.
We thank B. Shklovskii for helpful and constructive remarks. This work is supported by Laboratoires Europens Associs e e (LEA).

References
1. J. Israelachvili, Intermolecular and Surface Forces (Academic, London, 1992). 2. D.F. Evans, H. Wennerstrm, The Colloidal Domain o (Wiley-VCH, New York, 1999). 3. T.L. Hill, Statistical Mechanics (Addison-Wesley, Reading, Mass., 1960). 4. H. Wennerstrm, B. Jnsson, P. Linse, J. Chem. Phys. 76, o o 4665 (1982).