Weak Long-Ranged Casimir Attraction in Colloidal Crystals

Ajay Gopinathantex2html_wrap_inline^(1), Tong Zhoutex2html_wrap_inline^(2), S. N. Coppersmithtex2html_wrap_inline^(3) , L. P. Kadanofftex2html_wrap_inline^(1) and David G. Griertex2html_wrap_inline^(1)
tex2html_wrap_inline^(1) Dept. of Physics and James Franck Institute, The University of Chicago, 5640 S. Ellis Avenue, Chicago, IL 60637
tex2html_wrap_inline^(2) Physics Department, University of California, Santa Barbara, CA 93106-9530
tex2html_wrap_inline^(3) Physics Department, University of Wisconsin, 1150 University Avenue, Madison, WI 53705

Abstract:

We investigate the influence of geometric confinement on the free energy of an idealized model for charge-stabilized colloidal suspensions. The mean-field Poisson-Boltzmann formulation for this system predicts pure repulsion among macroionic colloidal spheres. Fluctuations in the simple ions' distribution provide a mechanism for the macroions to attract each other at large separations. Although this Casimir interaction is long-ranged, it is too weak to influence colloidal crystals' properties.

Experimental evidence collected over 20 years [1] suggests that similarly-charged colloidal spheres dispersed in water need not simply repel each other. Under some circumstances they instead experience an unexpected long-ranged attraction. For example, like-charge attractions are implicated in the cohesion of metastable superheated colloidal crystals [2,3] even though isolated pairs of the constituent spheres are observed to repel each other [4,5]. Comparable attractions have been measured for pairs of spheres confined by two [5,6] charged planar walls. Recent calculations [7,8] reveal that such confinement-induced attractions cannot be accounted for by local density theory nor by electrohydrodynamic coupling [9,10]. Such anomalous effects in charge-stabilized colloid therefore challenge our general understanding of interactions and dynamics in macroionic systems.

This Letter addresses fluctuations' contribution to the free energy of highly charged colloidal spheres surrounded by a neutralizing cloud of small singly-charged counterions. Highly symmetric monopolar fluctuations in the counterion distribution increase the system's free energy. We demonstrate that their suppression by boundary conditions at the spheres' surfaces introduces a long-range attraction into the crystal's free energy analogous to the Casimir force in quantum electrodynamics, but that it is too weak to account for anomalous behavior in charge-stabilized suspensions.

Our treatment is based on the Wigner-Seitz cell model introduced by Wennerström, Jönsson and Linse [11] which has been studied extensively [12] both theoretically and through Monte Carlo simulation. It consists of a single spherical macroion of radius $a$ carrying a uniformly distributed surface charge $-Ze$ and surrounded by a thermal cloud of $Z$ point-like counterions at temperature $T$, each carrying a single charge $e$. The macroion and counterions are confined by a concentric conducting spherical shell of radius $R$. This outer shell plays a role analogous to the Wigner-Seitz cell boundary in a colloidal crystal. More generally, it models the crowding or geometric confinement characteristic of colloidal crystals [13].

Previous investigations of this and related models [12,14] have found short-ranged correlation-driven attractions between the bounding surfaces under some conditions, particularly when the counterions are polyvalent. They have not found evidence for long-ranged attractions in monovalent electrolytes [13].

Our method for evaluating the counterions' partition function allows us to investigate much higher macroion charges than have been considered before. The outer boundary's suppression of counterion fluctuations induces a long-ranged Casimir-like attraction [14] between the macroion and its neighbors across the Wigner-Seitz cell boundary. Although this cell model is far too simple to describe the behavior of real charge-stabilized suspensions, it highlights a previously unexplored mechanism for long-ranged confinement-induced like-charge colloidal attractions.

We adopt the path integral formalism reviewed in [14] and write the counterions' canonical partition function as a functional integral over all possible counterion distributions, $n({\vec r})$:

$${\cal Q}= \sqrt{2\pi Z} \, Z^Z \, \int^\prime e^{-\beta f[n]} \, Dn,$$ (1)

where $\beta^{-1} = k_B T$ is the thermal energy scale at temperature $T$, the prime indicates that the number of charges is conserved ( $\int n \, dV = Z$), and [15] $f[n] = U[n] + k_B T \, \int n \, \ln n \, dV.$ The potential energy functional $U[n] = \frac{1}{2} \int ne \phi \, dV,$ describes the counterions' interaction with the local electric potential $\phi({\vec r})$. The system's Helmholtz free energy is then $F = -k_B T \ln {\cal Q}$.

One ionic distribution, $n_0({\vec r})$, minimizes $f$, and thus has the greatest statistical weight in ${\cal Q}$. We factor the partition function ${\cal Q}= {\cal Q}_0 {\cal Q}_{fl}$ into the saddle point contribution ${\cal Q}_0 = Z^Z \, e^{-\beta f_0},$ where $f_0 = f[n_0]$, and a term ${\cal Q}_{fl}$ accounting for fluctuations, $\delta n$, away from $n_0$. Expressing ${\cal Q}_{fl}$ as a series expansion in $\delta n$ yields as the lowest-order non-vanishing term

$${\cal Q}_{fl} \simeq \sqrt{2\pi Z} \, \int^\prime \, e^{-\beta \delta^2f} \; Dn,$$ (2)

where $\delta^2f$ is the second-order change in $f[n]$ due to $\delta n$.

Terminating this expansion at Gaussian order is justified if corrections at higher order in $\delta n$ contribute negligibly to ${\cal Q}_{fl}$. This condition is met if $Z$ is large [16] and $n_0(r)$ itself changes negligibly over the mean radial counterion separation:

$$\frac{1}{n_0} \, \frac{dn_0}{dr} \ll 4 \pi r^2 n_0,$$ (3)

and can be tested a posteriori once $n_0(r)$ is evaluated. Expanding around any other distribution [17,18] would not give such a criterion for establishing convergence.

It has been shown [16,18] that the saddle point corresponds to the mean-field result $n_0 = n_s e^{-\beta e \phi}$ which, combined with the Poisson equation $\nabla^2 \phi = - \frac{en}{\epsilon},$ yields the familiar Poisson-Boltzmann (PB) equation,

$$\nabla^2 \phi = - \frac{e n_s}{\epsilon} \, e^{-\beta e \phi}.$$ (4)

Selecting $n_s = 3Z/(4\pi a^3)$ sets the potential's reference point conveniently without loss of generality. Eq. (4) accounts for the solvent's influence in the so-called primitive model through its dielectric constant $\epsilon$.

Following conventional practice, we introduce the Bjerrum length $\lambda_B = \beta e^2 / (4 \pi \epsilon)$ and an effective screening length $\kappa^{-1} = (4 \pi n_s \lambda_B)^{-1/2}$. Taking the system's radial symmetry into account leads to

$$\frac{d^2\phi}{dr^2}+\frac{2}{r}\frac{d\phi}{dr} - \frac{\kappa^2}{e \beta} e^{- \beta e \phi}=0.$$ (5)

We solve Eq. (5) subject to two boundary conditions: Gauss' theorem at the macroion's surface gives

$$\left. \frac{d\phi}{dr} \right\vert _{r= a} = \frac{\kappa^2 a}{3\beta e},$$ (6)

and electroneutrality requires

$$\left. \frac{d\phi}{dr} \right\vert _{r = R} = 0.$$ (7)

Equations (5-7) can be solved numerically for $\phi(r)$ from which we can calculate the mean-field free energy,

$$f_0 = - \frac{Z}{2} \, \left[ e \phi(a) + \frac{3e}{a^3} \, \int_a^R \phi \, e^{-\beta e \phi} \, r^2 \, dr \right]$$ (8)

$$= \frac{Z}{2} \, k_B T \, \left[ \ln \left( \frac{n_0(a)}{n_s} \rig... ...\int_a^R \frac{n_0}{n_s} \, \ln \left(\frac{n_0}{n_s}\right) \, r^2 dr \right].$$ (9)

Figure 1 shows $f_0$ as a function of $R/a$, calculated numerically from the solution of the full PB equation. Since $f_0(R)$ decreases monotonically with $R$, the mean-field theory predicts pure repulsion.

Figure 1: Dependence of the mean-field free energy on the confining boundary's radius, $R$, from Eq. (9) with $Z=7300$, $a = 0.326 \mu$m and $\lambda _B = 0.7$ nm. The solid curve shows the solution of the full PB equation while the dashed line is from the linearized PB theory.

Figure: Numerically evaluated ratio $\vert \chi_a/\chi_r\vert$ maximised over $R$ for various values of the macroion radius $a$ and the macroion charge $Z$. A smooth surface is fitted to the data points to aid the eye. The plot clearly demonstrates that the attraction is four orders of magnitude too small to overcome the repulsion.

An equivalent result obtained by solving the linearized PB equation for the same system does have a local minimum for $R$ slightly bigger than $a$. In the context of colloidal interactions, the linearized mean-field description yields the electrostatic component of the pair potential due to Derjaguin, Verwey, Landau and Overbeek (DLVO) [19]. The local minimum in the free energy suggests an effective electrostatic attraction between like-charged colloidal spheres within the DLVO theory. Comparison with the full calculation, however, shows this to be an artifact of linearization.

In addition to confirming the absence of attractions in the mean-field description, our numerical results also satisfy the condition in Eq. (3). Thus, we are justified in using Eq. (2) to calculate ${\cal Q}_{fl}$. The second-order change in $U[n]$ due to $\delta n$ is

$$\delta^2U = \frac{e^2}{8\pi\epsilon } \, \int d^3r_1\int d^3r_2 \, \frac{\delta n(\vec{r}_1)\delta n(\vec{r}_2)}{\vert\vec{r}_1-\vec{r}_2\vert}.$$ (10)

So, from Eq. (2),

$${\cal Q}_{fl} = \sqrt{2\pi Z} \, \int^\prime \exp \left[-\beta\delta^2U-\int\frac{(\delta n)^2}{2n_0} dV \right] \, Dn.$$ (11)

We evaluate the functional integral by partitioning the system into $N$ concentric shells,

$${\cal Q}_{fl} = \sqrt{2\pi Z} \, \int^\prime \prod_{i=1}^N \frac{... ...nt_{V_j} \int_{V_k} \frac{d^3r_jd^3r_k}{\vert\vec{r}_j-\vec{r}_k\vert} \right],$$ (12)

where $V_i$ is the volume of cell $i$, $n_i$ the number density of counterions therein, and $n_{0i}$ the equivalent number density for the mean-field case. This highly symmetric partition of the system is appropriate for radial or monopolar density fluctuations. A more complete description including multipole fluctuations is not necessary for our purposes because recent calculations reveal that these yield only short-ranged attractions [20]. Cross terms cancel at the Gaussian level of approximation, so we can consider monopole and multipole contributions separately.

We define $x_i = \delta n_i/\sqrt{2n_{0i}/V_i}$, and divide the system so that every cell has the same number of counterions in the mean-field distribution, i.e. $n_{0i}V_i = Z/N$. Rescaled in this way,

$${\cal Q}_{fl} = \sqrt{2\pi Z} \, \int^\prime \prod_{i=1}^N \frac{... ...{V_k} \frac{d^3r_j d^3r_k}{\left\vert \vec r_j - \vec r_k \right\vert} \right].$$ (13)

Then, numbering the cells from the center outwards so that $r_i < r_j$ if $i<j$, we obtain

$${\cal Q}_{fl} = \sqrt{2\pi Z}\int^\prime \prod_{i=1}^N \frac{dx_i}{\sqrt{\pi}} \e... ...}{r_j}\right)x_j^2 -\frac{Z \lambda_B}{N}\sum_{j<k} \frac{2}{r_k}x_jx_k \right]$$ (14)

$$= \sqrt{2\pi Z}\int^\prime \prod_{i=1}^N\frac{dx_i}{\sqrt{\pi}} \exp\left[-\vec{x}\cdot {\mathbf{A}} \cdot\vec{x}\right].$$ (15)

The prime indicates that the integral is constrained to maintain electroneutrality: $\sum_{i=1}^N x_i=0$. Because of this condition, we can subtract a constant from each element of $\mathbf{A}$ without changing the value of ${\mathbf{A}}\cdot\vec{x}$. Let us subtract $Z \lambda_B/(N R)$ from each element of $\mathbf{A}$ to obtain matrix $\mathbf{\tilde A}$. All of the elements in the last row and column of $\mathbf{\tilde A}$ vanish except ${\mathbf{\tilde A}}_{N,N}$ which is $1$. Separating mode $x_N$ from the other $N-1$ modes in this way allows us to remove the explicit constraint:

$$\int^\prime \prod_{i=1}^N\frac{dx_i}{\sqrt{\pi}} \, e^{-{x_N}^2} = \int\prod_{i=1}^N \frac{dx_i}{\sqrt{\pi}} \, e^{-{x_N}^2} \, \delta\left(\sqrt{\frac{2Z}{N}} \sum_{j=1}^N x_j\right)$$ (16)

$$= \sqrt{\frac{N}{2\pi Z}} \, \int \prod_{i=1}^{N-1}\frac{dx_i}{\sqrt{\pi}} \; \exp\left[-\left(\sum_{j=1}^{N-1}x_j\right)^2\right].$$ (17)

Consequently,

$${\cal Q}_{fl} = \sqrt{N} \, \int \prod_{i=1}^{N-1}\frac{dx_i}{\sq... ...-\vec{x}\cdot {\mathbf{B}} \cdot\vec{x}} = \sqrt{\frac{N}{\det {\mathbf{B}} }},$$ (18)

where $\vec{x}$ is now a vector of dimension $N-1$, and ${\mathbf{B}}$ is a matrix of dimension $N-1$ obtained by adding $1$ to each element in the first $N-1$ rows and columns of matrix ${\mathbf{\tilde A}}$.

$\mathbf{B}$ may be expressed as the sum of two matrices, $\mathbf{C}$ and $\mathbf{D}$, whose components are ${\mathbf{C}}_{ij} = 1 + \delta_{ij}$ and

$${\mathbf{D}}_{ij} = \frac{Z}{N} \, \lambda_B \, \left( \frac{1}{r_p} - \frac{1}{R} \right),$$ (19)

where $p$ is the greater of $i$ and $j$. $\det {\mathbf{C}} = N$, so that $\det {\mathbf{B}} = N \, \det {{\mathbf{I}} + {\mathbf{C}}^{-1}{\mathbf{D}}}$ where $\mathbf{I}$ is the identity matrix.

Evaluating $\det \mathbf{B}$ is greatly facilitated if the components of ${\mathbf{C}}^{-1}{\mathbf{D}}$ are much smaller than 1. ${\mathbf{C}}_{ij}^{-1} = \delta_{ij} - 1/N$ differs little from the identity matrix. The components of $\mathbf{D}$, on the other hand, are bounded above by ${\mathbf{D}}_{ij} < (Z/N) \, (\lambda_B/a)$. We previously assumed $Z/N \gg 1$ in deriving Eq.  (2). But $\lambda_B/a \ll 1$ for the micron-sized spheres in experimental observations, so that we may reasonably assume ${\mathbf{D}}_{ij} < 1$. Even if this were not the case, we would be justified in formally taking the limit $Z/N \ll 1$ at this point because the final result cannot depend on $N$. Consequently,

$$\det {\mathbf{B}} = N \, \exp\left[ \mbox{Tr}\ln \left({\mathbf{I... ...\right] \approx N \, \exp \left( \mbox{Tr}{\mathbf{C}}^{-1}{\mathbf{D}}\right).$$ (20)

In this approximation, the fluctuation contribution to the free energy is

$$\delta F = \frac{1}{2} \, k_B T \, {\mathbf{C}}^{-1}{\mathbf{D}} = \frac{Z}{2} \, k_B T \, \frac{1}{... ... \sum_{j=1}^k \left( \frac{\lambda_B}{r_j} - \frac{\lambda_B}{r_{k+1}} \right).$$ (21)

Rewriting the sums over shell indices as integrals over radii, we obtain

$$\delta F \approx \frac{Z}{2} \, k_B T \, \left( \frac{4\pi}{Z} \r... ...0(r^\prime) \, \left( \frac{\lambda_B}{r^\prime} - \frac{\lambda_B}{r} \right).$$ (22)

Most of the boundaries, $r_i$, between cells in the mean-field distribution are clustered near $a$. Consequently,

$$\delta F \approx 2 \pi k_B T \, \int_a^R n_0(r) \, w_0(r) \, \left( \frac{\lambda_B}{a} - \frac{\lambda_B}{r} \right) \, r^2 dr,$$ (23)

where $w_0(r) = (4\pi/Z) \, \int_a^r n_0(r^\prime) \, {r^\prime}^2 \, dr^\prime$ is the fraction of counterions within radius $r$, in the mean field approximation. Integrating by parts then yields

$$\delta F \approx \frac{Z}{2} \, k_B T \, \lambda_B \, \int_a^R \frac{1 - w_0^2(r)}{r^2} \, dr.$$ (24)

Unlike $f_0(R)$, $\delta F(R)$ decreases with decreasing $R$ because the outer boundary condition suppresses fluctuations as $R$ approaches $a$. The resulting attraction therefore is reminiscent of the Casimir attractions previously identified in confined electrolytes as well as other systems [14]. It is interesting to note that this mechanism does not yield an attractive contribution in all geometries; the second term in the exponent of Eq. (13) vanishes for unbounded systems such as parallel plates and concentric cylinders. This is consistent with the absence of long-ranged like-charge attractions in measurements [21], theoretical treatments [12,14], and simulations [12,14] of unbounded systems. Such attractions, therefore, are peculiar to closed systems, such as the Wigner-Seitz cells of colloidal crystals. If there is an $R$ at which the attractive force

$$\chi_a \equiv - \frac {\partial \delta F(R)}{\partial R} = \frac{... ...{B}}{2} \int_a^R \frac{2 w_{0}(r)}{r^2} \frac{\partial w_{0}(r)}{\partial R} dr$$ (25)

has larger absolute value than the mean-field repulsive force

$$\chi_r \equiv -\frac{\partial f_0(R)}{\partial R}$$

$$= -\frac{Zk_{B}T}{2} \left[ n_{0}(R) \ln \left(\frac{n_{0}(R)}{n_{s... ...}{\partial R} \ln \left(\frac{n_{0}(r)}{n_{s}}\right ) \right ) r^2 dr \right],$$ (26)

then $F(R) = f_0(R) + \delta F(R)$ would have a minimum at that $R$. Such a minimum would correspond to a fluctuation-mediated bound state for a colloidal crystal of nearest neighbor spacing $2R$. However, numerical results shown in Fig. 2 reveal $\vert \chi_a/\chi_r\vert < 10^{-4}$ over the entire range of conditions studied experimentally. On this basis, we conclude that monopolar fluctuations are not responsible for the strong and long-ranged attractions reported in measurements on charge-stabilized colloid [2,3,5,6].

To understand why the ratio $\vert \chi_a/\chi_r\vert$ is small we estimate it analytically by assuming that most of the counterions are clustered close to the macroion surface and that only a small fraction $p(Z,a,R)$ of counterions are affected by a change in the radius $R$ of the confining shell, and that this fraction is uniformly distributed in the volume. We then evaluate $\chi_a(R)$ and $\chi_r(R)$ to ${\cal O}\left(p\right)$. Truncating to ${\cal O}\left(p\right)$ is justified by numerical investigation of the mean-field solution which indicates $p \leq 0.05$ in the region of interest. In this approximation

$$\vert \chi_a\vert = \frac{Zk_B T \lambda_B}{2a^2} \frac{3p}{2} \f... ...rt \chi_r \vert} = \frac{\lambda_B}{a} \frac{(x^2+x-2)}{(x^2+x+1)} \frac {1}{x}$$ (27)

where $x = R/a$, is independent of $p$ and has a maximium value $0.286~\lambda_B/a$ at $R = 2.067 a$. The location and magnitude of the maximum value agrees with the full numerical solution to within factors of 1.2 and 4, respectively. This value of $R$ would correspond to half the nearest-neighbor separation in a fluctuation-stabilized colloidal crystal; it also is consistent with the lattice constants observed in superheated metastable colloidal crystallites [2,3]. However, water has a dielectric constant $\epsilon = 80$ which corresponds to a Bjerrum length $\lambda _B = 0.7$ nm at room temperature. Consequently, $\vert \chi_a / \chi_r \vert \ll 1$ for all reasonable colloidal radii. Changing solvents would not affect this conclusion because of the limited range of accessible values of $\epsilon$.

We also investigated the possible influence of multivalent counterions carrying charge $qe$. Ignoring the relatively weak $Z$ dependence of $\chi_a/\chi_r$ yields $\chi_a/\chi_r \sim q^2$. Despite the relative strength of the attraction being larger for higher valency counterions, the effect is still too weak to induce measurable attractions for physically plausible values of $q$.

We have demonstrated that suppression of monopolar ionic fluctuations by surfaces induces a long ranged attraction remniscent of Casimir attractions [14]. This interaction is distinct from and complementary to attractions driven by multipolar fluctuations which have been studied elsewhere [20]. Neither mechanism, however, accounts for the strong and long ranged attractions observed experimentally between highly charged colloidal spheres.

The long-ranged like-charge attractions observed in confined colloid are not consistent with mean field theories for electrolyte structure. Possible explanations must incorporate mechanisms such as fluctuations and high-order correlations not captured by mean field theory. While multipole fluctuations in the distribution of simple ions induce strong attractions [20], they are short ranged. The present study demonstrates that radially symmetric fluctuations can induce long-ranged attractions, but that they are too weak to influence colloidal behavior. Consequently, the explanation must lie in another mechanism not yet considered and thus remains an important outstanding challenge.

We are grateful to Tom Witten, Stuart Rice, Phil Pincus, Mehran Kardar, Ramin Golestanian, An-Chang Shi and Adrian Parsegian for enlightening conversations. This work was supported by the MRSEC Program of the National Science Foundation under Award Number DMR-9808595.