Charge Renormalization, Thermodynamics, and Structure of Deionized Colloidal Suspensions
aa r X i v : . [ c ond - m a t . s o f t ] D ec Charge Renormalization, Thermodynamics, andStructure of Deionized Colloidal Suspensions
Ben Lu and Alan R. Denton ∗ Dept. of Physics, North Dakota State University, Fargo, ND, U.S.A. 58108-6050
Abstract.
In charge-stabilized colloidal suspensions, highly charged macroions, dressedby strongly correlated counterions, carry an effective charge that can be substantiallyreduced (renormalized) from the bare charge. Interactions between dressed macroionsare screened by weakly correlated counterions and salt ions. Thermodynamic andstructural properties of colloidal suspensions depend sensitively on the magnitudesof both the effective charge and the effective screening constant. Combining a chargerenormalization theory of effective electrostatic interactions with Monte Carlo simula-tions of a one-component model, we compute osmotic pressures and pair distributionfunctions of deionized colloidal suspensions. This computationally practical approachyields close agreement with corresponding results from large-scale simulations of theprimitive model up to modest electrostatic coupling strengths.
PACS (2006) : 82.70.Dd, 83.70.Hq, 05.20.Jj, 05.70.-a
Charge-stabilized colloidal suspensions [1–3] exhibit rich thermodynamic phase behav-ior and tunable materials properties (e.g., thermal, mechanical, optical, rheological) thatare the basis of many industrial and technological applications. Thermally excited (Brow-nian) motion of nm- µ m-sized particles dispersed in a fluid medium drives the self-assemblyof ordered phases. Important examples are nanoscale structures [4] and colloidal crys-tals, whose diverse crystalline symmetries and variable lattice constants can convenientlytemplate photonic band-gap materials [5–7].Interparticle interactions and correlations determine the distribution of microions(counterions and salt ions) around the colloidal macroions and thereby govern microion-mediated electrostatic interactions between macroions. As explained by the landmarktheory of Derjaguin, Landau, Verwey, and Overbeek (DLVO) [8, 9], repulsive interactionsbetween like-charged macroions can stabilize a suspension against aggregation inducedby van der Waals attractive interactions [10, 11]. Equilibrium and dynamical properties ∗ Corresponding author.
Email address: [email protected] of charged colloids depend sensitively on electrostatic interactions, which can be widelytuned by adjusting system parameters: size, charge, and concentration of ion species; pHand dielectric constant of solvent.A powerful approach to modeling charged colloids is molecular simulation of eitherthe primitive model, which includes all ions explicitly, or all-atom models, which includeeven individual solvent (e.g., water) molecules. At such microscopic resolution, simu-lations can potentially yield valuable insights into the thermodynamic, structural, anddynamical properties of real suspensions. Currently available processors and algorithmsare limited, however, to relatively small systems and low ion size and charge ratios. Fur-thermore, simulations of microscopically detailed models do not necessarily elucidatethe physical mechanisms underlying complex cooperative behavior.An alternative to brute-force modeling is simulation of a one-component model, de-rived by pre-averaging over the microion coordinates to obtain effective interactions be-tween macroions, screened by implicitly modeled microions. This coarse-grained strat-egy finesses the computational challenges that plague more explicit models, but relies onpractical and accurate approximations for the effective interactions. A previous simula-tion study [12] validated the one-component model — implemented with linear-responseand mean-field approximations — by direct comparison with pressure data from primi-tive model simulations at modest electrostatic couplings. The eventual breakdown of themodel at higher couplings was attributed to failure of linear-screening approximations toaccount for nonlinear counterion association near macroions.A recently proposed charge renormalization theory of effective interactions in chargedcolloids [13] addresses limitations of linear-screening approximations by incorporatingan effective (renormalized) macroion charge into the one-component model. Calcula-tions based on a variational approximation for the free energy indicated the potentialof the theory to accurately predict thermodynamic properties of suspensions well intothe nonlinear-screening parameter regime. The present paper describes complementaryMonte Carlo simulations designed to test predictions for both thermodynamics (osmoticpressure) and structure (pair distribution function) of deionized suspensions of highlycharged colloids.The remainder of the paper is organized as follows. Section 2 first reviews the prim-itive and effective one-component models of charged colloids. Section 3 outlines thecharge renormalization theory, which predicts renormalized system parameters: effec-tive macroion charge, volume fraction, and screening constant. Section 4 describes ourMonte Carlo simulations, which take as input the renormalized effective interactions.Section 5 compares our results for the pressure and pair distribution function of deion-ized suspensions with corresponding data from simulations of the primitive model andwith predictions of a variational free energy theory. Excellent agreement is obtained, withminimal computational effort, over ranges of macroion charge, volume fraction, and elec-trostatic coupling strength. Finally, Sec. 6 closes with a summary and conclusions.
The primitive model idealizes the colloidal macroions as charged hard spheres of monodis-perse radius a and bare valence Z (charge − Z e ), the microions as point charges, and thesolvent as a continuum of uniform dielectric constant ǫ . We limit present considerationto monovalent microions, for which microion-microion correlations are relatively weak,and further assume index-matching of macroions and solvent, justifying neglect of imagecharges and polarization effects. At absolute temperature T , a characteristic length scaleis the Bjerrum length λ B = e / ( ǫ k B T ) , defined as the distance between a microion pair atwhich the electrostatic interaction energy rivals the average thermal energy. In a closedvolume V , the suspension comprises N m macroions and N ± positive/negative microions,of which N c are counterions and N s are salt ion pairs. Global electroneutrality relates theion numbers via N c = N + − N − = Z N m . Alternatively, the suspension may be in Donnanequilibrium (e.g., across a semi-permeable membrane) with an electrolyte reservoir thatfixes the salt chemical potential.While the primitive model has the virtue of explicitly representing all ions, simu-lations face severe computational challenges posed by long-range Coulomb interparti-cle interactions and large size and charge asymmetries between ion species. Sophisti-cated numerical algorithms, such as cluster moves [14, 15] and Ewald or Lekner sum-mation [16, 17], can significantly broaden the range of accessible length and time scales.Nevertheless, currently feasible studies are still limited to relatively small systems andmodest asymmetries.An alternative approach to modeling charged colloids is derived from the Hamilto-nian H of the primitive model by formally tracing over the degrees of freedom of themicroions ( µ ) in the partition function of the multi-component mixture: Z = hh exp ( − H / k B T ) i µ i m = h exp ( − H eff / k B T ) i m , (2.1)leaving an explicit trace over only the macroion ( m ) degrees of freedom. The resultingone-component model (OCM) is governed by an effective Hamiltonian H eff = E + N m ∑ i = j = v eff ( | r i − r j | )+ ··· , (2.2)a function of the macroion coordinates r i , which comprises a one-body volume energy E ,an effective macroion-macroion pair potential v eff ( r ) , and higher-order terms involvingsums over effective many-body potentials. The volume energy, although independentof macroion coordinates, depends on the average densities of macroions and salt ions.The effective pair potential arises from screening of the bare Coulomb potential by theimplicitly modeled microions.The computational benefit of reducing the number of components and the interac-tion range comes at a cost of increased complexity in the effective interactions. Linear-screening and mean-field approximations, which are justifiable in the case of weakly + + +++ + + ++++ + +− ++++ + + + + Ζ + − − +++ +− ++− − +−+++ δ a Ζ Figure 1: Model of charged colloidal suspension: spherical macroions of radius a and point microions dispersedin a dielectric continuum. Strongly associated counterions in a spherical shell of thickness δ renormalize thebare macroion valence Z to an effective (lower) valence Z . correlated monovalent microions, yield analytical forms for the volume energy and ef-fective pair potential [18–20]. Moreover, many-body effective interactions are usuallyweak [21, 22]. The range of validity of the conventional OCM is nevertheless limitedby the inevitable onset of nonlinear screening with strengthening macroion-counterioncorrelations. Previous studies [12,22] establish a parameter threshold Z λ B / a ≃ Z ≤ Z .Counterions bind to a macroion at a distance at which the electrostatic energy of at-traction is comparable to the average thermal energy per counterion. Denoting by φ ( r ) the electrostatic potential at distance r from a macroion center, the width of the associa-tion shell δ can be defined via e | φ ( a + δ ) | = Ck B T , (2.3)where C is a dimensionless parameter of order unity. In general, the reference potential should be chosen as the Donnan (average) potential φ D of the suspension, C then beinga function of salt concentration. Counterions within an association shell ( a < r < a + δ )are presumed bound by the potential well of the respective macroion, while more distantcounterions roam free and screen the dressed macroions. In other studies, a similar ther-mal criterion has been applied to the electrostatic potential [26–28] and the effective pairpotential [29]. The effective interaction theory proposed in ref. [13] combines a charge renormalizationtheory for the effective macroion charge with a mean-field, linear-response theory of theOCM [18–20]. Here we briefly outline the theory, whose predictions for effective interac-tions are subsequently used as input to our simulations (Sec. 4).By invoking a mean-field random-phase approximation for the response functions(partial static structure factors) of a microion plasma, linear-response theory yields ana-lytical expressions for effective electrostatic interactions. This approach proves formallyequivalent to density-functional [30–34], extended Debye-H ¨uckel [35], and distribution-function [36, 37] formulations of effective interaction theory [4]. Linear-response theoryis also equivalent to linearized Poisson-Boltzmann (PB) theory, based on a consistent lin-earization of the PB equation and expansion of the ideal-gas free energy functional toquadratic order in the microion density profiles [38].Linearized theories predict the reduced electrostatic potential generated by a singlebare macroion as ψ ( r ) ≡ βφ ( r ) = − Z λ B e κ a + κ a e − κ r r , r ≥ a , (3.1)where β = k B T , κ = p πλ B ( n + + n − ) denotes the bare Debye screening constant, and n ± = N ± / [ V ( − η )] are the mean number densities of microions in the free volume, i.e.,the total volume V reduced by the fraction η occupied by the macroion hard cores. Theconstraint of global electroneutrality ensures that κ depends implicitly on the averagedensities of both macroions and salt ions.Adapting Eq. (3.1) to the charge renormalization model (Fig. 1), the potential arounda dressed macroion, of effective valence Z and effective radius a + δ , becomes˜ ψ ( r ) = − Z λ B e ˜ κ ( a + δ ) + ˜ κ ( a + δ ) e − ˜ κ r r , r ≥ a + δ , (3.2)where ˜ κ = p πλ B ( ˜ n + + ˜ n − ) now represents the effective (renormalized) screening con-stant, ˜ n ± = ˜ N ± / [ V ( − ˜ η )] and ˜ N ± are the mean number densities and numbers of freemicroions, and ˜ η = η ( + δ / a ) is the effective volume fraction of the dressed macroions.Combining Eqs. (2.3) and (3.2) yields the transcendental equation Z λ B [ + ˜ κ ( a + δ )]( a + δ ) = C , (3.3) which may be solved for the association shell thickness δ for given values of Z , η , and C ,noting that ˜ κ depends self-consistently on δ .The distinction between free and bound microions implies a corresponding separa-tion of the total free energy F = F free + F bound + F m , (3.4)into three terms representing, respectively, contributions from free and bound microionsand from macroion effective interactions. The free microions are, by construction, suf-ficiently weakly correlated with the macroions to be well described by linear-responsetheory. Their free energy (per macroion) is therefore accurately approximated by thelinear-screening prediction for the volume energy (with renormalized parameters): β f free = ∑ i = ± ˜ x i [ ln ( ˜ n i Λ ) − ] − Z κλ B + ˜ κ ( a + δ ) − Z n m ˜ n + + ˜ n − , (3.5)where ˜ x ± = ˜ N ± / N m are the free microion concentrations, Λ is the microion thermal deBroglie wavelength, and n m = N m / V is the mean macroion number density. The threeterms on the right side have physical interpretations as the ideal-gas (entropic) free en-ergy of the free microions, the self-energy of the dressed macroions, and the Donnan(average) potential energy of the microions. Although the bound counterions, being rela-tively strongly correlated with the macroions, do not yield to a linear-response treatment,their free energy (per macroion) can be reasonably approximated by β f bound = ( Z − Z ) (cid:20) ln (cid:18) Z − Zv s Λ (cid:19) − (cid:21) + Z λ B a , (3.6)where the first term on the right side is an ideal-gas contribution, ( Z − Z ) / v s being themean counterion density in the association shell of volume v s = ( π /3 )[( a + δ ) − a ] , andthe second term is the electrostatic energy, assuming tight localization of the counterionsnear r = a .The effective valence Z is now prescribed, for a given bare valence Z , by the varia-tional ansatz [39–41] (cid:18) ∂∂ Z ( f free + f bound ) (cid:19) T , n ± = Z and δ are connected by Eq. (3.3). The effective valence and shellthickness then determine the effective screening constant ˜ κ .The macroion free energy F m in Eq. (3.4) depends on correlations and effective inter-actions between dressed macroions. Incorporating renormalized system parameters intolinear-response theory [19, 20], the effective pair potential is given by β v eff ( r ) = Z λ B (cid:18) e ˜ κ a + ˜ κ a (cid:19) e − ˜ κ r r , r ≥ ( a + δ ) , (3.8) whose screened-Coulomb form is identical to the long-range limit of the DLVO pair po-tential [9]. The renormalized effective interactions and system parameters can be inputinto computer simulations or liquid-state theories [42] of the OCM to compute thermo-dynamic and structural properties of bulk suspensions of charged colloids. In the nextsection, we compare our simulation results with predictions of a variational approxima-tion [22, 30] for the macroion free energy (per macroion) based on first-order thermody-namic perturbation theory with a hard-sphere reference system: f m ( n m , ˜ n ± ) = min ( d ) (cid:26) f HS ( n m , ˜ n ± ; d )+ π n m Z ∞ d d rr g HS ( r , n m ; d ) ˜ v eff ( r , n m , ˜ n ± ) (cid:27) . (3.9)Here f HS and g HS are, respectively, the excess free energy density and pair distributionfunction of the HS fluid, computed from the near-exact Carnahan-Starling and Verlet-Weis expressions [42]. Minimization with respect to the effective hard-sphere diameter d yields a least upper bound to the free energy. In practice, the renormalized system param-eters ( Z , δ , ˜ κ ) must be held fixed in this minimization and in all partial thermodynamicderivatives. The corresponding prediction for the thermodynamic pressure is β p = n m + ˜ n + + ˜ n − − Z ( ˜ n + − ˜ n − ) ˜ κλ B [ + ˜ κ ( a + δ )] + n m β (cid:18) ∂ f m ∂ n m (cid:19) T , N s / N m , (3.10)where the first three terms on the right side are ideal-gas contributions from macroionsand free microions, the fourth term arises from density dependence of the self-energy ofthe dressed macroions, and the final term is generated by effective interactions betweenpairs of dressed macroions. While this simple variational theory yields the pressure andother thermodynamic properties, computer simulations or integral-equation theory arerequired to determine structural properties. Working within the canonical (constant-
NVT ) ensemble, we consider a one-componentfluid of macroions in a cubic box, subject to periodic boundary conditions, at fixed tem-perature, volume, macroion number, and mean salt concentration. The macroions in-teract with one another via effective electrostatic interactions that depend on the meandensities of both macroions and (implicitly modeled) salt. According to the standardMetropolis algorithm [16, 17], trial particle displacements are accepted with probability P = min { exp ( − β ∆ U ) ,1 } , (4.1)where ∆ U = U ( n ) − U ( o ) is the change in pair potential energy between the new (n) andold (o) states and U = N m ∑ i < j = v eff ( | r i − r j | ) (4.2) is a sum of hard-sphere-repulsive-Yukawa (screened-Coulomb) pair potentials. Note thatthe acceptance probability for trial displacements does not involve the volume energy,since the mean density is fixed. Consequently, the volume energy does not affect themacroion structure, although it does contribute to the pressure (see below). To achievehigh numerical precision, pair interactions were cut off at a distance r c ≃
20/ ˜ κ , i.e., 20effective screening lengths. The cut-off radius determined the minimum box side length, L = r c , necessary to avoid interactions of a particle with its own periodic images. Fora given volume fraction, the requisite number of macroions was therefore prescribed as N m ≃ ( η /4 π )(
40/ ˜ κ a ) .We performed a series of simulations, using renormalized effective interaction pa-rameters (effective macroion valence, volume fraction, and screening constant), startingfrom face-centered cubic crystal configurations. Trial moves were executed by randomlydisplacing macroions with a step size adjusted to yield an acceptance rate of about 50%.Following an equilibration phase of 10 cycles, statistics were accumulated for pressuresand pair distribution functions over the next 10 cycles ( N m × particle displacements).Test simulations for larger systems confirmed finite-size effects to be negligible. Rela-tively modest computing resources were required, with typical runs on a single IntelXeon-HT processor lasting 30, 120, and 200 hours for N m = NVT ensemble, the total pressure may be expressed as β p = n m + ˜ n + + ˜ n − − Z ( ˜ n + − ˜ n − ) ˜ κλ B [ + ˜ κ ( a + δ )] + β p pair , (4.3)where p pair is generated by effective interactions between pairs of macroions and corre-sponds to the last term on the right side of Eq. (3.10). The pair pressure is computed fromthe virial expression for a density-dependent pair potential [43]: p pair = (cid:28) V int V (cid:29) − *(cid:18) ∂ U ∂ V (cid:19) ˜ N s / N m + + p tail , (4.4)where V int is the internal virial, the volume derivative term accounts for the density de-pendence of the effective pair potential, angular brackets denote an ensemble averageover configurations in the canonical ensemble, and p tail corrects for cutting off the long-range tail of the pair potential. The internal virial is given by V int = N m ∑ i = r i · f i = N m ∑ i < j = ( + ˜ κ r ij ) v eff ( r ij ) , (4.5)where f i = − ∑ j = i v ′ eff ( r ij ) is the effective force exerted on macroion i by all neighboringmacroions within a sphere of radius r c . The second term on the right side of Eq. (4.4) iscomputed as the ensemble average of (cid:18) ∂ U ∂ V (cid:19) ˜ N s / N m = − n m V N m ∑ i < j = (cid:18) ∂ v eff ( r ij ) ∂ n m (cid:19) ˜ N s / N m = V ( − ˜ η ) N m ∑ i < j = (cid:18) ˜ κ r ij − ˜ κ a + ˜ κ a (cid:19) v eff ( r ij ) . (4.6) Finally, neglecting pair correlations for r > r c , the tail pressure is approximated by p tail = − π n m Z ∞ r c d rr v ′ eff ( r ) = π n m (cid:18) ˜ κ r c + κ r c + κ (cid:19) r c v eff ( r c ) . (4.7)The macroion structure of the suspension is characterized by the macroion-macroionpair distribution function g ( r ) , defined such that 4 π r g ( r ) d r equals the average numberof macroions in a spherical shell of radius r and thickness d r centered on a macroion.With each particle regarded in turn as the central particle in a given configuration, neigh-boring particles were assigned to concentric spherical shells (bins), of thickness ∆ r = a ,according to their radial distance r from the central particle. Following equilibration, g ( r ) was computed in the range 2 ( a + δ ) < r < L /2 by accumulating statistics of particlenumbers in radial bins and averaging over all configurations. The raw distributions werefinally smoothed by averaging each bin together with its neighbors in a moving averagealgorithm. The renormalized effective interactions (Sec. 3) previously provided the basis for vari-ational theory calculations of the pressure of deionized suspensions [13]. The sameeffective interactions are here input into Monte Carlo simulations (Sec. 4) of the one-component model to compute both the pressure and macroion-macroion pair distribu-tion function. Comparing our results with corresponding data from primitive modelsimulations allows testing the accuracy of the effective interactions and the variationalfree energy approximation.All results presented below are for the case of monovalent counterions, zero salt con-centration, and an aqueous solvent at room temperature ( λ B = C =
3, a value shown in ref. [13]to give satisfactory agreement with pressure data from primitive model simulations. Forsalt-free suspensions, this choice of C corresponds to e | φ ( a + δ ) − φ D | = k B T . In passing,we note that the parameter C is analogous to the adjustable cell radius b in the PB celltheory of Zoetekouw and van Roij [44].Figure 2 illustrates the distinction between the bare and effective macroion valences.For a sufficiently small bare valence, Eq. (3.3) has no real-valued solution for any nonzeroassociation shell thickness. In this case, all counterions are free ( δ = v s =
0) and the freeenergy is minimized by Z = Z . Beyond a threshold bare valence, however, the associationshell appears and rapidly thickens with increasing Z . Correspondingly, the free energyminimum shifts to Z < Z . With increasing Z , the effective valence grows logarithmically,in contrast to the saturation observed for polyelectrolytes [24, 25] and predicted for col-loidal suspensions by PB cell theories [23,45,46] and Debye-H ¨uckel theories [39–41,47,48].We first test the capacity of the charge renormalization theory [13] to predict themacroion structure. Figures 3 and 4 compare results for the pair distribution function Bare Valence Z E ff e c t i v e V a l en c e Z Z δ / a Shell Thickness
Z=Z Figure 2: Effective valence Z vs. bare valence Z for a deionized suspension ( c s ≃ ) of macroions of radius a = nm and volume fraction η = . Inset: counterion association shell emerges and thickens beyond threshold Z . g ( r ) from our simulations of the OCM, using charge-renormalized effective interactionsas input, and corresponding data from extensive Monte Carlo simulations of the primi-tive model (PM) [49] for salt-free suspensions over ranges of bare valence, volume frac-tion, and electrostatic coupling strength Γ = λ B / a . Well beyond the threshold for chargerenormalization ( Z Γ ≃ Z Γ >
15, the OCM consistently overpredicting the macroion structure.For reference, Fig. 3 also includes OCM results obtained with unrenormalized effectiveinteractions, illustrating the impact of charge renormalization on strongly coupled sus-pensions.Finally, we test predictions of the charge renormalization theory [13] for thermo-dynamics. Figures 5 and 6 directly compare the pressures computed from our OCMsimulations [via Eqs. (4.3)-(4.7)] with corresponding data from primitive model simu-lations [49] of salt-free suspensions over ranges of bare valence, volume fraction, andelectrostatic coupling strength. The OCM and primitive model results are seen to be inexcellent agreement up to the highest couplings for which primitive model data are avail-able ( Z Γ ≃ r/ a g (r) OCM2OCM1PM Z =40 η =0.01(a) r/ a g (r) OCM2PM Z =80=0.01 η (b) r/ a g (r) OCM2PM Z =160=0.01 η (c) Figure 3: Macroion-macroion pair distribution function g ( r ) vs. radial distance r (units of macroion radius a )of salt-free suspensions computed from Monte Carlo simulations of the effective one-component model, with(OCM2, dashed curves) and without (OCM1, dotted curves) charge renormalization, and of the primitive model(PM) [49] with explicit counterions (solid curves) for volume fraction η = , electrostatic coupling parameters Γ = λ B / a = Z = (a), 80 (b), and 160 (c). For clarity, curves are vertically offset in steps of 0.5. r/ a g (r) OCM2PM Z =40 Γ =0.1779(a) r/ a g (r) OCM2PM Z =40 Γ =0.3558(b) r/ a g (r) OCM2PM Z =80 Γ =0.1779(c) r/ a g (r) OCM2PM Z =80 Γ =0.3558(d) Figure 4: Macroion-macroion pair distribution function g ( r ) vs. radial distance r (units of macroion radius a ) of salt-free suspensions computed from Monte Carlo simulations of the charge-renormalized effective one-component model (OCM2, dashed curves) and of the primitive model [49] (PM, solid curves) for various baremacroion valences Z , electrostatic coupling parameters Γ = λ B / a , and volume fractions η = Macroion Volume Fraction η P r e ss u r e β p / n t o t Figure 5: Total reduced pressure β p / n tot vs. macroion volume fraction η , where n tot = ( Z + ) n m (totalion density), of salt-free suspensions with bare macroion valence Z = and electrostatic coupling constants(from top to bottom) Γ = Γ = . The dashed curve for Γ = is off-scale, the pressure being negative. Coupling Constant
Γ=λ B / a P r e ss u r e β p / n t o t Z/Z Bare Valence Z Figure 6: Total reduced pressure β p / n tot vs. electrostatic coupling parameter Γ of salt-free suspensions withfixed volume fraction η = and bare macroion valence (top to bottom) Z =
10, 20, 40, 80. Open symbols:Monte Carlo simulations of the primitive model [49]. Filled symbols: Monte Carlo simulations of the effectiveone-component model. (Symbol sizes exceed error bars.) Curves: charge-renormalized variational theory. Inset:Ratio of effective to bare macroion valence Z / Z vs. Γ . Summarizing, we have performed Monte Carlo simulations of the one-component modelof charged colloids, using as input effective electrostatic interactions predicted by a re-cently proposed charge renormalization theory [13]. Structural and thermodynamic prop-erties of salt-free suspensions are computed and directly compared with simulations ofthe primitive model. For bare valences Z and electrostatic coupling strengths Γ consid-erably exceeding the renormalization threshold, the pair distribution function and pres-sure are in close agreement with corresponding results from simulations of the primitivemodel. Significant discrepancies between the OCM and PM simulations are observedin g ( r ) for Z Γ >
15, and in the pressure only for Z Γ >
30. The level of agreement ap-pears to be comparable to the charge-renormalization scheme of ref. [14], which is basedon a PB cell model [23]. The computationally practical approach described here mayprovide a useful modeling alternative, for some applications, to molecular-scale simu-lations of charged colloids. Further comparisons with primitive model simulation data,particularly at nonzero salt concentrations, would help to chart the validity range of thecharge-renormalized OCM.The threshold for charge renormalization closely coincides with the onset of a spin-odal phase instability predicted by linearized effective-interaction theories [22, 30, 32, 35]for deionized (but not salt-free) suspensions of highly charged colloids. Although obser-vations consistent with bulk phase separation have been reported [50–57], the experimen-tal picture of deionized suspensions remains cloudy. Furthermore, the relevant param-eter regime is not yet accessible to primitive model simulations. Recent studies basedon the PB cell model have concluded that the predicted instability may be an artifactof linearization approximations [58–60]. Preliminary calculations based on the presentapproach indicate, however, that the phase instability may survive charge renormaliza-tion [61]. Future work will continue to explore the remarkable phase behavior of deion-ized suspensions.
Acknowledgments
This work was supported by the National Science Foundation (DMR-0204020) and thePetroleum Research Fund (PRF 44365-AC7). We thank Per Linse for helpful correspon-dence and for providing primitive model simulation data.
References [1] D. F. Evans and H. Wennerstr ¨om,
The Colloidal Domain (2nd ed.) (Wiley-VCH, New York,1999).[2] K. S. Schmitz,
Macroions in Solution and Colloidal Suspension (VCH, New York, 1993).[3] P. N. Pusey, in
Liquids, Freezing and Glass Transition , session 51, ed. J.-P. Hansen, D. Levesque,and J. Zinn-Justin (North-Holland, Amsterdam, 1991). [4] A. R. Denton, in Nanostructured Soft Matter: Experiment, Theory, Simulation and Perspectives ,ed. A. V. Zvelindovsky (Springer, Dordrecht, 2007).[5] C. M. Soukoulis,
Photonic Band Gap Materials (Kluwer, Dordrecht, 1996).[6] J. E. G. J. Wijnhoven and W. L. Vos, Science , 802 (1998).[7] J. D. Joannopoulos, Nature , 257, (2001); Y. A. Vlasov, X. Z. Bo, J. C. Sturm, and D. J. Nor-ris, Nature , 289, (2001).[8] B. V. Derjaguin and L. Landau, Acta Physicochimica (URSS) , 633 (1941).[9] E. J. W. Verwey and J. T. G. Overbeek, Theory of the Stability of Lyophobic Colloids (Elsevier,Amsterdam, 1948).[10] J. Israelachvili,
Intermolecular and Surface Forces (Academic, London, 1992).[11] For reviews of colloidal interactions, see C. N. Likos, Phys. Rep. , 267 (2001);L. Belloni, J. Phys.: Condens. Matter , R549 (2000); J. P. Hansen and H. L ¨owen,Annu. Rev. Phys. Chem. , 209 (2000).[12] B. Lu and A. R. Denton, Phys. Rev. E , 061403 (2007).[13] A. R. Denton, J. Phys.: Condens. Matter , 494230 (2008).[14] V. Lobaskin and P. Linse, J. Chem. Phys. , 4300 (1999).[15] J. Liu and E. Luijten, Phys. Rev. Lett. , 035504 (2004); Phys. Rev. E , 066701 (2005).[16] D. Frenkel and B. Smit, Understanding Molecular Simulation (Academic, London, 2001).[17] M. P. Allen and D. J. Tildesley,
Computer Simulation of Liquids (Oxford, Oxford, 1987).[18] M. J. Grimson and M. Silbert, Mol. Phys. , 397 (1991).[19] A. R. Denton, J. Phys.: Condens. Matter , 10061 (1999).[20] A. R. Denton, Phys. Rev. E , 3855 (2000).[21] A. R. Denton, Phys. Rev. E , 31404 (2004).[22] A. R. Denton, Phys. Rev. E , 41407 (2006).[23] S. Alexander, P. M. Chaikin, P. Grant, G. J. Morales, and P. Pincus, J. Chem. Phys. , 5776(1984).[24] G. S. Manning, J. Chem. Phys. , 924 (1969).[25] F. Oosawa, Polyelectrolytes (Dekker, New York, 1971).[26] K. S. Schmitz, Langmuir , 4093 (1999).[27] V. Sanghiran and K. S. Schmitz, Langmuir , 7566 (2000).[28] A. K. Mukherjee, K. S. Schmitz, and L. B. Bhuiyan, Langmuir , 4210 (2002).[29] D. B. Lukatsky and S. A. Safran, Phys. Rev. E , 011405 (2000).[30] R. van Roij and J.-P. Hansen, Phys. Rev. Lett. , 3082 (1997).[31] H. Graf and H. L ¨owen, Phys. Rev. E , 5744 (1998).[32] R. van Roij, M. Dijkstra, and J.-P. Hansen, Phys. Rev. E , 2010 (1999).[33] R. van Roij and R. Evans, J. Phys.: Condens. Matter , 10047 (1999).[34] B. Zoetekouw and R. van Roij, Phys. Rev. E , 21403 (2006).[35] P. B. Warren, J. Chem. Phys. , 4683 (2000); J. Phys.: Condens. Matter , S3467 (2003);Phys. Rev. E , 011411 (2006).[36] B. Beresford-Smith, D. Y. C. Chan, and D. J. Mitchell, J. Coll. Int. Sci. , 216 (1985).[37] D. Y. C. Chan, P. Linse, and S. N. Petris, Langmuir , 4202 (2001).[38] A. R. Denton, Phys. Rev. E , 051401 (2007).[39] Y. Levin, M. C. Barbosa, and M. N. Tamashiro, Europhys. Lett. , 123 (1998).[40] M. N. Tamashiro, Y. Levin, and M. C. Barbosa, Eur. Phys. J. B , 337 (1998).[41] A. Diehl, M. C. Barbosa, and Y. Levin, Europhys. Lett. , 86 (2001).[42] J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids , 2 nd edition (Academic, London,1986). [43] A. A. Louis, J. Phys.: Condens. Matter , 9187 (2002).[44] B. Zoetekouw and R. van Roij, Phys. Rev. Lett. , 258302 (2006).[45] L. Belloni, Colloids Surf. A , 227 (1998).[46] L. Bocquet, E. Trizac, and M. Aubouy, J. Chem. Phys. , 8138 (2002).[47] Y. Levin, E. Trizac, and L. Bocquet, J. Phys.: Condens. Matter , S3523 (2003).[48] E. Trizac and Y. Levin, Phys. Rev. E , 031403 (2004).[49] P. Linse, J. Chem. Phys. , 4359 (2000). Reduced pressures in Table III of this paper shouldbe β p / n tot = Z = Γ = η = Z = Γ = η = Z = Γ = η = , 3778 (1992).[51] K. Ito, H. Yoshida, and N. Ise, Science , 66 (1994).[52] N. Ise and H. Yoshida,
Acc. Chem. Res. , 3 (1996).[53] B. V. R. Tata, E. Yamahara, P. V. Rajamani, and N. Ise, Phys. Rev. Lett. , 2660 (1997).[54] N. Ise, T. Konishi, and B. V. R. Tata, Langmuir , 4176 (1999).[55] H. Matsuoka, T. Harada, and H. Yamaoka, Langmuir , 4423 (1994); H. Matsuoka, T. Harada,K. Kago, and H. Yamaoka, ibid , 5588 (1996); T. Harada, H. Matsuoka, T. Ikeda, and H. Ya-maoka, ibid , 573 (1999).[56] A. E. Larsen and D. G. Grier, Nature , 230 (1997).[57] F. Gr ¨ohn and M. Antonietti,
Macromolecules , 5938 (2000).[58] R. Klein and H. H. von Gr ¨unberg, Pure Appl. Chem. , 1705 (2001).[59] M. Deserno and H. H. von Gr ¨unberg, Phys. Rev. E , 011401 (2002).[60] M. N. Tamashiro and H. Schiessel, J. Chem. Phys.119