Surface Charge Induced Reentrant Crystalline−Liquid Transition in

8028

Langmuir 2001, 17, 8028-8039

Surface Charge Induced Reentrant Crystalline-Liquid Transition in Colloidal Systems: The Role of the Microion Disposition Kenneth S. Schmitz Department of Chemistry, University of MissourisKansas City, Kansas City, Missouri 64110 Received May 15, 2001. In Final Form: September 21, 2001 The liquid f crystalline f liquid transition as a function of surface charge density (σ) for colloidal dispersions reported by Yamanaka, Yoshida, Koga, Ise, and Hashimoto [Phys. Rev. Lett. 1998, 80, 5806] is examined in detail by Brownian dynamics simulations of the microion distribution. In the absence of added electrolyte, two different distribution functions are determined: Gpc(r) centered at the center of the computation cell and gpc(r) centered at a selected macroion. The “thermal radii” rGtherm and rgtherm are determined for the respective distribution functions at the value exp(1), i.e., when the total electrostatic energy of the system, 〈Esys〉, is equal to the thermal energy kBT. The methods of juxtaposition of potential fields [Langmuir 1997, 13, 5849] and Brownian dynamics simulations are used to describe the disposition of the counterions, or an “orbital” picture of colloidal structures. It is suggested from this orbital model that the crystalline regime obtains when the thermal radii extend far from the parent macroions as this leads to maximum “sharing” of the counterions, and that the reentry transition occurs because of a drastic decrease in rgtherm with high surface charge densities.

1. Introduction The structure of colloidal suspensions and polyelectrolyte solutions is an interesting topic with a polarizing effect on the scientific community in regard to the origin of the forces responsible for observations that seem to defy the conventional wisdom of the mean field theories. It is observed that a “crystalline” structure (highly ordered distribution of macroions) coexists with either a “liquid” structure (disordered or random distribution of particles)1-4 or a “void” structure (absence of colloidal particles).5-8 These “two-state” structures are not an artifact of data analysis as these structures can be directly viewed by digital video microscopic (DVM) methods1-4 and confocal laser scanning microscopy (CLSM) methods.5-8 Several attempts have been made in the literature to explain the two-state structure in terms of a “mean field pair potential” between two colloidal particles. The analytical expressions for the Coulombic contribution to the effective pair potential can be classified as repulsive only (R) or repulsive with a long-range attractive tail (RA). The general form of the R category pair potential is

βURelec )

Zp2F2λB exp(-yr) ap r

(1)

where β ) 1/kBT with kB being the Boltzmann constant and T the absolute temperature, r ) R/ap is a reduced distance of separation normalized by the particle radius (1) Ise, N.; Okubo, T.; Ito, K.; Dosho, S.; Sogami, I. Langmuir 1985, 1, 176. (2) Ito, K.; Nakamura, H.; Yoshida, H.; Ise, N. J. Am. Chem. Soc. 1988, 110, 6955. (3) Larsen, A. E.; Grier, D. G. Nature 1997, 385, 230. (4) Carbajal-Tinoco, M. D.; Castro-Roma´n, F.; Arauz-Lara, J. L. Phys. Rev. E 1996, 53, 3745. (5) Ise, N.; Matsuoka, H.; Ito, K.; Yoshida, H. Faraday Discuss. Chem. Soc. 1990, 90, 13. (6) Ito, K.; Yoshida, H.; Ise, N. Chem. Lett. 1992, 1992, 2081. (7) Ito, K.; Yoshida, H.; Ise, N. Science 1994, 263, 66. (8) Yoshida, H.; Ise, N.; Hashimoto, T. J. Chem. Phys. 1995, 103, 10146.

ap, F is a model-dependent parameter usually associated with the charge Zp to define a “model-dependent charge” ZM ) ZpF, λB ) qe2/4πorkT is the Bjerrum length with qe being the magnitude of the electron charge, o is the permittivity of the vacuum, r is the relative permittivity of the bulk medium, and y ) κap is the dimensionless screening parameter where the screening parameter κ is defined by the equation

κ2 ) 4πλB(

∑j njZj2 + |Zp|npZc2) ) κDH2 + κc2

(2)

where nj is the number concentration of added salt ions with charge (with sign) Zj, np is the number concentration of the macroions of charge (with sign) Zp, and Zc is the charge (with sign) of the counterions which have a number concentration nc ) |Zp|np. We partitioned the screening parameter into two contributions: that from the added salt and denoted as κDH, the Debye-Hu¨ckel screening parameter, and the contribution from the counterions released from the macroions, κc. The DLVO potential9 obtains by setting F ) exp(y)/(1 + y) and the Yukawa potential with F ) 1. Within the set of potentials in the RA classification is the Sogami-Ise (SI) potential,10 and the more recent “volume term” (VT) theories11-15 which employ the DLVO potential but whose formalism leads to the SI form through the calculation of the chemical potentials as recently shown.16 The RA category of colloid-colloid pair potentials has the general form

βURAelec(r) )

Zp2F2λB yr exp(-yr) Bap 2 r

(

)

(3)

(9) Verwey, E. J.; Overbeek, J. Th. G. Theory of the Stability of Lyophobic Colloids; Elsevier Publishing Co., Inc.: New York, 1948. (10) Sogami, I.; Ise, N. J. Chem. Phys. 1984, 81, 6320. (11) van Roij, R.; Hansen, J.-P. Phys. Rev. Lett. 1984, 81, 6320. (12) van Roij, R.; Hansen, J.-P. Prog. Colloid Polym. Sci. 1998, 110, 50.

10.1021/la010725w CCC: $20.00 © 2001 American Chemical Society Published on Web 11/13/2001

Crystalline-Liquid Transition in Colloidal Systems

Langmuir, Vol. 17, No. 26, 2001 8029

For the SI potential10 one sets F ) sinh(y)/y. A feature of the RA category potentials is that a minimum exists in the pair potential profile as a function of the separation distance. The location of the minimum is given, in reduced coordinates rmin ) Rmin/ap, by

rmin )

B + [B(B + 2)]1/2 y

(4)

where

B ) 1 + y coth(y)

(5)

That the form of the SI potential is implicit in the formalism of the VT theories is readily shown by the example of the Yukawa (Y) potential, which obtains by setting F ) 1 ) F. The total mean field reduced macroion pair interaction energy expressed per unit volume, as employed in the VT expressions,11-13 thus may be written as

np2 Zp2λB exp(-yr) np2 βAppelec elec ) φpp ) ) βUYelec V 2 ap r 2 (6) The chemical potential for the jth component, µj, is then calculated from the unit volume Helmholtz free energies by the standard thermodynamic relationship µj ) ∂φ/∂nj, and the chemical potentials for the components in each phase are equated.11-15 Thus, for the macroion-rich phase the total contribution of the macroion and the microion that arises from the pair interaction energy in eq 6 is

np

( ) ( ) ∂φppelec ∂np

+ ns

ns

∂φppelec ∂ns

np

)

[ ]

np2 yr 2βUYelec (7) 2 2

The preexponential multiplicative term in eq 7 exhibits an attraction component exactly as expressed in eq 3. Thus, the free energy contribution in the macroion-rich phase calculated by means of the VT approach has a long-range attractive tail to the free energy pair contribution in a form equivalent to that of the SI theory. While either the R or RA potential adequately describes most of the experimental observations, there are studies reviewed herein to which neither analytical expression applies. Among these are the ultra-small-angle X-ray (USAXS) studies of Matsuoka et al.,17-19 the surface charge density induced reentry transition reported by Yamanaka and co-workers,20 and the light scattering study by Gro¨hn and Antonietti21 indicating a macroscopic phase separation of microgels. The focus of this paper is on the reentry transition as reported by Yamanaka and co-workers.20 The problem is defined by the scattering profiles for the three transition regions liquid (region I) f solid (region II) f reentry liquid (region III) given in Figure 1. (13) van Roij, R.; Dijkstra, M.; Hansen, J.-P. Phys. Rev. E 1999, 59, 2010. (14) Denton, A. R. J. Phys. Condens. Matter 1999, 11, 10061. (15) Warren, P. B. J. Chem. Phys. 2000, 112, 4683. (16) Schmitz, K. S.; Bhyuiyan, L. B. Phys. Rev. E 2001, 63, 011503. (17) Matsuoka, H.; Harada, T.; Yamaoka, H. Langmuir 1994, 10, 4423. (18) Matsuoka, H.; Harada, T.; Kago, K.; Yamaoka, H. Langmuir 1996, 12, 5588. (19) Harada, T.; Matsuoka, H.; Ikeda, T.; Yamaoka, H. Langmuir 1999, 15, 573. (20) Yamanaka, J.; Yoshida, H.; Koga, T.; Ise, N.; Hashimoto, T. Phys. Rev. Lett. 1998, 80, 5806. (21) Gro¨hn, F.; Antonietti, M. Macromolecules 2000, 33, 5938.

Figure 1. USAXS data of Yamanaka et al. on the surface charge induced reentry transition. USAXS intensity (arbitrary units) at different surface charge densities with Cs ) 10 µM: region I, liquid (σe ) 0.07 µC/cm2); region II, solid (σe ) 0.36 µC/cm2); region III, liquid (σe ) 0.72 µC/cm2).

By the very nature of the DLVO theory, increasing the surface charge density of the macroion must result in stronger repulsion interactions between the macroions and therefore an increasingly more rigid crystalline structure. Structural stability in the RA potential is measured by the depth of the potential well. Hence, it is not obvious with the RA potentials as to whether a reentry transition occurs with an increase in the surface charge of the macroions. However, it is shown in the following section that the mean field RA potentials likewise fail to account for a reentry transition. We explore in the present study an “orbital approach” to the description of colloidal suspensions as results from the juxtaposition of potential fields (JPF) approach22 as reviewed in section 3. Brownian dynamics (BD) simulations are given in section 4 for two types of microion distribution functions: Gpc(r) centered at the origin of the spherical computation cell and gpc(r) centered at a “probe” macroion. The JPF and BD results for selected geometries are reported in section 5, and discussed in terms of the orbital description of colloidal clusters. The orbital description of the counterion dispositions is then applied as a novel explanation of the reentry transition in section 6. 2. Failure of the DLVO and SI Potentials To Explain the Surface Charge Induced Reentry Diagram of Yamanaka et al. Yamanaka and co-workers20 reported USAXS studies on the structure of colloidal silica and latex dispersions as a function of surface charge density (σ), molar salt concentration (Cs), and particle volume fraction (φp). In their analysis Yamanaka et al. made a distinction between three definitions of σ: the analytic value σa as determined by titration (bare surface charge), the “effective” value σe determined from conductivity measurements extrapolated to zero colloid concentration, and the value σ* computed from the “charge renormalization” theory of Alexander et al.23 The inequalities σa > σe > σ* as a function of σa were reported, where σ* was constant. Except for the lowest values of φp, it was found that the transitions liquid f crystalline f liquid obtained for increasing values of σe (see Figure 1). They compared the numerical simulations (22) Schmitz, K. S. Langmuir 1997, 13, 5849. (23) Alexander, S.; Chaikin, P. M.; Grant, P.; Morales, G. J.; Pincus, P. J. Chem. Phys. 1984, 80, 5776.

8030

Langmuir, Vol. 17, No. 26, 2001

of Kremer and co-workers24,25 on the phase transition of spheres interacting through the DLVO potential using renormalized charge values with their experimentally determined profile of the phase boundary as a function of Cs and σa. The simulations showed a monotonic increase to a constant plateau, whereas the experimental profile showed a maximum in the vicinity of φp ) 0.005, with a slowly decreasing tail. Yamanaka et al.20 thus concluded that the DLVO theory and charge renormalization methods were inadequate for the interpretation of the reentrant phase diagram. Yamanaka et al. did no further calculations but did suggest the possibility that the reentrant behavior from solid to liquid was governed by an increase in the value of κ such that the magnitude of the pair potential was sufficiently screened as to return to the liquid state. Tata and Ise26,27 reported Monte Carlo simulations on the surface induced reentry phase diagram using the SI pair potential. They focused on surface charge densities of σ ) 2 × 10-7, 4 × 10-7, and 6.8 × 10-7 C/cm2 to represent reentry diagram. [Note that these values correspond to the effective surface charge densities reported by Yamanaka et al.20 for the liquid (σe ) 0.07 µC/cm2), solid (σe ) 0.36 µC/cm2), and reentry liquid (σe ) 0.72 µC/cm2) regimes.] For the σ ) 2 × 10-7 C/cm2 calculation the interparticle spacing obtained from the SI potential was equal to that of the uniform distribution, i.e., 2Dsim/2Do ) 1, where the subscript “sim” refers to “simulation”. For larger values of σ the ratio was found to be 2Ds/2Do < 1, indicating an attraction between the colloidal particles. Their projection diagrams of the colloidal particle positions for surface charge densities showed a highly ordered crystalline structure for σ ) 2 × 10-7 C/cm2, a two-state system with one void region for 4 × 10-7 C/cm2, and a two-state system with two void regions for 6.8 × 10-7 C/cm2. In all cases the pair distribution function exhibited more than one peak. Nonetheless, they26 concluded that, in regard to the 6.8 × 10-7 C/cm2 calculation, “The presence of only a few broad peaks in the case of 6.8 × 10-7 C/cm2 suggests that it is disordered.” This conclusion that a few broad peaks suggests disorder is in contrast with an earlier statement by Ito and Ise28 on their two-dimensional Fourier analysis of colloidal crystals, viz., “The study suggests that the single, broad scattering peak often observed for solutions of macroions and suspensions of latex particles reflects a structural correlation between the solutes.” If one broad peak indicates the coexistence of ordered and disordered regions, then a few broad peaks must also indicate the coexistence of the two structures with a more highly ordered crystalline regime. Furthermore, the USAXS data of Yamanaka and co-workers indicate no difference between, in our notation, regions I and III (cf. Figure 1). Hence, the projection calculation for σ ) 2 × 10-7 C/cm2 (region I) should be indistinguishable from that for σ ) 6.8 × 10-7 C/cm2 (region III). It may thus be concluded that the SI potential does not result in a completely homogeneous structure with a continued increase in the surface charge density. This is in contrast to the Monte Carlo simulations with the SI potential which does give a salt-induced reentry diagram.29-31 Coordinate (24) Kremer, K.; Robbins, M. O.; Grest, G. S. Phys. Rev. Lett. 1986, 57, 2694. (25) Robbins, M. O.; Kremer, K.; Grest, G. S. J. Chem. Phys. 1988, 88, 3286. (26) Tata, B. V. R.; Ise, N. Phys. Rev. B 1996, 54, 6050. (27) Tata, B. V. R.; Ise, N. Phys. Rev. E 1998, 58, 2237. (28) Ito, K.; Ise, N. J. Chem. Phys. 1987, 15, 4176.

Schmitz

Figure 2. Salt dependence of the Sogami-Ise pair potential. The reduced SI potential βUSIelec(r) (cf. eq 4) is plotted as a function of the reduced separation distance r ) R/ap where ap ) 600 Å, Zp ) 1000, and φp ) 0.03. The three curves are for the added electrolyte concentrations: Cs ) 10-6 M (s); Cs ) 10-5 M (---); Cs ) 10-4 M (- - -).

projections of the Monte Carlo (MC) calculations in Figure 7 of Tata, Arora, and Valsakumar30 indicate a completely homogeneous distribution of the colloidal particles for high salt (ns ) 1000npZp) and for low salt (ns ) 2.6npZp), and a well-defined colloidal cluster for intermediate salt (ns ) 78npZp). To explain the differences in the two behaviors of the projection diagrams of the salt and charge variations using the SI potential, we first review the principles set forth by Tata and co-workers32 in regard to the salt-induced “phase separations” and the SI potential. They use the analogy of the vapor-liquid transition in atomic systems where “...as a result of lowering the temperature, when the thermal energy kBT becomes lower than the depth of the potential minimum, this transition takes place”. In other words the two-state structure obtains when |USI(rmin)/kBT| ≡ WSI > 1, where WSI is the magnitude of the depth of the well at the location of the minimum in the pair potential. Using this criterion, it necessarily follows that a signature of the reentry phase transition is that WSI must exhibit a minimum when plotted against the appropriate variable. Shown in Figure 2 are the reduced potentials βUSI(r) for three different salt concentrations (Cs ) 10-4, 10-5, and 10-6 M) with φp ) 0.03, Zp ) 1000, and ap ) 600 Å. Clearly both rmin and WSI are sensitive to changes in Cs. Shown in Figure 3 is the dependence of WSI for Zp ) 700 and ap ) 600 Å on Cs for the volume fractions φp ) 0.01, 0.02, and 0.04. The curves in Figure 3 thus indicate a salt-induced reentry transition can be explained in terms of the SI potential in accordance with the criterion of Tata, Rajalakshmi, and Arora.32 At low salt concentrations the minimum WSI is too shallow (or nonexistent) to hold the particles in a lattice position. An initial increase in salt results first in a greater depth of the potential well. Since rmin is less than the average spacing in a uniform distribution, the colloidal particles “phase separate” into colloidal clusters and void regions as previously noted in Monte Carlo simulations.26,27,29-32 A further increase in the added salt again results in a shallow potential well and thus reentry into the liquid configuration. From the preceding paragraph we take as the signature of a reentry “phase” diagram the presence of a minimum (29) Arora, A. K.; Tata, B. V. R.; Sood, A. K.; Kesavamoorthy, R. Phys. Rev. Lett. 1988, 60, 2438. (30) Tata, B. V. R.; Arora, A. K.; Valsakumar, M. C. Phys. Rev. E 1993, 47, 3404. (31) Tata, B. V. R.; Yamahara, E.; Rajamani, P. V.; Ise, N. Phys. Rev. Lett. 1997, 78, 2660. (32) Tata, A. K.; Rajalakshimi, M.; Arora, A. K. Phys. Rev. Lett. 1992, 69, 3778.

Crystalline-Liquid Transition in Colloidal Systems

Figure 3. Salt dependence of the depth of the well of the SI pair potential. The reduced depth of the SI potential WSI ) βUSIelec(rmin) is plotted as a function of the added electrolyte concentration Cs where ap ) 600 Å and Zp ) 700. The three curves are for the volume fractions: φp ) 0.01 (s); φp ) 0.02 (---); φp ) 0.04 (- - -). The presence of a minimum is indicative of a salt induced re-entrant transition.

Figure 4. Surface charge dependence of the depth of the well of the SI pair potential. The reduced depth of the SI potential WSI ) βUelec SI (rmin) is plotted as a function of the surface charge Zp where ap ) 600 Å and Cs ) 2 × 10-6 M. The three curves are for the volume fractions: φp ) 0.001 (s); φp ) 0.005 (---); φp ) 0.01 (- - -). The absence of a minimum indicates that there is no surface charge induced reentrant transition.

in the depth of the potential well, WSI, where the homogeneous structure exists on either side of the minimum with WSI < 1. We now show that these characteristics are not present with the SI potential when the surface charge is varied over the experimental range of interest. To mimic the experimental data in Figure 3 of Yamanaka et al.,20 we set ap ) 600 Å and Cs ) 2 µM and vary the effective surface charge over the range 0 < Zp < 2000. Shown in Figure 4 is WSI for the three volume fractions φp ) 0.001, 0.005, and 0.01. The lack of a minimum in the profile of WSI vs |Zp| indicates that the SI potential does not describe fully the experimental observations over the range of experimental parameters. 3. Juxtaposition of Potential Fields: “Orbital” Model of Colloidal Clusters It is an experimental observation of the two-state structure in colloidal suspensions that any particular crystalline configuration persists for several minutes, as monitored by DVM1-4 and CLSM5-8 methods. One may therefore focus on the distribution of the microions in the presence of a fixed configuration of macroions over this experimental time frame. We may thus view these colloidal clusters in a manner similar to that for molecules and crystals. One such method for the determination of atomic participation in chemical bonding is the gradient vector field approach of Bader and co-workers.33,34 The more massive components, i.e., the atoms, are held in fixed positions, and the electrons are then allowed to distribute (33) Bader, R. F.; Slee, T. S.; Cremer, D.; Krata, E. J. Am. Chem. Soc. 1983, 105, 5061. (34) Bader, R. F. Acc. Chem. Res. 1985, 18, 9.

Langmuir, Vol. 17, No. 26, 2001 8031

themselves in accordance to the electric field set up by the atoms. The characteristics of the gradient of the vector field are thus used to determine which atoms interact with other atoms and which atoms participate in the formation of a chemical bond. Identification of interacting atoms is manifested as a “corridor of vectors” connecting the atoms, whereas chemical bonding is inferred by the presence of closed surfaces defined by tangent vectors about the participating atoms. The gradient vector field method was adapted to the description of colloidal clusters as the JPF22 method, where the macroions function in the same manner as the atoms and the microions as the electrons. Thus, the connotation of an “orbital” model of colloidal clusters. The JPF approach views the system as being composed of a “solution phase” containing the solvent and solute particles, i.e., the microions, but not the macroions. The solvent functions solely as a continuum dielectric medium. The solution phase is further subdivided into subregions that are sufficiently large to contain several microions but viewed as having a constant chemical potential throughout the subregion. For the purpose of illustration of the principles, we cite the standard textbook development related to sedimentation equilibrium, where the local gravitational field establishes the local concentration of macromolecules. In this approximation the cumulative electric field due to all of the macroions as evaluated in the subregion plays the same role as the local gravitational field in sedimentation equilibrium studies. Similarly the electrical charge on the microions plays the same role as the mass in sedimentation equilibrium, but with one important difference. The microions can be positively or negatively charged. Hence, the increase or decrease in the microion concentration with a change in ΨM(r) is dependent upon whether the charge on the microions has the same sign as the macroion charge. With this understanding of the parallels in the two systems, the reduced form of the chemical potential of the subregion located at r takes on the form

µions(r) µ ) ΨM(r) + kBT kBT

(8)

where the reduced potential due to all the colloidal particles is given by the sum M

ΨM(r) ) λB

∑

Zm

m)1am|rm|

(9)

|rm| ) |Rm|/am is the reduced distance from the center of the mth macroion of diameter am and charge Zm, and the chemical potential due to the microions in the region is Js

µions(r) )

∑ j)1

Js

µj(r) )

[µjo + kBT ln(γj(r) Cj(r))] ∑ j)1

(10)

where Js is the number of ion types, Cj(r) is the molar concentration of the jth ion type, and γj(r) is the activity coefficient that reflects the electrical interactions of the microion. At this level of illustration we assume that the interactions between the microions can be viewed as a perturbation on the local chemical potential set up by eq 8. These microion-microion interactions, as well as hard sphere interactions, are of course included in the simulations as described in the following section. According to eq 9 the configuration of the macroions determines the value of ΨM(r) in each subregion. The

8032

Langmuir, Vol. 17, No. 26, 2001

Schmitz

lations were done with a spherical cell model whose volume reflects the volume fraction of the macroions. Hence, the reduced radius of the computational cell for Np colloidal particles, denoted by rcell, is given by

rcell ) (Np/φp)1/3

(11)

There are two types of microion distribution functions to be calculated: Gpc(r), which is centered at the origin of the coordinate system of the computation cell, and gpc(r), in which the coordinate system is centered on one of the macroions in the cluster not located at the cell center. Both Gpc(r) and gpc(r) consisted of 500 equally spaced bins over the range 0 < r < rb, where rb ) rcell - b is the reduced radius of the appropriate spherical cell region, and b is the radial distance of the center location of the distribution function from the cell origin. If there is only one particle in the calculation, then b ) 0, r0 ) rcell, and gpc(r) ) Gpc(r). In the multi-macroion calculations rb ) r0 ) rcell for Gpc(r) and rb ) rcell - b for gpc(r). The initial positions of the Nc counterions were randomly distributed within the spherical computational cell. The position of the generalized coordinate of the jth microion for the k + 1 iteration is given formally by the evolution expression Figure 5. Macroion potential profiles in the x-y plane. The potential ΨM along the x axis in the x-y plane is shown above for the macroion configuration shown in the inset. The charge is Zp ) 10, and all distances are normalized to the macroion radius, viz., ap ) 1. Top: isolated macroion. Bottom: cluster of seven macroions with a separation distance d ) 3. Notice that the potential near the surface of the macroions is greater than the surface potential of each of the macroions in the cluster due to the superposition of the fields of the other macroions.

microion concentration in eq 10 thus adjusts itself in such a way as to maintain the constant chemical potential in the subregion as given by eq 8. To illustrate, we compare the potential field profile for an isolated macroion with a cluster of seven macroions. In this comparison we employ the macroion charge Zp ) 10 and a reduced macroion radius of ap ) 1. For the cluster we place the probe macroion at the origin and the remaining six macroions along each axis in a symmetric way at reduced distances d ) (3. Shown in Figure 5 are the potential profiles for these two systems taken along the x axis in the x-y plane at z ) 0. The horizontal lines interior to the macroion locations are the reduced surface potential for that particular macroion. In the case of the isolated macroion the potential shows a monotonic decrease from the macroion surface. However, for the cluster of seven macroions the potential immediately from the surface of the centrally located macroion exceeds that of its interior due to the cumulative effect of the surrounding macroions. Likewise a similar relationship obtains for the surfaces facing the interior of the cluster for the other macroions in this cluster, but not for the sides facing away from the cluster origin. The implication is that the “effective charge” of the macroions interior to the cluster should be less than that anticipated from extrapolation of the data to infinite dilute conditions. This deduction from the JPF profiles is verified by BD simulations described in the following section. 4. Brownian Dynamics Simulations The BD simulations follow standard methods.35 For simplification and computational purposes the BD simu(35) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon Press, Oxford Science Publications: Oxford, England, 1993.

qjk+1 ) qjk + βDj∆tFj + ω(Dj∆t)1/2

(12)

where Dj is the diffusion coefficient, ∆t is the time interval for the move, -1 e ω e + 1 is a random number, and Fj is the electrical force exerted on the jth microion by all of the other particles present, including the macroion. In practice the following substitution is made for the simulations: 2D∆t ) 〈S2〉, where S is a step size whose value is proportional to the radius ac of the counterion, viz., S ) Mac, where M is a multiplicative factor. All of the distances employed in these calculations were scaled to the radius of the macroion. The step size in the simulations was thus sr ) S/ap. We performed calculations in which all of the parameters are scaled to the radius of the macroions, leading to the reduced parameters used in these simulations: apr ) 1 and acr ) ac/ap. The major obstacle with simulations involving macroions is the number of particles. Clearly the calculations are prohibitive for cluster simulations for large values of Zp and Np. We therefore take advantage of the tabulated data of Roberts et al36 in which a universal curve was reported of the ratio of the effective charge to bare charge, Zeff/Zp, as a function of the surface potential, |Zp|/ap. Hence, simulations of a macroion of large charge and size may be mimicked by using smaller values of Zp and ap as long as the ratio |Zp|/ap is the same. The BD simulations of the microions are performed in two stages: a preliminary calculation to “equilibrate” the system, usually consisting of 106 iterations for each microion in the system, and a “final” calculation generally in the range of 106 to 108 iterations, in which case the positions of each microion are tabulated for the calculation of the microion distribution function. Hard sphere overlap with other particles, i.e., interpenetrating, was checked after each move, and the particle was returned to its initial position if the test was positive. The partition function for the system, Q, was calculated from the accepted moves for all of the microions in the system for a particular configuration J, and summed over all configurations J

Q)

∑J exp(-βEJ)

(13)

Crystalline-Liquid Transition in Colloidal Systems

Langmuir, Vol. 17, No. 26, 2001 8033

and the reduced average internal energy β〈Esys〉 was calculated from the standard expression in statistical mechanics

β〈Esys〉 ) β

∑J EJ exp(-βEJ) Q

(14)

According to Verwey and Overbeek9 the interaction between double layers involves only the counterions if the solution is salt-free, viz. (p 58 of ref 9), “The very simple result is that we find the total free energy of the double layer if only we calculate the electric work necessary to discharge stepwise all ions of the solution.” They also claim this extends to added salt solutions,9 “A further simplification is that we need only to consider those ions which are present in excess in the solution and carry the liquid charge of the double layer, since for the mutual neutralization of all other ions, present in equal amounts in each point in space, the net work is zero.” In other words, the Verwey-Overbeek formulation of the problem focuses only on one aspect of the electrical system, namely, the interaction of counterions centered around one macroion, forming a double layer, with a corresponding double layer of the second macroion. To describe the complete system, one must included electrical interactions of all origins. Limiting these contributions to pairwise interactions, the component parts of the electrical systems are the macroion-macroion interaction, β〈Epp〉, the macroion-counterion interaction, β〈Epc〉, the counterion-counterion interaction, β〈Ecc〉, the macroion-added electrolyte interaction, β〈Eps〉, the counterion-added electrolyte interaction, β〈Ecs〉, and the interaction between the added electrolyte particles, β〈Ess〉. In this notation the subscript “s” refers to both the cations and the anions of the added electrolyte. While it is our belief that the added electrolyte is a significant contributor to the phase stability of the system, our present focus is on “salt-free” systems to emphasize the primary features of “sharing” counterions between macroions comprising a cluster. Omission of the added salt contributions is partially justified by the fact that the added electrolyte in the experimental systems of interest is in the micromolar range and that under these conditions the counterions released from the macroions are more than an order of magnitude larger in concentration than the added electrolyte. The total energy of the salt-free system in the present study is therefore given by its component pairwise interactions

β〈Esys〉 ) β〈Epp〉 + β〈Epc〉 + β〈Ecc〉

(15)

where the average counterion-counterion interaction was calculated using the expression Nc-1 Nc

β〈Ecc〉 ) β

∑J ( ∑ ∑Eij)J exp(-βEJ) i)1 j> i Q

(16)

where the double sum is over the counterion pairs for the system state J. Since the macroions are in fixed positions, the average pairwise macroion-macroion interaction was simply calculated from the unscreened Coulomb expres(36) Roberts, J. M.; O’Dea, J. J.; Osteryoung, J. G. Anal. Chem. 1998, 70, 3667.

sion

β〈Epp〉 )

Zp2λBNp-1 Np 1

∑ ∑ i)1 j>i r

ap

(17)

ij

where the Bjerrum length was chosen to be λB ) 7.18 Å in these calculations. The macroion-counterion interaction was indirectly calculated from eqs 14-17, viz.

β〈Epc〉 ) β〈Esys〉 - β〈Epp〉 - β〈Ecc〉

(18)

We employ two measures of the degree of distortion of the counterion cloud for macroions within a cluster: the “thermal radius” and an “effective charge”. Both of these concepts are based on the pair distribution function for the counterions in the presence of the macroions, Gpc(r) and gpc(r). The thermal radius is defined as that distance from the appropriate origin for which the total interaction energy is equal to the thermal energy. Hence, rGtherm and rgtherm are defined, respectively, at the distance for which Gpc(r)rGtherm) ) exp(1) and rgpc(r)rgtherm) ) exp(1). The thermal radius is therefore a measure of the range of the electrostatic influence of a macroion on the counterions in the system in the presence of all pairwise interactions. There are two types of effective charges calculated in the present study: the “thermal charge” Ztherm and the “surface charge” Zsurf. The thermal charge for a single macroion is composed of the bare macroion charge plus the charges of all the counterions within the thermal radius (i.e., charge neutralization counterions), as previously defined22,37

∫1rg

Ztherm ) Zp + ZcNc

Fpc(r) dr

therm

(19)

where Fpc(r) is a number density as measured from the center of the chosen colloidal particle. This definition of Ztherm was shown successful37 in reproducing the trend of Zeff/Zp versus |Zp|/ap in the Roberts et al.36 summary, and also the empirical equation of Yamanaka and co-workers.20 In the present study of a cluster of macroions the definition of the effective charge is slightly modified by the normalization of the number density to that fraction of counterions found within the cell radius rb ) rcell - b, viz.

∫1r Fpc(r) dr ) fb b

(20)

where fb is the fraction of the counterions within the radius rb as monitored during the BD simulations. To distinguish the effective charges obtained from the two different types of distribution functions, we employ the notations Zgtherm and ZGtherm. We also define a surface charge Zsurf as being the bare charge plus the neutralizing counterions within two counterion diameters from the surface of the colloid particle, viz.

∫11+4a Fpc(r) dr

Zsurf ) Zp + ZcNc

r

c

(21)

All calculations were performed on the DEC Alpha AXP2100/M500 at the University of Missouri Computing Facilities. The results were then expressed in graphics form using Mathematica. 5. JPF and BD Results: The Orbital Description of Colloidal Clusters Present day computers are incapable of handling multibodied microion effects on multibodied macroion (37) Sanghiran, V.; Schmitz, K. S. Langmuir 2000, 16, 7566.

8034

Langmuir, Vol. 17, No. 26, 2001

Schmitz

Table 1. BD Results for Clusters of One, Two, and Three Charged Spheresa run Np 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 c

1 1 1 1b 1c 2d 2e 2e 2f 2f 3d 3d 3e 3e,g 3e 3f

〈Nav〉 1× 5 × 106 1 × 106 1 × 106 1 × 106 1 × 106 5 × 107 5 × 106 1 × 107 1 × 106 1 × 106 1 × 107 1 × 107 1 × 107 1 × 106 1 × 106 107

rgtherm Zgtherm/Zp 1.482 1.509 1.482 1.643 1.429 1.317 1.363 1.363 1.597 1.597 1.333 1.323 1.366 1.366 1.366 1.624

0.2413 0.2186 0.2440 0.1868 0.1292 0.3033 0.2273 0.2094 0.1168 0.1209 0.2236 0.2595 0.2372 0.2358 0.2408 0.0836

rsurf

Zsurf/Zp

β〈Esys〉

β〈Ecc〉

1.04 1.04 1.04 1.80 1.004 1.04 1.04 1.04 1.04 1.04 1.04 1.04 1.04 1.04 1.04 1.04

0.7405 0.7364 0.7364 0.1619 0.8617 0.7442 0.7601 0.7556 0.7345 0.7329 0.7438 0.7553 0.7241 0.7264 0.7263 0.7272

-41.46 -39.27 -38.86 -35.47 -42.68 -82.37 -83.18 -78.26 -83.07 -81.13 -122.4 -125.4 -125.10 -125.20 -122.30 -123.90

19.98 18.21 17.55 17.69 20.35 58.97 63.82 58.03 66.53 63.59 118.60 120.5 127.60 129.20 123.90 134.90

a Z ) 10, a ) 1 nm, a r ) 0.01, sr ) 0.01, φ ) 0.01. b a r ) 0.2. p p c p c acr ) 0.001. d d ) 3.5. e d ) 3.0. f d ) 2.5. g sr ) 0.2.

systems if one is required to use actual charge and radius parameters. As previously noted, the data of Roberts et al.36 give rise to the suggestion that a suitable parameter is the ratio |Zp|/ap rather than their actual values. Recent BD simulations37 indicate that this hypothesis is reasonable, where the charge fraction Ztherm/Zp as calculated from eq 19 was found to be a function of |Zp|/ap over a wide range of individual values of Zp and ap. We therefore use as our standard system the ratio |Zp|/ap ) 10 nm-1 and a volume fraction of φp ) 0.01, where the operational variables are Zp ) 10, ap ) 1 nm, and acr ) 0.01. In the discussion of Figure 5 it was pointed out that the total potential due only to the macroions was enhanced near the surface of a central macroion by the presence of all the other macroions in the system. This behavior is a primary source of the nonlinearity in these complex systems, and thus invalidates the superposition approximation as contained in linearized theories such as the DLVO potential. Failure of the DLVO potential to describe the distribution of counterions about two isolated macroions was previously reported by Sa´nchez-Sa´nchez and Lozada-Cassou38 in their exact numerical solution to the integral Poisson-Boltzmann equation, and also by Lo¨wen, Hansen, and Madden39 in their Monte Carlo simulations. We therefore examined dimer and trimer clusters using BD simulations for the purpose of comparison with these studies, and also as reference points for larger clusters. Summarized in Table 1 are the results of BD simulations for the standard system with a relative step size of sr ) 0.01. The average number of moves, 〈Nav〉, for each counterion particle is in the range 106 to 107. We also include in Table 1 for comparison purposes simulations using different values of acr and sr. We first focus on run 3 as the reference simulation, and compare with runs 4 and 5. In all three of these simulations the relative step size is sr ) 0.01. In run 4 the relative size of the counterion is increased by a factor of 20 to the value acr ) 0.2, whereas in run 5 it is decreased by a factor of 10 to the value acr ) 0.001. In run 4 the values of rgtherm and rsurf are increased with a concomitant decrease in the values of Zgtherm/Zp and Zsurf/Zp. These variations are interpreted as an excluded volume effect since the radius of the counterion is comparable to that of the macroion. Quite clearly the increase in rsurf is a direct result of the definition of rsurf. Less obvious is the increase in rtherm. As the counterions (38) Sa´nchez-Sa´nchez, J. E.; Lozada-Cassou, M. Chem. Phys. Lett. 1992, 190, 202. (39) Lo¨wen, H.; Hansen, J.-P.; Madden, P. A. J. Chem. Phys. 1993, 98, 3275.

are attracted to the macroion, their larger volumes begin to play an increasing role in excluding other counterions, and the process now becomes a problem in packing of the counterions. Since fewer counterions can approach the surface, the total energy of the region is larger due to the presence of the macroion charge. Therefore, to achieve a total interaction energy equal to the thermal energy, one must go further from the macroion surface than for the smaller counterion size. As a result of the larger values of rtherm and rsurf the volume of the surrounding “shell” associated with charge neutralizing counterions is larger, and therefore, the ratios Zgtherm/Zp and Zsurf/Zp are smaller. The situation is different in run 5, where the value of rtherm is comparable to that in run 3 but with a decrease in the value of Zgtherm/Zp. In this case the electrostatic effects dictate the value of rtherm, and the relative size of the counterions allows more to be within the neutralizing charge volume. This conclusion is partly based on the observation that β〈Ecc〉 is larger for run 5 than for either run 3 or 4. Hence, the ratio Zgtherm/Zp decreases in value. However, rsurf is again subject to its definition in regard to counterion size. Therefore, the ratio Zsurf/Zp increases relative to the value in run 3. Simulation run 14 was performed with the step size sr ) 0.2 and is compared with run 13 with step size sr ) 0.01. Even though there is a difference by a factor of 20 in step size, there is no significant difference in the calculated parameters. The above comparisons emphasize both the importance and the limitations of using ratios in the simulations. If one is to examine the effects of surface charge density on the effective charge of the macroion, the simulations must be carried out with variable Zp/ap ratios and a fixed as/ap ratio to minimize these “excluded volume” effects of the counterion size. For example, if one chooses an ion radius of ac ) 0.2 nm, as usually employed in simulations, with the ratio ac/ap ) 0.01, then the corresponding macroion parameters are ap ) 20 nm and Zp ) 200 for Zp/ap ) 10. Likewise one may apply these results to a system of a more complex and larger counterion. If we choose ac ) 0.4 nm, then the calculations at the ratios Zp/ap ) 10 and as/ap ) 0.01 are for macroions of charge Zp ) 400 and radius ap ) 40 nm. It is for this reason that the simulations are restricted to fixed values of Zp/ap ) 10 and as/ap ) 0.01. It is also noted that the choices of step length sr and number of cycles 〈Nav〉 have small effects in the reported simulations. We next focus on the dimer simulations. The thermal radius rgtherm tends to increase in magnitude as the separation distance d is decreased. Likewise the effective charge ratio Ztherm/Zp decreases from 0.30 to 0.12, which reflects this distortion of the counterion distribution as implied in the rgtherm values. The distribution function Gpc(r) for two macroions separated by a distance d ) 2.5 is compared with the 3-dimensional representation of the results of Sa´nchez- Sa´nchez and Lozada-Cassou38 in Figure 6. The parameters for the Sa´nchez-Sa´nchez and LozadaCassou simulations are for macroions of radius ap ) 7 Å and surface charge density σa ) 0.018 C/m2 separated by a distance d ) 3.214, with a counterion radius ac ) 2.125 Å and a salt concentration of 0.01 M. Notice that in the 3-dimensional representation the two macroions cannot be discerned by the envelope of the counterions. In contrast, separating these two macroions to a distance d ) 5.03 results in distinguishable macroion locations in the Sa´nchez- Sa´nchez and Lozada-Cassou simulations.38,40 Hence, both Gpc(r) and the Sa´nchez-Sa´nchez and LozadaCassou results clearly show an extensive accumulation of (40) Stigter, D. J. Phys. Chem. 1960, 64, 838.

Crystalline-Liquid Transition in Colloidal Systems

Langmuir, Vol. 17, No. 26, 2001 8035

Figure 6. Comparison of counterion distributions about two macroions. The distribution function Gpc(r) is shown for Zp ) 10 and ap ) 1 for two spheres separated by a distance d ) 2.5 with 〈Nav〉 ) 1 × 107. The location of rGtherm ) 2.371 indicated by the “+” extends into the region occupied by both charged spheres. Also shown is the 3-dimensional plot of the results of Sa´nchez-Sa´nchez and Lozado-Cassou38 of the counterion distribution about two macroions with a surface charge density σo ) 0.018 C/m2 and a salt concentration of 0.01 M. The radius of each sphere is 7 Å, and the counterion radius is 2.125 Å. The reduced separation distance is d ) 3.214. Notice that the two macroions are enclosed by the counterion cloud, as also obtains for the Gpc(r) shown above. (The photograph was taken by Jorge Lodogiani and kindly provided by Marcelo Lozado-Cassou.) Table 2. BD Results for Clusters of Seven Charged Spheresa run 17f 18b,f 19f 20g 21g 22c,g 23d,g 24d,g

〈Nav〉

rGtherm ZGtherm/Zp rgtherm Zgtherm/Zp

1× 1 × 107 1 × 106 1 × 107 1 × 106 1 × 106 1 × 106 3 × 106 107

1.632 1.614 1.668 4.080 4.080 4.116 1.331 1.348

0.1163 0.1263 0.0836 0.1437e 0.1413e 0.2655e 0.0619 0.0486

1.267 1.278 1.267 1.315 1.327 1.374 1.209 1.209

0.292 0.298 0.300 0.267 0.267 0.531 0.120 0.120

β〈Esys〉

β〈Ecc〉

-290.5 454.8 -288.7 451.8 -283.7 439.3 -289.8 500.4 -284.7 485.5 -67.6 108.0 -1128.0 2138.0 -1146.0 2190.0

a Z ) 10, a ) 1 nm, a r ) 0.01, sr ) 0.01, φ ) 0.01. rG p p c p therm and ZGtherm relate to the sphere located at the coordinate origin (0,0,0), and rgtherm and Zgtherm to a sphere at (d, 0, 0). b acr ) 0.2, sr) 0.2. c Z ) 5. d Z ) 20. e The bare charge being neutralized by the p p counterions is 70 since the value of rGtherm < d includes all of the macroions in the cluster. f d ) 3.5. g d ) 3.0.

Table 3. BD Results for Clusters of 27 Charged Spheresa macroion designation

coordinates

environment degeneracy

b

central face-centered edge-centered corner

(0, 0, 0) (3.5, 0, 0) (3.5, 3.5, 0) (3.5, 3.5, 3.5)

1 8 12 8

0 3.5 4.95 6.06

rgtherm Zgtherm/Zp 1.724 1.414 1.307 1.27

0.0858 0.2062 0.2593 0.2631

a Run 25: Z ) 10, a ) 1 nm, a r ) 0.01, sr ) 0.01, φ ) 0.01, p p c p d ) 3.5, 〈Nav〉 ) 1 × 106, β〈Esys〉 ) -1053, β〈Ess〉 ) 4482, β〈Epp〉 ) 4189.03, β〈Epc〉 ) -9724.03.

counterions between the two macroions at sufficiently close distances of separation. The trimer simulations in Table 1 further indicate an excess accumulation of counterions within the boundaries defined by their locations. Note, however, that compared to the dimer results the ratio Ztherm/Zp for the trimer is considerably smaller for the same macroion separation distance of d ) 2.5. To explore the “orbital model” in more detail, we turn from the two-dimensional structures of equivalent macroions in Table 1 to the asymmetric three-dimensional structures in Table 2 for Np ) 7 and Table 3 for Np ) 27. The structures of 7 and 27 macroions were chosen so that Gpc(r) now describes the counterion distribution about the centrally located macroion in a symmetric environment in regard to the distribution of macroions whereas gpc(r) represents the counterion distribution about the macroions in an asymmetric macroion environment. It is clear that

Figure 7. Macroion potential field for a cubic cluster of 27 spheres. The above potential profile is in the x-y plane at z ) 0 for a cubic lattice of three spheres on an edge with the lattice constant d ) 3.5. The charge is Zp ) 10 with λB/ap ) 0.718.

the centrally located macroion is more “highly shielded” than the externally situated macroions as inferred from the ratio Zgtherm/Zp. It is further noted that for Np ) 7, d ) 3, and Zp ) 10 (runs 20 and 21) and Zp ) 5 (run 22) the value rGtherm includes all seven macroions. Hence, the “bare charge” to be neutralized in these examples is 70, i.e., that of all seven macroions in the clusters. However, for run 23 for which Zp ) 20 the value of rGtherm is less than d, and therefore the bare charge that is neutralized in this system is that of the single macroion, viz., Zp ) 10. These observations are significant in the interpretation of the reentry transition. In the case of the cubic array with Np ) 27 there are four different macroion environments tabulated in Table 3: the “central” macroion having the coordinates (0, 0, 0) for which b ) 0, the “face-centered” macroion with coordinates (3.5, 0, 0) for which b ) 3.5, the “edge-centered” macroion with coordinates (3.5, 3.5, 0) with b ) 4.95, and the “corner” macroion at (3.5, 3.5, 3.5) and b ) 6.06. The parameter pairs (Zgtherm/Zp, degeneracy of the environment) for these four macroions are central (0.086, 1), facecentered (0.2062, 8), edge-centered (0.2593, 12), and corner (0.2631, 8). In the context of the JPF approach, the potential field of all the macroions is maximum at the center of the array, and greater for the face center location than at the edge center or corner locations. The values of the effective charges correlate well with the local potential fields set up by the macroion distribution in accordance with the JPF view. Shown in Figure 7 is the threedimensional plot of ΨM in the x-y plane at z ) 0 for the 27 macroion cubic array used in the BD simulations. This potential plot supports the effective charge simulations: the largest peak is at the central location, where it is anticipated that the counterion concentration is the greatest and should therefore provide the largest screening of the macroion charge, the next largest peak is at the face-centered location, and the least screened charge is at the edge-centered location of the cluster. Not shown in Figure 7 is the potential about the corner macroions as they do not lie in the x-y plane. Shown in Figure 8 is the distribution function Gpc(r) about the central macroion, where the four environments are readily discerned. The relative magnitude of the macroion peaks is deceptive since one must take into consideration the degeneracy of the environment (cf. Table 3) and the relative amount of

8036

Langmuir, Vol. 17, No. 26, 2001

Schmitz

Figure 8. Distribution function Gpc(r) for the central macroion in a 27 macroion simple cubic cluster. The distribution function Gpc(r) was computed for the cubic lattice of 27 spheres with d ) 3.5 as the lattice constant. The peaks in the distribution function correspond to the central macroion at (0, 0, 0) for which b ) 0, the face-centered macroion with coordinates (3.5, 0, 0) for which b ) 3.5, the edge-centered macroion with coordinates (3.5, 3.5, 0) with b ) 4.95, and the corner macroion at (3.5, 3.5, 3.5) and b ) 6.06. The parameter pairs (Zgtherm/Zp, degeneracy of the environment) for these four macroions are central (0.086, 1), face-centered (0.2062, 8), edge-centered (0.2593, 12), and corner (0.2631, 8).

Figure 10. Vector gradient fields of the cubic macroion cluster. Shown above are the gradient vector fields -Zc∇ΨM ) ∇ΨM in the x-y plane at z ) 0 for the cubic cluster described in Figures 7-9. Top: immediate vicinity of the cluster. Bottom: zoomed away by a factor of 2. See the text for details.

Figure 9. Constant potential contours of the 27 macroion simple cubic cluster. Shown above are the constant potential contours for the 27 macroions in the simple cubic array with d ) 3.5. There are three distinct characteristics in this plot: (1) constant contour lines that encircle only one macroion (inner orbitals), (2) constant contour lines that encompass two or more macroions (shared orbitals), and (3) constant contour lines that enclose the entire cluster of macroions (outer orbitals).

“empty space” in the respective radial shells wherein the macroions reside. We now draw analogies of the counterion distributions in the colloidal cluster with electrons in molecules. Shown in Figure 9 are the constant potential contours in the x-y plane at z ) 0, where the lighter the shading the larger the value of ΨM. The fact that the locations of the macroions are dark emphasizes the fact that the sum of the potential fields in the vicinity of any one macroion in the cluster is greater than at the surface potential of the isolated macroion. From the constant potential contours we may identify three distinct categories: (1) constant contour lines that encircle only one macroion, referred to as “inner orbitals”, (2) constant contour lines that encompass two or more macroions, referred to as “shared orbitals”, and (3) constant contour lines that enclose the entire cluster of macroions, referred to as “outer orbitals”. Notice that

the inner orbitals about the central macroion are represented by concentric circles whereas the inner orbitals about the other macroions are distorted in the direction of the coordinate origin. This exaggeration is manifested in the spherical averaged values of rgtherm. Shown in Figure 10 are the gradient vector fields -Zc∇ΨM ) ∇ΨM, which are also identified with the force vectors acting on the counterions by the positively charged macroions, where the top panel is in the immediate vicinity of the cluster and the bottom panel is zoomed away from the cluster by a factor of 2. In comparing Figure 9 with the top of Figure 10, we note that the force within the boundaries of the inner orbitals is directed toward the encompassed macroion. The counterions within these regions play a parallel role as inner electrons in atoms comprising the molecule. That is, the counterions within this domain shield the charge of the macroion in a manner similar to that of inner electrons which shield the nuclear charge of the atoms. We next focus on the “corridors” connecting the macroions in the top of Figure 10 as defined by the parallel vectors joining these macroions. These corridors are indicative of interactions between these macroions just as similar corridors observed in the vector gradient field of molecules indicate direct interactions between the atoms.33,34 However, one cannot state that the counterions within this potential region are associated, in our analogy, with the electrons forming a chemical bond. In the Bader model the chemical bond between two atoms is characterized by a “closed surface” defined by vectors tangent to the surface.33,34 The potential field in the shared orbital region in Figure 10 is relatively flat, or shallow. Hence, the counterions within this region are “shared” by the

Crystalline-Liquid Transition in Colloidal Systems Table 4. Partitioning of the Average Pairwise Energiesa run

Np

β〈Esys〉

β〈Upp〉

β〈Ecc〉

β〈Epc〉

3 6 11 19 25

1 2 3 7 27

-38.86 -82.37 -122.4 -283.7 -1053

0 20.51 61.54 327.93 4189.03

17.55 58.97 118.60 439.3 4482.0

-56.41 -161.85 -302.54 -1050.93 -9724.03

Zp ) 10, ap ) 1 nm, acr ) 0.01, sr ) 0.01, φp ) 0.01, d ) 3.5, 〈Nav〉 ) 1 × 106. a

Table 5. Average Unit Energiesa run

Np

β〈Esys〉/Np

β〈Upp〉/Np

β〈Ecc〉/Np

β〈Epc〉/Np

3 6 11 19 25

1 2 3 7 27

-38.86 -41.85 -40.80 -40.52 -39.0

0 10.25 20.51 46.85 155.15

17.55 29.49 39.53 62.76 166.0

-56.41 -80.93 -100.85 -150.13 -360.15

a Z ) 10, a ) 1 nm, a r ) 0.01, sr ) 0.01, φ ) 0.01, d ) 3.5, p p c p 〈Nav〉 ) 1 × 106.

Table 6. Pairwise Average Interaction Energiesa run

Np

β〈Esys〉

β〈Upp〉pp

β〈Ucc〉cc

β〈Ucc〉pc

3 6 11 19 25

1 2 3 7 27

-38.86 -82.37 -122.4 -283.7 -1053

0 20.51 20.51 15.62 11.93

0.39 0.31 0.27 0.18 0.12

-5.64 -4.05 -3.36 -2.14 -1.33

Zp ) 10, ap ) 1 nm, acr ) 0.01, sr ) 0.01, φp ) 0.01, d ) 3.5, 〈Nav〉 ) 1 × 106. 〈Upp〉pp ) 2〈Upp〉/Np(Np - 1), 〈Ucc〉cc ) 2〈Ucc〉/Nc(Nc - 1), 〈Upc〉pc ) 〈Upc〉/NpNc. a

enclosed macroions since there is no strong force acting on these counterions to direct them to specific macroions. That is, this region is topologically similar to a valley or meadow between neighboring mountains. Turning now to the bottom of Figure 10, we note that there is a region about the entire cluster for which the directed field gradients are along the tangent of a larger circle encompassing the entire cluster. In the present context this boundary means that the counterions are constrained to remain within the domain defined by the cluster. The average energies for systems of 2, 3, 7, and 27 macroions with a separation distance d ) 3.5 are summarized in Table 4 along with the isolated macroion. To test the “additivity” principle, the average energies of these systems are also reported in terms of their “unit cell values” by division by Np in Table 5. It is noted that the unit cell system energy remains relatively constant, being in the range -41 < β〈Esys〉/Np < -39. It may be concluded that β〈Esys〉 does indeed follow the additivity principle. Paradoxically the unit cell values of the pairwise interactions comprising β〈Esys〉 do not obey the additivity principle. We therefore express the components β〈Epp〉, β〈Ecc〉, and β〈Epc〉 as the average interaction per pair by dividing by the respective degeneracy in Table 6. The monotonic decrease in the average pair interactions in Table 6 is a general consequence of the fact that for larger cluster sizes there are more distant pair terms with diminished contributions to the average. We interpret the numbers in Tables 5 and 6 as follows. If we use the mean field approach, then the total energy of the system is given by the average quantities

β〈Esys〉 )

Np(Np - 1) Nc(Nc - 1) 〈Epp〉pp + 〈Ecc〉cc + 2 2 NpNc〈Epc〉pc (22)

If we now make the approximations Np . 1 and Nc . 1

Langmuir, Vol. 17, No. 26, 2001 8037

and the substitution Nc ) |Zp|Np, then we can rewrite eq 22 as

[

]

〈A〉pp Zc2〈A〉cc β〈Esys〉 ≈ (NpZp) + - |Zc|〈A〉pc 2 2 2

(23)

where 〈A〉ij are functions of the average separation distances between the ith and jth species. It is now clear that for β〈Esys〉 to be negative the value of 〈A〉pc must be larger than the sum of the other two terms given univalent counterions as in the present study. In other words, the counterions must statistically remain in the vicinity of the macroions. This is consistent with the restriction imposed on the inner orbital counterions about each macroion. Having established that it is the macroion-counterion interaction that results in a negative energy for the system, we now present an argument that it is nonlinear effects that are responsible for the stability of colloidal clusters. To do this, we focus on the failure of the nonadditivity rule for the trimer at an interparticle separation distance of d ) 3.5. If additivity obtains, then the interaction energy for the trimer, E1,2,3, is equal to the sum of the dimer interactions or 3 times that of the dimer at d ) 3.5. The mathematical expression thus may be written as

βE1,2,3 ) βE1,2 + βE1,3 + βE2,3 ) 3E1,2

(24)

According to the data in Table 4 it is readily apparent that β〈Epp〉(trimer) = 3β〈Epp〉(dimer). It is also evident from the data in this table that β〈Ecc〉(trimer) < 3β〈Ecc〉(dimer) and β〈Epc〉(trimer) < 3β〈Epc〉(dimer). That is, the additivity rule fails for the pairwise interactions involving the counterions. The reason for this failure is found in the values of rgtherm and Zgtherm/Zp in Table 1. In runs 6 and 11 the values of rgtherm are 1.317 and 1.333, respectively, with corresponding values of Zgtherm/Zp of 0.3033 and 0.2236. This means that for the trimer more counterions are found in the inner orbital than for the dimer, but their relative distances are further apart because of the larger value of rgtherm. As a result the average values of β〈Ecc〉 and β〈Epc〉 do not increase in proportion to the number of macroions in the cluster. 6. Orbital Interpretation of the Surface Charge Induced Reentry Transition It is now argued that the reentry from the two-state structure to the homogeneous liquid condition is consistent with the variation of the orbital structure of colloidal suspensions. We refer to runs 21, 22, and 23 in Table 2. For low surface charge densities of run 22 with Zp/ap ) 5 nm-1 the range of interaction inferred from rGtherm encompasses all of the macroions in the seven macroion cluster. Likewise in run 21 with Zp/ap ) 10 nm-1. However, in run 22 the total interaction energy of the system is β〈Esys〉 ) -67.6, whereas in run 21 it is β〈Esys〉 ) -284.7, or slightly larger than 4 times more attractive. This additional “attraction” is also reflected in the relative values of the reduced charge ratios at the two different environments of the seven macroion cluster. For run 22, for the smaller charge density, we have from Table 2 the ratio Zgtherm/Zp ) 0.2655, whereas for run 21 we have Zgtherm/Zp ) 0.1413. This means that more counterions are drawn into the interior of the cluster of run 21 than run 22. If we now focus on the values of Zgtherm/Zp, we find that this ratio is 0.267 for run 21 and 0.531 for run 22. We interpret this to mean that the counterions drawn into the cluster in run 22 are not as “tightly” associated

8038

Langmuir, Vol. 17, No. 26, 2001

Schmitz

obtains from current VT theories,11-15 and thus one might conclude that the region II f region III transition likewise cannot be explained by the VT approach. However, the VT theories draw attention to the thermodynamics of the system as a whole to explain phase separation behavior rather than focus on the effective pair potential. From this holistic point of view one must consider the variation of the multibodied interactions with the variation of the macroion configuration. Such variations must include the counterion disposition as well, as dictated by the well-known Gibbs-Duhem relationship, evaluated at the point r in the medium

dµc(r) ) -

Figure 11. The three stages of the liquid f crystalline f liquid reentrant transitions. Stage I: The macroions have a low surface charge density and exhibit a liquidlike structural configuration. This represents the “liquid” state in the phase diagram of Yamanaka et al.20 Stage II: This corresponds to the coexistent “crystalline-liquid” region of Yamanaka et al. The surface charge density is sufficiently large that counterions contribute to the effective charge of the macroions, as indicated by the outer shell about each macroion. Furthermore, this region represents the “overlap” of the neighboring ion clouds with extensive sharing of the counterions between macroions within the cluster. The counterions and solvent within the cluster can exchange with those in the solution region, as indicated by the arrows. Stage III: The surface charge density is very high, and there are a large number of counterions within the distance rgtherm from each macroion, as indicated by the darker shell surrounding each macroion. This large reduction in the effective charge (Zgtherm/Zp < 0.1) results in a dissolution of the clusters in stage II, and the homogeneous liquid state of the system results.

with the more distant macroions and hence this cluster is less stable than that in run 10. Hence, the colloidal structure in run 22 represents the homogeneous state of the system, whereas in run 21 it represents the more ordered crystalline state. If one now increases the surface charge density to that in run 23 (Zp/ap ) 20), the value of rGtherm now collapses to encompass only the central macroion in the cluster. The value of rGtherm ) 1.331 is comparable to the thermal radius about the exterior macroions, viz., rgtherm ) 1.209, and the macroion system again exhibits homogeneous behavior. A schematic of the proposed liquid f crystalline f liquid transitions in terms of the orbital model of colloidal systems is given in Figure 11. 7. Discussion Although there is little difference in the magnitude of the USAXS intensity for the three surface charge densities shown in Figure 1, these angle profiles clearly indicate that only the middle value of the surface charge density exhibits any degree of ordering of the macroions in the suspension. If one uses the depth of the minimum in the pair potential as a measure of the degree of order in the system, then as we have shown the SI potential does not explain the observation of a reentry transition from region II to region III. We have also shown16 that the SI potential

np(r) nc(r)

dµp(r)

(25)

The mean field approach is to reduce the multicomponent system problem to a one-component system by the integration over all of the microion positions prior to movement of the macroions to simplify the calculations but at the expense of counterion-counterion correlations and nonlinear effects. The shortcomings of the mean field theories, therefore, are the unproven assumptions that (1) the superposition approximation is valid and (2) the correlations between the counterions are negligible. The mean field theories, therefore, have highly restrictive applications. Sa´nchez-Sa´nchez and Lozada-Cassou38 and also Lo¨wen, Hansen, and Madden39 have shown that the counterions accumulate between two macroions to a larger degree than predicted in the linearized theories. Stigter40 pointed out that the linearized form of the Poisson-Boltzmann equation was strictly valid only in the “far field” region defined as separation distances larger than that for which the ratio of the DLVO potential was less than kBT. More recently Stevens, Falk, and Robbins41 pointed out that the Yukawa form of the pair potential and nonlinear Poisson-Boltzmann solutions may be valid for regions between the macroions but break downs near the surface of the macroion. In the present analysis we used a variation of the theme by defining a thermal radius, denoted by rGtherm and rgtherm, as being the distance from the macroion for which the total electrostatic interaction energy is equal to kBT. In view of the fact that rGtherm may encompass all of the macroions in the cluster (cf. runs 21, 22, and 23 in Table 2), the mean field expressions based on a linearized PB equation are clearly inadequate to characterize these systems. Although the proposed orbital model for colloidal systems does not rely on the concept of an “effective charge” or “thermodynamic radius”, these concepts nonetheless are valuable in understanding the physical properties of colloidal systems. The use of an effective radius to interpret colloidal properties is not new. It has been shown, for example, that by equating the DLVO potential to the thermal energy and taking the limit y(1 - r) , 1 in the present notation, the following form obtains for the thermal radius in nonreduced form:42

Rtherm ≈ ap + 1/κ

(26)

The relationship given by eq 26 has been used to define an “effective hard sphere radius” to interpret colloidal structures.43-46 Rather than use rGtherm and rgtherm to define (41) Stevens. M. J.; Falk, M. L.; Robbins, M. O. J. Chem. Phys. 1996, 104, 5209. (42) Schmitz, K. S. Phys. Chem. Chem. Phys. 1999, 1, 2109. (43) Doty, P.; R. F. Steiner, R. F. J. Chem. Phys. 1949, 17, 743. (44) Stigter, D. Recl. Trav. Chim. Pays-Bas 1954, 73, 593. (45) Okubo, T. J. Am. Chem. Soc. 1990, 112, 5420.

Crystalline-Liquid Transition in Colloidal Systems

excluded volume regions about the parent macroion, these thermal radii are used in our analysis as a measure of the “sphere of influence” of the electrical field about the central macroion. Overlapping of these spheres of influence denotes regions in which the counterions are shared by the participating macroions. By using rGtherm and rgtherm as a measure of the range of electrostatic interactions along with ZGtherm and Zgtherm as screened charges, an analogy is drawn between colloidal structures and molecular structures. This comparison is effected through the JPF and BD methods of characterization of colloidal clusters, which provides an orbital model for the description of the physical properties of colloidal systems. There are three types of “orbitals” identified in our simulations: (1) an “inner” orbital in which the counterions are strongly associated with one macroion, (2) a “shared” orbital in which the counterions act under the influence of two or more macroions in the cluster, and (3) “outer” orbitals in which the counterions are free to roam throughout the entire cluster. It must be emphasized that these orbitals are defined by the constant potential contours and vector gradient fields as proposed by Bader and co-workers33,34 and not by the thermal radii. Our simulations indicate that the electrical energy of the system, β〈Esys〉, is negative (cf. Tables 1-3) and therefore represents an overall attraction in the holistic characterization of the system. A similar result is reported by Lobaskin and Linse47 for 80 macroions with a macroion: microion charge asymmetry of 60:1. These authors did not partition the total interaction into its component parts β〈Epp〉, β〈Ecc〉, and β〈Epc〉. From our analysis we conclude (46) Okubo, T. Prog. Polym. Sci. 1993, 18, 481. (47) Lobaskin, V.; Linse, P. J. Chem. Phys. 1999, 111, 4300.

Langmuir, Vol. 17, No. 26, 2001 8039

that the dominance of β〈Epc〉 over the other two contributions is attributed to the inner orbitals surrounding each macroion. Finally, the transitions liquid f crystalline f liquid are interpreted in terms of the range of the electrostatic interaction for a cluster of macroions, as measured by the thermal radii rGtherm and rgtherm. Since the JPF and BD simulations used the “bare” macroion charge, we have thus avoided the introduction of “effective charges” and, therefore, “effective screening lengths” that are in common use in the mean field theories. One of the reviewers of this paper suggested that the results of the BD simulations for a system of a single macroion be compared with Monte Carlo and PB simulations. This comparison has in fact been done since submission of this paper, and the results will be reported in a later publication.48 The comparison indicates that for monovalent counterions all three methods give comparable results. However, when the counterions are multivalent, the PB simulations deviate the most from the other two methods, due to the omission of counterion-counterion correlations in the PB expression. Acknowledgment. I acknowledge the Kyoto University Foundation for providing partial support during my visit with Professor Matsuoka in the Department of Polymer Chemistry of Kyoto University, and the many stimulating discussions with Professor Matsuoka regarding these systems. LA010725W (48) Mukherjee, A. K.; Schmitz, K. S.; Bhuiyan, L. B. Submitted for publication.