Introduction to statistical mechanics of charged systems

Introduction to statistical mechanics of charged systems

Yan Levin

Instituto de Física, Universidade Federal do Rio Grande do Sul, Caixa Postal 15051, CEP 91501-970, Porto Alegre, RS, Brazil

ABSTRACT

The paper is the summary of lectures given in São Carlos, Brazil during the 2004 Summer School on Statistical Mechanics. My objective was to provide the students with some basic tools necessary to study the thermodynamics of Coulomb systems. I have restricted myself to simple models and techniques, which nevertheless, when used correctly can give a clear insight into the fundamental physics behind various complex phenomena that appear when the interactions between the system's constituents are dominated by the long ranged Coulomb force.

1 Introduction

Electrostatic interactions are ubiquitous. Yet, our understanding of the thermodynamics of these systems is far from complete. Even such seemingly simple question as whether or not a symmetric electrolyte or plasma can undergo a liquid-gas phase separation has been conclusively resolved only very recently [1,2,3]. Even though, the universality class of the phase transition still remains a source of ongoing debate [4]. For strongly asymmetric system such an aqueous colloidal suspension with monovalent salt even the existence of the liquid-gas phase transition still remains unsettled [5]. For two dimensional plasma, in which ions interact through a logarithmic potential, in addition to the liquid-gas phase separation, one also finds a metal-insulator transition, commonly designated as the Kosterlitz-Thouless transition [6]. The importance of this lies in the fact that the phase transitions in many two dimensional models can be mapped directly onto the metal-insulator transition of the two dimensional Coulomb gas. An example of these are: the roughening transition of a crystal interface [7], superfluid ⁴He films, two dimensional crystalline solids, etc [8].

Besides the question of thermodynamic stability, statistical mechanics of Coulomb systems can lead to a number of surprising conclusions [5]. Thus it is actually possible for two like-charged colloidal particles inside a suspension containing electrolyte to attract one another. The mechanism of attraction is purely electrostatic and is not due to some other, yet unknown, force. Another curious finding, is that the electrophoretic mobility of a highly charged colloidal particle can actually become reversed, if the suspension contains multivalent counterions. Thus, if the electric field is applied to the suspension, the particle will drift in the direction opposite to the one expected based on their bare charge.

The goal of the mini-course was to provide the students with the basic tools needed to understand the thermodynamics of charged systems. Because of the time constraint, I have restricted the presentation to simple models and techniques, only mentioning in passage the more advanced approaches, such as integral equations and field theories.

During the preparation of the mini-course I have relied heavily on the recent review which I wrote for the Institute of Physics Publishing, entitled Electrostatic correlations: from plasma to biology, Reports on Progress in Physics, 65 1577-1632, (2002).

2 Ideal gas

Lets begin by considering the gas of N non-interacting particles confined to a box with dimensions L × L × L. The Hamiltonian for this system is

The wave function y(x₁,x₂...x_N) satisfies the Schroedinger equation

the solution to which can be written in the form of a product

Quantum mechanical particles are indistinguishable, so that the wave function should be symmetric (for bosons) or antisymmetric (for fermions) under the permutation of indices. Eq. (3) must, therefore, be suitably symmetrized/antisymmetrized following the usual quantum mechanical procedure.

Substituting Eq.(3) into the Schroedinger equation (35) we find that functions f_j(x) satisfy the Helmholtz equation

where the eigenvalue k is determined from the boundary conditions. If we suppose that the probability of finding a particle outside the box is zero, then the solution to Eq.(4) has the form

where to simplify the notation we have dropped the particle index j. The eigenvalue , where

with i = {x,y,z} and {n_i} are the integer labels of the quantum states of a particle. For a given distribution of particles among the quantum states, the energy in the box is

The partition function

where b = 1/k_BT and Tr is the sum over all the possible quantum states of the particles in the box. The partition function can be rewritten as

where we first perform the trace over all the accessible quantum states of a particle j, and only then perform the product over all the particles. Furthermore, since all the N particles are indistinguishable and have the same accessible quantum states

where

with possible values of k_i given by Eq.(6). The trace over the quantum states can be done by transforming the sum into an integral,

The integration can now be easily performed by going to the spherical coordinates yielding,

where V = L³ is the volume of the box and

is the thermal de Broglie wavelength. The canonical partition function for N non-interacting particles is then

The Helmholtz free energy for the ideal gas is

which with the help of the Stirling approximation reduces to

where r is the density of the particles in the box. Given the free energy, all the thermodynamic functions can be easily calculated. For example the internal energy is

and the pressure is

For interacting particles the canonical partition function can be written as

where

and H is the potential energy of interaction between the particles.

3 Symmetric electrolyte

Consider the simplest model of an electrolyte solution: N ions idealized as hard spheres of diameter a carrying charge ±q at their center confined to volume V. The charge neutrality of solution requires that N₊ = N_– = N/2. The solvent will be modeled as a continuum of dielectric constant Î.

The Hamiltonian for this system is

where x_i,j = |x_i – x_j|, and v(x_i,j) is the potential of hardcore repulsion between particles i and j. Following the common nomenclature, we refer to the ''Hamiltonian'' where in reality we mean only the potential part of the total energy. The notation is a form of shorthand, since the integration over momenta in the partition function completely decouples from the integration over the positions of the particles, and can be performed explicitly. Thus, the free energy of an interacting system can always be written as a combination of the free energy of an ideal gas, plus the excess contribution coming from the interactions between the particles.

For a bulk, charge neutral system, the electrostatic potential at any point x inside the electrolyte is constant and can be set to zero. This means that the mean-field contribution to the free energy vanishes and the excess free energy is due to positional correlations between the ions of the electrolyte. The electrostatic free energy, F^el = F – F₀, (F₀ is the free energy when all the electrostatic interactions are turned off) follows from,

where

Lets define

as the potential that ion j feels due to interaction with all other ions located at {x_i}. The average electrostatic potential felt by ion j is

and we see that

Using the definition of the Helmholtz free energy, Eq. (27) can be written as

Integrating Eq. (28), the electrostatic free energy inside the electrolyte is

This equation is known as the Debye charging process [9].

To calculate the electrostatic contribution to the Helmholtz free energy, let us fix one ion of charge +q at the origin r = 0 and see how the other ions distribute around it, see Fig. 1. Inside the region 0 < r < a there are no other charges except for the one fixed at the origin, and the electrostatic potential f(r) satisfies the Laplace equation,

For r > a the electrostatic potential satisfies the Poisson equation

where the charge density can be expressed in terms of the charge-charge correlation functions g₊₊(r) = g_––(r) and g_+–(r) = g_–+(r)

The average densities of positive and negative ions are r₊ = N₊/V, r_– = N_–/V; r₊ = r_– = r/2.

The correlation functions can be written in terms of the potential of mean force w_ij

where b = 1/k_BT. The w_ij(r) is the work required to bring ions i and j from infinity to separation r inside the electrolyte solution. In their paper Debye and Hückel [9] made an implicit approximation of replacing the potential of mean force by the electrostatic potential

where q_j is the charge of j'th ion and f_i(r) is the electrostatic potential at distance r from the ion i fixed at the origin r = 0. With this approximation, Eq.(31) reduces to the non-linear Poisson-Boltzmann equation (PB),

Debye and Hückel proceeded to linearize this equation. Technically, linearization is only valid if bqf 1, however, being practically minded Debye and Hückel linearized first and worried about the consequences later. As was noted later by Onsager, linearization of Eq.(35) is a necessary step in order to produce a self-consistent theory[10]. The linearized Poisson-Boltzmann equation reduces to the Helmholtz equation

where the inverse Debye length is

The Laplace equation (30) for r < a and the Helmholtz equation (36) for r > a must be integrated, subject to the boundary condition of continuity of the electrostatic potential and the electric field across the boundary surface r = a. For r < a the electrostatic potential is found to be

while for r > a,

Equation (39) shows that the electrostatic potential produced by the central charge is exponentially screened by the surrounding ionic cloud. Because of the hardcore repulsion the screening, however, appears only at distances larger than r = a. This accounts for the presence of geometric factor q(ka) in Eq. (39). The screening of electrostatic interactions inside the electrolyte solutions and plasmas is responsible for the existence of thermodynamic limit in these systems with extremely long range forces.

The electrostatic potential f _< (r), Eq. (38), consists of two terms: the potential produced by the central ion q/Îr, and the electrostatic potential induced by the surrounding ionic cloud,

The electrostatic free energy can now be obtained using the Debye charging process Eq. (29). While performing the charging, it is important to remember that k(lq) = lk(q) . Defining the free energy density as f = F/V, the integral in Eq. (29) can be performed explicitly yielding

For large dilutions Eq. (41) reduces to the famous Debye limiting law,

Given the free energy, the limiting laws for the osmotic pressure and activity can be easily found [2].

The free energy is not analytic at r = 0. The singularity at r = 0 is a consequence of long-range Coulomb interactions, which also manifest themselves in the divergence of the standard virial expansion [11]. The total free energy of the electrolyte, F, is the sum of electrostatic Eq. (41), and entropic contributions. The entropic contribution to the free energy arises from the integration over the momentum degrees of freedom in the partition function Eq. 20, and is equivalent to the free energy of an ideal gas,

where the de Broglie thermal wavelength is given by Eq.(14).

The osmotic pressure of the electrolyte is

which can also be expressed in terms of the Legendre transform of the negative free energy density –f [2],

where the chemical potential is

It is a simple matter to see that below the critical temperature T_c the total free energy F = F^ent + F^el fails to be a convex function of the electrolyte concentration. This implies the presence of a phase transition. Alternatively the phase separation can be observed from the appearance of a van der Waals loop in the osmotic pressure Eq. (45), below the critical temperature T_c. The critical parameters are determined from

The coexistence curve can be obtained using the standard Maxwell construction. It is convenient to define the reduced temperature and density as T^* = k_BTaÎ/q² and r^* = ra³. The critical point of the plasma, within the DH theory, is found to be located at [1,2]

and

It is interesting to note that at criticality k = 1/a. This means that in spite of a very low concentration of electrolyte at the critical point, the screening remains very strong. We also observe that the reduced critical temperature for the electrolyte is almost an order of magnitude lower, than for systems in which the particles interact by the short-ranged isotropic potentials. Since the critical point within the DH theory occurs at extremely low density, we are justified in neglecting the excluded volume contribution to the total free energy.

Phase separation of an electrolyte or of a two component plasma is the result of an electrostatic instability arising from the strong positional correlations between the oppositely charged ions. This mechanism is very different from the one driving the phase separation in systems dominated by the short ranged isotropic forces. In that case the thermodynamic instability is a consequence of the competition between the interparticle attraction and the hardcore repulsion.

The reduced temperature can be written as T^* = a/l_B, where l_B = q²/k_BTÎ is the Bjerrum length. For water at room temperature l_B » 7 Å. This means that one would need ions of size less than 0.4 Å, in order to observe phase separation at room temperature. This is clearly impossible since the minimum hydrated ionic size is about 2 – 4 Å. Therefore, in order to see phase separation, one is required to look for solutions with l_B on the order of 40 Å or more. For water such large values of l_B correspond to temperatures well below the freezing. An alternative is to work with organic solvents which have dielectric constants significantly lower than water. This was the strategy adopted by K.S. Pitzer in his studies of ionic criticality [12,13,14]. Pitzer used liquid salt triethyl-n-hexylammonium triethyl-n-hexylboride (N₂₂₂₆B₂₂₂₆) in the diphenyl ether. With this he was able to observe the critical point at room temperature. Pitzer's work has provoked a lot of stimulating controversy because his measurements suggested that the Coulombic criticality belonged to a new universality class [15]. At first sight this might not seem very surprising, after all the Coulomb force is extremely long ranged. On further reflection the situation is not so clear. Although the bare interaction potential between any two ions is long ranged, inside the electrolyte solution it is screened by the surrounding particles, as is seen from Eq. (39). The effective interaction potential, therefore, is short ranged, which should place the ionic criticality firmly in the Ising universality class. In fact all the theoretical arguments lead to this conclusion, which seems to be contradicted by the Pitzer's experiments. In principle, it is possible that one has to be very close to the critical point before the Ising behavior becomes apparent. However, even this conclusion is hard to justify theoretically. Estimates of the Ginzburg criterion suggest that the width of the critical region for the Coulombic criticality should be comparable to that of systems with short ranged isotropic interactions [1,16]. The situation remains unclear.

An alternative to working with electrolyte solutions is to study molten salts, which are classical two component plasmas. In this case the dielectric constant can be taken to be that of vacuum, and ions are no longer hydrated. The reduced critical temperature T^* = 1/16 and the characteristic ionic diameter of about 2 Å, imply that at criticality l_B » 30, which means that the critical point for a molten salt is located at about 5000K. It is, indeed, very hard to study critical phenomena at such high temperatures! It seems, therefore, that we are stuck with the low dielectric solvents. An alternative is the computer simulations, which are becoming sufficiently accurate to allow measurements of the critical exponents, at least for symmetric 1:1 electrolytes. Indeed the most recent simulations suggest that the Coulombic criticality belongs to the Ising universality class [4].

4 The Bjerrum association

The DH theory presented in the previous section was based on the linearization of the Poisson-Boltzmann equation. In view of the strong screening and the rapid decrease of the electrostatic potential away from the central ion, such a linearization can be justified at intermediate and long distances. It is clear, however, that the linearization strongly diminishes the weight of configurations in which two oppositely charged ions are in a close proximity. Linearization underestimates the strength of electrostatic correlations which result in dipole-like structures. At low reduced temperatures characteristic of the critical point, these configurations should be quite important and must be taken into account. One way of doing this, while preserving the linearity of the theory, is to postulate existence of dipoles whose concentration is governed by the law of mass action. In the leading-order approximation the dipoles can be treated as ideal non-interacting specie [17,18,19]. The total number of particles N = rV is then subdivided into monopoles N₁ = r₁V and dipoles N₂ = r₂V. The particle conservation requires that, N = N₁+2N₂. The free energy of the mixture is F = + + F^el, where F^el and are the entropic and the electrostatic free energies of the monopoles, given by the Eqs. (41) and (43), but with N ® N₁ and r ® r₁. The entropic free energy of dipoles is,

where the internal partition function of a dipole is,

At low temperatures, the precise value of the cutoff R at which the two ions can be considered to be associated is not very important. Following the original suggestion of Bjerrum[17] we can take this value to be the inflection point of the integral in Eq. (52), R_Bj = l_B/2. This choice corresponds to the minimum of integrand in Eq. (52), which in turn can be interpreted as the probability of finding two oppositely charged ions at the separation r. The minimum then correspond to a liminal between bound and unbound configurations. A much more careful analysis of the dipolar partition function has been carried out by Falkenhagen and Ebeling based on the resummed virial expansion [19]. They found that that the low temperature expansion of the Bjerrum equilibrium constant is identical to the equilibrium constant which can be constructed on the basis of the resummed virial expansion. Since we are interested in the low temperature regime where the critical point is located, the Bjerrum equilibrium constant, z₂º z₂(R_Bj), will be sufficient.

It is important to keep in mind that at this level of approximation the electrostatic free energy is only a function of the density of free unassociated ions r₁, since the dipoles are treated as ideal non-interacting specie. The concentration of dipoles is obtained from the law of mass action

where the chemical potential of a specie s is

Substituting the expression for the total free energy into the law of mass action leads to

where the excess chemical potential is m^ex = ¶f^el/¶r₁. The critical point can be located from the study of the convexity of the total free energy as a function of ion concentration r. There is, however, a simpler way [2]. We observe that at Bjerrum level of approximation, dipoles are ideal non-interacting specie. This means that they are only present as spectators and do not interact with the monopoles in any way. This implies that only the monopoles can drive the phase separation. Thus, at the critical point the temperature must still be = 1/16 and the density of monopoles must still remain = 1/64p, as in the case of the pure DH theory. The corresponding density of dipoles at criticality is then given by Eq. (55), with = 1/16 = 0.0625 and = 1/64p = 0.00497. We find that at the critical point the density of dipoles is » 0.02. In the vicinity of the critical point there are many more dipoles than monopoles, / » 4. Within the Bjerrum approximation the non-linear correlations, in the form of dipoles, do not affect the critical temperature, but strongly modify the critical density, = + 2 = 0.045. In spite of the crudeness of approximations, the location of the critical point agrees reasonably well with the Monte Carlo simulations [20,3,21], = 0.051 and = 0.079. The coexistence curve, however, is found to have an unrealistic ''banana'' shape [2]. To correct this deficiency one must go beyond the ''ideal'' dipole approximation and allow for the dipole-ion interaction [1,2]. Most of the fundamental physics of electrostatic correlations, however, is already captured at the level of the Bjerrum approximation.

5 Plasma in d-dimensions

The theory presented above can be easily extended to arbitrary dimensions [22]. Specifically in d-dimensions the Poisson equation becomes,

where

is the surface area of sphere in d-dimensions: C₂ = 2p, C₃ = 4p, C₄ = 2p², etc.

As before, we shall approximate the potential of mean force by the electrostatic potential and then linearize the Boltzmann factor. Fixing one ion at r = 0, the electrostatic potential for distances r < a satisfies the Laplace equation

while for r > a the potential satisfies the Helmholtz equation Ѳf = k²f, with the inverse Debye length now given by

where r = r₊ + r_–. These equations must be solved subject to the boundary conditions of the continuity of the electrostatic potential and the electric field across the excluded volume region. For r < a, Eq. (59) can be easily integrated yielding

For r > a the electrostatic potential is

where G(r) is the solution of

The integration constants y and A are determined from the boundary conditions. Taking the Fourier transform of Eq. (62) we obtain,

To perform the integration we rewrite G as

Interchanging the limits of integration,

The integral over q can now be done easily since it involves only integration of decoupled Gaussians,

The above integral can be performed using the modified Bessel functions, yielding

The boundary conditions determine the integration constants to be,

and

The free energy can now be obtained using the Debye charging process

Not withstanding the apparent complexity of the Eq. (70), the integration can be done explicitly using the identities relating Bessel functions of different orders. The electrostatic free energy density for a d-dimensional plasma is found to be [22]

In d = 3, C₃ = 4p,

and the free energy density reduces to the one found earlier, Eq. (41).

It is now possible to study the thermodynamic stability of a general d-dimensional plasma against a liquid-gas phase separation [22]. A particularly interesting case is a plasma in two dimensions, to which we shall now turn our attention.

6 Two-dimensional plasma and the Kosterlitz-Thouless transition

The 2d plasma has attracted much attention over the years because various important physical systems can be mapped directly onto it. Examples include superfluid ⁴He films, two-dimensional crystalline solids, and XY magnets [8]. Although a continuous symmetry can not be broken in two dimensions [23], if the Hamiltonian of a system is invariant under an Abelian group, a finite temperature phase transition is possible. This transition occurs as the result of unbinding of the topological defects or ''charges''. The defect-mediated phase transitions belong to the universality class of the two-dimensional plasma.

Thirty years ago Kosterlitz and Thouless (KT) have presented a renormalization group study of the 2d plasma [6]. They concluded that at sufficiently low temperature, the 2d plasma becomes an insulator. All the positive and negative ions pair-up forming dipoles. The metal-insulator transition was found to be of infinite order, characterized by an essential singularities in the thermodynamic functions. The KT analysis, however, was restricted to the low ionic densities and it is not clear what happens when the concentration of charged particles is increased. To study this we can apply to the 2d plasma the theory developed above. Lets define d = d – 2, then in the limit d ® 0 the free energy in Eq. (71) becomes,

The second term of Eq. (73) diverges in the limit that d ® 0. This divergence of the free energy, however, does not influence the pressure which remains well defined. This is because (ka)² ~ r, which means that the contribution to the electrostatic free energy coming from this term is volume independent, and will vanish when derived with respect to volume. In fact it is quite straightforward to see the origin of the divergence appearing in Eq. (73). Recall that the bare Coulomb potential between two ions in the d-dimensional electrolyte is,

see Eq. (60). In the limit d ® 0, Eq. (74) can be expanded in powers of d yielding

Therefore, the divergence appearing in Eq.(73) can be traced back to the divergence found in Eq. (75). This divergence is easily renormalized away by redefining the point of zero potential. Thus if in d-dimensions the bare electrostatic potential is redefined as [22]

the limit d ® 0 is well defined, and the electrostatic potential in 2d reduces to

Within the Debye-Hückel approximation the free energy for particles interacting through this potential is finite and is given by

The electrostatic contribution to the chemical potential and the pressure can now be calculated. We find

and

where the reduced density and temperature in 2d are defined as: r^* = ra² and T^* = k_BTÎ/q², respectively. In the limit of large dilution r ® 0, the pressure inside the plasma is

We observe that for T^* < 1/4, the pressure at low density of electrolyte is negative, signifying presence of a thermodynamic instability. The critical point is located at = 0 and = 1/4. For < 1/4 the 2d plasma phase separates into two coexisting phases one of which has zero density [22]. This is clearly an artifact of the approximations made.

While the Debye-Hückel theory is sufficient to account for electrostatic correlations on long length scales, due to its linear nature it fails for short distances. On short length scales, the correlations manifest themselves in the form of dipolar pairs of positively and negatively charged ions. As in the case of 3d electrolyte we should, therefore, separate the total density of hard discs into the density of monopoles r₁, and dipoles r₂, so that r = r₁ + 2r₂. At the level of Bjerrum approximation the dipoles are treated as ideal non-interacting particles. The electrostatic free energy then comes only from the interactions between free un-associated ions, and is given by Eq. (78) with

The total free energy density is , where

and

The internal partition function for a 2d dipole is

We note that for low temperatures, T^*< 1/2, the integral in Eq (85) converges uniformly as R ® ¥. In this regime it is possible, therefore, to define the internal partition function of dipole as

The thermodynamic equilibrium requires that for fixed volume and number of particles the Helmholtz free energy be minimum. This is equivalent to the law of mass action Eq. (53), which upon the substitution of free energy simplifies to Eq. (55). In the limit of small concentrations, the excess chemical potential can be expanded in powers of r₁ yielding

where g_E is the Euler constant. Substituting Eq. (87) into Eq. (55), we find that the concentration of dipoles in the limit r ® 0 scales as

where

For T^* < 1/4, the exponent q(T^*) < 0, and in the limit r₁® 0 the law of mass action can not be satisfied. This means that in the temperature density plane (T^*,r^*), for sufficiently small densities, the line T^* = 1/4 corresponds to the critical locus of metal-insulator transitions. Below this line, and for sufficiently small ionic concentrations, no free monopoles can exist. All the ions are paired up into neutral dipolar pairs. The critical line terminates at the tricritical point located at = 1/4 and

For T^* < 1/4 and r^* > there is a phase separation between an insulating vapor and a conducting liquid phases, Fig. 2.

As the critical line is approached from high temperatures, the Debye length diverges as

where

and

The critical exponent is n = 1. The Kosterlitz-Thouless (KT) renormalization group theory [6] predict the same behavior for x_D except that n = 1/2. The KT theory, however, leaves unanswered the question of what happens to the metal-insulator transition for higher plasma concentrations. The current theory, on the other hand, shows that the critical line terminates in a tricritical point, after which the metal-insulator transition becomes first order [22]. This topology is also consistent with the findings of Monte Carlo simulations [24,25]. A more sophisticated theory introduced by Minhagen [26], leads to a very similar phase diagram, except that the tricritical point is replaced by a critical end-point.

We see that the electrostatic correlations are even more important in 2d than in 3d. While in three dimensions the electrolyte phase separates into the coexisting liquid and gas phases, both of which contain monopoles and dipoles, in two dimensions the low density vapor phase does not contain any free charges and is an insulator.

7 Confined one component plasma

In 1971 Crandall and Williams suggested that electrons trapped on the surface of liquid helium ⁴He can crystallize, forming a two dimensional Wigner crystal [27]. Eight years later this order-disorder transition was observed experimentally by Grimes and Adams [28]. In this system electrons obey the classical mechanics, since the Fermi energy is much smaller than k_BT. Similar crystallization can occur in the inversion layer near the surface of a semiconductor, however, in this case the quantum effects are important and the electrons form a degenerate quantum gas [29].

The trapped electrons above the liquid ⁴He can be modeled as a confined quasi-two-dimensional plasma of particles interacting by 1/r potential. This model is also appropriate for the study of correlations between the condensed counterions on the surface of colloidal particles.

The average spacing between the confined electrons is d = (ps)^–1/2, where s is the average surface density, s = N/A. The dimensionless quantity parameterizing the strength of electrostatic interactions is G = q²/Îk_BTd. For an infinitesimally thin layer separating two mediums of dielectric constants Î₁ and Î₂, the important parameter is the average dielectric constant Î = (Î₁ + Î₂)/2. It has been observed in computer simulations [30] that the 2d OCP crystallizes into triangular Wigner crystal for G > 130. This value is also in a close agreement with the experiments of Grimes and Adams.

We can gain much insight into thermodynamics of 2d OCP using the, now familiar, Debye-Hückel theory. Our model consists of a plasma of point particles of charge q and of a neutralizing background, confined to an interface located at z = 0 between the two dielectric half-spaces. For z < 0 the dielectric constant is Î₁ and for z > 0 the dielectric constant is Î₂. Since the half-spaces do not contain any free charges, the electrostatic potential everywhere satisfies the Laplace equation Ѳf = 0. The electrostatic free energy is obtained by fixing one particle and calculating the induced potential resulting from the redistributions of other ions in the z = 0 plane. It is convenient to adopt the cylindrical coordinate system, (, j, z), so that the fixed ion is located at = 0, z = 0. Using the azimuthal symmetry and the fact that the electrostatic potential vanishes at infinity, the solution to Laplace equation can be written as [31]

and

where J₀(x) is the Bessel function of order zero.

The functions A₁(k) and A₂(k) can be determined from the boundary conditions, which are: continuity of the electrostatic potential,

and discontinuity of the displacement field across the z = 0 plane. The discontinuity results from the inhomogeneous distribution of interfacial charge induced by the fixed ion,

From charge neutrality the average interfacial charge is zero, so that s_q () is the result of ionic correlations,

The first term of Eq. (98) is the surface charge density of the fixed ion, the second term is due to the uniform negative background, while the last term is the surface charge density of ions confined to the interface. We have, once again, approximated the potential of mean force by the electrostatic potential. In the spirit of Debye-Hückel theory we shall now linearize the Boltzmann factor. The surface charge density becomes

where

is the Gouy-Chapman length.

The continuity of the electrostatic potential requires that A₁(k) = A₂(k). Substituting Eqs. (94) and (95) into Eq. (97) and using Eq. (99), we find the electrostatic potential over the full range –¥ < z < ¥ to be

For z = 0 the integral can be performed explicitly yielding

where the functions t_n(x) are defined as [32],

with H_n(x) and N_n(x) being the Struve and the Bessel functions of order n, respectively. For large values of x, t₀» 1/x², so that asymptotically,

We conclude that in the case of a confined plasma there is no exponential screening, instead the electrostatic potential is purely algebraic and has the form of a dipole-dipole interaction [33].

There is a well known argument in the condensed matter physics going back all the way to Bloch [34], Peierls [35] and Landau [36] in the 1930's, which states that a continuous symmetry can not be broken in two dimensions. This means that there can not exist a true two dimensional crystalline order, since it requires breaking translational symmetry. The argument was made rigorous by Mermin, who proved it for particles interacting by short-ranged potentials [23]. It is quite simple to see how this conclusion arises. Suppose that there is a 2d crystal, one can then calculate the mean-square displacement d² of one particle from its equilibrium position due to thermal fluctuations. It is found that d² ~ T ln L, where L is the characteristic crystal size. For L ® ¥, the mean square displacement diverges for any finite temperature, implying that in thermodynamic limit a 2d crystal is unstable to thermal fluctuations. Although there is no true long-range order in two dimensions for systems with short-range forces, there exists a pseudo-long-range order characterized by an algebraically decaying correlation functions. It is not clear, however, to what extent this conclusion applies to the 2d OCP, whose particles interact by a long-ranged 1/r potential. Certainly in this case Mermin's proof is no longer valid. However, since the effective interaction potential inside a 2d OCP decays as 1/r³, which is short-ranged in two dimensions, suggests that there should not be any long-range order. Whether there is a true long-range order or a pseudo-long-range order for a 2d OCP remains uncertain. Simulations find that for G » 130 there is a crystallization transition. It is, however, difficult to say whether the crystalline state has a true long-range order or a pseudo-long-range order [30]. It is also unclear if the transition is of first order or continuous, belonging to the Kosterlitz-Thouless universality class [6,8]. Existence of the thermodynamic limit for confined 2d plasmas can also be attributed to the effective renormalization of the interaction potential from a non-integrable 1/r to integrable (in two dimensions) 1/r³ form.

The Helmholtz free energy of a 2d plasma can be obtained directly from Eq. (102). We need to know the induced potential felt by the central ion due to other particles. In the limit ® 0, the electrostatic potential reduces to

The first term of this expression is the potential produced by the central ion, while the second term is the induced potential felt by the fixed ion. We note that the induced potential is actually divergent in the limit ® 0. This is a consequence of the failure of linearization of the Poisson-Boltzmann equation. The linear equation allows for a much closer mutual approach between two ions, than is actually possible. A suitable cutoff must, therefore, be introduced into the linear theory to to account for the omitted non-linearity of the Boltzmann factor. It is reasonable to postulate that the distance of closest approach h between the two ions is such that their electrostatic energy of repulsion is comparable to their characteristic kinetic (thermal) energy,

This means that h » l_B. We can use this value as the short-distance cutoff in the calculation of free energy. The induced potential then becomes

The free energy is obtained through the usual Debye charging process, Eq. (29). Recalling that l_B(lq) = l²l_B(q) and l_GC(lq) = l_GC(q)/l², where l is the charging parameter, in the limit of high temperatures G ® 0, the reduced free energy per particle is found to be

Eq. (108) is precisely the leading order term of the resumed virial expansion obtained by Totsuji [37,38].

For low temperatures, the OCP crystallizes into a triangular lattice. The Madelung energy of this lattice is,

This equation provides a surprisingly good fit not only for the free energy of solid, but also for the free energy of fluid at sufficiently high values of G. Comparing to the results of the Monte Carlo simulations [30] we find that for G = 5 the error accrued from using Eq. (109) to calculate the total electrostatic free energy is about 30%. For G = 20 this error drops to 11% and for G = 50 it goes down to 6%. Recalling that the crystallization transition occurs at G » 130, we see that the Eq. (109) works well into the fluid phase. It is reasonable, therefore, to approximate the electrostatic free energy of a fluid for G > 5 by

The reason why the electrostatic free energy of a fluid is so well approximated by the free energy of the crystal, is a consequence of strong electrostatic correlations.

8 Asymmetric systems

Up to now we have considered only symmetric plasmas and electrolytes. In practice, however, it is unlikely that both cations and anions have exactly the same size and magnitude of charge. It is, therefore, important to explore the thermodynamics of a general Z : 1 electrolyte in which cations have charge Zq and diameter a_c, while anions have charge –q and diameter a_a. Unfortunately, as soon as the asymmetry is introduced, the internal inconsistency enters into the Poisson-Boltzmann equation [10]. Recall that the cation-anion correlation function can be expressed in terms of the potential of mean force w_+–

The w_+–(r) is the work needed to bring cation and anion from infinity to separation r inside an electrolyte. Clearly this work is invariant under the permutation of particle labels w_+–(r) = w_–+(r). This means that

The Poisson-Boltzmann equation, which serves as the basis, for the Debye-Huckel theory, approximates the potential of mean force by w_+–(r) = q_–f₊(r). The self consistency condition, Eq. (112), then requires that

Because of the non-linear nature of the PB equation this condition can not be satisfied except for symmetric electrolytes. The linearization prescription intrinsic to the Debye-Hückel theory allows Eq. (113) to hold for ions of different valence, but with the same ionic diameter, a_c = a_a.

We see that as soon as the symmetry between the cations and anions is broken the physics and the mathematics of the problem becomes significantly more complex. In the limit of very large asymmetries, Z ® ¥ and a_c a_a a new simplification, however, enters into the game.

9 Colloidal suspensions

A typical colloidal suspension often studied experimentally consists of polystyrene sulphonate spheres of diameter 10nm – 1mm and 10³ – 10⁴ ionizable surface groups. Because of the large surface charge, the colloidal particles tend to repel each other, forming crystals, even at fairly low volume fraction of less than 10%. Using the periodic structure of the lattice, the thermodynamics of a colloidal crystal can be studied fairly straightforwardly. Each colloidal particle can be thought to be confined to a Wigner-Seitz (WS) polyhedral cell. A further approximation replaces the polyhedral WS cell by a sphere [39].

9.1 Colloidal lattices

We shall model the colloidal particles as hard spheres of radius a carrying Z ionizable groups of charge –q distributed uniformly on the surface. The counterions will be idealized as point particles of charge +q. The suspension of N_p = r_pV polyions and N_c = ZN_p = r_cV counterions is confined to a volume V. As usual, the solvent will be treated as a uniform continuum of dielectric constant Î. For sufficiently large polyion concentrations colloidal suspension crystallizes. Using the lattice symmetry, we restrict our attention to one colloidal particle and its counterions inside a spherical WS cell of radius R such that

The osmotic pressure inside the WS cell is

where the canonical partition function is

Using the chain rule and the fact that all the particles are identical,

On the other hand, we recall that the concentration of counterions at the WS cell boundary is

Substituting Eqs. (117) and (118) into Eq. (115) we obtain [40],

The thermodynamics of a crystalline colloidal suspension now reduces to the calculation of the distribution of counterions inside the WS cell. This can be achieved using a simple mean-field picture. The electrostatic potential inside the WS cell satisfies the Poisson equation (31), with the counterion charge density approximated by the normalized spherically symmetric Boltzmann distribution,

The non-linear Poisson-Boltzmann equation can be solved numerically to yield the electrostatic potential and the distribution of counterions inside the cell. In practice it is more convenient to work with the electric field

The Poisson equation can then be rewritten as

where q_p(r) is the polyion charge density,

To simplify the calculations we have uniformly smeared the charge of the polyion over its surface. Integrating both sides of Eq. (122) and taking advantage of the divergence theorem, the electric field at distance r from the polyion is

where

is the counterion charge inside a sphere of radius r centered on the colloidal particle. Using the gauge in which f(a) = 0 the electrostatic potential is

and the Poisson-Boltzmann equation reduces to an integral equation for the electric field. Note that Eq. (124) naturally incorporates the boundary conditions

and

Eq. (124) can be solved iteratively to yield the counterion density profile from which all other thermodynamic functions are straightforwardly determined.

For aqueous colloidal lattices with monovalent counterions, the Poisson-Boltzmann equation is in excellent agreement with the experiments and simulation. The PB equation, however, does not account for the correlations between the counterions and breaks down for low dielectric solvents or for aqueous suspensions with multivalent ions. Fortunately, in the case of colloidal lattices, it is fairly straightforward to account for these effects using the density functional theory.

At lower concentrations, when the crystalline structure has melted the situation, unfortunately, is no longer so clear cut. In this case a simple picture based on the Wigner-Seitz cell is not sufficient and new methods must be developed [41,42]. Unfortunately the standard techniques of the liquid state theory based on integral equations are powerless in the case of highly asymmetric colloidal systems. The field theoretic methods also fail when applied to this difficult problem [5]. Furthermore, even the experimental situation is far from clear. Ise and coworkers claim to have seen stable clusters of colloidal particles in highly deionized colloidal suspensions. Tata et al. even report an observation of a full equilibrium vapor-liquid-like phase separation [43]. These experiments, however, have been challenged by Palberg and Würth, who demonstrated that the phase separation observed by Tata et al. was the result of non-equilibrium salt gradients produced by the ion exchange resin [44,45]. In the colloidal science community the possibility of a liquid-vapor phase separation in highly deionized colloidal suspensions has met with a large amount of scepticism. The usual argument against the phase transition is based on the Derjaguin-Landau-Verwey-Overbeek (DLVO) colloidal pair potential [46,47].

It is easy to understand the nature of the DLVO potential based on the Debye-Hückel theory. If the size of colloidal particles is shrunk to zero, a ® 0, then due to screening by counterions, the interaction energy between two ''point'' colloids would be of a Yukawa form,

where the inverse Debye length is

Now, consider the electrostatic potential outside the fixed colloidal particle of radius a and charge –Zq, Eq. (39)

Evidently the factor q(ka) accounts for the fact that screening starts only outside the cavity, r > a. We also can think of Eq. (131) as the potential of a point particle with an effective charge Q_p = Zqq(ka). An advantage of this alternative point of view is that the interaction energy for two ''point'' particles is simply given by Eq. (129) with Zq ® Q_p. This leads directly to the famous DLVO potential

This potential is purely repulsive [48], which naively suggests that a charged colloidal suspension is stable against a liquid-gas phase separation. Sogami and Ise, therefore, have argued that the DLVO potential must be incorrect, since it cannot account for the inhomogeneities observed experimentally [49]. In its stead, they proposed a different interaction potential derived on the basis of the Gibbs free energy. The potential found by Sogami and Ise contains a minimum [49], which implies that at short enough separations the two like-charged colloidal particles attract! What is most surprising is that the attraction appears even for monovalent counterions, i.e in the absence of strong correlations between the colloidal double layers. Furthermore, water is an incompressible fluid so that it is difficult to see how a change of paradigm from Helmholtz to the Gibbs free energy can lead to such a profound modification of the interaction potential. Inconsistency in the results based on the Helmholtz and the Gibbs free energies has been carefully reexamined by Overbeek, who has traced the discrepancy to a flaw in the Sogami and Ise's calculations [50].

It is important to stress that the repulsive two-body interactions do not, in general, preclude the possibility of a liquid-gas phase separation in a multicomponent fluid. In fact van Roij and Hansen found, within the linearized density functional theory, that it is possible for a colloidal suspension with polyions interacting by the repulsive DLVO potential to phase separate into coexisting liquid and gas phases [51]. Before entering into the discussion of colloidal fluids it is, however, important to introduce a new fundamental concept — the colloidal charge renormalization.

9.2 Charge Renormalization

Although the non-linear Poisson-Boltzmann equation can not be solved analytically for a spherical geometry, the numerical solution indicates that the electrostatic potential far from colloidal particle saturates as a function of the bare colloidal charge [39]. This suggests that the thermodynamics of a highly charged colloidal systems can be based on the linearized PB equation but with the bare colloidal charge replaced by an effective renormalized charge. The original concept of colloidal charge renormalization is due to Alexander et al., but is well predated in the polyelectrolyte literature, where the phenomenon is known as the Manning counterion condensation [52,53,54].

To understand better colloidal charge renormalization, let us first consider a uniformly charged plane at fixed potential y_s inside a salt solution of concentration c. The electrostatic potential at distance x from the plane satisfies the PB equation,

Since at the moment we are considering aqueous suspensions containing only monovalent ions, the electrostatic correlations are insignificant and the mean-field Poisson-Boltzmann approximation is sufficient. Multiplying both sides of Eq. (133) by df/dx allows us to perform the first integration. Since the potential vanishes in the limit x ® ¥, we find

The second integration yields [55]

where the inverse Debye length is,

In the limit of large surface potentials this expression simplifies to

For separations from the plane larger than the Debye length, Eq. (137) becomes

An important observation is that for large surface potentials, bqy_s/4 1, the electrostatics away from the plane is completely insensitive to the surface charge density.

Now, let us consider a highly charged colloidal particle of valence Z and radius a inside a symmetric 1 : 1 electrolyte of concentration c. The electrostatic potential at distance r from the center of a colloidal particle satisfies the PB equation (35). For distances r > a + x_D the electrostatic potential is small and the PB equation can be safely linearized leading to the Helmholtz equation (36). This can be easily integrated yielding the electrostatic potential,

where A is the integration constant. To find its value, lets restrict our attention to suspensions in which the x_D

a. In practice this is not a very strong restriction. For salt solutions at physiological concentrations x_D » 8 Å while the characteristic colloidal size is on the order of 1000 Å. Even for solutions with very low salt content, in the mM range, the Debye length is on the order of 100 Å. Under these conditions all the curvature effects associated with the spherical geometry of colloidal particle are effectively screened at separations a + x_D < r < 2a, and the electrostatic potential is well approximated by that of a uniformly charged plane, Eq. (138). Comparing Eqs. (138) and Eq. (139) the value of the integration constant follows directly, and the electrostatic potential at distance r > a + x_D from the center of colloidal particle is

This is the asymptotic solution of the full non-linear PB equation for ka 1. Comparing this to the solution of linearized PB, Eq. (131), it is evident that the two are identical as long as the bare colloidal charge is replaced by the renormalized charge. For highly charged particles, Eq. (140) shows that the renormalized charge saturates at [56]

While the previous analysis was carried out for one colloidal particle inside an electrolyte solution, the concept of charge renormalization is quite general and can be applied to colloidal suspensions under various conditions [57,57,58,59]. The difficulty of defining the effective charge for suspensions at non-zero concentrations resides in the complexity of accounting for the consequences of colloidal interactions. The standard practice is to to study one colloidal particle inside a spherical Wigner-Seitz cell whose radius is determined by the volume fraction of colloids [39]. While this procedure is fully justified for colloidal lattices, its foundation is less certain for fluidized suspensions. To find the renormalized charge one numerically solves the full non-linear PB equation and matches the electrostatic potential to the solution of the linearized equation at the cell boundary. Alternatively, the osmotic pressures inside the WS cell calculated using the non-linear and linear equations are matched in order to define the effective charge. One should remember, however, that while at the level of non-linear PB equation the osmotic pressure is directly proportional to the concentration of ions at the cell boundary, Eq. (119), this is not the case for the linearized PB equation. The various procedures lead to similar values of the renormalized charge. In the case of salt-free suspensions, the effective charge is found to saturate at [39]

where c is an approximately linearly increasing function of colloidal concentration for suspensions with volume fraction larger than 1%. For suspensions with colloidal volume fraction between 1% and 10% the value of c varies from around 9 to 15 [60,61].

9.3 Colloidal Fluid

In this section we will apply the insights gained from the study of one and two component plasmas to the exploration of stability of charged colloidal suspensions against a gas-liquid phase separations. We note that the large size asymmetry between colloids and counterions leads to very different equilibration time scales. On the time scale of polyion motion, the counterions are always equilibrated. This suggests that the calculation of free energy should be done in two stages [62]. First, we shall trace out the counterion degrees of freedom, leading to effective many-body interactions between the colloidal particles. Then we will use these effective interactions to calculate the colloid-colloid contribution to the total free energy. The procedure is similar to the one used in McMillan-Mayer theory of solutions [63].

We shall first calculate the contribution to the total free energy arising from the polyion-counterion interactions [41,42]. Consider a suspension in thermal equilibrium. While the colloidal particles are more or less uniformly distributed throughout the solution, the positions of counterions are strongly correlated with the positions of polyions. As a leading order approximation we can, therefore, take the polyion-polyion correlation function to be

and the polyion-counterion correlation function to be

Choosing the coordinate system in such a way that it is centered on top of one of the colloidal particles, the electrostatic potential at distance r < a satisfies the Laplace equation, while for distances r > a it satisfies the Poisson equation Eq. (31). Based on Eqs. (143) and (144) the charge density in the region r > a can be approximated by,

In the spirit of the Debye-Hückel theory we shall linearize the exponential [9,2]. The distribution of charge around the colloid reduces to

For r > a the electrostatic potential, therefore, satisfies the Helmholtz equation (36) with k given by Eq. (130). The solution to this equation is

while the solution to the Laplace equation for r < a is

The electrostatic energy due to polyion-counterion interaction is

where r_q(r) is the charge density of counterions given by Eq. (146), and q_p(r) is the charge density of a polyion,

Performing the integration we find

The electrostatic free energy of a polyion inside the suspension is obtained using the Debye charging process [64],

Note that this free energy is the sum of the polyion self energy

and the solvation energy

which the polyion gains from being inside the ''ionic sea''. The electrostatic free energy due to interaction between all the polyions and counterions is

We have effectively integrated out the counterion degrees of freedom. This, however, leaves us with the effective many-body potentials of interaction between the colloidal particles. For dilute suspensions, the pairwise interaction potential should be the dominant one. The two-body interaction potential can be obtained from the solution of Helmholtz equation for two colloidal particles [65,66]. At large separations this leads directly to the DLVO interaction potential, Eq. (132). This potential has been extensively tested experimentally and found to work very well for bulk colloidal suspensions [67]. Since the DLVO potential is short ranged, the contribution to the total free energy arising from the colloid-colloid interaction can be calculated in the spirit of the traditional van der Waals theory, through the second virial term. A more sophisticated calculation of the colloid-colloid free energy relies on the Gibbs-Bogoliubov variational bound,

where the reference system is taken to be the fluid of hard spheres, whose diameter plays the role of a variational parameter. The free energy resulting from the polyion-polyion interaction, F^pp, can be approximated by the lowest variational bound of Eq. (156). The calculation is somewhat involved, so we refer the interested reader to the original papers [68,69,51,70].

The entropic mixing free energy of colloids and their counterions is simply that of an ideal gas,

where L_c and L_p are the de Broglie thermal wavelengths of counterions and polyions, respectively.

The total free energy of colloidal suspension is the free energy needed to solvate colloids in the sea of other polyions and counterions F^pc + F^pp, and the free energy of mixing F^ent,

The osmotic pressure is

It is found that for suspensions with

the pressure is not a convex function of the colloidal concentration, implying existence of a thermodynamic instability. At criticality the colloidal volume fraction is around 1%. The crucial question is whether this result is reliable? In order to calculate the electrostatic free energies, we were forced to linearize the Boltzmann factor. While this is a reasonable approximation at large separations away from the polyions, linearization is clearly invalid in the vicinity of colloidal surface. There, the strong electrostatic interactions result in an accumulation of counterions and the effective polyion charge renormalization. Therefore, the linear theory can be used only if the bare colloidal charge is replaced by the effective renormalized charge, Z ® Z_eff, in all the expressions. It was found, however, that the bare charge does not increase without limit but saturates at the value given by the Eq. (142). Substituting Z ® Z_eff, into the definition of C Eq. (160), we see that C < 15 for all the values of the bare charge Z in the critical region. The critical threshold, therefore, can not be reached [71], meaning that a deionized aqueous suspensions with monovalent counterions is stable against a liquid-gas phase separation for all colloidal charges and sizes. This conclusion has also been confirmed by more detailed calculations and simulations [41,72,73,74,75].

The result that the non-linear terms omitted within the Debye-Hückel approximation stabilize a deionized colloidal suspensions against a liquid-vapor phase separation has also been obtained by von Grünberg et al. [76,77,78] and Tamashiro and Schiessel [79] based on the analysis of the full non-linear Poisson-Boltzmann equation inside a Wigner-Seitz cell. The numerical integration of the non-linear PB shows that the osmotic pressure is a monotonically increasing function of colloidal concentration. This means that at the level of WS approximation suspension is thermodynamically stable. Von Grünberg et al. and Tamashiro and Schiessel, however, demonstrate that the linearized PB equation leads to the negative compressibility and the osmotic pressures for highly charged colloidal particles. This erroneously suggests presence of a thermodynamic instability. Clearly the instability is an artifact of the linearization. Furthermore, our calculations show that any linear theory, which does not take into account the colloidal charge renormalization, is likely to lead to an incorrect prediction of a liquid-vapor phase separation [51,80] in deionized aqueous suspensions with monovalent counterions.

It is curious that the ''linear'' correlations between the colloids and the counterions, responsible for the screening of electrostatic interactions, are also the ones driving the suspension towards the phase separation. On the other hand, the ''non-linear'' correlations responsible for the counterion condensation and the colloidal charge renormalization, stabilize the suspension against the phase transition.

10 Multivalent counterions

Up to now we have been concentrating our attention on aqueous solutions with monovalent counterions. It was already mentioned that in this case the correlations between the condensed counterions can be neglected. To understand why, let us compare the characteristic electrostatic energy of a counterion-counterion interaction to the characteristic thermal energy k_BT,

where a is the counterion valence and d is the average separation between the n condensed counterions on the surface of a colloidal particle of radius a. Since np(d/2)² = 4pa²,

and the coupling strength becomes,

Now, lets consider highly charged latex particles with Z = 7000 and a = 1000 Å, in water at room temperature. From Eq. (142), taking c = 15, = 2100, which means that 4900 monovalent (a = 1) counterions are condensed onto the particle. The coupling strength of the counterion-counterion interaction is then G » 0.13, which clearly shows that the electrostatic interactions between the condensed counterions are very weak. We can make this observation even more general. The high surface charge concentration s_m encountered in nature is on the order of one elementary charge per 100 Ų. Lets suppose that suspension consists of highly charged colloidal particles with surface charge density s_m. Clearly this means that there will be a lot of counterion condensation. For a salt-free colloidal suspension containing multivalent counterions, the number of condensed counterions will be approximately

The radius of a colloidal particle can be expressed as

Substituting Eqs. (164) and (165) into Eq. (163) we find that the maximum counterion-counterion coupling strength is,

For monovalent counterions G_max » 0.65, for divalent counterions G_max » 1.8, and for trivalent counterions G_max » 3.4. Although G_max is an overestimate, it clearly shows that for highly charged colloidal particles, correlations between the condensed multivalent counterions cannot be ignored.

10.1 Overcharging

One consequence of strong electrostatic correlations is the phenomenon known as the ''overcharging'' [81-94,5]. Overcharging occurs as the result of highly favorable gain in electrostatic free energy due to strong positional correlations between the condensed counterions.

To understand better how the overcharging of colloidal particles comes about let us consider a simple case of one colloidal particle with a uniform surface charge –Zq and radius a, at zero temperature [95,96]. The question that we would like to answer is how many a-valent counterions should be placed on top of the colloidal particle in order to minimize the electrostatic energy of the resultant polyion-counterion complex? Naively we might suppose that the number of condensed counterions should be such as to neutralize completely the colloidal charge. This, indeed, would be the case if the charge of counterions was uniformly smeared over the surface of colloid. In reality, the counterions are discrete entities and can gain favorable energy by maximizing their separation from one another. Lets calculate the electrostatic energy of the polyion-counterion complex,

The first term is the self energy of a polyion, the second term is the electrostatic energy of interaction between the polyion and n condensed a-ions, and the last term is the electrostatic energy of repulsion between the condensed counterions. Now, consider a one component plasma of n a-ions on the surface of a sphere of radius a but with a uniform neutralizing background charge –anq. The electrostatic energy of this OCP can be expressed as the sum of contributions arising from the counterion-counterion interaction, counterion-background interaction, and the self energy of the background,

Substituting this expression into Eq. (167), the electrostatic energy of the polyion-counterion complex simplifies to

For low temperatures, the condensed counterions try to maximize their separation from one another. In the planar geometry the ground state corresponds to a triangular Wigner crystal. A similar arrangement of counterions will also be found on the surface of a spherical colloidal particle, up to some topological defects. The electrostatic energy of a planar OCP has been discussed in Section 7. For a spherical OCP the electrostatic energy at zero temperature is

where M is the Madelung constant. At this point it is interesting to make a historical aside. The energy in Eq. (169) is related to the very old and famous problem in the history of modern physics. The question was first posed by J.J. Thomson in the context of his exploration of what is ''atom'' [97]? After his discovery of electron in 1897, Thomson was trying to understand the structure of the periodic table and thought that he could do this if he could figure out the location of the electrons inside atoms. The only problem was that at that time proton was still not discovered, so in order to keep his atom neutral, Thomson proposed that the ''counter-charge'' to electron is uniformly distributed inside the atom. The point-like electrons were allowed to move freely in the interior of the atom, with the ground state determined by the minimum of the electrostatic energy. This became know as the ''plum pudding'' model of an atom: electrons confined to an interior of a sphere with a uniform neutralizing positive background.

The general solution to the Thomson problem is still unknown. In the limit that the background charge vanishes, the atom turns into conductor and all the electrons move to the surface of the sphere. Even this simplified version of the Thomson problem has no general solution, except for a small number of particles — this is one hundred years after its original formulation! The difficulty in determining the ground state configuration of charges on the surface of a sphere is a consequence of an exponentially large number of metastable states, with energies very close to that of the true ground state [98]. Nevertheless, although in general the exact ground state is not available, its energy is very well approximated by the simple expression given by Eq. (169) with Z = 0 [99,100].

Because of a topological difference between a plane and the surface of a sphere, we expect that the Madelung constant will not be exactly the same in the two cases. The difference, however, is found to be very small M_s = 1.102 for the sphere, as opposed to M_p = 1.106 for the planar OCP.

Returning to our charge inversion problem, the effective charge of the polyion-counterion complex in units of –q is

where n^* is the number of condensed a-ions which minimize the electrostatic energy,

The effective charge is found to be

where g is,

We see that the effective charge of the polyion-counterion complex is reversed compared to the bare charge Z of the colloidal particle. The complex is overcharged. For highly charged colloids, the effective charge scales as the square root of the bare charge [101,96],

The analysis above was conducted for one colloidal particle at zero temperature. The charges were placed on top of a sphere without taking into account the energy cost for the charge transfer. This is clear unrealistic, the transfer of charge from a reservoir to the colloidal surface requires an expenditure of work, which in general can be quite large. This will certainly strongly influence the degree of overcharging. Nevertheless, the simple calculation presented above provides us with an insight into how the electrostatic correlations can lead to a charge reversal in a colloidal suspension.

11 Conclusions

The goal of the mini-course was to provide the students with some basic tools necessary to understand the thermodynamics of Coulomb systems. We saw that much of the qualitative and often semi-quantitative insight into these complex systems can be gained from studying simple models, using some fairly simple theories. One finds, however, that it is often difficult to go beyond the ''simple'' theories. As one tries to improve on these theories, the results often become worse instead of better [5]. Thus, what might look like an improvement, from the theoretical stand point, often fails when compared to experiment or Monte Carlo simulations. The advantage of ''simple'' theories, however, is that their structure often allows to correct them in a physically transparent sort of way. A nice example of this is the Bjerrum modification of the Debye-Hkel theory (DHBj), Section 4. The drawback, on the other hand, is that the approximations are uncontrolled. Thus, there is no simple perturbative parameter which allows one to clearly state the ''order'' to which the theory is correct. Nevertheless, the physically based theories are found, in many cases, to work better than the more rigorous field theoretic perturbative expansions [5]. One possible explanation for this is that the DHBj types of theories are intrinsically non-perturbative. This can already be seen from the form of the Debye-Hückel free energy, Eq. (41). The free energy contains terms to all orders in ka. The absence of a good resummation prescription, in the case of perturbative field theories, might be partially responsible for their failure in application to Coulombic criticality.

Received on 10 May, 2004

• [1] M. E. Fisher and Y. Levin, Phys. Rev. Lett. 71, 3826 (1993).
• [2] Y. Levin and M. E. Fisher, Physica A 225, 164 (1996).
• [3] G. Orkoulas and A. Z. Panagiotopoulos, J. Chem. Phys. 110, 1581 (1999).
• [4] E. Luijten, M. E. Fisher, and A. Z. Panagiotopoulos, Phys. Rev. Lett. 88, 185701 (2002).
• [5] Y. Levin, Rep. Prog. Phys. 65, 1577 (2002).
• [6] J. Kosterlitz and D. Thouless, J. Phys. C 6, 1181 (1973).
• [7] S. T. Chui and J. D. Weeks, Phys. Rev. B 14, 4978 (1976).
• [8] D. Nelson, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic Press, NY, 1983), p. 1.
• [9] P. W. Debye and E. Hückel, Phys. Z. 24, 185 (1923).
• [10] L. Onsager, Chem. Revs. 13, 1933 (1933).
• [11] M. E. Fisher, X.-J. Li, and Y. Levin, J. Stat. Phys. 79, 1 (1995).
• [12] K. S. Pitzer, J. Phys. Chem. 99, 13070 (1995).
• [13] T. Narayanan and K. S. Pitzer, J. Chem. Phys. 102, 8118 (1995).
• [14] T. Narayanan and K. S. Pitzer, Phys. Rev. Lett. 73, 3002 (1994).
• [15] S. Wiegand and et al., J. Chem. Phys. 109, 9038 (1998).
• [16] M. E. Fisher and B. P. Lee, Phys. Rev. Lett. 77, 3561 (1996).
• [17] N. Bjerrum, Kgl. Dan. Vidensk. Selsk. Mat.-fys. Medd. 7, 1 (1926).
• [18] W. Ebeling, Z. Phys. Chem.(Leipzig) 238, 400 (1968).
• [19] H. Falkenhagen and W. Ebeling, in Ionic Interactions, edited by S. Petrucci (Academic Press, NY, 1971).
• [20] J. M. Caillol, D. Levewque, and J. J. Weis, J. Chem. Phys. 107, 1565 (1997).
• [21] Q. Yan and J. J. de Pablo, J.Chem.Phys. 111, 9509 (1999).
• [22] Y. Levin, X.-J. Li, and M. E. Fisher, Phys. Rev. Lett. 73, 2716 (1994).
• [23] N. D. Mermin, Phys. Rev. 171, 272 (1968).
• [24] J.-M. Caillol, J. Chem. Phys. 100, 2161 (1994).
• [25] G. Orkoulas and A. Z. Panagiotoupolos, J. Chem. Phys. 104, 7205 (1996).
• [26] P. Minnhagen, Rev. Mod. Phys. 59, 1001 (1987).
• [27] R. S. Crandall and R. Williams, Phys. Letl. A 34, 404 (1971).
• [28] C. C. Grimes and G. Adams, Phys. Rev. Lett. 42, 795 (1979).
• [29] A. V. Chaplik, Sov. Phys. JETP 35, 395 (1971).
• [30] R. C. Gann, S. Chakravarty, and G. V. Chester, Phys. Rev. B 20, 326 (1979).
• [31] J. E. Flores-Mena, M. C. Barbosa, and Y. Levin, Phys. Rev. E 63, 066104 (2001).
• [32] E. S. Velazquez and L. Blum, Physica A 244, 453 (1997).
• [33] M. Baus and J. P. Hansen, Physics Reports 59, 1 (1980).
• [34] F. Bloch, Z. Phys. 61, 206 (1930).
• [35] R. E. Peierls, Ann. Inst. Henri Poincare 5, 177 (1935).
• [36] L. D. Landau, Phys. Z. Sowjetunion II, 26 (1937).
• [37] H. Totsuji, J. Phys. Soc. Jpn. 39, 253 (1975).
• [38] H. Totsuji, J. Phys. Soc. Jpn. 40, 857 (1976).
• [39] S. Alexander et al., J. Chem. Phys. 80, 5776 (1984).
• [40] H. Wennerströn, B. Jönsson, and P. Linse, J. Chem. Phys. 76, 4665 (1982).
• [41] Y. Levin, M. C. Barbosa, and M. N. Tamashiro, Europhys. Lett. 41, 123 (1998).
• [42] M. N. Tamashiro, Y. Levin, and M. C. Barbosa, Physica A 258, 341 (1998).
• [43] B. V. R. Tata, M. Rajalakshmi, and A. Arora, Phys. Rev. Lett. 69, 3778 (1992).
• [44] T. Palberg and M. Wurth, Phys. Rev. Lett. 72, 786 (1994).
• [45] B. V. R. Tata and A. K. Arora, Phys. Rev. Lett. 72, 787 (1994).
• [46] B. V. Derjaguin and L. Landau, Acta Physicochimica (USSR) 14, 1941 .
• [47] E. J. W. Verwey and J. T. G. Overbeek, Theory of the Stability of Lyophobic Colloids (Elsevier, Amsterdam, 1948).
• [48] M. O. Robbins, K. Kremer, and G. S. Grest, J. Chem. Phys. 88, 3286 (1988).
• [49] I. Sogami and N. Ise, J. Chem. Phys. 81, 6320 (1984).
• [50] J. T. G. Overbeek, J. Chem. Phys. 87, 4406 (1987).
• [51] R. van Roij and J. P. Hansen, Phys. Rev. Lett. 79, 3082 (1997).
• [52] G. S. Manning, J. Chem. Phys. 51, 924 (1969).
• [53] F. Oosawa, Polyelectrolytes (Marcel Dekker, NY, 1971).
• [54] G. S. Manning, Q. Rev. Biophys. II 2, 179 (1978).
• [55] W. B. Russel, D. A. Saville, and W. R. Schowalter, in Colloidal Dispersions, edited by G. K. Batchelor (Cambridge University Press, Cambridge, 1989).
• [56] E. Trizac, L. Bocquet, and M. Aubouy, Phys. Rev. Lett. 89, 248301 (2002).
• [57] M. Quesada-Pérez, J. Callejas-Fernández, and R. Hidalgo-Álvarez, Phys. Rev. E 61, 574 (2000).
• [58] A. Fernández-Nieves, A. Fernández-Barbero, and F. J. de las Nieves, Langmuir 16, 4090 (2000).
• [59] A. Fernandez-Nieves, A. Fernandez-Barbero, and F. J. de las Nieves, Phys. Rev. E 63, 041404 (2001).
• [60] M. J. Stevens, M. L. Falk, and M. O. Robbins, J.Chem.Phys. 104, 5209 (1995).
• [61] L. Belloni, Colloids and Surfaces A 140, 227 (1998).
• [62] H. Löwen, J. P. Hansen, and P. A. Madden, J. Chem. Phys. 98, 3275 (1993).
• [63] W. G. McMillan and J. E. Mayer, J. Chem. Phys. 13, 276 (1945).
• [64] R. A. Marcus, J. Chem. Phys. 23, 1057 (1955).
• [65] X.-J. Li, Y. Levin, and M. E. Fisher, Europhys. Lett. 26, 683 (1994).
• [66] M. E. Fisher, Y. Levin, and X.-J. Li, J. Chem. Phys. 101, 2273 (1994).
• [67] J. M. Crocker and D. G. Grier, Phys. Rev. Lett. 73, 352 (1994).
• [68] G. A. Mansoori and F. B. Canfield, J. Chem. Phys. 51, 4958 (1969).
• [69] B. Firey and N. W. Ashcroft, Phys. Rev. A 15, 2072 (1977).
• [70] P. B. Warren, J. Chem. Phys. 112, 4683 (2000).
• [71] Y. Levin, E. Trizac, and L. Bocquet, J. Phys.: Condens. Matter 15, S3523 (2003).
• [72] A. Diehl, M. C. Barbosa, and Y. Levin, Europhys. Lett. 53, 86 (2001).
• [73] P. Linse and V. Lobaskin, J. Chem. Phys. 112, 3917 (2000).
• [74] P. Linse and V. Lobaskin, Phys. Rev. Lett. 83, 4208 (1999).
• [75] P. Linse, J. Chem. Phys. 113, 4359 (2000).
• [76] H. H. von Grünberg, R. van Roij, and G. Klein, Europhys. Lett. 55, 580 (2001).
• [77] M. Deserno and H. H. Grünberg, Phys. Rev. E 66, 011401 (2002).
• [78] R. Klein and H. H. von Grünberg, Pure Appl. Chem. 73, 1705 (2001).
• [79] M. N. Tamashiro and H. Schiessel, J. Chem. Phys. 119, 1855 (2003).
• [80] D. Y. C. Chan, P. Linse, and S. N. Petris, Langmuir 17, 4202 (2001).
• [81] G. M. Torrie and J. P. Valleau, J. Chem. Phys. 73, 5807 (1980).
• [82] E. Gonzales-Tovar, M. Lozada-Cassou, and D. Henderson, J. Chem. Phys. 83, 361 (1985).
• [83] M. Lozada-Cassou, R. Saavedra-Barrera, and D. Henderson, J. Chem. Phys. 77, 5150 (1982).
• [84] R. Messina, E. G. Tovar, M. Lozada-Cassou, and C. Holm, Europhys. Lett. 60, 383 (2002).
• [85] T. T. Nguyen and B. I. Shklovskii, Physica A 293, 324 (2001).
• [86] T. T. Nguyen, A. Y. Grosberg, and B. I. Shklovskii, Phys. Rev. Lett. 85, 1568 (2000).
• [87] A. Y. Grosberg, T. T. Nguyen, and B. I. Shklovskii, Rev. Mod. Phys. 74, 329 (2002).
• [88] T. T. Nguyen, A. Y. Grosberg, and B. I. Shklovskii, Phys. Rev. Lett. 85, 1568 (2000).
• [89] T. T. Nguyen, A. Y. Grosberg, and B. I. Shklovskii, J. Chem. Phys. 113, 1110 (2000).
• [90] M. Tanaka and A. Y. Grosberg, J. Chem. Phys. 115, 567 (2001).
• [91] Y. Levin and J. J. Arenzon, J. Phys. A: Math. Gen 36, 5857 (2003).
• [92] A. Martin-Molina, M. Quesada-Perez, F. Galisteo-Gonzalez, and R. Hidalgo-Alvarez, J. Chem. Phys. 118, 4183 (2003).
• [93] A. Martin-Molina, M. Quesada-Perez, F. Galisteo-Gonzalez, and R. Hidalgo-Alvarez, J. Phys. Condensed. Mat. 15, S3475 (2003).
• [94] M. Tanaka, Phys. Rev. E 68, 061501 (2003).
• [95] R. Messina, C. Holm, and K. Kremer, Phys. Rev. Lett. 85, 872 (2000).
• [96] R. Messina, C. Holm, and K. Kremer, Phys. Rev. E 64, 021405 (2001).
• [97] J. J. Thomson, Philos. Mag. 7, 237 (1904).
• [98] T. Erber and G. M. Hockney, Phys. Rev. Lett. 74, 1482 (1995).
• [99] T. Erber and G. M. Hockney, J. Phys. A 24, L1369 (1991).
• [100] Y. Levin and J. J. Arenzon, Europhys. Lett. 63, 415 (2003).
• [101] B. I. Shklovskii, Phys. Rev. E 60, 5802 (1999).

Publication Dates

• Publication in this collection
  05 Nov 2004
• Date of issue
  Sept 2004

History

• Received
  10 May 2004