# Effective charges and virial pressure of concentrated macroion solutions

^{a}Department of Materials Science and Engineering, Northwestern University, Evanston, IL 60208;^{b}Instituto de Física, Universidad Autónoma de San Luis Potosí, 78000 San Luis Potosí, Mexico;^{c}Institute for Theoretical Physics, Center for Extreme Matter and Emergent Phenomena, Utrecht University, 3584 CE Utrecht, The Netherlands

See allHide authors and affiliations

Contributed by Monica Olvera de la Cruz, June 18, 2015 (sent for review April 29, 2015; reviewed by Roland R. Netz)

## Significance

Colloids constitute the basic components of many everyday products and are integrated into the fabric of modern society. Understanding their assembly is key for nanotechnological and biotechnological advances. At the single-particle level, colloids commonly possess electric charge. Consequently, the structure to which they conform is strongly influenced by electrostatic interactions. In solution, these interactions are modified by the presence of ions. We have developed a model for computing the corresponding effective electrostatic interactions as well as the osmotic pressure. Our model extends the applicability of Derjaguin−Landau−Verwey−Overbeek theory to dense systems in which many-body effects are crucial. This will allow previously impossible, mesoscale studies of colloidal assembly to be performed analytically or by simulation with implicit ions models.

## Abstract

The stability of colloidal suspensions is crucial in a wide variety of processes, including the fabrication of photonic materials and scaffolds for biological assemblies. The ionic strength of the electrolyte that suspends charged colloids is widely used to control the physical properties of colloidal suspensions. The extensively used two-body Derjaguin−Landau−Verwey−Overbeek (DLVO) approach allows for a quantitative analysis of the effective electrostatic forces between colloidal particles. DLVO relates the ionic double layers, which enclose the particles, to their effective electrostatic repulsion. Nevertheless, the double layer is distorted at high macroion volume fractions. Therefore, DLVO cannot describe the many-body effects that arise in concentrated suspensions. We show that this problem can be largely resolved by identifying effective point charges for the macroions using cell theory. This extrapolated point charge (EPC) method assigns effective point charges in a consistent way, taking into account the excluded volume of highly charged macroions at any concentration, and thereby naturally accounting for high volume fractions in both salt-free and added-salt conditions. We provide an analytical expression for the effective pair potential and validate the EPC method by comparing molecular dynamics simulations of macroions and monovalent microions that interact via Coulombic potentials to simulations of macroions interacting via the derived EPC effective potential. The simulations reproduce the macroion−macroion spatial correlation and the virial pressure obtained with the EPC model. Our findings provide a route to relate the physical properties such as pressure in systems of screened Coulomb particles to experimental measurements.

Coulombic interactions between ionized species affect colloidal suspensions at the microscopic level and have an indirect, yet crucial, impact on the observable macroscopic characteristics of the system (1). The Derjaguin–Landau–Verwey–Overbeek (DLVO) theory (2, 3), proposed in the 1940s, has been crucial for understanding like-charged colloidal dispersions in a wide variety of experimental conditions. In this theory, the effective pair potential between two equally charged macroions immersed in an electrolyte is expressed as the sum of three terms: a hard-core potential that takes into account the excluded volume of macroions (preventing their overlap), an attractive potential due to short-range (van der Waals) interactions, and an electrostatic screened Coulomb or Yukawa potential resulting from the linearized Poisson–Boltzmann theory, which is the Debye–Hückel approximation. Many additions and modifications to the original theory have been proposed, including polarization effects, patchiness, or charge regulation, just to mention a few. Special care should be taken for nonaqueous solvents, divalent ions, or high salt concentrations, since, in these regimes, ion correlations are usually important (4⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–16). Generally speaking, modifications to the DLVO theory have been pivotal for systems in which the electrostatics are not well described by the linearized Poisson–Boltzmann (PB) theory. Although the DLVO theory has been used extensively to model colloidal dispersions (17), this approach is not exact within the context of the underlying Debye–Hückel approximation. In the precise analysis of the force between two charged spheres in an electrolytic solution, Verwey and Overbeek encountered additional terms that can be considered cross terms, resulting from the exclusion of the ionic double layer surrounding the first sphere by the hard core of the second sphere (3). Numerical methods exist to quantify such effects (18⇓–20), yet these approaches explicitly deal with two particles in an otherwise empty system. While in dilute solutions of macroions the resulting correction to DLVO is typically small and can usually be safely neglected, in dense macroion systems, the deviation from DLVO can become very significant due to the overlap between each macroion electrical double layer with the hard cores of all neighboring particles. As a result, the performance of the classical DLVO equation is limited to the description of dilute systems of macroions (21⇓⇓⇓–25), while many colloidal processes such as crystallization or glass formation predominantly occur in dense systems where many-body effects prevail. Marcelja et al. recognized the importance of many-body effects at high colloidal densities and low electrolyte concentration, and they described a method that uses cell theory to project a charged colloidal dispersion to a system of Coulomb particles (26). This enabled use of Monte Carlo simulations of the Wigner lattice to study the crystallization of latex suspensions. That method has recently been rediscovered and extended to charge-regulating particles to describe reentrant melting on addition of a charging agent to a colloidal suspension (27). Such an effective Coulomb representation is, however, not suitable to describe the structure of charged suspensions, particularly in the crystalline phase. This is also true of methods comprising the repulsive forces among macroions via hard-sphere interactions with effective hard-sphere radii (28⇓⇓–31). Different approaches such as the (renormalized) Jellium model (32) and methods that calculate the osmotic pressure within a Wigner–Seitz cell (33, 34) have been proposed. However, they do not yield information on the spatial configuration of the macroions and consequently are limited in describing dense macroion systems.

## Model

In this work, we introduce a method to calculate the effective electrostatic pair interaction between macroions in dense systems through the identification of their corresponding effective point charges. We verify the corresponding accuracy by comparing the resulting radial distribution functions and pressures to the primitive model (PM). To begin, we consider spherical and impenetrable macroions of valence *Z* and radius *a* immersed in a 1:1 electrolyte with bulk concentration *e* is the elementary charge, and *ϵ* is the relative dielectric permittivity. For sufficiently small charges, the PB equation can be linearized by using *D*,

Apart from being restricted to dilute systems, the DLVO equation above cannot directly be applied to strongly charged macroions, since the Debye–Hückel approximation no longer holds for these systems, which, strictly speaking, leads to nonpairwise additive interaction potentials (36, 37). Alexander et al. (33), however, showed that nonlinear ion behavior close to the macroion surface can be embodied in an effective linear screening model by calculating a renormalized surface charge *F*. Regarding a system of macroions at a concentration *C*. In this spherical geometry, the nonlinear PB equation and the associated boundary conditions can be written as*r*. The boundary conditions follow from Gauss’ law and include the global electroneutrality condition of the whole system. This set of equations is typically solved numerically, as no general analytical solution is known. Once the numerical solution is determined, one proceeds by linearizing Eq. **2a** around the obtained potential at the cell boundary, which can be regarded as the Donnan potential, *F*), and one finds *κ* and

The accuracy and simplicity of the previous cell model approach can, however, be improved by calculating an effective point charge *Q* directly through identification of a point charge at *G*), yielding the form *κ* and *Q* can then be used to approximate the effective electrostatic interactions in the original macroion system by those of point charges, using Eq. **1** to find the pairwise interaction energy. Although high macroion volume fractions render DLVO-based approaches inaccurate (21⇓⇓⇓–25), the effective system of point charges has no hard-core volumes that will overlap with ionic double layers. We therefore expect that Eq. **3** in combination with Eq. **1** will remain accurate even in dense macroion systems. Note that the hard-core repulsions for

In the regime where *Z* is small and the resulting potential profile is sufficiently flat throughout the cell, **3** will become exact on the entire space between the cell boundary and the macroion surface. As a consequence, *Z* and *Q* follows (Fig. 1*D*). Tantalizingly, inserting this *Q* into Eq. **1** yields a pair potential similar to the DLVO equation,*κ* that enters Eq. **4** reduces to the reservoir value *η*-independent) DLVO theory is reobtained for dilute suspensions, which is the limit **4** reduces to the Coulombic form *Z* due to the expulsion of the ionic background from the hard core (26, 27, 44⇓–46).

To verify our proposed prescription, molecular dynamic (MD) simulations of macroion/microion mixtures with particle diameters *L* under periodic boundary conditions. In the PM representation that we applied here, ionic species are represented by repulsive-core spheres with point charges in their centers immersed in a continuous solvent (48⇓⇓–51). The pairwise forces among all particles have a short-range repulsive-core potential component, *i* and *j*, respectively. These interactions are handled properly, using the particle mesh Ewald technique (52). We model the repulsive-core pair potential between a particle of species *i* and a particle of species *j*, separated by a distance *D*, as an impenetrable hard-core *σ* regulates the hardness of the repulsive-core interactions. To mimic the hard-core interaction characteristic of the PM, *σ* is set equal to 0.1 nm. We use

## Results

In Fig. 2, we compare radial distributions from computationally expensive PM MD simulations (circles) to much faster and more economic effective-model descriptions using MD simulations (solid lines), and integral equations (Fig. 2*B*, dashed lines). In the PM approach, we use a cubic simulation box of length *e*), and 2,880 small monovalent coions (+*e*). In the effective-model approach, microions are included implicitly in the Yukawa interactions between macroions with an effective charge *Q* and inverse screening length *κ*. The charges associated to the macroion profiles shown in Fig. 2 are *A*. In contrast, surface charge renormalization in combination with the DLVO theory deviates significantly from PM simulation results, as expected at this volume fraction. The use of integral equations theory allows for an even faster numerical calculation of the radial distribution functions within the effective model. The Rogers–Young (RY) closure (53), which is known for its superb accuracy for hard-core Yukawa systems (46, 54⇓–56), has good qualitative agreement compared with PM results, as seen in Fig. 2*B*. Note that PB techniques are grand canonical and therefore require a reservoir ion density

The total microion/macroion pressure resulting from the PM, *N* sums all particles, in the PM, or only the macroions, in the effective model, and *i* and *j*. However, to relate *SI Text* for this. Note that Eq. **5** implies that the macroion osmotic pressure, which is the pressure with respect to the ion reservoir, *κ* approaches **5** reduce to

Fig. 3 shows the PM pressure **5** applying the RY closure. Even though the agreement between the radial distribution functions obtained from the EPC method using the RY closure and the PM simulations is not as accurate as that obtained from the EPC effective screened Coulomb simulations (see, e.g., Fig. 2), we have observed that the virial pressure obtained from the RY closure and the effective screened Coulomb simulations using the EPC approach agreed with the pressure *Z* (Fig. 3*A*), and in particular for high *η* (Fig. 3*B*). While the dashed lines in Fig. 3 represent added salt cases, the full lines correspond to a system without coions, i.e., a salt-free “reservoir”: **2a** and **2b** in principle cannot be solved. Our cell calculations, however, show that the salt-free limit *Q*, and *κ* converge to their values in a salt-free environment. Salt-free cases therefore do not require a root-finding procedure with respect to *κ* is in accordance with a homogeneous distribution of neutralizing counterions, **4**, yields a description at any density of the pair interactions in salt-free systems depending on *Z*, *a*, and *η* only. The pressure correction in Eq. **5** reduces to *κ* on *η*, is particularly important in other geometries, such as the charged two-plate system, from what we study here. In the latter, the distance between the plates also sets the system’s volume and therefore *κ*, yielding a nonexponential form for the pair potential (63). However, for macroion suspensions this is not the case, as *η* is a fixed parameter independent of the configuration.

## Conclusions

In this work, we have introduced a method that extends the capabilities of the DLVO theory to high valences and volume fractions of colloidal macroions using no additional assumptions besides the underlying PB theory. Akin to the original theory, this approach is mainly applicable to systems with monovalent microions for which ion correlations are unimportant, although approximate extensions to systems with correlated multivalent counterions might be obtainable. We also propose a route to relate the pressure in the effective system of macroions to the osmotic pressure that can be measured experimentally in colloidal systems, for example, in sedimentation profiles (64, 65). Our method demonstrates accuracy with respect to acquiring the measurable properties of charged colloidal suspensions, and can therefore be applied to guide and interpret experiments on related systems.

## SI Text

Here we will demonstrate how the pressure difference between the effective Yukawa model and the PM, as stated in Eq. **5**, can be derived. We consider a very general charge distribution **S1** with respect to the ion densities **S1** can be approximated as a quadratic function of *V* the system’s volume. For a specific macroion/microion mixture, we determine *Model*. The effective point charge *Q* yields a charge distribution *M* macroions, and we choose **S2** may be transformed into*V*, choosing to keep **S3** for infinitesimal changes of *V*. We obtain

## Acknowledgments

We thank the support of the Center for Bio-Inspired Energy Science, which is an Energy Frontier Research Center funded by the US Department of Energy, Office of Science, Basic Energy Sciences under Award DESC0000989. G.I.G.-G. acknowledges the Mexican National Council of Science and Technology (CONACYT) for the financial support via the program “Catedras CONACYT.” The computer cluster where part of the simulations were performed was funded by the Office of the Director of Defense Research and Engineering and the Air Force Office of Scientific Research under Award FA9550-10-1-0167. This work is part of the Delta Institute for Theoretical Physics (D-ITP) consortium, a program of The Netherlands Organisation for Scientific Research, that is funded by the Dutch Ministry of Education, Culture and Science.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: m-olvera{at}northwestern.edu.

Author contributions: N.B., G.I.G.-G., R.v.R., and M.O.d.l.C. designed research; N.B. and G.I.G.-G. performed research; N.B., G.I.G.-G., R.v.R., and M.O.d.l.C. analyzed data; and N.B., G.I.G.-G., R.v.R., and M.O.d.l.C. wrote the paper.

Reviewers included: R.R.N., Free University of Berlin.

The authors declare no conflict of interest.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1511798112/-/DCSupplemental.

## References

- ↵
- ↵.
- Derjaguin B,
- Landau L

- ↵.
- Verwey EJW,
- Overbeek JTG

- ↵
- ↵
- ↵
- ↵
- ↵.
- Valleau J,
- Ivkov R,
- Torrie G

- ↵.
- Oosawa F

- ↵
- ↵
- ↵
- ↵
- ↵.
- Kjellander R

- ↵.
- Wernersson E,
- Kjellander R,
- Lyklema J

_{4}interface: Toward the solution of a long-standing issue. J Phys Chem C 114(4):1849–1866 - ↵.
- Zwanikken JW,
- Olvera de la Cruz M

- ↵
- ↵.
- Glendinning A,
- Russel W

- ↵
- ↵
- ↵
- ↵
- ↵.
- Warren PB

- ↵
- ↵
- ↵
- ↵
- ↵.
- Stigter D

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Griffiths D

- ↵
- ↵
- ↵.
- Belloni L

- ↵
- ↵
- ↵
- ↵
- ↵.
- Zoetekouw B

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Zhou S,
- Zhang X

- ↵
- ↵.
- Graf H,
- Löwen H

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences