Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Forces between like-charged walls immersed in electrolyte solution Otto, Frank 2000

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2000-565998.pdf [ 8.25MB ]
Metadata
JSON: 831-1.0061490.json
JSON-LD: 831-1.0061490-ld.json
RDF/XML (Pretty): 831-1.0061490-rdf.xml
RDF/JSON: 831-1.0061490-rdf.json
Turtle: 831-1.0061490-turtle.txt
N-Triples: 831-1.0061490-rdf-ntriples.txt
Original Record: 831-1.0061490-source.json
Full Text
831-1.0061490-fulltext.txt
Citation
831-1.0061490.ris

Full Text

FORCES BETWEEN LIKE-CHARGED WALLS IMMERSED IN ELECTROLYTE SOLUTION By Frank Otto Diplom (Physics) Free University Berlin (1994)  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in  THE FACULTY OF GRADUATE STUDIES (Department of Chemistry) We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA July 2000 © Frank Otto, 2000  In presenting degree  this  at the  thesis  in  partial fulfilment  of  University of  British Columbia,  I agree  freely available for reference copying  of  department  this or  publication of  and study.  this  his  or  her  Department The University of British Columbia Vancouver, Canada  DE-6 (2/88)  that the  representatives.  may be It  thesis for financial gain shall not  permission.  requirements  1 further agree  thesis for scholarly purposes by  the  is  that  an  advanced  Library shall make it  permission for extensive  granted by the understood  be  for  allowed  that without  head  of  my  copying  or  my written  Abstract  The effects of a molecular solvent on the forces between two infinite like-charged walls immersed in aqueous electrolyte solution are investigated. The solution models are chosen such that the numerical solution of accurate anisotropic integral equation theories is possible, and the anisotropic hypernetted-chain (AHNC) approximation is used in the present calculations. The resulting average particle densities and the pair correlation functions are used to analyze the force acting between the walls in detail. Mixtures of neutral hard-sphere solvent particles and counterions are employed to investigate the influence of the finite size of the solvent molecules. At wall separations of a few solvent diameters, even at relatively high surface charge and moderate solvent density, the ionic contribution to the force tends to be dominated by the hard-core or packing component. If the ions and solvent particles are of equal size, then the net pressure between the walls can be reasonably well approximated by adding the pressures of pure one-component ionic and solvent systems. However, if the ion and solvent diameters are significantly different, the pressure curve is more complex, and the hard-core component must be evaluated for a mixture of neutral hard spheres. The interaction of the walls is also investigated at the McMillan-Mayer (MM) level of description. In these models, the solvent is not represented by discrete particles but exerts its influence through solvent-averaged ion-ion potentials of mean force which serve as effective potentials. The approximations involved in applying M M theory to the inhomogeneous slit system are discussed. The wall-wall interactions obtained can differ dramatically from the primitive model (dielectric continuum solvent) case. Most interestingly, at the M M level, the force between like-charged  ii  walls at small separations and with realistic surface charges can be attractive for monovalent counterions, which is due to solvent effects on the effective counterion-counterion interaction. Several ion-ion potentials of mean force for "realistic" models for N a and +  C I in water are available in the literature. We find that the forces strongly depend on -  the model employed and different models for the same ion can give qualitatively different results. To construct a self-consistent model, ion-wall potentials of mean force must also be considered. Typical model ion-wall potentials show an additional counterion attraction towards the walls. Incorporation of these wall potentials yields an overall increase of the attraction between the walls. The possible relevance of our observations in the interpretation of experimental force measurements is briefly discussed. '  iii  Table of Contents  Abstract  ii  List of Figures  vii  List of Tables  xv  List of Important Symbols and Abbreviations Acknowledgement  xvi xx  1 Introduction  1  2 Equilibrium Statistical Mechanics  6  2.1  Introduction  6  2.2  Basic Definitions  7  2.3  Integral Equation Theories  10  2.4  Closure Approximations  12  2.5  General Model  14  2.6  A H N C Theory for an Ionic Fluid Between Planar Walls  17  2.7  Pressure Calculation  22  3 Mixtures of Neutral and Charged Hard-Spheres - Packing Effects  25  3.1  Introduction  25  3.2  Model and Remarks on the Numerical Method  26  3.3  Results  28 iv  3.4 4  29  3.3.2  Larger Ions  39  3.3.3  Monovalent Counterions  44  Summary  49 51  4.1  Introduction  51  4.2  Pairwise Additive McMillan-Mayer Theory  52  4.3  The Model  55  4.4  Results  57  4.4.1  Systems in Equilibrium with a Bulk at Infinite Dilution  58  4.4.2  Systems at Finite Bulk Concentration  65  4.4.3  Pair Distribution Functions  76  Summary  83  Further Investigations at the M c M i l l a n - M a y e r Level  84  5.1  McMillan-Mayer Results for Several Electrolyte Models  84  5.1.1  Models  85  5.1.2  Results  88  5.1.3  Summary  94  5.2  5.3 6  Particles of Equal Size  Molecular Solvent Effects at the M c M i l l a n - M a y e r Level  4.5 5  3.3.1  Ion-Wall Potentials of Mean Force  97  5.2.1  The Model  98  5.2.2  Results  102  5.2.3  Summary  104  Remarks on Some Experimental Results  Summary and Conclusions  106 111  v  Bibliography  117  Appendices  <  A Numerical Algorithm A.l  122  General Method  122  A. 2 MM Systems with Strong Anion-Cation Attraction  B Formulas for P i s  B. l  122  128  130  i t  Midplane Formulation  130  B.2 Calculation of  ft  A  B.3 Remarks on the Numerics  133 134  C P B and D L V O Theory  136  D McMillan-Mayer Theory  139  vi  List of Figures  2.1  The frame of reference used in the calculations. Note that this a twodimensional representation of the three-dimensional system. For clarity, only the case where all particles have the same diameter d = d is considered. 15 a  2.2  Schematic view of the walls immersed in the bulk fluid  16  2.3  Comparison of the electrostatic energy of a charge q far away from a charged macroparticle with a screening double layer around it for two macroparticle geometries. In (a) the macroparticle is an infinite planar wall with a surface charge density —a. The charge distribution of the double layer is centered at a distance A from the surface and cancels the wall charge exactly. The Coulombic contribution to the energy of a charge q far away (in the bulk) is independent of the distance from the wall and given by v = qa/A/(2ee ) el  0  • The corresponding case for a spherical macroparticle  in b) yields v = 0. Here, the total surface charge — Q is the negative of the el  charge associated with the double layer (but the surface charge densities of both charged spheres in (b) are different) 3.1  21  Density profiles for the mixture of ions and hard spheres (PM+HS), the primitive model (PM), and the pure hard-sphere model (HS). In all cases dion =  4 s and Z  ion  = +2. In (a) d n = 2.6d s and in (b) wa  h  vii  d  w a l l  = 3.1 d  h s  -  • • •  30  3.2  Contact densities for (a) hard-sphere and (b) ion components for different models. The model labeled U C P M + H S is identical to P M + H S but the charges are switched off as described in the text. The remaining models are as in Fig. 3.1. In all cases d[ = 4s and zT = +2 0n  3.3  32  ion  Densities at the midplane for (a) hard-sphere and (b) ion components for different models. The models are as in Fig. 3.2. In all cases d  i0Ti  = 4s and  zT = +2  34  ion  3-4  #ionion(7")  parallel to the walls in the midplane and contact plane for the  (rf = 4s) with d  P M and P M + H S 3.5  ion  w a l l  = 2.05 d  hs  and Z  ion  35  = +2  The net pressure acting between the walls obtained for the different models and with different superposition approximations. In all cases 4>n —4s and Z\on = +2- The labels are as in Fig. 3.2 and as discussed in the text. . . .  3.6  The electrostatic component of the pressure for the P M and P M + H S with ^ion = +2. Results for 4)n = 4s and for 4, = 1.52 4s are shown  38  n  3.7  37  Density profiles for the mixture of ions and hard spheres (PM+HS), the P M , and the pure HS model. In all cases 4m = 1-52 4s and Z\ = +2. In oxl  (a) 4aii = 2.5 4s and in (b) 4aii = 2.9 4s  3.8  40  Contact densities at the wall for (a) hard-sphere and (b) ion components for different models. The models are as in Fig. 3.2. In all cases d-  l0n  1.52 4  3.9  and Z  s  ion  =  = +2  42  The net pressure acting between the walls obtained for the different models and with different superposition approximations. In all cases d  10n  —1.52 4s  and zTj = +2. The labels are as in Fig. 3.2 and as discussed in the text. 0n  43  3.10 Density profiles for the mixture of ions and hard spheres (PM+HS), the primitive model (PM), and the pure hard-sphere model (HS). In all cases dion =  4s and Z  ion  =+l.  In (a)  4aii = 2.6 4s and in (b) 4aii = 3.0 4s- • • • viii  45  3.11 Contour plots of p i o n i o n ^ , 0, z) for the P M and the P M + H S model with  ^ion =  and Z\  on  = +1 for d n = 2.6 d . r and z are in units of d . wa  ns  n s  The black horizontal bars represent the walls and the legend on the right explains the colour scale, with the number between two squares indicating the boundary value  46  3.12 Contour plots of gionion(r, —z , , z) for the P M and the P M + H S model ma  with 3-13  G?ion =  ^ionion(^)  x  ^hs and zT; = + l for d n = 2.6dh - r and z are in units of d swa  on  s  n  47  parallel to the walls in the midplane and contact plane for the  P M and P M + H S (d  ion  = d ) with d hs  Vfan  = 2.0d  and Z  hs  i o n  = +l  48  4.1 Ion-ion pair potentials of mean force at infinite dilution for monovalent hard-sphere ions with T = 298 i f , e = 78.7, d = 2.8 A . Both M M and P M curves are shown 4.2 P  n e t  56  and its components for systems including only monovalent counterions  with a = -0.1307 e/d = -0.267 C / m . The Poisson-Boltzmann (PB) curve 2  2  is shown for comparison (see Appendix C), where d n[PB] is redefined to wa  give the same accessible space (i.e., the same z models. P e[MM] is less than 2 x 1 0  -5  k T/d  3  B  cor  m a x  ) as the finite size ion  throughout and is therefore  not shown. For comparison: 0.020fc T/d = 1.5 M P T = 3.7 x 1 0 N / m , 3  6  2  B  where R is the gas constant  59  4.3 P t and its components for systems including only monovalent counterions ne  with a= -0.0653 e/d = -0.134 C / m . The Poisson-Boltzmann (PB) curve 2  2  is shown for comparison (see Appendix C). P 10  -6  core  [ M M ] is less than 2 x  A;sT/d throughout and is therefore not shown 3  ix  61  4.4  P  and its components for systems including only monovalent coun-  n e t  terions with a = -0.1634 e/d  2  4 x 10 4.5  - 5  = -0.334 C / m .  P  2  core  [ M M ] is less than  k T/d? throughout and is therefore not shown  62  B  The a dependence of P  n e t  and P ff for systems including only monovalent e  counterions at fixed d n = 4 d  63  wa  4.6  and its components for systems including only divalent counterions  P t n e  with cr = -0.1307e/d  2  = -0.267 C / m . 2  The Poisson-Boltzmann (PB)  curve is shown for comparison (see Appendix C). P 10~  17  k T/d  3  B  and P  [ P M ] is less than 4x 10~ k T/d 4  core  core  3  B  [ M M ] is less than throughout, there-  fore these curves are not shown 4.7  64  Comparison of P u calculated via the contact theorem and via the mids  t  plane formula at a = —0.1307 e/d . Results are shown for a 1:1 electrolyte 2  solution in equilibrium with a bulk at 1.0 M , and for a system where only monovalent counterions are included. For the former case the derivative of the total grand potential (per unit surface area) with respect to d n is wa  also shown 4.8  66  Density profiles for 1:1 electrolyte solutions between uncharged walls, p is LOOM = 0.0132/d and 1.03M = 0.0136/d for the P M and the M M 3  3  model, respectively, p is indicated by horizontal lines. Results,are shown for d n = 1 d, 2 d, and 8 d. The profiles for cations and anions are identical. 69 wa  4.9  Density profiles for 1:1 electrolyte solutions between charged walls with a = —0.1307 e/d  2  = —0.267 C / m  2  for (a) counterions and (b) coions. p  is LOOM = 0.0132/d and 1.03M = 0.0136/d for the P M and the M M 3  3  model, respectively, p is indicated by horizontal lines. Results are shown for d n = 2.5 d, 4.0 d, and 6.0 d  70  wa  x  4.10 Midplane densities of monovalent counterions and coions at and  CTI  o w  n  —0.1307e/d  g  sities p are as in Fig. 4.9 and are indicated by horizontal dotted lines. . . 4.11 P  for 1:1 electrolyte solutions at p = 1.0M for cr j  n e t  h  ciow =  4.12 P  n e t  Chigh/2,  2  for (a) the M M system and (b) the P M . The bulk den-  = a i h/2 h  <7 igh =  7  = —0.1307e/d , 2  gh  and for uncharged walls  7  for M M systems in equilibrium with 1:1 electrolyte solutions at various  bulk concentrations and o — —0.1307 e/d = —0.267 C / m . The curves for 2  2  p = 0.5 M and 1.0 M are indistinguishable within numerical error and therefore only the latter case is shown 4.13 The dependence of P  7  on a at fixed d n — Ad. Results are shown for  n e t  wa  1:1 electrolyte solutions at p = 1.0 M and for the the infinite dilution (counterions only) case  7  4.14 $07( ) r  a  t the midplane parallel to the walls (z\ — z^ — 0) for 1:1 electrolyte  a = —0.1307e/d .  =  2  is shown in (a) and g  ++  and g  ("4-" denotes  counterion and "—" the coion) are shown in (b) 4.15 Contour plots of g++(r, — z  m a x  7  , z) for systems where only monovalent coun-  terions are present with d u = 3.6d and cr= —0.1307 e/d . 2  wa  r and z are in  units of d. The black horizontal bars represent the walls and the legend on the right explains the colour scale, with the number between two squares indicating the boundary value 4.16 Contour plots of g++(r, —z  7 , z) for systems where only divalent counte-  mSiX  rions are present with rf 4.17 Contour plots of p (z)g^. +  = 2.8rf and a = —0.1307 e/d  8  2  wall  (r, 0, z)d , the average counterion density around z  +  a midplane coion. The ions are monovalent and p = 1 . 0 M , d n = 4d, and wa  o — —0.1307 e/d  8  2  xi  4.18 Contour plots of p-.(z)g -(r,  — z , z)d , the average coion density around 3  +  max  a counterion in contact with one wall. The system parameters are as in Fig. 4.17 5.1  82  u ( P ) components of the ion-ion potentials of mean force at infinite dilueff  tion. The curves labeled [KP], [GRP] and [PR] are from Refs [52], [74,75] and [76], respectively 5.2  86  The dependence of P  and its components on the wall separation for  net  Na  +  [GRP] and N a  +  [PR] counterions.  —0.334 C / m = — le/48A . 2  2  The wall charge density cr =  M represents moles/litre and R is the gas  constant 5.3  89  The dependence of P  net  and its components on the wall separation for C l ~  [GRP] and CI" [PR] counterions. The wall charge density a = +0.334 C / m = 2  + l e / 4 8 A . The notation is as in Fig. 5.2  90  2  5.4  Density profiles for the systems shown in Figs. 5.2 and 5.3 at the wall separation D = 7.25A. The curve for N a [PR] is within ±0.5 M of the +  result for monovalent charged hard spheres (PM) with d  10n  = 2.94A. The  notation is as in Fig. 5.2 5.5  92  Ion-ion pair distribution functions at the wall plane (z\ = z = D/2) for 2  the systems shown in Figs. 5.2 and 5.3 at the wall separation D = 7.25A. Curves obtained for E q [KP] and monovalent charged hard spheres with +  d  i o n  = 2.94A (PM) at the same wall charge are also shown. The C l ~ [PR]  curve reaches a maximum of 8.3 5.6  P  n e t  93  for the NaCl [GRP] model [74,75] at different salt concentrations,  p, and cr = —0.334C/m = - l e / 4 8 A . The inset shows the tails on an 2  2  expanded scale. The notation is as in Fig. 5.2  xii  95  5.7  Comparison of P  for four systems with monovalent counterions, includ-  n e t  ing the P M with d  = 2.94.A, and a = - 0 . 3 3 4 C / m = - l e / 4 8 A . 2  ion  The  2  curves for Na+ [GRP] and N a [PR] are as in Fig. 5.2 and the Eq+ [KP] +  result is from Fig. 4.4 . The notation is as in Fig. 5.2 5.8  Two model ion-wall potentials for counterions, V+>  a s  96 functions of the  distance to the planes of closest approach to the walls, £. For C= (z  max  z)>ld  and for coions at all distances tp = 0. For x = C/d<l a  ±  the curves  are given by the polynomials 2.75a; — 7.75a; + 7.25a; — 2.25 (-01) and 3  2  71.5a; - 145.9a; + 91.0a; - 9.3a; - 7.7a; + 0.4 (^2) 5  5.9  4  3  101  2  P t for monovalent counterions and a wall charge o = —0.1307e/d , with 2  ne  curves for the M M model (see Fig. 4.1) and the P M shown. The ionwall potentials of mean force are either given by the ip in Fig. 5.8 or no additional short-range wall potential is used (-0 = 0). The triangles show the results for the M M model with wall potential ipl in equilibrium with a bulk of salt concentration of 1.0 M . The size of the triangles roughly indicates the magnitude of the numerical error due to the finite grid size.  103  5.10 Comparison of the pressure components for two systems with monovalent counterions interacting via the M M ion-ion potential of mean force with a = —0.1307 e/d . 2  The lines and symbols indicate the results for effective  ion-wall potential ipl and for the case with no additional short-range ionwall potential (•0 = 0), respectively  xiii  105  5.11 Schematic view of the pressure curves obtained in some SFA experiments. The curve is thought to result from three (additive?) contributions: a monotonic repulsive component, an oscillatory component due to packing of the solvent, and an attractive component with a minimum around the wall separation corresponding to point " E " . Only the solid lines represent values accessible for measurement. Once point " A " is reached, a further increase of the external force results in a sudden change of the wall separation corresponding to a point between " B " and " C " . However, an increase of the external force at point " C " yields a jump to a wall separation corresponding to a point between "F" and " G " and the local maximum at  A.l  " E " cannot be detected experimentally  110  Iteration scheme for the algorithm  126  xiv  List of Tables  5.1  F i t parameters for the Pettitt-Rossky [76] ion-ion potentials of mean force. The units are chosen to yield u  eS  in fc T if R is given in A B  xv  87  List of Important Symbols and Abbreviations  Symbol  Description  a  activity  b  bridge function  c  direct correlation function  D  width of the accessible space for particles in the slit  d  hard-sphere diameter  d n  wall separation  e  elementary charge  g  pair distribution function  h  pair correlation function  h  Planck constant  k  coordinate parallel to the walls after Hankel transformation  I, I  labels of the walls  ks  Boltzmann constant  m  mass  N  number of particles  wa  7Y  number of grid points parallel to the walls  N  total number of layers (grid points perpendicular to the walls)  P  pressure  Pbuik  bulk pressure  Pnet  pressure difference between the slit and the bulk  r  z  xvi  P iit  pressure in the slit  q  charge of a particle  R  position vector  R  particle separation  r  cylindrical coordinate parallel to the walls  s  r  m a x  cut-off length parallel to the walls  T  temperature  UN  potential energy of n particles due to particle interactions  u  potential energy of a pair of particles  u  effective potential component of w due to molecular solvent  u  Coulomb pair potential  eS  el  M  H S  hard-sphere pair potential  u  short-range component of pair potential  V  volume  VAT  potential energy of n particles due to the external potential  v  potential energy due to the external field (wall)  sh  u  el  Coulomb interaction with the walls  v  HS  hard-sphere interaction with the walls  v  short-range wall potential  W  one-particle potential of mean force  w  pair potential of mean force  w^>  n-particle potential of mean force  x, y  Cartesian coordinates parallel to the walls  y  correlation function h — c  Z  ion valence  sh  xvii  z  coordinate perpendicular to the walls z-coordinates of planes of closest approach to the walls  i-^max  a,  7,  5  indices for particle species  7  activity coefficient  Az  layer width  e  dielectric constant  eo  vacuum permittivity  A  thermal wavelength chemical potential excess chemical potential grand potential short-range wall potential measured from z  p  one-particle density  p  {n)  n-particle density function  a  surface charge density  m a x  grand partition function  c  distance from plane of closest approach to a wall  V  gradient operator  /  Hankel transform of f  AHNC  anisotropic hypernetted-chain  BGY  Born-Green-Yvon  HNC  hypernetted-chain  HS  hard-sphere  M  moles/litre  —*  xviii  MD  molecular dynamics  MM  McMillan-Mayer  MSA  mean-spherical approximation  PB  Poisson-Boltzmann  PM  primitive model  RHNC  reference hypernetted-chain"  RHS  right-hand side  SFA  surface force apparatus  xix  Acknowledgement  It is my pleasure to thank my supervisor, Dr. Gren Patey, for the excellent academic guidance, for his patience, and for making the past five years a great learning experience. Lots of postdocs and grad students were (or still are) part of the group and have helped to create an enjoyable workplace: John, Mark, Gary, Natalie, Claudio, Richard, Tyler, Liam, Hans, Phil, Vladimir, Sabine, and Chris; with special thanks to those who bore with all the German coffee-klatsch for the last two years. I am grateful for two years of U G F funding by the University of British Columbia and for financial support by the Department of Chemistry.  My time here in Canada has also been a fascinating cultural adventure, and many people contributed to make me feel "at home" in Vancouver. It is hard to think of a way to thank for all the emotional support and love I have got from that special someone in my life for the past three years - Jon, you've been the sunshine of my life. A big merci to Jean-Marc for experiencing together the peculiarities of Canadian culture. And a special danke to my sister Silke, my mum Monika, and my dad'Wilfried for the long-distance support and many "Care packets", and to all my other friends back in Berlin, who demonstrated that it is possible to live far away and still be in constant touch.  xx  Chapter 1 Introduction  The interaction of macroparticles immersed in electrolyte solution is influenced by many factors. Properties of the macroparticles themselves (charge, shape, material), of the solution (both ions and solvent), and of the surface (roughness, interaction with the electrolyte) all play a significant role [1,2]. A n electrical double layer close to the surface of a macroparticle in an aqueous environment can form if two ionic species, having charges of opposite sign, have different affinities for physical or chemical binding with the surface. Both ion dissociation from an originally neutral surface and preferential ion absorption from an electrolye solution onto the surface yield a net surface charge density. A certain number of counterions, sufficient to neutralize the surface charge, remains close to the macroparticle but dispersed in the solution, balancing electrostatic attraction and entropy. The average force acting between macroparticles in electrolyte solutions plays a fundamental role in determining the behaviour of colloidal systems. A n understanding of how to control the transition from a dispersed state to a coagulated state of the colloidal particles is important for many technical applications such as the manufacturing of paints and inks. In soil science, one seeks insight into the mechanisms behind the swelling of clays. Important implications for biological systems arise when considering the formation of biological membranes, the interaction of macromolecules, and even protein folding. A typical simplified system used to examine these features theoretically consists of a model  1  Chapter 1. Introduction  2  electrolyte solution between two structureless planar walls with a constant uniform surface charge density in an overall electroneutral setup. The forces acting between the walls depend on details of the interacting double layers as the surfaces are brought together, and are very sensitive to the underlying model and the approximations applied. Sometimes the results can be rather unexpected. A good review of many theoretical aspects of bulk electrolyte solution models, their confined counterparts, and the electrical double layer has been given by Attard [3]. Continuum solvent solution models, where the ionic interaction is reduced by a separation independent dielectric constant, have received considerable attention and their influence on wall-wall forces is now reasonably well understood [4]. The classical PoissonBoltzmann (PB) theory treats point ions at a mean-field level and always yields a repulsion between like-charged macroparticles, which decays exponentially with increasing wall separation. B y dropping the mean-field approximation and considering finite size ions, one arrives at the primitive model (PM), for which the exact results at small slit widths can be contrary to the predictions of the P B theory. For example, it is now well established [4] that like-charged walls immersed in a P M solution with divalent counterions can experience an attractive interaction at small wall separations. The attraction arises because of correlated charge fluctuations in the electroneutral slit system, which are ignored by the P B mean-field approximation. The possible importance of this mechanism was pointed out by Oosawa in 1971 [5], and it was first demonstrated by Patey in 1980 [6] that the magnitude of the effect might be big enough to outweigh the repulsive interaction component. The latter calculation overestimated the attraction, as an approximate integral equation theory was employed. However, major improvements of numerical methodology and computational power in the last two decades have helped to confirm the importance of Oosawa's insight. Related interesting developments for models that employ P M electrolyte solutions  Chapter 1.  Introduction  3  are, for example, the attractions found between charged spherical macroparticles at finite macroparticle concentrations [7], and additional macroparticle attractions for models that allow for a variation of the wall charge via ion absorption potentials [8]. Recent experimental results of attractions between two confined macroparticles and the failure to explain this phenomenon theoretically have been summarized by Grier [9]. Note that the effect of attractions between charged macroparticles (due to correlated charge fluctuations) is of different origin than the long-range attraction found experimentally between hydrophobic surfaces [10,11]. In the latter case, the liquid between the macroparticles is destabilized with respect to the vapour phase below a certain macroparticle separation, due to the unfavourable contribution to the free energy from the surfaces (i.e., the surface tension) [12]. The mechanism of the attraction between hydrophobic walls can be attributed to a slow decrease of the average fluid density in the slit as the walls are brought together until sudden cavitation occurs [13]. As significant as the effects found for P M solutions between two walls are, it is thought that the influence of a discrete molecular solvent will be at least as important.  The  concentration of water molecules at room temperature is about 55 moles/litre (M), which is much higher than the salt concentrations typically considered (^ 1 M ) . Therefore, continuum solvent models ignore the influence of particle packing on the force between the walls. Moreover, the restrictions the presence of a dense component places on the distribution of the ions might become considerable. Another important aspect of the molecular nature of the solvent is that at small ion-ion separations the assumption of a distance independent dielectric constant breaks down. The relative lack of dielectric screening and preferential solvent ordering around the ions cause deviations from the asymptotic Coulombic part, which are significant for separations of a few ion "diameters". Similarly, the interaction between an ion and the wall (or macroion) is also altered by the solvent molecules, but now, close to a wall, the structure of the solvent itself (density,  Chapter 1. Introduction  4  orientation, ...) is different from that of the corresponding bulk solvent. Clearly, such effects require more attention if all factors contributing to macroparticle interactions are to be understood. The present thesis is a contribution towards that end and addresses some aspects of discrete solvent effects. The determination of the equilibrium structure of a confined electrolyte consisting of "realistic" models for the ions and solvent particles remains a difficult computational task. Solving the problem exactly via computer simulation methods is not feasible for the systems considered here, especially if one wishes to explore the influence of the different model parameters in detail. Both the technical difficulties due to the strong interactions between ions and water molecules [14,15] and the different orders of magnitude for the numbers of ions and solvent particles render direct simulation impractical. Some effects of the discrete solvent can be studied by employing McMillan-Mayer solution theory, where the ions interact via "effective" potentials. These potentials of mean force include averaged solvation effects due to the molecular solvent and circumvent the need for explicitly considering a solvent species in the numerical calculations. Although simulation of such a system is much less demanding than the corresponding full model calculation, it remains very advantageous to use an approximate but accurate theoretical treatment for our investigations. The choice of an appropriate approximate method poses a problem common in the theory of fluids, in that a good compromise between accuracy and computational efficiency must be found. We employ a sophisticated integral equation theory to calculate particle densities and pair correlations accurately enough to obtain reliable results for the forces between the walls. The remainder of this thesis is divided into 5 chapters.  In Chapter 2, we give an  overview of the theoretical foundation of our method, which relies on the principles of statistical mechanics.  For a general electrolyte solution model and slit geometry, we  present the basic set of equations to be solved numerically and discuss the calculation of  Chapter 1. Introduction  5  the pressure (force per unit surface area) between the walls, as this quantity will receive the most attention in our analysis. Chapter 3 describes a study of how the finite size of the solvent molecules can influence the force between the walls via packing effects. The ion-ion pair interaction is not directly altered, but ions have to share the slit with a dense solvent component and are therefore excluded from a certain portion of the space. In Chapters 4 and 5, we investigate the slit pressure with the help of McMillanMayer solution theory, where "effective" ion-ion potentials are employed to represent the influence of the molecular solvent. Possible consequences of the approximations made in translating the McMillan-Mayer ideas into a theory for the slit system are discussed. The introduction of "effective" ion-wall potentials poses particular problems. In Chapter 4, we present results for a relatively simple solution model, ignoring effects due to the effective ion-wall potentials. The fluid structure is investigated in detail and the dependence of the net pressure on the wall charge and on the bulk salt concentration is discussed. We employ several more "realistic" solution models in Chapter 5, and the effect of the solvent-induced wall potentials is also considered. The relevance of our calculations for the interpretation of experimental results, obtained, for example, with the surface force apparatus, is briefly discussed. A summary of our findings is given in Chapter 6 together with some general conclusions.  Chapter 2  Equilibrium Statistical Mechanics  2.1  Introduction  Statistical mechanics provides the molecular motivation for the thermodynamic behaviour of a system. For over a century, it has successfully been employed as a microscopic theory for the properties of macroscopic objects.  The word "mechanics" points to the  fundamental belief that certain laws (either classical or quantum in nature) determine the microscopic evolution of the system, whereas "statistical" puts the impossibility of actually calculating the deterministic variables (either the positions, orientations, and momenta for all particles or the overall wavefunction) in perspective. Due to the large number of particles involved, the statistical point of view helps to simplify the description of the system, concentrating on average quantities and their fluctuations. The distinct states of matter lend themselves to different theoretical approaches. For instance,the low-temperature solid is well represented by the solid lattice (determined by the overall potential energy) with deviations from this ideal state described by phonons. The dilute gas, which is dominated by the kinetic contribution, can be successfully treated by perturbation techniques starting from the ideal gas. For liquids and,fluids however, 1  for which the magnitudes of the average kinetic energy and the average potential energy are comparable, no simple ideal model is available and a variety of more refined theoretical treatments is employed [16,17]. The formulation of perturbation theories, for example, is Above the critical temperature no distinction between gas and liquid states can be made and we refer to the state of matter as "fluid". 1  6  Chapter 2. Equilibrium Statistical Mechanics  7  still possible (but more complicated than for the dilute gas) starting from the observation that the structure of a liquid is usually dominated by core repulsion. In some cases an analytically solvable model, often using the mean-spherical approximation (MSA), plays the role of a convenient starting point. This concept can also be applied to ionic fluids, which are of special interest here. A more direct way of dealing with the multi-body problem is often preferred, namely, the use of "computer experiments" to solve the problem at hand "exactly". Two methods for classical systems exist [18]: molecular dynamics (MD) simulations, involving the solution of Newton's equations of motion for a set of N particles, and Monte-Carlo (MC) simulations, where phase-space configurations are sampled according to their probability of occurrence (without creation of a "real" time evolution). "Exact solution" in this context means within statistical error and subject to certain boundary conditions (usually periodic), as M D and M C methods are necessarily restricted to simulate finite systems. Exploiting the symmetry of the system directly, integral equation theories or density functional methods yield only approximate results, but are in general much faster to compute than M D / M C simulations for identical models. Approximate theoretical approaches still have to be complemented by simulation results, in order to assess the severity of the simplifications made. The results presented in this thesis have been calculated using a sophisticated integral equation method, which is known to be accurate for ionic fluids.  2.2  B a s i c Definitions  We assume in the following that only two particle species, A and B, are present, that classical mechanics is appropriate to describe the behaviour of the system, and that only pairwise particle interactions have to be considered, which do not depend on the orientation of the particles. The definitions are straightforward adaptations from standard  Chapter 2. Equilibrium Statistical Mechanics  8  formulas given, for example, in [17,19]. The Hamiltonian, ri , of a system consisting of N  —*  ./V particles, with coordinates Ri, momenta pi, particle types a,, and masses m , is given ai  by H ({Ri},{pi})  =  N  £„({&})+  V ({&})+U ({R }) N  N  i  = £ ^ +£ ^ ) +£ £ W ^ ) > i=l '' ai Z  (2- - ) 2  1  i=l j>i  i=l  i  where /CAT, VAT, and UN are the contributions due to the kinetic energy, an external potential (v (R)),  and the pair particle interactions (u (Ri,Rj)),  a  grand partition function,  E(p ,  p ,V,T),  A  A  + NB(3pB) i N  =  £  N  N +N =0 A  A  b  ^  m  ^A-^B-n  R  r  nj B\ (t  a  em-^N({Ri},{pi})\  j  B  l  B  (2.2.3)  a  B  a  (2-2.2)  and the activity for species a is given by a = e^/Al.  V, T, k , h, and p  r-ivi  N  a  B  P  d  e  where N + N = N, j3 = (k T)~ , A  SN  d  M~PV ({P^}) - PUA{^})\ ,  dRN  A  r,^f,  H3  ^TJ^J a  The  is defined by  B  exp(N pp -  respectively.  ai  are the volume, the temperature, the Boltzmann constant, the  Planck constant, and the chemical potential of species a, respectively. A = Q  .  h  ,  (2.2.4)  is called the thermal wavelength and has to be small compared to typical particle separations to justify the restriction to classical mechanics. All thermodynamic information about the system can be obtained by employing the appropriate link between statistical mechanics and thermodynamics. Here, the link is given by relating the grand potential Sl{l*>A,V>B,V,T) and the grand partition function via Q(p ,p ,V,T) A  B  = -k T\nE(p ,p ,V,T) B  A  B  .  (2.2.5)  Chapter 2. Equilibrium Statistical Mechanics  9  W i t h this relationship, the system's pressure, P, for example, is given by  -{dv)  = M  p==  -  kBT  T  --  (2 2 6)  T  The outlined route is usually not practical though, as often H cannot be calculated with sufficient accuracy for complex systems. For M C methods, the basic idea is to use the exponential term in Eq. (2.2.2) as a statistical weight to evaluate averages of variables, according to the ensemble of interest, "exactly" (if statistical noise and problems due to boundary conditions are disregarded). density functions, p^[„ (Ri,  A different approach employs the n-particle  Rn), which are obtained by integration over (N — n)  an  particle positions. The underlying idea for this procedure is that the low-order functions are generally sufficient to describe the fluid structure accurately enough for the calculation of thermodynamic properties. Only the functions for n= 1, 2 will be considered in this thesis. Assuming particle 1 being of species A, one has for the 1-particle density 1 PA(RI)  oo  oo  4E  E  N  N  A  B  Z  A (  N  ,  N  N  dR ...dR 2  e  N  .  -W"W  (2.2.7)  The AA and AB 2-particle densities are given by i  p!&(£i,&)  oo  N  oo  =, E E  *  - N =2n =0 a  i  fikRiX)  ("a  b  oo  E  (  The pair distribution functions, g , ay  3  n  d  ( 2  .2.  8 )  R  1  dh\...dR  V  N  p {Ri)p {R ) 1  e~W»W  .(2.2.9)  J  J^S^ll^L —-=ra  (oo-in\ (2.2.10)  •  2  9a-y(Ri, ^2)^7(^2) is the particle density of species 7 at R  2  a particle of species a is located at Ri.  a  B  B  tufr. -iy.{N -iy. are defined via  2  e-flv*-"*)  N  N  N  n  gn (K(P ,H ) P \ -l  dR ...dR  B  B  ai  , NA  - N =I N =I A  B  ~ 2j!/V ! j  oo  4E  N  A  a  under the condition that  If \Ri — R \ is large compared to the typical 2  length over which particle interactions are significant (and the fluid is not close to a  Chapter 2. Equilibrium Statistical Mechanics  10  critical point) then ^ ( P i , P ) ~ 1 and the particle densities at R\ and R a 7  2  2  are not  correlated. Pair distribution functions can often be measured experimentally via X-ray or neutron scattering, at least for bulk fluids, where g y{Ri, R2) = ^ ( l - R i — -R2I) is called a  the radial distribution function. One motivation to employ theories that calculate radial distribution functions directly is therefore the possible comparison with these data. This also facilitates the "fitting" of parameters for "realistic" models.  2.3  Integral Equation Theories  The objective of the integral equation scheme presented here is the (approximate) calculation of the functions p (R) a  and g ^(Ri, R2), whereas higher order particle correlations a  are not explicitly considered [but implicitly included in Eq. (2.3.3)]. For this purpose, it is convenient to first introduce the pair correlation function, h (R\, ai  ai  R ) = g y{Ri, R2) — 1 j 2  via the integral relationship  aj  h -y(Ri, R2) = c ~f(Ri, R2) + ^2  dR c s(R\,  a  5  (2.3.1)  a  and the direct correlation function, c ,  a  h ,  3  a  R3)ps(R3)h$ (R , 1  3  2  R2) •  (2.3.2)  J  The latter equation is called the Ornstein-Zernike (OZ) equation, which was introduced in 1914 to study critical density fluctuations in fluids [20,21], and further proved to be extremely useful for investigating fluid structures in general. Eq. (2.3.2) is called more precisely the inhomogeneous OZ equation, as the influence of an external field is accounted for. Without an external field (v (R) = 0) p is translationally invariant and a  a  the 2-particle correlation functions depend only on \R\ — Ri\. The resulting simplified version of Eq. (2.3.2) is referred to as the homogeneous (or singlet) OZ equation. 2  Note, that c ^(Ri,R2)  can be alternatively defined via the functional derivative of \n[p (Ri)/a ]  a  0v (Ri) a  a  with respect to p ( ^ ) 7  a  +  Chapter 2. Equilibrium  Statistical  In order to calculate h  ai  11  Mechanics  for a given set of p , another equation relating h a  ai  and c  a 7  is needed. Several authors studied these relationships for homogeneous systems in the late 1950's and 1960's, using diagramatic and functional differentiation techniques (for 3  references see [17]). Relevant contributions for the inhomogeneous case have been made by Morita and Hiroike [22] and Percus (II-3, Appendix D in [21]). As a general result one can write /* (Pi,i? ) + l = e x p [ - / ? ^ ( ^ ^ Q7  , (2.3.3)  2  where b  is called the bridge function.  ai  Although certain features of the bridge function  are known (it is short-ranged for hard-sphere (HS) systems [23]), the exact calculation of b ^(Ri, R ) is not possible, as an infinite sum of "bridge diagrams" is involved. Hence, all a  2  integral equation schemes have to replace Eq. (2.3.3) with an appropriate approximation, which will be discussed in the subsequent section. Note, the external field can in many cases be considered as arising due to a particle itself. Then, both the homogeneous and inhomogeneous approach are possible choices for theoretical levels of description. The inhomogeneous theory yields in general more accurate results, as the approximation for Eq. (2.3.3) is made at effectively a 3-particle level and more bridge diagrams are included in the description of the fluid structure. However, the numerical evaluation of the correlation functions is more demanding than for the homogeneous theory. The inhomogeneous case, where p {R) is not a constant, requires an additional equaa  tion for the average particle density. One option is the first member of the Born-GreenYvon (BGY) hierarchy, derived by taking the gradient with respect to Ri of Eq. (2.2.7) yielding V i l n p a ^ ) = -BVxVaiRx)  - 8J2 [dR p (R )g (R ^V.u^R,,R ) 2  7  y  2  aj  u  2  ,  (2.3.4)  J  Diagrams are a convenient graphical representation of the multi-dimensional integrals arising in expressions such as Eqs. (2.2.7), (2.2.8), or (2.2.9). 3  Chapter 2. Equilibrium  Statistical  Mechanics  12  although the Wertheim-Lovett-Mou-Buff ( W L M B ) equations [24], which are of similar integrodifferential type, are often preferred in numerical calculations [25]. A different route involves the calculation of the excess chemical potential, p = p}^ (R) l  a  + p *(R)  = k T[3\nA  e  B  + \n (R)}  a  defined via  X  + v (R)  Pa  p a{R) represents the change in £1(PA, PB,V,T)  —* /^ (PL),  + p™(R) .  a  (2.3.5)  due to particle interactions when a parti-  e  —*  cle of species a is added at R, and can be calculated via a coupling parameter integration [26]. The addition of the particle is described by the coupling parameter 0 < £ < 1, where £ = 0 means no interaction with the rest of the system, u (R\,  i?2|£ =  ai  stands for the fully "switched on" particle, u (Ri,  R-2\l) = u (R\,  ay  define g (Ri,R,21£) ai  ai  0) = 0, and £ = 1  R ).  We further  2  as the corresponding pair distribution function when particle a at  Ri is coupled to the system according to £. Increasing the coupling strength from £ to £ + d£ raises the interaction energy with particle 7 at R  2  by d^du (Ri,  R \^)/d^  ai  2  and  integration over the whole system for the total coupling process, £ = 0 —> £ = 1 , yields p *(Ri) e  such that = E /dR  Eq. (2.3.5) gives for p (R) a  where W (R) a  2  (R2) J\  d  u  Pl  « ^  m  )  R \0  a  2  •  (2.3.6)  p (R) a  = a exp{-(3[v {R) a  + p,f(R)]}  a  = a exp{-0W (R)} a  a  ,  (2.3.7)  is the potential of mean force acting on particle a located at R. Introduction  of the activity coefficient, j (R), a  via a = p (R)j (R) a  a  a  1*(R) = exp{/3W (R)} a  2.4  g ,(Ru  yields therefore  .  (2.3.8)  Closure Approximations  The choice of the approximation for Eq. (2.3.3) determines the accuracy of the resulting correlation functions, depending on the system at hand. Setting b -y(Ri,R.2) = 0 in Eq. a  Chapter 2. Equilibrium  Statistical  (2.3.3) gives the hypernetted-chain  (HNC) closure  hay(Ri, R ) + l = exp[-pu (R 2  ay  13  Mechanics  R ) + h {Ri,  u  2  R ) - c (Ri,  ai  2  R )} .  ai  (2.4.1)  2  This closure yields good results when the system is dominated by long-range forces [27,28] such as the 1/R  2  Coulomb forces, and is therefore selected for the calculations presented  here. Another advantage of Eq. (2.4.1) is that the coupling parameter integration for p™(R),  Eq. (2.3.6), can be performed explicitly [29,26]. One finds  PtC(Ri)  = E/dR p {R )[hl (R R ) 2  2  y  7  u  2  - c (R R )] aj  u  2  - ^[c^RuRJ  + 1] . (2.4.2)  The H N C approximation can be applied both at the homogeneous and inhomogeneous level. The latter level of theory has been named anisotropic hypernetted-chain  (AHNC)  theory by Kjellander and co-workers to distinguish it from schemes that calculate the inhomogeneous particle density via an equation similar to Eq. (2.3.4) but inserting the homogeneous correlation functions (at an appropriate density) [30]. More sophisticated closures employ approximations to the bridge functions of a suitable HS system, which can usually be calculated with higher accuracy than those of the original system. A n example of this strategy is the reference hypernetted-chain (RHNC) theory (see [25] for an application at the A H N C level). For ionic fluids this scheme works only as long as b (Ri,R ) ai  2  is mainly influenced by the core-core interaction between the ions, which is  true at high temperatures. For HS fluids and some other systems without long-range forces the Percus-Yevick (PY) closure c^(R R ) u  2  = ( l - vcp{(3u {R ai  u  R )}) 2  g (R ai  u  R) 2  (2.4.3)  is preferable over the H N C closure [31] (although fewer diagrams are taken into account) and yields excellent results for inhomogeneous systems [32,33].  Chapter 2. Equilibrium  Statistical  14  Mechanics  When particles interact via a hard-core potential ( u ( i ? i , R ) = oo, if \Ri — R \ Q7  2  <d ),  2  ai  the M S A closure is valuable as for some important models (including the restricted primitive model) an analytic solution can be found. This type of closure consists of setting c (Ri,  R ) = — (3u (Ri,  ai  2  ai  R ) outside the hard core and exploiting <? (i?i, R ) = 0 2  Q7  2  inside the hard core.  2.5  General M o d e l  We consider model fluids confined to a slit by two parallel, identical, and infinite walls. For the most part, the fluid represents an electrolyte solution in which the ions interact via an "effective potential", including effects due to the surrounding solvent to some degree, but without the need to model the solvent molecules explicitly. Only one model employs a simple neutral hard-sphere representation of the solvent particles.  Specific  details of the interaction potentials and in particular their motivation will be discussed later. Here, we give the most general form of the particle-particle and the particle-wall interactions encompassing all systems studied in this thesis. These are used to derive the basic formulas employed in the numerical calculations. Associated with a fluid particle of type a are the hard-sphere diameter d  a  and the  charge q = Z e, where e is the elementary charge. Particles of species a and 7 at a a  a  distance R interact via the pair potential u (R),  which can be expressed in the form  ai  u (R) ai  =  u™(R) + u%(R) + u%(R)  =  <  00 , «9L^  I AirtocR  if R < \(d  (2.5.1) a  +  U  sh  ( i 2 )  >  + d) , 7  otherwise,  '  where e is the dielectric constant of the solvent and eo is the vacuum permittivity. The terms u^(R) while u^(R)  and u^(R)  are the hard-sphere and Coulombic interactions, respectively,  is an additional short-range contribution. A n electrolyte model comprising  Chapter 2.  Equilibrium  Statistical  15  Mechanics  d/2 t T  Z  max  z =0  D  r  Z2  -z max  t  Figure 2.1: The frame of reference used in the calculations. Note that this is a two-dimensional representation of the three-dimensional system. For clarity, only the case where all particles have the same diameter d = d is considered. a  only the two former terms  (M^  h  7  (i?) =  0) is called the primitive model (PM), in which the  solvent is reduced to a background dielectric continuum. The confining walls are located at a separation  Hard-core particle-wall in-  d aiiW  teractions are assumed and therefore the space accessible to particle a has the width D = d \i — d . Each wall is charged with a surface charge density a and taken to have a  wa  a  the same dielectric constant as the solution, which circumvents the need for treating image charge effects, although this could be done with the anisotropic integral equations scheme in a rather straightforward manner [34]. In fact, for the P M it has been shown [35] that image effects are not very important at high surface charge densities ( < 7 > 0 . l C / m ) . 2  In our calculations we employ the reference frame shown in Fig. 2.1 with coordinates z and r = \/x  2  + y , where z = 0 and z = ±z 2  maiXiCt  (z ,x, = D /2) ma  a  a  designate the midplane  between the walls and the planes of closest approach to the walls, respectively. The interaction energy of a particle with both walls (I and II) is given by (2.5.2)  Chapter  2. Equilibrium  Statistical  16  Mechanics  ©  bulk ©  o  o  o  o  Figure 2.2: Schematic view of the walls immersed in the bulk fluid. where a constant from the electrostatic contribution has been dropped as it is canceled by an identical term from the rest of the electrolyte solution due to the imposition of overall electroneutrality of the slit system [34].  f  a  h , /  ( z ) =w  Q  h  '  7 /  is an additional short-range  (—z)  particle-wall interaction. This wall potential is more conveniently described by shifting the origin of the coordinate system to the contact planes, i.e., by introduction of the + z a x , | ) , and  function ip , with ^ (0) = ^ ( - ^ m a x . a ) = < ' (^max,a), vf' (z)=ip (\z h  a  va S  ,n  Q  (z)=if) (\z — z a  /7  I  a  \).  maXiQ  m  a  The assumptions implicitly made by writing this contribution  in the above way is that the effects from both (identical) walls are superimposable and that ip does not depend on d n itself. a  wa  The fluid between the walls is assumed to be in equilibrium with a bulk system, which is shown schematically in Fig. 2.2. This fixes the chemical potential, p , for each species, a  Chapter 2.  Equilibrium  StatisiicaUMechanics  17  or, in the case of an electrolyte solution, p±. Specifying the p is equivalent to fixing the a  bulk densities, p .  The walls are assumed to be wide and long enough, so that neither  a  interactions between slit and bulk particles nor boundary effects have to be considered.  2.6  A H N C Theory for an Ionic Fluid Between Planar Walls  We now apply the A H N C integral equation theory to the model and geometry presented above.  As the system is invariant under translations parallel to the walls, the den-  sity profiles, p (z),  depend only on the z-coordinate, whereas the correlation functions,  a  and c (r, 'z\, z ), are functions of three variables. The OZ equation, Eq.  h (r,Zi,Z2)  Q7  ay  2  (2.3.2), takes the form h (r,z ,z ) aj  1  =  2  c (r,z z ) ai  with r = \Jx + y and r 3  aj  2  3  z ) + 1 = exp[-/?u (i2) + h (r, z  u  with R = \jr  (2.6.1)  2  2S  h (r, z  +  2  = ^J(x — r) + y\. The H N C closure, Eq. (2.4.1), becomes  2  3  u  Q7  2  ai  u  z ) - c (r, Zi,z )] Q7  2  2  ,  (2.6.2)  + (z — Zi) , and the equation for the density profiles is given by combining 2  2  Eqs. (2.3.7) and (2.4.2), yielding p (zi) a  = a exp{ a  -8v (zi) a  - ^[c (0, z z ) + 1] QQ  " H Jdz p (z ) 27r  2  1  2  u  (2.6.3)  x  jdrr\^h {r, 2  ai  z z ) - c (r, z u  2  Q7  u  z )] } 2  7  We prefer this approach over the use of the B G Y or W L M B equations, as once the a  a  have been determined the system of equations can be solved for any wall separation, whereas the B G Y / W L M B route requires a slow step by step decrease of d i i to keep wa  the fluid between the walls in equilibrium with the bulk. In addition, we apply the  Chapter 2.  Equilibrium  electroneutrality  Statistical  18  Mechanics  condition to the slit system, which means that 2o + ^Qa a  fdzp (z) a  =0 .  (2.6.4)  J  Electroneutrality for the slit at small d n (i.e., a few ion diameters) can only be expected wa  if the walls are sufficiently thick, as has been shown by Lozada-Cassou and co-workers [36].  For P M electrolytes they found (employing a singlet integral equation theory) that  the required wall thickness ranges from around 40 A (2:2 electrolyte, p = 2 M [M represents moles/litre]) up to several hundred Angstroms (1:1 electrolyte, p = 0.01M) . Once the 4  walls are chosen to be effectively infinitely thick the electroneutrality condition follows from the fact that the electrostatic energy of an ion in a non-electroneutral slit tends to ±oo  as the slit becomes infinitely long . 5  If only a single ionic species is present between the walls, a  a  Eq.  follows directly from  (2.6.4), as the wall charges are exactly neutralized by an equivalent number of  counterions. This case corresponds to a system in equilibrium with a bulk at infinite dilution (salt concentration p—>0), where no coions are drawn into the cavity. For systems in equilibrium with a bulk electrolyte solution at finite concentration the procedure is slightly more involved, because of the fact that the individual ion activities are not independent.  Generalizing the ideas from Kjellander and Marcelja [39], we consider a z  z  strong electrolyte of the form M\ _\X ~, +  Z  cations of valency Z , and Z +  mean chemical potential,  +  Z  i.e., one unit of the salt dissociates into  anions of valency Z_.  We start from the definition for the  p±,  t*± =  1 Z + +  | _|(I -I^+ + Z  Z  '  (--) 2  6  5  A study for equivalent systems using density functional methods has been recently published by Henderson et al. [37]. A simulation study, done by Lee and Chan [38], where deviations from the slit electroneutrality have been achieved by adjusting the individual [i+ and u- is therefore misguided. 4  5  Chapter 2. Equilibrium  Statistical  Mechanics  19  which is satisfied by setting A -j± 3  p= +  p± + Z p + k Tln  where A =  +  B  and  p,_ = p± + Z_~p + k Tln  A  B  3  ,  (2.6.6)  is the mean thermal wavelength and p is a free constant  +  to be chosen to satisfy the electroneutrality condition. Combining Eqs. (2.2.3), (2.3.7), (2.3.8), (2.6.4), and (2.6.6) yields X +a z  with 7a(-2) — a /p (z) a  [-*$+ X -a„ J 7+(z)  1° = 0, exp{/fy }/A  z  +  a  J 7_(z)  (2.6.7)  3  ±  given by Eq. (2.6.3) and X = exp{(3p}. Eq. (2.6.7) is of order  (Z+ + \Z-\) in X and has to be solved together with Eqs. (2.6.1), (2.6.2), and (2.6.3). This way, the density profiles are finally obtained via 1  ""  e^  M±  "^)-w  --  (z) =xz  (2 6 8)  Details of the numerical algorithm used to solve this set of equations by iteration are given in Appendix A . Equation (2.3.5) for the bulk (v = 0) applied to Eq. (2.6.5) gives a  - ^ ^  =  7  ±  p  +  = 7 ±  ( j ^ i J  "-UrJ  •  (2  -'  6 9)  where the mean activity coefficient 7± is defined via ( z  7± =  +  + | z _ | )  /  | Z _ | Z  V7+ 7-  / o c i n \  +  ;  (2.6.10)  and expresses the idea that the deviation from nonideality is shared equally between the ions. 7-t for the bulk and slit electrolyte solutions must be adjusted in order to satisfy the condition that for large d n the midplane densities must equal the desired bulk wa  concentrations,  and p_. The advantage of writing the constant exp{/3p±}/A  3  P  +  as in  Eq. (2.6.9) lies in the possible use of bulk j± results for identical/similar models from simulations as a good starting point for this adjustment.  Chapter 2.  Equilibrium  Statistical  20  Mechanics  It is interesting to note that the exact location of the surface charge planes does not enter the problem, as long as the accessible space for the fluid particles, D , a  remains  unchanged. For example, relocating the charged planes further apart (while keeping D  a  constant) changes the constant d n —» d i in the expression for u® [see Eq. (A. 1.14)], 1  wa  wal  i.e., 7_|_ and 7_ modify their values accordingly. However, the resulting density profiles can be shown to be invariant under such a transformation due to a corresponding adjustment of X in Eq. (2.6.7) (i.e., because of the imposition of strict electroneutrality in the slit). A related observation is that because of Eq. (2.6.6), the chemical potentials of the individual ion species, for systems where an infinite wall is present, are different from the corresponding bulk values (even for a symmetric electrolyte). This is due to the formation of an infinite planar dipolar sheet at the surface (see Fig. 2.3) which gives an energetic contribution independent of the distance from the walls . The effect disappears when 6  spherical macroparticles are considered. We stress again that p, and p- never have to +  be specified explicitly, as their values follow from fixing p± and from the electroneutrality of the slit system. The accuracy of the A H N C theory can be expected to be good to excellent whenever core-core interactions are not significant [28,41,42]. Deviations occur when counterions are frequently in contact with each other (as for monovalent ions at high density [25]) or when coion-counterion pairing becomes important (for example, the low cr/low p case, which is not considered here). With the system parameters typically considered in this thesis the contact values of c/ _ are overestimated by the A H N C theory [43], as the +  corresponding bridge functions (ignored in the H N C closure) give a short-ranged repulsive contribution in Eq. (2.3.3). We note that an alleged failure of the A H N C scheme reported by Lozada-Cassou et al. [36] is due to a misinterpretation of simulation results. The chemical potential of the counterions close to a charged wall for a P M electrolyte solution has been studied by Svensson and Woodward [40]. 6  Chapter 2. Equilibrium  Statistical Mechanics  21  Figure 2.3: Comparison of the electrostatic energy of a charge q far away from a charged macroparticle with a screening double layer around it for two macroparticle geometries. In (a) the macroparticle is an infinite planar wall with a surface charge density —a. The charge distribution of the double layer is centered at a distance A from the surface and cancels the wall charge exactly. The Coulombic contribution to the energy of a charge q far away (in the bulk) is independent of the distance from the wall and given by v = 0crA/(2ee ). The corresponding case for a spherical macroparticle in b) yields v = 0. Here, the total surface charge — Q is the negative of the charge associated with the double layer (but the surface charge densities of both charged spheres in (b) are different). el  o  el  Chapter 2.  2.7  Equilibrium  Statistical  Mechanics  22  Pressure Calculation  A quantity that is of major interest is the average net force acting between macroparticles. For example, in colloidal systems this force determines whether the colloidal particles stay dispersed in the solution or aggregate. In the case of two planar walls we are interested in the force per unit surface area, i.e, the pressure, P. Note that unlike the bulk system case, here the pressure is not isotropic but is a tensor quantity having two components for the slit geometry: the pressure perpendicular to the walls, which represents the forces on the walls (and will be considered the pressure from now on) and a parallel pressure component . For the former component two contributions must be considered. One is 7  due to the attractive short-ranged van der Waals interaction which decays as d~l and will n  not be included in our discussions. The other, here called P n , arises from the interaction s  t  of the electrical double layers as the surfaces are brought together. Henderson et al. have suggested [45] and proven [46] the contact theorem for the pressure of a fluid in contact with a planar wall (I), starting from an equation of the form Psiit = - E / ° °  dzp (z)^^-  ,  a  (2.7.1)  i.e., the derivative of the average interaction energy of the fluid with the wall. Inserting () = a ()  vI  z  a  v  S,I  z  + t { ) + a''{ ) v  J  z  vS  z  i n t o  ^  E c  l - (2.7.1) and using Eq. (2.6.4) one obtains  (J  Pm = ^ T ] > > ( - z a  m a x  - .  2  , ) ~Q  £  0  dV (z) Sh,I  /-2max,Q  /  dzp {z) a  a  {  ]  .  (2.7.2)  The first term on the right-hand side (RHS) depends on the particle density at the plane of contact, whereas the second term represents the electrostatic interaction between a wall and the remainder of the system (i.e., the fluid and the other wall). Carnie and Chan [47] arrived at identical results employing a different method and remaining questions have The normal pressure component, P , is constant throughout the system, whereas the parallel component, pH (z), is not uniquely determined and depends on the position between the walls. For limd _^ one has, of course, PH(0) = P . and the quantity J cte[P- - P"(z)] is the surface tension [44]. 7  x  wall  x  2max  0  L  0O  Chapter  2. Equilibrium  Statistical  23  Mechanics  been cleared up by Henderson and Blum [48]. As the normal pressure can be evaluated at every plane parallel to the walls, Kjellander and Marcelja [35,39] have suggested that instead of direct application of Eq. (2.7.2), it is desirable to transform to an equation for the midplane with the help of the B G Y equation. In this way P ij is split into a s  t  term depending on the midplane density, Pkj , and contributions representing the various n  particle interactions [see Eq. (2.5.1)] across the midplane.  One finds (see Appendix B  for details) -fslit  -Pkin + Pcore  =  + P +P el  Pki„ =  +P  sh  , with  wall  (2.7.3)  A: T5> (0), Q  B  a  -Pcore  =  /c T2_,27r /  dzip (zi)  B  x a,7  fO  fZma.x,a  k  a  /  dz p (z ) 2  2  Q7  Q7  2  5<'(r,2i,Z )  x / 2-ivrdr—2XL 22 dzip (zi) Jo azi  2  ~E  /"'"""'^lPo^i) /  x /  2nrdr  a  i  v  hdz(r,z p (z, ) al  a  =  2  (zi — z ) <7 (r , zi, z ) , /•oo  Fsh  1  2  7  x  z) , 2  2  dz {z ) 2Pl  \  2  g (r,z z ) ai  u  , and  2  wall  The integrations over Z ! and z for P 2  (zi — z )  2  2  are to be carried out only if r  2  c o r e  = ^(d + d ) — 2  tJ  a  7  is positive, i.e., if two particles at z\ and z can be in contact. Note that 2  the calculation of P i involves h , and not g , as the second term on the R H S of Eq. e  ay  aj  (2.7.2) is absorbed into P i ; only deviations from the mean field case (i.e., correlated e  electrostatic fluctuations) contribute and in general one has P i < 0 . The use of h 8  e  ai  The pressure component P \ can be predicted quite accurately by assuming quasi-2D surface layers of counterions and considering the opposite effects of the Coulombic repulsion and the thermal motion 8  e  Chapter 2. Equilibrium  instead of g  ai  Statistical  24  Mechanics  in an equivalent expression for P .given by Eq. (3) of Ref. [50] is due to a sh  typographical error [51]. Both Eqs. (2.7.2) and (2.7.3) are exact but can yield different results when applied in numerical calculations because the correlation functions are only known approximately [32], which will be further investigated in Section 4.4.2 . A third route to the pressure has been suggested by Kjellander and Marcelja [34,35], exploiting the thermodynamic relation Eq. (2.2.6). Specifically, one calculates the appropriate total thermodynamic potential of the system between the walls, here the grand potential per unit surface area, fl , and obtains the pressure via A  ^=-(IR \Cawall/  • {  /  1  }  )  <™>  T  The H N C approximation yields an explicit formula for Q , which is given in Appendix B . A  The net pressure, P , is given by the difference between the internal (slit) pressure net  and the outside (bulk) pressure such that -Pnet = -fslit -fbulk •  (2.7.5)  —  For large wall separations no net force acts between the walls and one has P  bulk  =  lim P  (2.7.6)  slit)  "wall-* OO  which is exploited in our calculations to obtain Pbuik- Wall separations yielding P  net  =0  represent a fluid-wall (macroparticle) system in equilibrium (when no additional outside force is applied), although the equilibrium position is stable only if dP t/<9d ii < 0. ne  of the ions in a simple model [49].  wa  Chapter 3 Mixtures of Neutral and Charged Hard-Spheres - Packing Effects  3.1  Introduction  While it is now possible [52] to examine bulk electrolyte solutions using quite realistic solvent models and reasonably accurate theories, the same level of treatment is not yet feasible for inhomogeneous systems. Therefore, investigations of two interacting double layers have usually employed the primitive model, where the ions are represented by charged hard spheres and the solvent is a dielectric continuum without any inherent granularity. We would expect discrete solvent effects to alter the P M results through at least two mechanisms. These are, particle packing constraints, and the relative lack of dielectric screening at small ion-ion and ion-wall separations. The latter effect will be included in the description of the electrolyte solution by employing McMillan-Mayer level theory with effective ion-ion and ion-wall interactions in the following chapter. In the present chapter we focus upon solvent granularity and its effects on the ion distributions and on the force between the immersed charged walls. It has been proposed [53,54] that the net pressure between walls in an ion-solvent mixture might be given, at least approximately, by the superposition of separate ion and solvent contributions obtained from simpler one-component models. A n important objective of the present work is to examine the validity of this appealing suggestion. Some related earlier work has been reported by Davis and coworkers [55-57] and by Patra and Ghosh [58]. These authors considered ions in hard-sphere solvents using  25  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing  Effects  26  the methods of density functional theory. They found that at high density the pressure between walls exhibited an oscillatory structure due to the hard-core interactions, which parallels the behaviour for two walls [33,59] and for two spherical macroparticles [60] immersed in pure hard-sphere (HS) fluids. Here, we consider similar models but employ the A H N C technique. However, our main purpose is not to compare results obtained by different theoretical methods, but rather to complement the earlier findings with a detailed analysis of the fluid structure, of how ions and neutral particles mutually influence each other in the mixture, and, particularly, of the net pressure in terms of its component parts. Cases where the ions and solvent particles are of equal and unequal size are discussed.  3.2  M o d e l and Remarks on the Numerical Method  The interaction potentials for the systems considered in this chapter are given by the HS and Coulombic parts of Eqs. (2.5.1) and (2.5.2), u (R)  =  u™(R) + u ^(R)  v (z)  =  v**(z) + v*(z),  ai  a  e  and  (3.2.1) (3.2.2)  i.e., by setting u^^R) = v^(z) = 0. The fluid between the walls is a binary mixture of hard sphere counterions and neutral hard sphere solvent particles, labeled a, 7 = ion and hs, respectively. The effects of coions are therefore totally ignored here and the model represents a system in equilibrium with a bulk solution with a HS particle density p  h s  and the salt at infinite dilution. The advantage of this procedure lies in the restriction to a two-component mixture which is much easier to treat than the full three-component model. While this simplified model clearly involves some level of physical approximation, it provides a reasonable description of the situation for highly charged walls immersed in dilute electrolyte solution [42]. As the short-range ion-ion and ion-wall interactions are  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing  Effects  27  zero for the present model, the contact theorem and midplane formula, Eq. (2.7.3), for the slit pressure simplify to Psm. =  kT B  [p (z hs  max,hs ) + Pion( max,ion)] 2  •Pslit  =  -Pkin + -Pel + -Pcore  ~~ ~  >  ,  a n o  -  (3.2.3) (3.2.4)  respectively. The accuracy of the A H N C theory for confined ionic fluids is well documented [41,42]. On the other hand, it is known that the H N C approximation does not perform comparatively well for dense homogeneous liquids with short-range pair potentials, and this holds true for the inhomogeneous case as well. For a dense HS fluid between two walls the anisotropic P Y theory yields excellent results [33,59], but the P Y closure is not appropriate when long-range forces play an important role. Henderson and Pljschke compared the contact densities (equivalent to Pbuik/ksT)  for HS fluids in contact with one wall  for the anisotropic and singlet versions of the P Y and H N C approximations with exact results [31]. Surprisingly, they found the worst performance for the A H N C theory, with even the singlet H N C approximation yielding more accurate results. In contrast, our calculations imply that the A H N C theory works better than both singlet theories. The deviations of our results from those in [31] might be due to lower numerical accuracy of their calculations (e.g. smaller number of grid points). Remarks in this direction have been made by Henderson and Plischke themselves and by Kjellander and Sarman concerning the corresponding accuracy of the anisotropic P Y curves [33] (although for the latter case the deviations are less severe). For a density of p  h s  ~ O^/df^, as used in  our study, the deviation of the A H N C P ik from the exact value is « 9 % , where we use Du  wall separations of ~ 9 dh or larger to determine Pbuik with sufficient accuracy. OverS  all, limitations in the treatment of the hard-core interactions do not seriously influence our conclusions as much of the discussion is qualitative in nature, and all quantitative  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  28  comparisons are made between systems involving the same level of approximation.  3.3  Results  In the following discussion, we refer to results for models of four types. These include the primitive model where only counterions are present between the walls, hard-sphere systems where only hard spheres occupy the slit between uncharged walls, and mixtures of counterions and hard spheres (labeled PM+HS) of particular interest here. In addition, we consider systems which are equivalent to P M + H S models except that all charges are switched off while keeping the number of "ions" between the walls fixed . We label 1  systems of this type U C P M + H S (UC for "uncharged") and they prove very useful in our attempt to obtain and understand different superposition approximations for the force between the walls. In all calculations, the relevant state and interaction parameters were chosen to be consistent with d = 2.8 A, e = 78.7 and a temperature T = 298 K. Thus, the hard-sphere n s  solvent particles are roughly the size of water molecules and the background dielectric continuum has the dielectric constant of pure water at 298 K . We consider mainly divalent counterions ( Z ; = +2) with diameters dj = dh and d; = 4.25 A —1.52 d on  0n  s  on  2 h s  . The latter  value is frequently used in P M calculations and is sometimes considered to represent an "effective" ion diameter which includes some portion of a strongly bound solvation shell. Here we have selected this value simply because it allows us to check our P M results against earlier work, and it suffices to show the large effects which occur when the ions and solvent particles differ significantly in size. The restriction to divalent counterions is This can be achieved by scaling down all charges (a and gi ) by a factor 1000, while still applying the electroneutrality condition. The reduced temperatures, T* = 4ne edi T/q , for the systems investigated are 0.1 {Z\ = 2, 4> = 2.8 A), 0.15 (Z = 2, d =4.25 A), and 0.39 (Z = l, d = 2.8 A), but note that the significance of T* for the charged slit system is changed in comparison to the corresponding bulk. 1  0n  2  2  0  n  ion  ion  on  on  ion  ion  Chapter 3. Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  29  basically due to the more interesting behaviour of P n for the corresponding P M (a range s  t  of attractive wall-wall interaction). A few results concerning monovalent counterions are also presented. The results reported are for phs = hm<2 _>oo Phs(O) =0.492/d . To achieve this value, wal  we used a bulk activity coefficient 7  ns  ns  = 57.1. This activity co efficient ^ differs from the  corresponding "exact" (i.e., Carnahan-Starling [61]) value of 40.9 due to approximations inherent in the H N C closure. We note that our bulk density is a little lower than that of water under ambient conditions (i.e., ~ 0.7/o! ). However, the numerical solutions are ns  dramatically easier at the lower density and the system considered is sufficiently dense to illustrate the large influence of solvent granularity which is of primary interest here. Furthermore, fluids of hard spheres are considerably more structured than water-like models at the same density [62], so using a somewhat lower density offsets the "over structuring" to some extent. The presentation of the results is divided into three parts, and we first consider the simpler case where the solvent particles and divalent ions are of the same size.  3.3.1  Particles of Equal Size  In Fig. 3.1 we give two examples of the density profiles of hard spheres and divalent ions in the mixture (PM+HS) compared with the profiles for divalent ions (PM) and hard spheres (HS) alone. Note that p (z) a  is symmetrical around the midplane and only one  half of the profiles are shown. The most obvious feature is the tendency of the ions to be more likely found close to the wall when the dense hard-sphere fluid is present, due to entropic effects as the free space between the walls is reduced. The presence of the ions at the wall lowers the contact density of the neutral component in the mixture in comparison to the pure hard-sphere system, whereas the densities around the midplane are similar in magnitude.  The profile for d u = 3.1 d wa  h s  plotted in Fig. 3.1(b) shows  Chapter 3. Mixtures  of Neutral and Charged Hard-Spheres  -1.0  -0.8  -0.6  - Packing  -0.4  z/ d  30  Effects  0.0  -0.2  h s  Figure 3.1: Density profiles for the mixture of ions and hard spheres (PM+HS), the primitive model (PM), and the pure hard-sphere model (HS). In all cases d = d and i o n  Z  ion  — +2. In (a) d aii = 2.6dhs and in (b) d n = 3.1 d W  wa  hs  h s  Chapter 3. Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  31  the appearance of maxima and minima typical of dense hard-core systems indicating the build up of particle layers. It should be noted that with the moderate hard-sphere bulk density used here we do not find any significant layering of the ions in contrast to the situation for p = 0.7/d  3  hs  hs  [55,56].  A n interesting point is that the average total particle density between the walls,  p, tot  (calculated by integrating the density profiles of all particles for a given model) is higher for P M + H S then for HS. This fact, which is true for all wall separations, is not due to the structural changes because of the charges present. A comparison between the charged and uncharged mixtures reveals that p [PM tot  + HS] and p ot[UCPM + HS] are t  always within 0.5% of each other, whereas p t[HS] is up to 5% lower than the values to  for the binary mixtures. This indicates that the U C P M + H S model is not equivalent to HS even though all the species have identical interactions. The difference in p  tot  for the  pure and the mixed systems is rooted in the distinction between particles with a fixed chemical potential and particles that are located essentially by definition (due to the electroneutrality condition) between the walls. The latter species may be viewed as part of the wall-wall system and their chemical potential is generally different from the former species even if their interactions are identical. A more complete discussion of the trends in the density profiles is possible if we look at the contact plane and midplane densities as functions of d u for the different models. wa  The contact densities (Fig. 3.2) for all components in the dense systems (all except the PM) show to some degree maxima (around d ji = 2.05 dh , 3.1 dh ) and minima (around S  wa  dwaii = 1.6dh j2.6dhs)S  S  These are related to fluid structures that are more "efficiently  packed", as in Fig. 3.1(b), or more "loosely packed", as in Fig. 3.1(a), respectively. At wall separations slightly larger than a multiple of d s neighboring particle layers interact n  more strongly and particles in the contact layer are pushed towards the wall [59]. As hard-core interactions in the P M are of no importance the corresponding curve is not  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  32  2.0  1.0 cox: T5  (a)  c  0.0 a  N Q  PM+HS HS - UCPM+HS ••••  h  J  I  J  L  I  I  J  L  I  I  L  J  L  L  J  L  PM+HS  0.8  PM UCPM+HS  0.4  0.0  - (b) J  I  I  1  i  J  I  I  2  J  L  3 dwall /  ^hs  Figure 3.2: Contact densities for (a) hard-sphere and (b) ion components for different models. The model labeled U C P M + H S is identical to P M + H S but the charges are switched off as described in the text. The remaining models are as in Fig. 3.1. In all cases djon = d  n s  and Z-  l0Xl  = +2.  Chapter 3.  Mixtures  of Neutral  and Charged Hard-Spheres  - Packing  Effects  33  structured, but has one broad minimum at d n « 2.15 d s due to correlated ion-ion wa  n  density fluctuations. The midplane densities (Fig. 3.3) for the dense systems exhibit layering features similar to those of the contact values. For d ii ^ 2.2 d s) the densities at the midplane and contact plane for P M + H S and wa  n  U C P M + H S are very similar. This demonstrates that charge effects at these small wall separations are only of secondary importance compared with packing constraints. Features such as higher ion contact densities and a faster decrease of the ion midplane density with increasing d n for the P M + H S (compared to the P M ) come about because of the wa  hard-core interactions. For larger wall separations the curves for the charged and uncharged mixtures deviate significantly. The "ions" in the U C P M + H S spread more evenly over the accessible space and their contact density falls slowly to zero as d n increases, wa  whereas a constant nonzero contact density is reached for the P M + H S [see Fig. 3.2(b)]. The solution of the integral equations also yields the particle distribution functions 9a-y( , zi, 2)r  z  These functions depend on three coordinates and are therefore hard to  depict, especially if small differences between models are to be shown. Here (Fig. 3.4), we present only a plot of #i nion( ) parallel to the wall at the contact plane (z\ = z = z ) r  2  0  and midplane (zi = z = 0) for d n = 2.05 d . The differences between the c/j 2  wa  ns  0nion  majX  (r) for  the P M and the mixture are rather small (i.e., < 0.1). This means that although the average positioning of the ions between the walls (as indicated by the density profiles) changes considerably upon insertion of the hard spheres, the ion-ion structure in the fluid is not disturbed very much. The presence of the hard-sphere solvent induces small undulations with a "wavelength" slightly smaller than the solvent diameter, similar to the effects produced by a more realistic solvent model in bulk solution [52]. For the P M midplane function, the single maximum at r = 2.7d s> which constitutes the global n  maximum as well (as seen in a full 2D contour plot), is split into two maxima in the mixture. The oscillations in <7inion(?") in the contact plane are smaller and no maximum 0  Chapter  3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing  Effects  34  Figure 3.3: Densities at the midplane for (a) hard-sphere and (b) ion components for different models. The models are as in Fig. 3.2. In all cases d- = d„s and Z = +2. wn  i o n  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  35  Figure 3.4: # i n i o n ( 0 parallel to the walls in the midplane and contact plane for the P M and P M + H S (o?j = d s) with d = 2.05 d s and Z ; = +2. O  0n  n  w a n  n  on  Chapter 3. Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  36  exists; the global maximum (and the most probable spot for the next ion) can be found at the opposite wall. We now discuss the pressure perpendicular to the walls.  Because of the contact  theorem, Eq. (3.2.3), most of the points raised in our discussion of the behaviour of the contact densities hold for the pressure as well. Fig. 3.5(a) shows the oscillating pressure for the pure hard sphere system together with the smooth, slightly attractive curve for the P M and the electrostatic component, P i , which is always attractive as it e  arises because of correlated ion-ion fluctuations across the midplane. The net pressure for the system of interest, P M + H S , together with results given by two superposition approximations, is plotted in Fig. 3.5(b). If one wishes to estimate P  n e t  for P M + H S  with simple addition schemes involving more basic components, different possibilities are available. The simplest suggestion is to add the values of P  net  for the pure systems (i.e.  P t[HS] + P [ P M ] ) . From Fig. 3.5(b), we see that for the equal-sized case this gives ne  net  good agreement down to d aii~l-8dhsW  For small wall separations another procedure yields much better values, specifically, the superposition of P  n e t  [ U C P M + HS] and P i[PM]. The idea leading to this choice e  is that for small wall separations packing effects are crucial for the fluid structure and should be separated from the electrostatic contribution. P i on the other hand shows e  rather small changes when comparing ion models with and without hard spheres (as is evident from Fig. 3.6) despite significant differences in the density profiles. However, this superposition scheme breaks down at larger wall separations as the "ions" in U C P M + H S do not stay close to the walls (as do those in PM+HS) and the pressure is systematically overestimated. A couple of remarks are appropriate here. First, the comparisons made in the preceding paragraphs are of rather academic interest as we do not suggest that any procedure  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing  37  Effects  Figure 3.5: The net pressure acting between the walls obtained for the different models and with different superposition approximations. In all cases d- = d and Z = +2. The labels are as in Fig. 3.2 and as discussed in the text. lon  n s  ion  Chapter 3.  Mixtures  I  of Neutral and Charged Hard-Spheres  i  i  0  i i I  i  i  1 (d  2  3  ii dion) /  38  Effects  i i I i i i i L_i i —  w a  - Packing  i  i  I  4  dhs  Figure 3.6: The electrostatic component of the pressure for the P M and P M + H S with ^ion = +2. Results for d i = <4s and for d 0n  ion  = 1.52d are shown. hs  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  39  of adding pressure components of different models will yield exactly P [ P M + H S ] . Furnet  ther, because of the restriction to equal-sized particles, most schemes of combining the oscillating pressure found for the HS fluid with a slight attraction of electrostatic origin will likely give a fair qualitative picture. The case where the ions and solvent particles are of unequal size considered in the following subsection is much more interesting. Secondly, note that P [HS] and P net  net  [UCPM + HS] are not the same although the "ions" and hard  spheres are identical in U C P M + H S . This is due to the above mentioned peculiarity of the present model, where the bulk is treated as a pure fluid with the mixture existing only in the cavity between the walls.  3.3.2  Larger Ions  A set of calculations analogous to those described above have been performed for systems with identical diameters for the neutral HS component, but with ions which are 50% larger, or, more precisely, dj  0n  = 4.25 A = 1.52 d . With this choice, the ions are hs  considerably larger than the solvent particles and deviations from the equal-sized case are expected to be significant. Moreover, this value of d;  on  is often used in studies of the  P M and results for the pure ion system have been published [42]. We note that for both ion sizes considered here the P M yields essentially identical results (with the appropriate redefinition of d n); the divalent counterions stay sufficiently apart from each other that wa  ion-ion hard-core interactions are not important. See, for example, the P i plots given in e  Fig. 3.6. For larger values of dj  on  and especially for monovalent ions [42] density profiles  and pressure curves depend more strongly on the ion diameter. A qualitative comparison of Figs. 3.7 and 3.1 indicates that with the different particle sizes the density profiles differ more significantly from those of the pure systems. This comes about because an "in phase" layering of ions and hard-sphere solvent particles is no longer possible. The incommensurate diameters appear to bring on an even more  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing  Effects  40  Figure 3.7: Density profiles for the mixture of ions and hard spheres (PM+HS), the P M , and the pure HS model. In all cases dj = 1.52 d and Z — +2. In (a) d n = 2.5 <4 and in (b) d i i = 2.9tZhson  wa  n s  ion  wa  S  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing  41  Effects  pronounced tendency to push the ions towards the walls. The plot of the contact densities given in Fig. 3.8 confirms this observation; the ion contact densities are higher and those of the neutral species lower compared to the mixture with particles of the same size (Fig. 3.2). Another point of interest is that the first maximum of the ion contact density in the mixture (PM+HS) is at d n ~ 2 . 5 dh = dhs + ^ion, which means that ions are pushed wa  S  closer to one wall by hard spheres associated with the contact layer at the opposite wall. The pure HS fluid has accordingly its first maximum slightly higher than d  wall  = 2 d [Fig. hs  3.8(a)]. Since hard spheres in P M + H S and U C P M + H S feel the interactions with both species at the opposite wall we find two maxima at the corresponding wall separations for these mixtures. This more complex behaviour can also be found for the pressure between the walls shown in Fig. 3.9. The distances between maxima and minima in P  net  for P M + H S  are irregular and no correlation with the oscillations characteristic of the HS model is apparent. The curves shown in Fig. 3.9(b) also demonstrate that the simple addition approximation P  n e t  [ P M + HS] ?a P  n e t  [HS] + P [ P M ] is not good at any wall separation if net  the particle diameters are sufficiently different. On the other hand, the more complicated scheme, P of  c/waii •  n e t  [ P M + HS] ^ P  n e t  [ U C P M + HS] + P i[PM], works very well for smaller values e  Deviations do occur at larger wall-wall separations where the fluid structure in  U C P M + H S differs significantly from that of P M + H S . Nevertheless, these results show that to a good approximation the net pressure in the P M + H S can be viewed as a superposition of an attractive electrostatic part and a "quasi oscillatory" component arising from the hard-core interactions in the system.  Figure 3.8: Contact densities at the wall for (a) hard-sphere and (b) ion components for different models. The models are as in Fig. 3.2. In all cases dj = 1.52d and zT = +2. on  hs  ion  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  43  Figure 3.9: The net pressure acting between the walls obtained for the different models and with different superposition approximations. In all cases d; = 1.52 d s and Z = +2. The labels are as in Fig. 3.2 and as discussed in the text. on  h  i o n  Chapter 3.  3.3.3  Mixtures  of Neutral  and Charged Hard-Spheres  - Packing  Effects  44  Monovalent Counterions  If divalent ions are replaced by monovalent ions, while keeping all other model parameters fixed, the average density of ions in the slit doubles. Moreover, monovalent ions repel each other less strongly and the overall effect is that ions are much more likely to be in contact, compared to the divalent case. However, the density profiles (Fig. 3.10) qualitatively resemble those for the models with divalent ions (Fig. 3.1), although a slight tendency of ion "layering" around the midplane can be detected for the P M + H S model at d n = 3.0d wa  ns  [Fig- 3.10(b)]. Also, the conclusions reached for the slit pressure  and the appropriate approximation schemes still hold. It is more interesting to consider the ion-ion pair distribution functions. In Figs. 3.11 and 3.12 we show results for ions located at the midplane (zi = 0) and, in contact with one wall  (z\  = —z  ma  x,ion)  f°  r a  wall separation of  d ii wa  = 2.6d - The white region around ns  the particle is forbidden to other ions because of the hard-core interactions. The P M contour plots show that beyond R=^r  + (z — Zi) = d\  2  2  on  the likelihood of encountering  another ion is reduced through the electrostatic repulsion, c/i ; (i?) < 1 is monotonic 0n  on  in all directions, and only a small degree of anisotropy is introduced by the walls. In contrast, the corresponding c/  ionion  (i?) for the P M + H S model show oscillatory behaviour,  with stronger anisotropy and, in particular, regions where the contact value is larger than 1. The last observation indicates that the monovalent ions in the dense HS fluid experience an effective attraction when close to each other, as has been pointed out for bulk systems by Rescic et al. [63]. The distribution function parallel to the walls for  d ii wa  = 2.0dh shown in Fig. 3.13 demonstrates the same effective attraction for the S  midplane, whereas along the walls <7i nion(0 remains <1 (similar to the P M curves). O  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  -1.0  -0.8  -0.6 z/d  - Packing Effects  -0.4  -0.2  45  0.0  h s  Figure 3.10: Density profiles for the mixture of ions and hard spheres (PM+HS), the primitive model (PM), and the pure hard-sphere model (HS). In all cases d = d and Zion = + l . ( ) d i i = 2.6d s and in (b) d = 3.0 d i o n  I n  a  wa  h  w a l l  hs  h s  Chapter 3. Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  46  g(r,z) for counterion/counterion PM 1.2000  • • •  1.1000 -1.0-  1.0100 -6.0  -4.0  -2.0  0.0  +2.0  +4.0  +6.0  0.9900  r  I1  0.9000  g(r,z) for counterion/counterion PM+HS  0.8000  • • •  0.7000 0.6000 0.0000  •  Figure 3.11: Contour plots of g\ \ {r, 0, z) for the P M and the P M + H S model with dion = dhs and Z j = +1 for d n = 2.6 d - r and z are in units of d^ . The black horizontal bars represent the walls and the legend on the right explains the colour scale, with the number between two squares indicating the boundary value. on on  o n  wa  ns  s  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing  Effects  47  r Figure 3.12: Contour plots of # n i o n ( r , -z , max  0  4>n = 4s and Z  ion  = +l for d =2.6d . waa  ha  z) for the P M and the P M + H S model with  r and z are in units of 4s-  Chapter 3. Mixtures  of Neutral and Charged Hard-Spheres  1.6 \-  - Packing  Effects  48  dw«u=20 dhs  1.2 5G o  0.8  midplane  s  0.4 \0.0  contact plane  1  2  PM PM+HS PM PM+HS 3  r/d  h s  Figure 3.13: #i nion(?") parallel to the walls in the midplane and contact plane for the P M and P M + H S (dj = 4s) with d i i = 2.0dhs and Z\ = +l. 0  on  w a  on  Chapter 3.  3.4  Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  49  Summary  In this chapter, we have solved the A H N C approximation for mixtures of neutral hard spheres and counterions between charged hard walls. In agreement with earlier calculations [56,58], we observe that adding a neutral solvent to the P M modifies the ion density profiles considerably. The hard-core ion-solvent interactions tend to push the ions towards the wall resulting in higher ion densities at contact. This effect is amplified if the ions are larger than the solvent particles. The results for the ion-ion pair distribution functions, <7ionion> are qualitatively different for the models with monovalent and divalent ions, (/ionion for divalent ions is only slightly influenced by the presence of the dense HS fluid, whereas monovalent ions are less resistant to being pushed together by the dense component.  The corresponding curves exhibit stronger oscillations and regions where  <7ionion > 1 close to ion-ion contact. We find that, even at moderate solvent density and high surface charge density, the net pressure between the walls at separations of a few solvent diameters tends to be largely dominated by oscillations associated with the hard-core interactions. At these separations, electrostatic effects are only of secondary importance in these models. For the restrictive case where the ions and solvent particles are of the same size, the oscillations are regular and similar to those for a pure hard sphere fluid. For this system, simply adding the P M and HS pressures gives a reasonable approximation to the forces found for the mixture. If the ions and solvent particles are significantly different in size, one obtains a more complex pressure curve which cannot be approximated by simply adding the P M and hard sphere pressures. However, for wall separations up to several solvent diameters the net pressure can still be well described by a better superposition approximation. This consists of adding the net pressure for a corresponding uncharged hard-sphere mixture to the purely electrostatic component, P i , of the P M . This scheme e  Chapter 3.  Mixtures  of Neutral and Charged Hard-Spheres  - Packing Effects  50  works because the hard-core interactions remain the most important contribution and the P i component in the mixture differs little from that of the P M . e  Chapter 4  M o l e c u l a r Solvent Effects at the M c M i l l a n - M a y e r L e v e l  4.1  Introduction  Incorporating sophisticated models of solvent molecules into calculations of the interaction between macroparticles in electrolyte solution at a reliable level of theory is a task at (and sometimes beyond) the edge of present computational possibilities. Although direct simulations of selected systems have recently become possible [64], more extensive investigations for a range of systems and state conditions still require model simplification and approximate numerical methods. For example, in a recent paper Trokhymchuk et al. [65] examined the interaction of a pair of large spherical macroions immersed in a solution of point ions and dipolar hard spheres. The model is solved at the level of the mean spherical approximation with the homogeneous H N C equation used to calculate the macroion-macroion potential of mean force. A n important step was taken by Kinoshita et al. [66] who examined potentials of mean force between a pair of spherical macroparticles in a mixture of monovalent hard-sphere ions and water-like solvent particles. They studied in detail the dependence of the potential of mean force on macroparticle charge, on electrolyte concentration, and on ion size. However, the accuracy of the homogeneous reference hypernetted-chain (RHNC) approximation used in their calculations is difficult to judge for these systems; previous experience with the P M suggests that H N C results for potentials of mean force between macroions can be inaccurate [67], and the homogeneous integral equation approach must be treated with caution in the absence of other  51  Chapter 4.  Molecular Solvent Effects at the McMillan-Mayer  Level  52  evidence. In the present chapter, we compare our findings with the results of Kinoshita et al. [66] wherever possible. In a different approach, the complexity of the model electrolyte solution might be reduced at the outset by averaging over all solvent coordinates to obtain ion-ion potentials of mean force. These "solvent-averaged" potentials of mean force can then be used as effective ion-ion potentials in calculations of the force between walls. The theoretical calculations are of course greatly simplified because the solvent particles are not treated explicitly and exert their influence only indirectly through the ion-ion potentials. For bulk solutions, this scheme was first outlined by McMillan and Mayer [68] and we will refer to it as the McMillan-Mayer (MM) level of description. A summary of the important ideas involved is given in the following section. We note that Marcelja [69] has also recently considered some related formal aspects of M M theory for surface force problems. A t the M M level of description the solution of the A H N C theory becomes computationally possible because the ionic fluids are at low density and the effective interaction potentials are spherically symmetric, if the bulk functions are used. Some results obtained in this way have been given by Marcelja [50,70] for a model NaCl solution. Unfortunately, due to a calculation error [51], the pressures reported by Marcelja are incorrect and cannot be used for comparison. In the present chapter we employ hard-core models for both the ions and the water-like solvent particles. These models were previously developed for bulk aqueous electrolyte solutions [52] and are a natural extension of the P M approach.  4.2  Pairwise Additive McMillan-Mayer Theory  We now consider the case of two slit systems in equilibrium with the same bulk solution consisting of a dense solvent and the salt at infinite dilution. The first (or "reference") system, indicated by a superscript "*", has no charges on the walls (cr* = 0) and therefore  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer  Level  53  the fluid in the slit consists only of solvent molecules. The second system has a finite a, which implies that the fluid between the walls is a mixture of solvent particles and counterions (due to the restriction to an infinitely diluted bulk solution). In Appendix D, we derive the grand canonical M M formula for the ratio of the partition functions at two different sets of activities, S(a , a )/S(a*, a*), where the indices <; and s stand for the ?  s  solvent and a single solute (here ion) species, respectively. The solvent activities for both slit systems are identical to the bulk value, i.e., a^ — a*. The ion activity in the reference system is zero, a* = 0, whereas for the finite a case a ^0  for the counterions because of  s  the surface potential due to the charge distribution of the double layers (see Fig. 2.3). One finds ^  =  = J  / d R i - d i L exp [-BWn({&h  a*, 0)] ,  where n is the number of counterions fixed by the wall charges, and W ({R}; n  (4.2.1) a , a = 0) f  s  represents the "effective" potential energy of n counterions (including "effective" ion-wall potentials) at infinite dilution in the solvent at activity a [compare with Eq. (D.16)]. If q  one substitutes the n-particle potential energy (V„ + U ) for W , the RHS of Eq. (4.2.1) n  n  is formally equivalent to the partition function of a one-component system. Equation (4.2.1) therefore shows that the ratio E/E* can be calculated by only considering the ions interacting via these effective potentials determined in the (here inhomogeneous) reference system. This is the main point made by M M theory. Moreover, as the pressure is proportional to ln E, the total pressure of the slit with a finite a can be written as a sum of the contributions from the pure solvent slit system and the pressure of the charged slit evaluated at the M M level. The case with a finite bulk salt concentration is slightly more complicated but the main arguments remain the same. We note that Eq. (4.2.1) is an exact description of the mixture if many-body potentials of mean force are included. However, in practice the effective interactions are usually  Chapter  4.  Molecular Solvent Effects at the McMillan-Mayer  Level  54  truncated at the pairwise level and W is set to n  = it  W {R R\) n  lf  (Ki)  + E E  W  i—l  i=l  (Ri> Rj) >  w  .  (4-2.2)  j>i  where W(R) is the (one-particle) potential of mean force from Eq. (2.3.7), which acts as an effective ion-wall potential, and w(Ri, Rj) is the ion-ion pair potential of mean force. As a further approximation, we do not employ the w(Ri,Rj)  obtained from the  inhomogeneous solvent system but the corresponding bulk functions. These depend only on the ion separation, R, and are given by (now written for the case with more than one ionic species) w (R) a7  = -k T\ng (R) B  ay  ,  (4.2.3)  where g -y{R) is the ion-ion radial distribution function obtained in the molecular solvent. a  For bulk electrolyte solutions it has been demonstrated [71] that good agreement is obtained between M M theory and full discrete solvent calculations, provided that no small ions (e.g. N a ) are present in high concentration. For the inhomogeneous system +  the approximations are expected to be more severe. The high wall charges typically used in our calculations always imply a "high" counterion concentration close to the walls and the validity of truncation at the pairwise level is not clear. In particular, the higher-order potentials of mean force also include effects due to changes in the solvent structure upon insertion of the ions, and this point will be important when effective ion-wall potentials are discussed in Section 5.2 . Finally, by employing the bulk pair potentials of mean force, one neglects the influence of solvent ordering close to the walls on the w . ai  Both the  fluctuation of the average solvent density near the walls and solvent particle orientation [62] might have significant effects. The remainder of the chapter investigates the M M level for one specific set of w (R), ai  W (z). a  while ignoring the effects of a molecular solvent on  Chapter 4.  4.3  Molecular Solvent Effects at the McMillan-Mayer  55  Level  The Model  The only fluid particles treated explicitly are hard-sphere ions. According to M M theory truncated at pairwise interactions, two ions at a separation R interact via a solventaveraged pair potential of mean force which can be conveniently expressed in the form w ,(R) a  = u™(R) + < ( P ) + v*(R)  ,  (4.3.1)  where u ^(R) represents the contribution from the molecular solvent and is not present in e  the P M . This "effective potential" will depend on the solvent model employed. Here, we 1  use results obtained [52] for a hard-core water model with a diameter of 2.8 A (i.e., equal diameters for ions and solvent) decorated with a dipole moment, a tetrahedral quadrupole moment, and with nonadditive polarizability effects treated at a self-consistent mean field level. A t a density appropriate for water at T = 298 K , and with the ions at infinite dilution, the pair correlation functions for this model can be obtained by solving the homogeneous R H N C equations [52]. The required ion-ion potentials of mean force are then obtained via Eq. (4.2.3). We note that more detailed water-like models (e.g. with an added octupole moment) have also been considered [72], but that the structural results are qualitatively similar to the simplified model described above. In the present calculations we are interested in general trends rather than specific details, hence we adopt the simplified model throughout. The potentials of mean force w ^(R) for monovalent ions with equal diameters d — a  +  cL = d = 2.8 A are shown in Fig. 4.1. Note that with the simplified water model employed here, positive and negative ions of the same size are solvated symmetrically such that w (R) ++  = w  (R). The most obvious feature at very short distances is the  decreased dielectric screening resulting from the fact that close to ion-ion contact solvent We use the notation w(R) to distinguish the pair potential of mean force from the pair potential u(R), although w(R) is employed as the "effective pair potential" in the HNC closure, Eq. (2.6.2). 1  Chapter 4.. Molecular Solvent Effects at the McMillan-Mayer  10.0  -  \  \ -\.  5.0 -  ? 0.0  MM +-  -  MM ++&-PM + PM + + & - -  -  —  1  56  Level  —  —  •  h  - 1  1.0  1  1  1.5  i  1  •  2.0 2.5 R/d  1 , 1  3.0  i  3.5  Figure 4.1: Ion-ion pair potentials of mean force at infinite dilution for monovalent hard-sphere ions with T = 298 K, e = 78.7, d = 2.8 A. Both M M and P M curves are shown.  Chapter 4. Molecular  Solvent Effects at the McMillan-Mayer  57  Level  particles cannot efficiently shield the Coulombic interaction. Hence, at small separations like-charged ions repel each other and oppositely charged ions attract each other much more strongly than in a corresponding continuum solvent. At intermediate distances the presence of the dense solvent induces further undulations which die off rapidly as R increases (at R&5d,  contributes less than 2.5% of w )  and the w  ay  ai  approach the  limiting 1/R behaviour. The simplified water model employed here gives<a bulk dielectric constant that is somewhat larger than the experimental value for water at 298 K (78.7), therefore we have adjusted u ^ (R) to be consistent with the correct dielectric constant. e  y  A n important feature (discussed below) is the minimum in w  at RtzlAd.  aj  note that the solvent-induced contribution u^ (R) a  hence w  aa  We also  has a short-range attractive tail and  approaches the continuum limit from below.  The ion-wall interaction is identical to the one employed in the last chapter and consists of the hard-sphere and electrostatic contributions v^ (z) and v^(z), respectively. s  The midplane formula for the pressure now includes an additional contribution due to the potential component  (i.e., because of the molecular nature of the solvent) and  one has -fslit = -Pkin + -Pcore + Pel + Peff •  P ff is calculated by using u (R) eS  e  4.4  as u (R) sh  in the formula for P  (4.3.2)  sh  [see Eq. (2.7.3)].  Results  A l l results reported are for T = 298 i f , d = 2.8 A , and e = 78.7. As noted above, the ion-ion potentials of mean force (see Fig. 4.1) were obtained with a water-like solvent model at a number density of 0.73/d , and with ions and solvent of equal diameter d. 3  Throughout we compare M M results with those given by corresponding (i.e., the same solvent dielectric constant and ion diameters) P M calculations.  Chapter 4. Molecular  Solvent Effects at the McMillan-Mayer  58  Level  We consider only symmetrical electrolytes and the slit systems are considered to be in equilibrium with a bulk solution with ion concentrations p = p_ = p. In the limit p —> 0 +  (infinite dilution) no coions are drawn into the cavity, and the counterions compensate the wall charges exactly. The assumption that no coions are present in the cavity is also a good approximation for finite but low bulk salt concentrations, and hence this simplified case is worth investigating before considering higher concentrations. In general, it is found that the smaller d  wal  i and the higher the wall charge density the better the fluid structure  (and therefore P i ) is approximated by the infinite dilution result. However, it must s  it  be noted that increasing the bulk density increases Pbuik- Therefore, because P i; is less s  t  sensitive to p for high surface charge, the net wall-wall interaction at small separations tends to become less repulsive (or more attractive) as the electrolyte concentration is increased [see Eq. (2.7.5)].  4.4.1  Systems in Equilibrium with a Bulk at Infinite Dilution  If the fluid slab consists purely of counterions, only the like-ion interactions (see Fig. 4.1) come into play. At the MM level the interaction between like ions is less repulsive than the P M case, except when the ions are very close to contact ( P < 1 . 2 d ) . Consequently, at the MM level the ions are more likely to be found near the midplane and less likely to be near the wall than in the PM. This means that for this system we would expect P  net  [MM] to be less than P t[PM] (note that for the infinite dilution case Pbuik = 0 and ne  -Pnet — -fslit) •  We first discuss a system of monovalent counterions with the wall charge density cr = -0.1307e/d = -0.267 C / m = - l e / 6 0 A . The density profiles are not shown because 2  2  2  they are similar to those given below for a 1 M system [see Fig. 4.9(a)]. The pressure curves obtained are plotted in Fig. 4.2. We see that not only is the net pressure for the MM system smaller than the corresponding P M result, but that the wall-wall interaction  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer  59  Level  0.08 0.06  MM, P  NET  PM, P  NEL  0.04 0.02 E-  0/00  CO  0.10  PH  0.05 0.00  &eroe*%%  -0.05  0  MM, P  EFF  PM,  / o  -0.10  A  fcore  JL  6  8  dwall/ d Figure 4.2: P and its components for systems including only monovalent counterions with o - -0.1307 e/d = -0.267 C / m . The Poisson-Boltzmann (PB) curve is shown for comparison (see Appendix C), where d n[PB] is redefined to give the same accessible space (i.e., the same z ) as the finite size ion models. P [ M M ] is less than 2 x 10~ ksT/d throughout and is therefore not shown. For comparison: 0.020 k T/d = 1.5 MRT = 3.7 x 10 N / m , where R is the gas constant. n e t  2  2  wa  m a x  5  3  B  3  6  2  core  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer  60  Level  is attractive (P t < 0) for 3.4 d < d n < 5 d. A comparison of the different pressure ne  wa  components from Eq. (4.3.2) demonstrates that this behaviour is solely due to the new pressure component P ff with P e  k i n  and P i being basically unchanged from the P M case. e  One has P ff < 0 as long as ions interact across the midplane mainly through the attractive e  part of uf {R).  The minimum of P ff at d n « 3.6 d corresponds to a configuration in  +  e  wa  which ions on one side of the midplane can interact with ions from the layer at the opposite wall at the optimal distance of R & 1.4. d corresponding to the minimum of u^ . +  At small wall separations, where the repulsive part of the ion-ion potential is more  important (both wall layers are in contact), one finds P ff[MM] RS P e  core  [ P M ] and moreover  P [MM] « P [PM]. net  net  The corresponding picture for a system with only half the previous wall charge (a = —0.0653 e/d ) 2  is given in Fig. 4.3. We note that P s is still negative for most wall e  separations, but much smaller in magnitude such that now P rather high wall charge o — —0.1634e/d  2  P  n e t  n e t  >0 for all <i n. For the wa  = —0.334C/m (achievable on mica sheets), 2  [MM] has a "double well" structure (Fig. 4.4) demonstrating that the location of  the minimum depends to some extent on a.  In order to investigate the influence  of the wall charge in greater detail, calculations were performed at fixed d n = Ad. wa  The results are shown in Fig. 4.5 and reveal that P t[MM] becomes more and more ne  negative with increasing \o\, again due to the pressure component P ff. e  Interestingly,  P [ P M ] remains comparatively constant over a wide range of a. These observations net  suggest that at high surface charges (when the counterion density between the walls is high) the solvent-induced "attraction" in the like-ion potential has a significant effect on macroion-macroion interactions. Clearly, under such conditions the P M description might be qualitatively incorrect. The case of divalent ions is more difficult to treat because the homogeneous H N C theory is not reliable for small highly charged ions. Nevertheless, it is possible to construct  Chapter 4.  Molecular Solvent Effects at the McMillan-Mayer  61  Level  0.08 MM, P  0.06 \  PM, P  \ \  0.04  •\  \  nel  net  PB, P  net  v.  0.02 E-  0.00  CO  x  0.10  MM, r^ PM, P  0.05  MM, P PM, P  o  in  kin  el el  0.00  *  —  MM, P  eff  -0.05  PM p  i i»i, i core  -0.10  4  6  *  —  *  -  H  A — —  8  dwall/ d Figure 4.3: P and its components for systems including only monovalent counterions with o= -0.0653 e/d = -0.134 C / m . The Poisson-Boltzmann (PB) curve is shown for comparison (see Appendix C). P [ M M ] is less than 2 x 1 0 ksT/d throughout and is therefore not shown. n e t  2  2  - 6  core  3  Chapter 4. Molecular  Solvent Effects at the McMillan-Mayer  62  Level  0.06 MM, P  0.04  PM, P  net  net  0.02 0.00 03  -0.02  MM, ri - x  CO  CU  in  0.10  PM, P  KM  o  MM, P  EL  0.05  PM, P  EL  -  -  0.00 -0.05 -0.10  PAAAA  A  MM, P  /C  >  A  EFF  PM P i xu, i core  4 dwall/  6  8  d  Figure 4.4: P and its components for systems including only monovalent counterions with o = -0.1634 e/d? = -0.334 C / m . P [ M M ] is less than 4 x 10~ k T/d throughout and is therefore not shown. n e t  2  5  core  3  B  Chapter 4.  Molecular  Solvent Effects at the McMillan-Mayer  Level  63  0.04  0.00  0.05  0.10  0.15  0.20  0.25  \o\d /e s  Figure 4.5: The o dependence of P counterions at fixed d n = 4 d. wa  n e t  and P ff for systems including only monovalent e  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer  64  Level  0.10  0.00  -  « • • 9-^*  CO  MM, It,, PM, P MM, P PM, P  DH  x  kin  o  el  el  MM, P  eff  4  A  6  8  dwaii/ d Figure 4.6: P t and its components for systems including only divalent counterions with o = -0.1307 e/d = -0.267 C / m . The Poisson-Boltzmann (PB) curve is shown for comparison (see Appendix C). P [ M M ] is less than 10~ k T/d and P [ P M ] is less than 4 x 10~ k T/d throughout, therefore these curves are not shown. ne  2  2  17  core  4  3  B  3  B  core  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer  reasonable w (R)  65  Level  for divalent ions by recognizing that in a linearized version of the  aj  theory (the L H N C approximation) w (R) ai  scales to a very good approximation as  \Z Z \ a  7  [73]; this can be interpreted as the attainment of structural saturation. Then, since for monovalent ions the potentials of mean force given by the H N C and L H N C theories are very similar, we take  = 4w (R)  W2+2+{R)  ++  to be at least a qualitative approximation for  divalent ions, assuming of course ion and solvent models analogous to those considered in the monovalent case. Results for a system with divalent counterions and a = —0.1307 e/d  2  ave given in Fig. 4.6. We note that in this case the P M also gives an attractive wall-wall interaction (in contrast to the P B theory [4]), which is attributed to correlated charge fluctuations. P [ M M ] is even more negative and has an extended attractive range due net  to the attractive component P ff. Also, we note that unlike the monovalent case, ion-ion e  hard-core interactions are of no importance for the divalent system. This is because the ions remain well apart from each other and, additionally, in the divalent system with the same o there are only half as many ions present as in the monovalent case. The minimum in P ff is now at e  d \\Xi2.8d, wa  where the favorable interaction via u ± e  +  takes place directly  between the two layers adjacent to the walls.  4.4.2  Systems at Finite Bulk Concentration  When coions as well as counterions occupy the slit the numerical calculations require extra care, as close to ion-ion contact  and c _ are very narrow and problems with +  long-range tails in Fourier space arise. Details are given in Appendix A.2 . In addition, for this case pressure results obtained from the contact theorem and the midplane formula differ qualitatively, as shown in Fig. 4.7. A third route to the pressure is given by * calculating the total grand potential (per unit surface area), Q , see Appendix B.2 . A  One then obtains the pressure by taking its (negative) derivative with respect to d n. wa  This method has the conceptual advantage that the derivation of the expression for Q  A  Chapter 4. Molecular  Solvent Effects at the McMillan-Mayer  66  Level  contact plane p=1.0M midplane -aQ /ad „  0.10  n  A  wa  EPQ  0.05  only contact plane counterions midplane  x  ^—a—a—  A  0.00 h  CO  CO  -0.05 \-  a •  -0.10  _ O •  Ijl  n a  D  •  •  d  °  D  4  6  8  dwall/ d Figure 4.7: Comparison of P i calculated via the contact theorem and via the midplane formula at o = —0.1307 e/d . Results are shown for a 1:1 electrolyte solution in equilibrium with a bulk at 1.0 M , and for a system where only monovalent counterions are included. For the former case the derivative of the total grand potential (per unit surface area) with respect to d n is also shown. s  2  wa  it  Chapter  4.  Molecular Solvent Effects at the McMillan-Mayer  67  Level  exploits the H N C approximation, and therefore the level of approximation involved is consistent with that of the correlation functions. Moreover, f i depends "symmetrically" A  on the densities and correlation functions from the whole slab and is less likely to be sensitive to inaccuracies close to the wall or around the midplane. The results given in Fig. 4.7 show that the contact theorem, Eq. (2.7.2), agrees well with  —dft /dd ,\\. A  wa  Of course, this does not prove that the contact theorem delivers the correct pressure, but it does indicate a high level of internal consistency in the A H N C calculations. The problem with Eq. (4.3.2) can be understood by considering that P strongly on g  ai  c o r e  and _Pff depend e  at, and close to ion-ion contact, respectively. From P M calculations  using anisotropic integral equations it is known that the H N C closure overestimates the g - contact values [43]. For the M M model the cation-anion attractions are very strong, +  therefore contact value errors likely contribute to the problems encountered using the midplane formula. This also explains why for systems with only counterions between the walls, the midplane and contact plane routes yield essentially the same pressures (see Fig. 4.7). In this subsection, we present only results for P \\ calculated via the contact s  t  theorem. For finite bulk concentrations one has first to determine the activity coefficient y± for the model in question. This is done most conveniently for uncharged walls; in this case the density profiles are rather flat and, with increasing d n, the midplane concentrations wa  reach the bulk value p faster than for t r ^ O . For the 1:1 electrolyte, we find j± = 0.38, p = 0.0136/d = 1.03 M for the M M case, and j 3  = 0.58, p = 0.0132/d = 1.00 M for the 3  ±  P M . Some sample density profiles for these systems are shown in Fig. 4.8 (Note that due to the symmetry of the system p ( ) z  a  = Pa(~z) and therefore only one half of each profile  is displayed). Compared with the P M , in the M M system the ions exhibit a stronger tendency to be located away from the walls. This is likely because in the M M case the cation-anion contact interaction is more favorable than in the P M , and more such  Chapter 4.  Molecular Solvent Effects at the McMillan-Mayer  Level  68  interactions are possible when the ions are away from the walls. For the same reason, with decreasing d n we observe a higher degree of depletion of ions from the slit in the wa  M M model than in the P M . Corresponding density profiles for a wall charge density of —0.1307 e/d  2  lowing discussion we shall refer to this value as  (Thigh)  (in the fol-  are shown in Fig. 4.9. For this  system the counterions are the overwhelmingly predominant species present in the slit and the majority of these are found close to the walls which maximizes the average ionion separation. Nevertheless, the contact densities are smaller than for the P M because there is less repulsion among the counterions. Interestingly, in the M M system the coion densities are much higher than in the P M case and even exceed p around the midplane for certain wall separations. This behaviour is analyzed in detail in Fig. 4.10(a). At (Thigh the coion midplane density rises above p for d n }t 4 d but always remains lower than the wa  corresponding counterion density. For d \\> 10d we find p+(0) « P-(0) and this value wa  decreases slowly to p as the walls are pulled apart. For a system with half the wall charge (ciow)  this effect has basically disappeared and the p (0) curves are very similar to the Q  P M [Fig. 4.10(b)]. In both Figs. 4.9(b) and 4.10(a) oscillations in the M M coion densities indicate strong structuring in the fluid (due to the increased coion-counterion attraction), which we will encounter again below in the discussion of correlation functions. The net pressures for three different values of o are compared in Fig. 4.11. For high o we find a region of attractive wall-wall interaction as in the infinite dilution case. Now, however, within numerical accuracy P  n e t  remains negative for large  d aii> W  i.e.,  the "barrier" seen in Fig. 4.2(a) has vanished. Lowering the wall charge results again in a purely repulsive interaction nearly identical in size to the P M case, but note the opposite trends for M M and P M systems when changing o. Our observations are in qualitative accord with results presented by Kinoshita and coworkers (see Fig. 19 of Ref. [66]) although they considered spherical macroions and the level of theory is distinctly  Chapter 4.  Molecular  Solvent Effects at the McMillan-Mayer  Level  69  0.015  co  0.010  N  MM  S-O.OOS -  PM I  0.000  -3.0  i  -2.0  -1.0  0.0  z/d Figure 4.8: Density profiles for 1:1 electrolyte solutions between uncharged walls, p is 1.00 M = 0.0132/d and 1.03 M = 0.0136/d for the P M and the M M model, respectively. p is indicated by horizontal lines. Results are shown for d n = Id, 2d, and 8d. The profiles for cations and anions are identical. 3  3  wa  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level  70  0.3  0.2 -4 0.1 CO  Tl N  0.0  o Q.  0.015  h  0.010 \0.005 r0.000  -2.5  -2.0  1.5 -1.0  0.5  0.0  z/d Figure 4.9: Density profiles for 1:1 electrolyte solutions between charged walls with cr = —0.1307e/d = —0.267C/m for (a) counterions and (b) coions. p is LOOM = 0.0132/d and 1.03M = 0.0136/d for the P M and the M M model, respectively. p is indicated by horizontal lines. Results are shown for d \\ = 2.5 d, 4.0 d, and 6.0 d. 2  3  2  3  w&  71  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level  0.025  -(a)  •  A  A  A  • a  A  coions  A A  A  • • A A  0.015  ^  A  cr * low  n  D  ^high  counterions  A  0.020  .  A  <7high  a  1(  6  D  ow  .  a  0.010 _PJ_  1  (b)  J  \\  I  I  I  counterions  \ *  0.020  I  I  L  a h i g h  ^low  0.015  0.010  2  4  6  8  10  12  14  dwaii/ d Figure 4.10: Midplane densities of monovalent counterions and coions at fhigh = -0.1307 e/d and c r = a /2 for (a) the MM system and (b) the PM. The bulk densities p are as in Fig. 4.9 and are indicated by horizontal dotted lines. 2  l o w  high  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer  -0.02  1  1  2  '  « 4  '  ' 6  72  Level  '  '  1  8  1  10  dwan/ d Figure 4.11: P for 1:1 electrolyte solutions at p = 1.0M for o"iow = Chigh/2, and for uncharged walls. net  a i h h  g  = -0.1307e/d , 2  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer  Level  73  different from ours. In Ref. [66] the macroion charges Qu = 20 e and 40 e correspond to ourCTIOWand o"high- The change from negative to positive slope of the macroion-macroion potential of mean force as Q M is increased is equivalent to the switching from a repulsive to an attractive wall-wall interaction in the present calculations. The concentration dependence of P  net  for cr igh is shown in Fig. 4.12. Curiously, n  the curves for p= 1.0 M and 0.5 M are almost indistinguishable. In general, a higher bulk concentration produces a less repulsive (or more attractive) wall-wall interaction which is in agreement with results presented in Ref. [66], and also with the PM. For very high wall charges ( c r ^ 0.18 e/d ) this trend is reversed, as shown i n Fig. 4.13 where 2  we compare the variation of P  net  with a at fixed wall separation for p — 1.0 M and  at infinite dilution. Note that these high values of a are unlikely to be encountered in colloidal systems, e.g., mica sheets have up to a fa 0.16 e/d . 2  From Fig. 4.13 we  learn as well that the overall influence of the coions on the mechanism of the observed wall-wall interaction is rather minor. The same explanation for the observed wall-wall attraction holds for systems in equilibrium with bulk electrolytes at high concentration and at infinite dilution. Specifically, the slightly attractive correction to the counterioncounterion pair potential (due to the molecular solvent) tends on average to pull the ions away from the walls when the walls are brought together. The strong coion-counterion attraction at contact does not contribute to the lowering of the counterion density at the walls at high surface charge, but it does cause intriguing changes of the fluid structure (as described above and in the following subsection). Coions have their largest effect as dwaii—>co (i-e., they influence the value of Pbuik), whereas their overall concentration is low at high o and small d n. wa  Chapter 4.  Molecular Solvent Effects at the McMillan-Mayer  Level  74  Figure 4.12: P for M M systems in equilibrium with 1:1 electrolyte solutions at various bulk concentrations and cr = -0.1307 e/d = -0.267 C / m . The curves for p = 0.5 M and 1.0 M are indistinguishable within numerical error and therefore only the latter case is shown. n e t  2  2  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level  75  0.02  0.00  0.05  0.10  0.15  k|d /e  0.20  0.25  2  Figure 4.13: The dependence of P on o at fixed d n = Ad. Results are shown for 1:1 electrolyte solutions at p = 1.0 M and for the the infinite dilution (counterions only) case. net  wa  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level  4.4.3  76  Pair Distribution Functions  The anisotropic pair distribution functions, # (r, Zi, z ), contain information about the Q7  2  correlation between two ions in the sense that given particle 1, type a, is located at position (r = 0,2:1) the probability density of finding particle 2, type 7, at (r, z ) is 2  27rrp (z )<jr (r, z\, z ). As g 7  2  Q7  2  ai  depends on three variables it is difficult to depict in full  and here we select a few examples to illustrate some general features. <?a (r, 0,0) (i.e., at the midplane, parallel to the walls) for a 1.0 M , 1:1 electrolyte 7  in a slit with d n = 4 d is shown in Fig. 4.14. Whereas in the P M case all functions wa  (g+-, g++, g  ) reach the limiting value of one monotonically, the M M results exhibit  strong oscillations indicating a much more "structured" fluid. Especially striking is c/ _ +  which shows a "shell-like" layering of anions and cations in the M M system. The contact value for  is very large (~ 100) [see Fig. 4.14(a), inset] and this creates the numerical  difficulties discussed in the previous section. The height of the first peak in c/ _ does not +  change much if one looks in directions not parallel to the walls, and a quick calculation (integration from r —1.0 d to 1.1 d, using an average counterion density of 0.08/d around 3  a midplane coion) reveals that the average number of counterions in the first "shell" is about 0.6, i.e., approximately 60% of the coions are in close contact with a counterion. This number is much higher than the bulk value (~ 10%) and increases further with decreasing d n as fewer coions are surrounded by more counterions. The coion-coion wa  correlation function, g  , shows a heightened affinity of coions for each other in the range  1.3 d ; $ r < 2.0 d, followed by a slowly decaying oscillatory tail. A possible explanation is a tendency towards the formation of mixed "clusters". This could account for the pronounced structural features in g  , whereas g  ++  is not affected to the same degree as  there are many more counterions than coions and only a few counterions (~ 8%) are in the close vicinity of a coion. In fact, g  ++  at infinite dilution is very similar to the 1.0 M  Chapter  4.  Molecular Solvent Effects at the McMillan-Mayer  Level  77  r/d Figure 4.14: g {r) at the midplane parallel to the walls (zi = z = 0) for 1:1 electrolyte cr = —0.1307 e/d . g+_=g_ is shown in (a) and g and g ("+" denotes counterion and "—" the coion) are shown in (b). ai  2  2  +  ++  Chapter 4. Molecular  Solvent Effects at the McMillan-Mayer  78  Level  case, but with slightly lower values and a less structured tail. We now turn to a few examples of two-dimensional contour plots where the analysis is much more qualitative in nature. As a general trend, one finds that the g  aj  parallel and  perpendicular to the walls are surprisingly similar and this is especially true for cationanion pairs. Figures 4.15 and 4.16 show g  for monovalent and divalent counterions,  ++  respectively, at the wall, with c? n corresponding to the minima of P ff in Figs. 4.2 and wa  e  4.6. The white region around the particle is forbidden to other ions because of the hardcore interactions. For monovalent ions in the P M case, g  ++  nowhere exceeds one but  the wall introduces a degree of anisotropy. The M M system exhibits structure similar to that of the midplane g  ++  shown in Fig. 4.14; the structure is qualitatively similar in  all directions and is distinctly not monotonic. The ordering of ions is more pronounced in the divalent case for both models. However, whereas the most probable position for a neighboring ion in the P M is clearly at the opposite wall with r « ±3d,  the M M  system shows a more complicated structure with an increased tendency to find ions at a separation R—Jr  2  + (z + z  ) ~l. 2  m a x  In order to interpret the contour plots in terms of average densities, g  ai  multiplied by the density profile of ion 7, p {z). y  must be  The average density of counterions  about a midplane coion is shown in Fig. 4.17 and demonstrates again the shell-like anion-cation structure which is totally absent in the P M . Another interesting detail (at the selected wall separation) is the decrease in counterion density at the walls near the coion location. In comparison, the corresponding picture for the P M shows an increased P+G-Zmaxl)-  Note that, as discussed above, this effect is not the reason for the wall-wall  attraction at d n « Ad. Finally, Fig. 4.18 shows the ordering of coions at preferred wa  distances from a counterion in contact with a wall (i.e., the most probable counterion position). This ordering results in the oscillatory structure in the coion densities seen in Figs. 4.9(b) and 4.10(a).  Chapter 4.  Molecular Solvent Effects at the McMillan-Mayer  79  Level  g(r,z) for counterion/counterion MM 1.1500  m 1.0500  •  1.0100  -6.0  I  •4.0  I  I  -2.0  0.0 r  •  +4.0  +2.0  +6.0  0.9900 r~|  0.9000  •  g(r,z) for counterion/counterion PM  0.8000  • •  04000  I  -6.0  l  -4.0  -2.0  0.0 r  ^  ^  ^  +2.0  +4.0  ^  +6.0  • 0.0000 •  Figure 4.15: Contour plots of g++(r, -z ,z) for systems where only monovalent counterions are present with d u = 3.6d and o = -0.1307e/d . r and z are in units of d. The black horizontal bars represent the walls and the legend on the right explains the colour scale, with the number between two squares indicating the boundary value. max  2  wa  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer  Level  80  g(r,z) for counterion/counterion  •  MM  1.1500  ^^^^^^^^^^^^P^r^^^^^^^^^T^ -6.0  -4.0  -2.0  0.0  +2.0  +4.0  +6.0  0.0000  D  Figure 4.16: Contour plots of g++(r, -z , z) for systems where only divalent counterions are present with d n = 2.8d and o = -0.1307 e/d . max  2  wa  Chapter 4. Molecular  Solvent Effects at the McMillan-Mayer  Level  81  fg(r,z) dens(z)J for coion/counterion  r Figure 4.17: Contour plots of p (z)g_ (r, 0, z)d , the average counterion density around a midplane coion. The ions are monovalent and p = 1 . 0 M , d , = 4d, and a = -0.1307e/d . 3  +  +  2  wa  n  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level  82  [g(r,z) dens(z)J for counterion/coion  MM 0.0500  •  0.0200  •  0.0100  •  0.0080 0.0060 [g(r,z) dens(z)J for counterion/coion  PM  • • •  0.0040 0.0020 0.0010  • 0.0000 • Figure 4.18: Contour plots of p_(z)p+_(r, -z , z)d?, the average coion density around a counterion in contact with one wall. The system parameters are as in Fig. 4.17. max  Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 4.5  83  Summary  The force between two like-charged walls immersed in electrolyte solution can be calculated by adding two contributions.  One is given by the force for a corresponding  uncharged slit system, whereas the other is obtained from the charged system with the electrolyte solution described at the McMillan-Mayer level. At this level the solvent is not represented by discrete particles but exerts its influence through solvent-averaged ion-ion potentials of mean force which serve as effective potentials. The scheme outlined is exact only if n-particle potentials of mean force are included. However, usually only pair potentials of mean force are considered, with the additional approximation that they are obtained from bulk models. The results reported in this chapter show that the introduction of solvent effects through effective ion-ion pair potentials can add an attractive component to the pressure not present in P M treatments of electrolyte solutions between charged walls. This effect can be very important, leading to an attractive wall-wall force for a model 1:1 electrolyte if the wall charge density is sufficiently high. The wall-wall attraction increases in magnitude and in range with increasing wall charge and electrolyte concentration. We find qualitative agreement with the homogeneous R H N C calculations of Kinoshita et al. [66] for spherical macroions in electrolyte solutions with solvent particles explicitly included. This suggests that the much simpler homogeneous R H N C approach could prove more useful for ion-solvent-macroion systems than one might have expected. The effect of the presence of coions is basically limited to their influence on the value of Pbuik> because their overall concentration in the slit is low at high o and small d i i . Nevertheless, the wa  fluid structure around the coions that are present is drastically altered compared to the PM.  Chapter 5 Further Investigations at the McMillan-Mayer Level  Three different topics are discussed in this chapter. First, we investigate the physical significance of the results described in the preceding chapter by employing other model potentials of mean force and comparing the results.  A n analysis of the influence of  effective ion-wall potentials of mean force on the pressure is also given. This requires a more thorough discussion of the pairwise level approximation of M M theory. We close the chapter by remarking on the possible relevance of the present model calculations in the interpretation of some reported experimental results.  5.1 M c M i l l a n - M a y e r Results for Several Electrolyte Models In this section we apply M M theory to what are generally regarded as more realistic models for aqueous electrolyte solutions. Ion-ion potentials of mean force for N a C l and other salts are available in the literature [52,74-76,30], but, even for models that attempt to mimic the same physical system, the potentials of mean force can differ markedly. Moreover, in the absence of any direct experimental measurements, one has no way of selecting the model that best represents a particular physical system. Therefore, we do not attempt to discuss the merits and/or disadvantages of particular models. Rather, we simply view the different potentials of mean force as inputs for the M M theory that hopefully span a significant portion of the range of physical possibilities. In fact the interaction of charged walls turns out to be quite sensitive to details of the models employed. 84  Chapter 5. Further Investigations at the McMillan-Mayer  5.1.1  85  Level  Models  The ionic fluids we consider here are based on electrolyte models with "soft" ions, which do not have well-defined diameters. Nevertheless, to avoid specifying short-range ionwall interactions, we assume that planes of closest approach exist that bound the space accessible to the particle centers. The separation of these planes (i.e, the width of the accessible space) is denoted D. Note that for a model consisting of hard spherical particles of diameter d;  on  and hard walls at a separation d n (as in Chapter 4) Z) = d n — dj . wa  wa  on  Formally, the solvent-averaged pair potential of mean force can still be expressed as in Eq. (4.3.1), with u ^ , representing the short-ranged part of the ion-ion potential of mean force that is due to the molecular nature of the solvent. Figure 5.1 summarizes the different  selected from the literature that are con-  sidered in the present section. A l l models are meant to represent aqueous electrolyte solutions at infinite dilution and at a temperature T = 298K. Guardia, Rey, and Padro (GRP) [74,75] performed M D simulations employing a flexible version of the S P C (simple point charge) water model [77]. The curves labeled G R P are obtained by reading their data points, subtracting the long-range Coulombic contribution assuming a dielectric constant of 78, and applying a simple polynomial interpolation. The G R P u^  1  to zero for R>7.$A.  are set  The label P R denotes results obtained by Pettitt and Rossky [76]  using R I S M theory. M D simulations have shown that the R I S M approximation is quite accurate for the models they consider [78]. In the Pettitt-Rossky calculations the solvent molecules are modeled by a modified version of the TIPS (transferable intermolecular potential functions) water model [79]. Unlike the G R P case, this model does not include intramolecular vibrations. Although the dielectric constant given by the TIPS model is much lower than the experimental value for water, this is corrected in the R I S M calculations where e is set to 78. The RISM curves can be fitted nicely (in contrast to the G R P  Chapter 5. Further Investigations  at the McMillan-Mayer  Level  86  R/A Figure 5.1: u (R) components of the ion-ion potentials of mean force at infinite dilution. The curves labeled [KP], [GRP] and [PR] are from Refs [52], [74,75] and [76], respectively. eff  Chapter 5. Further Investigations at the McMillan-Mayer  a / (k TA) 6.95 • 10 770. B  Na+/Na+ ci-/cr  12  b / A~  c/  l  -8.685 -1.247  {k TA) 7.96 657. B  d/1-  Level  1  -0.3172 -0.9618  e/  87  A"  f  1  3.01 1.96  -2.71 -3.94  Table 5.1: F i t parameters for the Pettitt-Rossky [76] ion-ion potentials of mean force. The units are chosen to yield u in k T if R is given in A. eff  B  results) by the function u (R) = a exp(bR)/R + c exp(dR) cos(eR + f)/R ef!  ,  (5.1.1)  where the fit parameters are given in Table 5.1. A similar fit was used by Zhong and Friedman [80]. The curve labeled K P was obtained by Kusalik and Patey [52] and is identical to the like-ion potential of mean force employed before in Chapter 4. We note that the first minimum in the K P w for the N a / N a +  +  eff  is only slightly deeper than the corresponding one  [PR] model, and comparing results for both alternatives gives insight  into how small differences in the effective ion-ion potentials can influence the wall-wall interaction. We remark that from a numerical perspective soft-core models are easier to solve at finite concentration because the complications arising from the very strong contact peak for oppositely charged hard-core ions do not occur in the soft-core case. Also, for numerical reasons, a hard-core cutoff is imposed on the ion-ion potentials at separations sufficiently small that the results are physically unaffected; specifically, we set w (R) = o o ay  if w (R) ^ 10 k T. We stress that although the exact physical location of the walls, and aj  B  therefore the location of the surface charges, is not defined in our model (see above), no free parameter exists, as has been noted in Section 2.6 . Relocating the surface charge planes further apart (keeping D fixed) adjusts the chemical potentials of individual ion species, but leaves p± and hence the physical behaviour unchanged. This follows from the imposition of strict electroneutrality in the slit.  Chapter 5. Further Investigations at the McMillan-Mayer Level  5.1.2  88  Results  A l l results reported are for aqueous systems at 298 K.  Unless otherwise stated, the  solutions between the walls are taken to be in equilibrium with a bulk at infinite dilution, i.e., there are no coions in the slit. The magnitude of the surface charge density, \a\, is 0.334 C / m = l e / 4 8 A . The walls may be negatively or positively charged depending on 2  2  the sign of the species we wish to consider as counterions. We first examine a system with N a [GRP] serving as counterions. This system was +  previously considered by Marcelja [50,70], also employing the A H N C theory, but due to a numerical error [51] the pressure results reported are incorrect. We note that for this model u^ (R) +  has a pronounced first minimum of ~ —1.8k T at R = 3.6 A (see Fig. B  5.1). The net pressure and its components are shown in Fig. 5.2. We see that P  n e t  has  a negative region beginning at a separation of ~ 13 A . This indicates a net attraction between the walls and this is due to the attractive part of P ff. At small D, P e  net  becomes  . positive partly because the ion layers at the planes of closest approach to the walls interact via the repulsive part of w . Also, the component Pki grows with decreasing eff  n  slit width as the average ion density between the walls increases. Results obtained for Na  +  [PR] counterions are also shown in Fig. 5.2. For the P R model u + has a shallow e  +  first minimum of ~ —0.6 k T followed by a slowly decaying oscillatory tail. The resulting B  pressure component  P ff e  has a very weak attractive part and exhibits small undulations  obviously reflecting the oscillations in u  e f F  .  Overall, unlike the G R P case discussed above  and the K P model considered in the last chapter, the P R potentials of mean force do not give a net attractive interaction. This is somewhat surprising given the similarity of the K P and P R potentials of mean force and serves to demonstrate the high sensitivity of the net pressure to details of the effective counterion-counterion interactions. Results for C l ~ [GRP] and C l ~ [PR] counterions (i.e., positively charged walls) are  Figure 5.2: The dependence of P and its components on the wall separation for N a [GRP] and N a [PR] counterions. The wall charge density cr = - 0 . 3 3 4 C / m = - l e / 4 8 A . M represents moles/litre and-R is the gas constant. +  n e t  +  2  2  Chapter 5. Further Investigations  at the McMillan-Mayer  90  Level  Figure 5.3: The dependence of P and its components on the wall separation for C l [GRP] and C l ~ [PR] counterions. The wall charge density cr = +0.334 C / m = + l e / 4 8 A . The notation is as in Fig. 5.2. n e t  2  2  Chapter 5.  Further  Investigations  at the McMillan-Mayer  Level  91  given in Fig. 5.3. For the G R P model, u__ is positive throughout and becomes strongly repulsive for R<§A  (i.e., the ions have large "diameters"). Due to the strong repulsive  core interaction, this model forms very pronounced layers at the wall. This gives rise to a highly repulsive pressure component  (and hence a repulsive P )  P ff  as the walls  nei  e  are brought together. In contrast, the midplane density and P^n remain rather constant as D decreases. In the C l ~ [PR] case, «!__ has a strong attractive region with a miniff  mum at R « 3.5 A. Ions close to the midplane are best able to take advantage of this favorable interaction, and hence the corresponding density profile (Fig. 5.4) exhibits a local maximum at z = 0. The resulting pressure curves shown in Fig. 5.3 are distinctly different from the G R P model in that P  n e t  now indicates a strong wall-wall attraction  at intermediate separations. Again, this demonstrates that the net pressure is strongly model dependent. In order to illustrate the dependence of the fluid structure on the u  eS  employed,  density profiles and pair distribution functions along the wall for D = 7.25 A are shown in Figs. 5.4 and 5.5, respectively. Large counterions (Cl~ [GRP]) and those with rather shallow first attractive wells ( N a [PR]) or no attractive well at all ( P M ions) have a +  higher density at the wall, followed by a monotonic decay towards the midplane. Smaller ions with clearly defined attractive wells tend to yield lower ion concentrations at the wall and some "layering" around the midplane. As one would expect, the pair distribution functions display maxima at ion-ion separations corresponding to the minima in  u. eS  Strong oscillations in u^ produce more pronounced oscillations and higher peaks in a  9aa{f)-  It is interesting to compare the N a  +  [PR] and E q  +  [KP] curves.  For these  models the heights of the first maxima are nearly equal, but for E q [KP] the long-range +  behaviour of g (r) is less oscillatory and lies closer to the P M result. aa  The dependence of P  n e t  on the bulk salt concentration was discussed in Section 4.4.2  for the EqEq [KP] case. The conclusions drawn earlier hold equally here, and therefore  Figure 5.4: Density profiles for the systems shown in Figs. 5.2 and 5.3 at the wall separation £> = 7.25A. The curve for N a [ P R ] is within ±0.5 M of the result for monovalent charged hard spheres ( P M ) with d i = 2.94A. The notation is as in Fig. 5.2. +  o n  Chapter 5. Further Investigations at the McMillan-Mayer Level  93  r/A Figure 5.5: Ion-ion pair distribution functions at the wall plane (z\ = z = D/2) for the systems shown in Figs. 5.2 and 5.3 at the wall separation D = 7.25A. Curves obtained for E q [KP] and monovalent charged hard spheres with dj = 2.94A (PM) at the same wall charge are also shown. The C l [PR] curve reaches a maximum of 8.3. 2  +  0n  -  Chapter 5. Further Investigations at the McMillan-Mayer  the present discussion is restricted to one model. P  n e t  Level  94  curves for the N a C l [GRP] model at  three different salt concentrations are plotted in Fig. 5.6. Overall the effect of increasing the bulk salt concentration is minor at the high wall charges considered here. This is because very few coions are found between the walls at small D. We do note that as the bulk concentration is increased, P  n e t  is lowered slightly (except for 8 A < P > < 1 3 A at  1.0 M), and the repulsive tail found for £ > > 1 3 A in the infinite dilution case disappears due to ionic screening.  5.1.3  Summary  We have investigated several model systems consisting of monovalent counterions between charged walls, where w  eff  was extracted from solvent-averaged potentials of mean force  that have been reported in the literature. Unfortunately, even when the same solution such as NaCl in water is the objective, different model assumptions for the ion-ion, ion-solvent and solvent-solvent interactions can lead to markedly different results for u . These in turn can give rise to different ion "structure" between the walls and to eff  dramatically different wall-wall interactions. Given this state of affairs, and the fact that one has no way of deciding which model best represents a particular physical system, the main purpose of this investigation is to demonstrate the range of possibilities for models in the "physically realistic" range. Some present and previous results are summarized and compared in Fig. 5.7. The strong model dependence is apparent from the figure. In particular, N a [GRP] yields a +  net attraction at slit widths of a few ion "diameters", whereas for N a [PR] the interaction +  is always repulsive and, moreover, P  n e t  lies close to the simple primitive model result. It  is also interesting to note that an attractive well in u  eS  significantly. For example, E q but E q  +  +  [KP] and N a  +  does not necessarily reduce P  n e t  [PR] have comparable attractive wells,  [KP] leads to an attractive wall-wall interaction whereas N a i [PR] does not.  Chapter 5. Further Investigations  0  at the McMillan-Mayer  5  10 D/A  Level  95  15  20  Figure 5.6: P for the NaClJGRP] model [74,75] at different salt concentrations, p, and o = -0.334 C / m = - l e / 4 8 A . The inset shows the tails on an expanded scale. The notation is as in Fig. 5.2. n e t  2  2  Chapter 5. Further Investigations at the McMillan-Mayer  Level  96  Figure 5.7: Comparison of P for four systems with monovalent counterions, including the P M with d i o n = 2.94.A, and o = -0.334 C / m = - l e / 4 8 l . The curves for N a [GRP] and N a [PR] are as in Fig. 5.2 and the E q [KP] result is from Fig. 4.4 . The notation is as in Fig. 5.2. n e t  2  +  +  2  +  Chapter 5. Further Investigations at the McMillan-Mayer  Level  Even calculations with a toy model that has 1.8 times the N a [PR] w +  the first well is —1.0 ksT) still yield a positive P  n e t  97  eff  (the depth of  at all wall separations, although its  magnitude is lower (by approximately a factor 1/2 for 5 A < D < 10 A). This indicates that the longer ranged behaviour of u  eff  is also influential. It appears that the equilibrium  structure of the fluid, and hence the net pressure, depends on subtle energetic balances that make qualitative predictions difficult. Nevertheless, for at least some models of N a in water, u  ef[  +  does display a strongly attractive well and/or an attractive tail [52,74,30].  Thus, solvent-mediated counterion-counterion interactions could lead to attractive wallwall forces, and this possibility should be considered seriously in the interpretation of experimental observations for negatively charged mica sheets in aqueous N a C l solution. It should be noted that the surface charge density used in the present calculations is rather high, but not outside of the physically relevant range. For smaller surface charges, there are fewer counterions in the slit and the short-range part of the ion-ion potential becomes less important. This results in net pressures that are more like those obtained in the primitive model.  5.2  Ion-Wall Potentials of Mean Force  So far the influence of the walls on the fluid structure in the slit has simply been to restrict the accessible space of the particles. With respect to solvent order induced by the surfaces two effects have been ignored. As the solvent structure in the vicinity of the walls is different from that of the bulk, the effective ion-ion potentials in the slit are changed in comparison to the bulk potentials of mean force. Further, solvent-dependent effective ion-wall potentials (akin to the effective ion-ion potentials in the bulk) come into play. The primary effect of solvent order at the wall on the ion-ion potentials is a reduction of the dielectric constant close to the surface [81], but this point is ignored  Chapter 5. Further Investigations at the McMillan-Mayer  here.  Level  98  However, the inclusion of the additional component of the ion-wall interaction  seems necessary in order to construct a self-consistent model where ions interact with each other and with the walls at the same (solvent mediated) level. To obtain the ion-wall potentials of mean force it is necessary to solve for the structure of the full electrolyte solution model close to a macroparticle. This has been done in a series of papers by Torrie et al. [62,82-84] for hard-core particles and the solvent model described in Chapter 4. Another interesting study by Vossen and Forstmann [30] of a model NaCl solution employs the central force water model close to a planar wall (with rather low surface charges) utilizing a mean density variant of the R H N C theory.  5.2.1  The Model  We consider again the electrolyte solution model employed in Chapter 4, which is based on hard-sphere ions and solvent particles with equal diameters d =2.8 A . The resulting bulk ion-ion pair potentials of mean force at infinite dilution, w (R),  calculated by  ai  Kusalik and Patey [52], are shown in Fig. 4.1. Here, we want to include the effect of the molecular solvent on the ion-wall potentials of mean force, W (z). a  This "effective"  ion-wall interaction can be written in the form [see Eq. (2.5.2)] W (z) a  = V (z) + V*(z) US  +V  Q  ( < W  +  Z) + MZmzx  ~ z)  ,  (5.2.1)  where tp (C) is the solvent-induced short-range interaction of an ion a with one wall as a  a function of the distance to the plane of closest approach, £ > 0. B y using the ansatz in Eq. (5.2.1) we implicitly assume that both walls give identical superimposable contributions and that ip does not depend on the wall separation itself. These assumptions are a  reasonable for large wall separations but are clearly an approximation at short d n. wa  A good model for ip (() is difficult to choose. a  Marcelja suggested the use of the  corresponding solvent-averaged potential of mean force for a single ion at an uncharged  Chapter  5. Further Investigations  at the McMillan-Mayer  Level  99  wall [69]. A t first glance this proposition seems reasonable when considering Eqs. (4.2.1) and (4.2.2), where the ion-wall potentials of mean force for an ion in the uncharged slit (i.e., pure solvent) system are required in W„. This contribution is omitted so far in our calculations (apart from the electrostatic and hard-core contributions) and a rather straightforward implementation of results obtained by Torrie et al. [82] for an ion immersed in a molecular solvent close to a neutral macroparticle is feasible. However, one can arrive at a different point of view, if one considers possible effects of truncating the n-particle potential of mean force at the pairwise level. The M M scheme at the pairwise level relies on the fact that the solvent structure is approximately independent of the ion concentration. Therefore, infinite dilution solventaveraged potentials of mean force can be used at higher ion concentrations, as has been demonstrated for the bulk case up to salt concentrations of several moles/litre [71]. In contrast, the solvent structure close to a single charged macroion can change dramatically due to polarization effects. Counterions at finite concentration placed between walls tend to yield dense ion wall layers, as the ions try to maximize their mutual distance, with typical ion densities at the wall around 25 M (for wall charges and models investigated in the preceeding chapter). The additional electric field due to the inhomogeneous ion distribution can polarize the solvent and the resulting equilibrium solvent structure is consequently distinctly different from the corresponding infinite dilution case [83]. M M theory including multibody ion potentials of mean force is exact and the changes in solvent structure are therefore described correctly, whereas at the pairwise level the effects due to solvent restructuring are lost. The use of ion-wall potentials of mean force calculated for a charged macroparticle and a finite counterion concentration therefore seems reasonable, as for such systems the solvent structure is similar to the system of interest. In this way, some of the errors introduced by neglecting higher order effective ion-ion interactions might be compensated.  Chapter 5. Further Investigations at the McMillan-Mayer  Level  100  Note that for small to medium wall separations none of the choices for ip(() discussed here might be appropriate as the solvent structure/polarization depends critically on the ion distribution, which changes with d ii- In addition, the superposition ansatz made in wa  Eq. (5.2.1) is no longer valid. Figure 5.8 shows the two models for ^>+(C) used here, both are guided by results for macroion-ion-solvent mixtures solved at the R H N C level [82,83]. The spherical macroions employed in those studies had a diameter of 30 d. Macroion curvature effects and problems with the planar wall limit are discussed in Ref. [62]. For convenience, the potentials are chosen to be nonzero only for a wall layer of thickness « 1 d. The curve labeled "•01" has been obtained by first extrapolating the results for a finite bulk salt concentration (0.1 M) [83] to the same macroion surface charge density as used in our calculations (le/60A ). Then, only the part close to contact with the wall has been retained, rep2  resented by a simple polynomial, and divided by two, as the full ip+(() did not yield a convergent numerical result. The hard-sphere counterions employed by Torrie et al. have a slightly larger diameter (1.08 d), but this will not influence our conclusions significantly. Note that for the salt concentration range, 0.1 M — 1.0M, investigated in [83], the short-range part of the ion-wall potential of mean force is nearly independent of p itself. The results for a single ion close to the macroion from the same study cannot be used here, as the solvent structure for this case is dominated by the unscreened field from the spherical macroion, which induces long-range order in the solvent. The curve labeled "ip2" is modelled after the case for a single ion close to an uncharged macroparticle [82]. However, both ip+(C)  a r e  better viewed as "test potentials" to study the general influence  of effective ion-wall potentials. The important common feature of both models is the (on average) added ion attraction towards the walls. The slit pressure now includes a term due to the wall potential, P n [see Eq. (2.7.3)], wa  Chapter 5.  Further  0.5 0.0 E-  Investigations  101  Level  \ . \  f  0.5 - \  m  at the McMillan-Mayer  l  ^  '  /  /  /  /  1.0 1.5  +  2.0 2.5  1  0.0  \  0.2  1  <  1  0.4  1  <  i  0.6  ,  0.8  1.0  Figure 5.8: Two model ion-wall potentials for counterions, ip , as functions of the distance to the planes of closest approach to the walls, (. For ( = (z ax ± z) > 1 d and for coions at all distances ijj = 0. For x = (,/d < 1 the curves are given by the polynomials 2.75a; - 7.75a; + 7.25x - 2.25 (ipl) and 71.bx - 145.9a; + 91.0a; - 9.3a; - 7.7x + 0.4 +  m  a  3  2  5  4  3  2  Chapter 5. Further Investigations  at the McMillan-Mayer  102  Level  and one has •Pslit = -fkin + -Pcore + Pel + -Peff + Pwall  Note that  Pwaii^O  •  (5.2.2)  if 0+(£) does not act across the midplane. This is true for our model  potentials if d n > 3 d. wa  5.2.2  Results  We present results for systems with monovalent ions at T = 298 K, e = 78.8, and a = —0.1307 e/d  2  = —0.267 C / m in Fig. 5.9. 2  Results indicated by lines denote systems  with only counterions between the walls, which correspond to bulk solutions at infinite dilution. For both the M M ion model, (a), and the P M , (b), the net pressure between the walls is reduced when the additional ion-wall potentials are considered. Apart from very small values of d n, the potential "02" yields more attractive (less repulsive) values. wa  Differences between the " ^ 1 " and ip2" cases come about because ions experiencing u  potential ipl are mostly confined to a rather thin wall layer, whereas ions interacting with the walls via potential ip2 are more spread out over a region of thickness 1 d, and the maximum of the corresponding density profiles lies at a distance of « 0 . 2 5 d from the walls. The more pronounced structure in the M M curves is due to the competing effects of the ion-ion pair potential of mean force, w (R), ++  which supports a preferred ion-ion  distance of « 1 . 4 d, and the ion-wall potential of mean force, W+(z), which tends to locate ions close to the walls. The triangles in Fig. 5.9(a) show results for a system equivalent to the M M / . ^ l model, but in equilibrium with a bulk of salt concentration 1.0 M . Overall the effect of finite bulk concentration is rather small. At large wall separations P  nei  is  more attractive than for the infinite dilution case (compare with Fig. 4.12), while at reduced d n the pressure oscillations are slightly increased in magnitude. The extremely wa  strong wall-wall attraction at small wall separations (the minimum for both ipl and ip2  Chapter 5. Further Investigations at the McMillan-Mayer Level  103  EDQ  CO  0)  e  0.00  H-  1 Pt  6  7  dwan/ d  Figure 5.9: for monovalent counterions and a wall charge a = —0.1307 e/d , with ne curves for the M M model (see Fig. 4.1) and the P M shown. The ion-wall potentials of mean force are either given by the ip in Fig. 5.8 or no additional short-range wall potential is used (ip = 0). The triangles show the results for the M M model with wall potential ipl in equilibrium with a bulk of salt concentration of 1.0 M . The size of the triangles roughly indicates the magnitude of the numerical error due to the finite grid size. 2  Chapter 5. Further Investigations at the McMillan-Mayer  is « —0.5d /k T) 3  B  Level  104  comes about as the potentials from both walls overlap and this is  therefore a feature dependent on the superposition ansatz. To further investigate the observed behaviour we compare the pressure components according to Eq. (5.2.2) for the M M models with wall potential tpl (lines) and no additional wall potential (symbols) in Fig. 5.10. We note that Pkinh/'l] <-Pkin[V  ,=  0] as ions  are drawn towards the walls and the midplane density is reduced, and P \[ipl]> P \[tp = Qi\ e  e  because the ions are more localized close to the walls and therefore the correlated charge fluctuations are less important. Pefff^l] is similar in shape to P ff['0 = O] but the minimum e  is shifted to smaller rf u, and there is an additional oscillatory tail. As long as P i i ~ 0 wa  w a  the biggest change upon introducing the short-range wall potential appears to stem from the decrease of the midplane density, which lowers the overall net pressure as well. A n analysis of the pressure components for the model with the wall potential ip2 (not shown) yields the same conclusion.  5.2.3  Summary  At the pairwise M M level solvent ordering and polarization effects due to the additional electric field from the ion distribution are neglected.  Molecular solvent effects might  therefore be more accurately described via effective ion-wall potentials evaluated at the wall charge of interest and at a finite ion concentration, instead of for the uncharged slit. However, the above discussion also shows that due to the approximations involved, the application of M M theory (truncated at the pairwise level) to the determination of the forces between charged walls is limited to a qualitative analysis. Nevertheless, calculations employing two different model ion-wall potentials show that a lower net pressure can be expected when solvent-induced ion attraction to the walls is considered. If the wall potential locates the ions in a thin wall layer (potential of mean force obtained with finite wall charge) an oscillatory structure in the pressure is observed.  Chapter  5. Further Investigations  1.0  at the McMillan-Mayer  2.0  3.0  105  Level  4.0  dwan/ d Figure 5.10: Comparison of the pressure components for two systems with monovalent counterions interacting via the M M ion-ion potential of mean force with cr = —0.1307 e/d . The lines and symbols indicate the results for effective ion-wall potential tpl and for the case with no additional short-range ion-wall potential (ip = 0), respectively. 2  Chapter 5. Further Investigations at the McMillan-Mayer Level  106  Our findings are in qualitative agreement with results published by Marcelja [54], who employed the A H N C theory to investigate how the absorption of divalent counterions onto the walls influences the pressure between two walls immersed in a mixture of 1:1 and 2:1 P M electrolyte solutions. He concluded that the wall-wall repulsion is drastically reduced (and can be attractive) if the ions are attracted to the walls. We note that solution models which do not include the solvent explicitly will in general overestimate the ion densities close to the walls. Without the dense solvent, only the ion-ion core and electrostatic repulsions can prevent the ions from overcrowding at the wall, whereas for a full solvent-ion model ions have to compete with the solvent particles for space at the walls. For the ion-wall potential ipl the contact ion densities reach values around 0.9/d . 3  This might be the reason why for even stronger ion-wall attractions a  numerical solution to the A H N C equations could not be found.  5.3  Remarks on Some Experimental Results  The objective of theoretical science is not necessarily to produce a "best fit" to experimentally obtained data and the investigations of the theoretical models presented in this thesis are interesting and important on their own. Nevertheless, the model parameters employed have been obviously chosen to mimic certain experimentally accessible systems (simple aqueous solutions at room temperature between smooth walls or macroparticles). Therefore, we give a brief overview of the experimental side and how our findings compare with the available measurements. This discussion must remain qualitative due to the model simplifications made and problems related to the reproducibility and the interpretation of the experimental data. A system that seems at least roughly comparable with our theoretical model is that of two mica sheets in aqueous electrolyte solution. The mica surfaces are molecularly  Chapter 5. Further Investigations  at the McMillan-Mayer  107  Level  smooth and bind negative ions to such an extent that wall charge densities similar to those used in our calculations ( R i l e / 4 8 A ) are achieved. W i t h increasing bulk concentration 2  0.01 M) these "wall charges" are partly compensated by absorption of the positively charged counterions, but the magnitude of the counterion absorption also depends on the wall separation. One way to measure the forces between the surfaces experimentally is with the help of the surface force apparatus (SFA), introduced by Israelachvili and Adams in 1978 [85]. The mica sheets are glued onto curved silica disks, which can be incrementally pushed together (but not incrementally pulled apart) with a mechanical spring system. Therefore, the difference in force between the slit and the outside bulk solution can be determined directly. The distance of the disks (not surfaces!) is measured via an optical interference technique and one problem with the interpretation of the experimental results is the determination of the surface separation. The results obtained are often reported as force per mean radius of curvature of the surfaces and are not directly comparable to theoretically calculated pressure results for the geometry of two flat walls . Only for 1  the case of high repulsive forces do the surfaces flatten and pressures can be estimated. Some of the problems encountered with the SFA can be circumvented by employing the more recent atomic force microscope technique [86]. Typical results from surface force measurements [87,88,53] agree well with the predictions from D L V O theory (see Appendix C) for large surface separations and/or low ( < 1 0 ~ M ) salt concentrations. For surface separations ^$20 A and medium to high salt 5  concentrations the experiments show an overall repulsive interaction between the sheets that has been termed the "hydration force" in the literature. Neither the magnitude nor Schoen et al. considered the force per radius of curvature, F(h)/R, for a fluid between a spherical macroparticle and a planar wall at a distance h. They concluded, by employing the Derjaguin approximation, that this quantity is related to the net pressure between a corresponding system of two planar walls at a separation h via -27rP (/i) = (d/dh)[F(h)/R}. 1  net  Chapter  5. Further Investigations  at the McMillan-Mayer  Level  108  the trend of increased repulsion for more concentrated solutions can be explained by P M or M M level calculations. For wall separations smaller than « 2 0 A strong oscillations also appear in the experimental pressure data. These oscillations have a periodicity of about the water diameter and can therefore be attributed to some kind of solvent layering, probably aided by the ions. In Chapter 3 we have investigated the validity of superposition schemes that combine hard-sphere and ionic contributions to the pressure for mixtures of P M ions and dense neutral hard spheres. The magnitude of the pressure oscillations for dense hard-sphere systems grows to a good approximation exponentially with decreasing dwaii,  whereas the growth observed in experiments is usually much slower; indeed, the  magnitude often seems rather constant over a wide range of separations. Trokhymchuk et al. [65] examined the interaction of two spherical macroions immersed in a solution of point ions and dipolar hard spheres. For their model and theory (MSA with the homogeneous H N C equation used to obtain the macroion-macroion potential of mean force), the authors show that solvent polarization in the vicinity of the macroions can significantly influence the results, particularly at higher salt concentrations . Furthermore, their ob2  servations are qualitatively consistent with experiment, suggesting that solvent ordering near the surface might be an important feature. Kjellander and coworkers [53,89] attempted to rationalize missing oscillations in the force between mica sheets immersed in C a C ^ solution with the help of calculations for a 2:1 P M . Their argument is that the attractive component due to correlated ion fluctuations will be present in the system and, if the attraction is sufficiently strong, it could outweigh the increase of the repulsive component as the walls are brought together. If this were true, then, due to the experimental procedure, one or more oscillations could be effectively "skipped over" and not observed (see Fig. 5.11). Hence, a missing oscillation In an earlier article.Henderson and Lozada-Cassou [60] concluded without employing a dense dipolar solvent, that simply assuming a lower dielectric constant in a layer close to the walls produces a similar effect. 2  Chapter 5. Further Investigations at the McMillan-Mayer  Level  109  might be an indication of an attractive pressure contribution at the corresponding wall separation. Similar experiments have been done for NaCl [87] and KC1 [88] solutions. In the KC1 case there are no missing oscillations but in the NaCl system the expected second to last oscillation before contact is absent. Unlike the divalent case, the missing oscillation observed for NaCl solution cannot be explained by correlated charge fluctuations because for 1:1 electrolytes this effect is not large enough. However, our M M level calculations might provide an explanation. At least for some models [52,76], the minima in the ion-ion pair potentials of mean force imply stronger attractions if N a rather than +  K  +  is the counterion. This could possibly explain the experimental observations for the  1:1 electrolyte solutions. Recently Kekicheff et al. [86] reported results for the force between mica surfaces in Ca(N03)2 solution, where the strong repulsive component at higher salt concentrations is absent. The explanation given is that the presence of halide ions causes the degradation of the mica sheets within a few hours and solvent molecules can easily get trapped between the mica sheets and the supporting disks. This, they claim, increased the repulsive pressure background in earlier measurements.  In the same article the authors find to  some degree qualitative agreement with A H N C calculations for a P M with additional counterion absorption to the walls (estimated to have a magnitude of « 2 . 5 k T). B  The  model was introduced and investigated before by Marcelja [54]. Note that an additional attraction of the counterions to the wall is consistent with "effective" ion-wall potentials due to the molecular nature of the solvent. These potentials depend on the wall charge but are relatively independent of the salt concentration [83], as discussed in the previous section.  Chapter 5. Further Investigations  at the McMillan-Mayer  Level  110  wall separation Figure 5.11: Schematic view of the pressure curves obtained in some SFA experiments. The curve is thought to result from three (additive?) contributions: a monotonic repulsive component, an oscillatory component due to packing of the solvent, and an attractive component with a minimum around the wall separation corresponding to point " E " . Only the solid lines represent values accessible for measurement. Once point " A " is reached, a further increase of the external force results in a sudden change of the wall separation corresponding to a point between " B " and " C " . However, an increase of the external force at point " C " yields a jump to a wall separation corresponding to a point between "F" and " G " and the local maximum at " E " cannot be detected experimentally.  Chapter 6  Summary and Conclusions  This thesis is concerned with exploring the effects of a molecular solvent on the forces between two infinite planar like-charged walls immersed in aqueous electrolyte solution. Investigating full models of ions and "realistic" solvent particles between walls is still beyond the computational power of present-day workstations, especially if a broad range of model parameters is to be studied. Here, the emphasis is put on comparing the results for a solution model that reduces the solvent to a dielectric continuum, called the primitive model, and models that include some (but not all) aspects of the presence of a molecular solvent. Our calculations employ numerical solution of anisotropic integral equations using the anisotropic hypernetted-chain approximation, which is known to give reliable results for such model systems. These calculations yield the average particle densities and the particle pair correlation functions for the fluid between the walls. With this information the pressure in the slit (force on the wall per unit surface area), P ij , can be s  t  obtained either from the contact theorem or via an equivalent "midplane formula". The latter route divides P n into several pressure components, representing particle interacs  t  tions on both sides of the midplane, and therefore facilitates the physical interpretation of the results. The fluid between the walls is assumed to be in equilibrium with a bulk system, which fixes the chemical potential of the individual species. In the limit of an infinitely diluted bulk we would expect no coions (to the wall charge) to be drawn into the cavity. The number of counterions (per unit surface area) is therefore dictated by applying overall  111  Chapter  6. Summary and  112  Conclusions  electroneutrality in the slit. This limiting case is easier to handle numerically (as only one ionic species is present), and it is a good approximation even for higher bulk salt concentrations if high surface charges and small wall separations are considered. The net pressure between the walls, P outside bulk pressure. P  n e t  n e t  , is given by the difference between P  s l i t  and the  is a quantity of major interest, as it indicates whether the  walls (macroions) show a net attraction ( P  net  < 0) or net repulsion ( P  n e t  > 0), which is  important for understanding the behaviour of colloidal systems and their technological applications. A simple tractable model useful in the investigation of the influence of the finite size of the solvent molecules is obtained by representing the solvent particles as neutral hard spheres. The ions interact via the usual Coulombic law, reduced'by a separation independent dielectric constant (as in the P M ) . In a mixture of neutral hard spheres and hard-sphere ions between charged walls the ions have to compete for space with the dense neutral component and hard-core interactions are important, in contrast with the P M case. In general, the ions are pushed further towards the walls and, especially for monovalent counterions, the ion-solvent hard-core interactions induce oscillations in the ion-ion correlations. At wall separations of a few solvent diameters, it is shown that even at relatively high surface charge and moderate solvent density, the ionic contribution to the pressure tends to be dominated by the hard-core or packing component. If the ions and solvent particles are of equal size, then the net pressure between the walls can be reasonably well approximated by adding the pressures of pure one-component ionic and solvent systems (see Fig 3.5). However, if the ion and solvent diameters are significantly different the pressure curve is more complex, and simple superposition of the ionic and solvent pressures is no longer sufficient. For this case, to a good approximation, it is still possible to divide the pressure into electrostatic and hard-core components, but now the appropriate hard-core system must itself be a mixture of neutral hard spheres (see  Chapter 6. Summary and Conclusions  113  Fig 3.9). Note that these superposition approximations yield reasonable results since the pressure component due to correlated charge fluctuations across the midplane, P i , e  changes only slightly upon addition of the neutral component to the P M . Apart from the finite particle size, the presence of a molecular solvent alters the ion potentials of mean force, which represent the "effective" ion interaction energy after averaging over the positions and orientations of the solvent molecules. For a continuum solvent model the ion pair potential of mean force (at infinite dilution), w (R), a7  is given  by adding the usual Coulombic and hard-core contributions. In a molecular solvent both the preferential orientation of the molecules around the ions and the lack of dielectric screening at short ion separations, P , contribute to deviations from the continuum case. In the formally exact McMillan-Mayer theory it is shown that the partition function for the full solution model (ions and solvent molecules) is given by the product of the partition functions for a pure solvent system and for a purely ionic system, where the ions interact via the potentials of mean force. The latter system can be viewed as a straightforward extension of the P M and the focus of much of this thesis is to obtain the ionic fluid structure and to evaluate the resulting pressure between the walls for such M M models. Although exact if n-particle potentials of mean force to all orders are included, for practical computations M M theory is usually truncated at the pairwise level. This is a good approximation only if the solvent structure in the vicinity of an ion pair is not much affected by the presence of other ions, e.g., for bulk solutions the salt concentration must not be too high (where "not too high" can mean up to several moles/litre). The situation for a solution between two walls charged with a high surface charge density is more complicated.  For such a system, the counterions form rather dense surface  layers, and the validity of the truncation at the pairwise level is not clear, but this is not the only (and perhaps not the most serious) approximation made. For the M M slit  Chapter 6. Summary and  114  Conclusions  calculations we employ the ion pair potentials of mean force for an ion pair immersed in a bulk solvent, instead of evaluating the potentials for a corresponding inhomogeneous (uncharged) solvent system, thereby ignoring the influence of solvent order near the walls on the potentials of mean force. The severity of the latter approximation is hard to assess, but some general consequences can be predicted, such as the reduction of the dielectric constant close to the walls due to the solvent order. In addition, the application of M M theory to the inhomogeneous slit system requires the introduction of "effective" ionwall potentials of mean force. These ion-wall potentials describe the solvent-mediated interaction of the ions with the walls. It is argued that the use of effective ion-wall potentials obtained from calculations for a single charged macroion in a model solution explicitly including the solvent molecules at a finite bulk salt concentration might cancel some of the errors introduced by the pairwise approximation. This suggestion is based on the fact that for this choice the solvent structure close to the macroion is similar to that in the charged slit system. We first considered a model solution consisting of hard-sphere ions and hard-sphere water molecules, where the solvent particles are decorated with appropriate dipole and quadrupole moments. For this model the ion-ion w (R) aj  were calculated before by Kusa-  lik and Patey. Initially the solvent-mediated wall potentials are ignored in order to investigate the effect of the pair potentials of mean force alone. It is found that for such systems the M M results can differ significantly from those of the P M . A new attractive component to the pressure, P s, arises, which is due to solvent effects on the effective e  counterion-counterion interaction and not to the correlated charge fluctuations that give rise to the overall wall-wall attraction found for divalent counterions in the P M case. The influence of P g can be very important, leading to an attractive net wall-wall force e  at small separations (;$ 15 A), even for monovalent counterions if the wall charge density is sufficiently high (see Fig. 4.5). With rising salt concentration in the bulk the  Chapter 6. Summary and  115  Conclusions  attraction increases in magnitude and in range (see Fig. 4.12). It is promising that calculations performed by Kinoshita et al. at the homogeneous R H N C level for a model of spherical macroions in electrolyte solutions with solvent particles explicitly included show qualitative agreement with our results. The homogeneous R H N C theory is easier to solve numerically than the A H N C approach, especially for the hard-core models employed here at finite salt concentrations. These systems pose severe numerical problems at the A H N C level due to the strong anion-cation attraction at ion-ion contact. It is found that this coion-counterion contact attraction causes striking changes in the fluid structure around the coions compared to the P M . Calculations were also performed employing "realistic" models for N a and C l in +  -  water, for which ion pair potentials of mean force are available in the literature. For some models attractive wall-wall forces are observed, similar to those obtained from the solution model discussed above. On the other hand, the M M theory is found to be rather sensitive to details of the counterion-counterion potential of mean force, and different models for the same counterion can give qualitatively different results (see Fig. 5.7). Therefore, if the M M approach is to enhance our understanding of real systems, physically accurate ion-ion potentials of mean force are required. Unfortunately, at present one has little means of evaluating this aspect of the various models proposed in the literature. The influence of the effective ion-wall potentials on the pressure between the walls was studied for the electrolyte solution model of Kusalik and Patey. For this model solution and a spherical macroion Torrie et al. have reported ion-wall potentials of mean force. Guided by their results we suggest two model ion-wall potentials for both a neutral and a highly charged macroparticle. The common feature of the two potentials is an additional counterion attraction towards the walls. Regardless of the ion-wall potential employed, one finds an increased wall-wall attraction for monovalent counterions (see Fig. 5.9). Therefore, it appears that the general trend of a short-range attractive force  Chapter 6. Summary and  116  Conclusions  between the walls, obtained for the M M models with and without additional ion-wall potentials, is not due to the pairwise level approximation made in the inhomogeneous M M theory. Again, the results are qualitatively different from those obtained for the 1:1 P M electrolyte solution case. These effects should be considered when experimental results for corresponding systems (e.g., mica sheets in salt solution) are to be interpreted. However, from the experimental data presently available it is difficult to determine if the attractions described here are important for "real" systems.  More reliable (i.e.,  reproducible) measurements for monovalent counterions are required, especially for small wall separations. On the theoretical side further investigations of the approximations involved in the inhomogeneous M M theory are necessary. This could be done for simple systems, such as mixtures of neutral particles and ions, with the integral equation theory employed here. More interesting would be a detailed analysis of the M M level with the help of simulation results for model solutions explicitly including "realistic" water models, once this kind of calculation becomes computationally feasible for inhomogeneous systems.  Bibliography  [1] J . N . Israelachvili. Intermolecular  & Surface Forces. Academic Press, 1991.  [2] J . Israelachvili and H . Wennerstrom. Nature, 379:219, 1996. [3] P. Attard. Adv. Chem. Phys., 92:1, 1996.  [4] R. Kjellander. Ber. Bunsenges.  Phys. Chem., 100:894, 1996.  [5] F . Oosawa. Poly electrolytes. Marcel Dekker, New York, 1971. [6] G . N . Patey. J. Chem. Phys., 72:5763, 1980. [7] P. Linse and V . Lobaskin. Phys. Rev. Lett., 83:4208, 1999. [8] L . Belloni and O. Spalla. J. Chem. Phys., 107:465, 1997. [9] D . G . Grier. J. Phys.:Condens.  Matter, 12:A85, 2000.  [10] J . Wood and R. Sharma. Langmuir, 11:4797, 1995. [11] A . Carambassis, L . C. Jonker, P. Attard, and M . W . Rutland.  Phys. Rev. Lett.,  80:5357, 1998. [12] K . Lum, D . Chandler, and J . D. Weeks. J. Phys. Chem. B, 103:4570, 1999. [13] D . R. Berard, P. Attard, and G . Patey. J. Chem. Phys., 98:7236, 1993. [14] J . C . Shelley and G . N . Patey. J. Chem. Phys., 100:8265, 1994. [15] J . C . Shelley and G . N . Patey. J. Chem. Phys., 102:7656, 1995.  [16] J . A . Barker and D. Henderson. Rev. Mod. Phys., 48:587, 1976. [17] J.-P. Hansen and I. R. McDonald.  Theory of Simple Liquids.  Academic Press,  second edition, 1991. [18] M . P. Allen and D. J . Tildesly. Computer Simulation  of Liquids. Claredon Press,  1989. [19] T. L . Hill. Statistical Mechanics. McGraw-Hill, 1956. [20] L . S. Ornstein and F . Zernicke. Proc. Akad. Sci. (Amsterdam),  117  17:793, 1914.  118  Bibliography  [21] H . L . Frisch and J . L . Lebowitz, editors. Fluids. W . A . Benjamin, Inc., 1964.  The Equilibrium  Theory of Classical  [22 T. Morita and K . Hiroike. Prog. Theor. Phys., 25:537, 1961. [23;  P. Attard and G . N . Patey. J. Chem. Phys., 92:4970, 1990.  [24;  J. S. Rowlinson and B . Widom. Molecular Theory of Capillarity. 1989.  Clarendon Press,  [25 H. Greberg and R. Kjellander. Mol. Phys., 83:789, 1994. [26; R. Kjellander and S. Sarman. J. Chem. Phys., 90:2768, 1989. [27 J. C. Rasaiah, D . N . Card, and J . P. Valleau. J. Chem. Phys., 56:248, 1972. [28 M . Plischke and D. Henderson. J. Chem. Phys., 88:2712, 1988. [29 R. Kjellander. J. Chem. Phys., 88:7129, 1988. Erratum - 89:7649.  [30 M . Vossen and F . Forstmann. Mol. Phys., 86:1493, 1995. [31 M . Plischke and D. Henderson. Proc. R. Soc. Lond. A, 404:323, 1986. [32 R. Kjellander and S. Sarman. Mol. Phys., 70:215, 1990. [33 R. Kjellander and S. Sarman. Chem. Phys. Lett, 149:102, 1988. [34 R. Kjellander and S. Marcelja. J. Chem. Phys., 82:2122, 1985. [35; R. Kjellander and S. Marcelja. 114:124.  Chem. Phys. Lett., 112:49, 1984. Erratum -  [36 M . Lozada-Cassou, W . Olivares, and B . Sulbaran. Phys. Rev. E, 53:522, 1996. [37; D. Henderson, P. Bryk, S. Sokolowski, and D . T. Wasan. Phys. Rev. E, 61:3896, 2000. [38 M . Lee and K . - Y . Chan. Chem. Phys. Lett:, 275:56, 1997. [39 R. Kjellander and S. Marcelja. Chem. Phys. Lett, 127:402, 1986. [40 B. R. Svensson and C. E . Woodward. Mol. Phys., 64:247, 1988. [41 J. P. Valleau, R. Ivkov, and G. M . Torrie. J. Chem. Phys., 95:520, 1991. [42; R. Kjellander, T. Akesson, B . Jonsson, and S. Marcelja. J. Chem. Phys., 97:1424, 1992.  Bibliography  [43] H . Greberg, R. Kjellander,  119  and T. Akesson. Mol. Phys.,  87:407, 1996.  [44] J . S. Rowlinson. J. Chem. Thermodynamics, 25:449, 1993. [45] D . Henderson and L . Blum. J. Chem. Phys., 69:5441, 1978. [46] D. Henderson, L . Blum, and J . L. Lebowitz. J. Electroanal. Chem., 102:315, 1979. [47] S. L . Carnie and D. Y . C. Chan. J. Chem. Phys., 74:1293, 1981. [48] D . Henderson and L . Blum. J. Chem. Phys., 75:2025, 1981. [49] I. Rouzina and V . A . Bloomfield. J. Phys. Chem., 100:9977, 1996. [50] S. Marcelja. Colloids Surfaces A, 129-130:321, 1997. [51] S. Marcelja, 1999. Private communication. [52] P. G . Kusalik and G . N . Patey. J. Chem. Phys., 88:7715, 1988. [53] R. Kjellander, S. Marcelja, R. M . Pashley, and J . P. Quirk. 92:6489, 1988.  J. Phys. Chem.,  [54] S. Marcelja. Biophys. J., 61:1117, 1992. [55] L . Zhang, H . T. Davis, and H . S. White. J. Chem. Phys., 98:5793, 1993. [56] Z. Tang, L . E . Scriven, and H . T. Davis. J. Chem. Phys., 97:494, 1992. [57] Z. Tang, L . E . Scriven, and H . T. Davis. J. Chem. Phys., 100:4527, 1994. [58] C. N . Patra and S. K . Ghosh. Phys. Rev. E, 49:2826, 1994. [59] R. Kjellander and S. Sarman. J. Chem. Soc. Faraday Trans., 87:1869, 1991. [60] D . Henderson and M . Lozada-Cassou. J. Colloid Interface Sci., 114:180, 1986. [61] D . J . Adams. Mol. Phys., 28:1241, 1974. [62] G . M . Torrie, P. G . Kusalik, and G . N . Patey. J. Chem. Phys., 88:7826, 1988. [63] J . Rescic, V . Vlachy, L . B . Bhuiyan, and C. W . Outhaite. 107:3611, 1997.  J. Chem. Phys.,  [64] R. M . Shroll and D . E . Smith. J. Chem. Phys., 111:9025, 1999. [65] A . Trokhymchuk, D . Henderson, and D . T. Wasan. 210:320, 1999.  J. Colloid Interface Sci.,  120  Bibliography  [66] M . Kinoshita, S. Iba, and M . Harada. J. Chem. Phys., 105:2487, 1996. [67] G . N . Patey. Ber. Bunsenges. Phys. Chem., 100:885, 1996. [68] W . G . McMillan and J. E . Mayer. J. Chem. Phys., 13:276, 1945. [69] S. Marcelja, 1999. Preprint. [70] S. Marcelja. Nature, 385:689, 1997. [71] P. G . Kusalik and G . N . Patey. J. Chem. Phys., 89:7478, 1988. [72] P. G . Kusalik and G. N . Patey. J. Chem. Phys., 89:5843, 1988. [73] P. G . Kusalik and G . N . Patey. J. Chem. Phys., 79:4468, 1983. [74] E . Guardia, R. Rey, and J . A . Padro. J. Chem. Phys., 95:2823, 1991. [75] E . Guardia, R. Rey, and J . A . Padro. Chem. Phys., 155:187, 1991. [76] B . M . Pettitt and P. J . Rossky. J. Chem. Phys., 84:5836, 1986.  ,  [77] K . Toukan and A . Rahman. Phys. Rev. B, 31:2643, 1985. [78] L. X . Dang and B . M . Pettitt. J. Phys. Chem., 94:4303, 1990. [79] W . L . Jorgensen. J. Am. Chem. Soc., 103:335, 1981. [80] E . C. Zhong and H . L. Friedman. J. Phys. Chem., 92:1685, 1988. [81] E . Di'az-Herrera and F . Forstmann. J. Chem. Phys., 102:9005, 1995. [82] G . M . Torrie, P. G . Kusalik, and G . N . Patey. J. Chem. Phys., 89:3285, 1988. [83] G . M . Torrie, P. G . Kusalik, and G . N . Patey. J. Chem. Phys., 90:4513, 1989. [84] G . M . Torrie, P. G . Kusalik, and G . N . Patey. J. Chem. Phys., 91:6367, 1989. [85] J . Israelachvili and G. E . Adams. J. Chem. Soc. Faraday Trans. 1, 74:975, 1978. [86] P. Kekicheff, S. Marcelja, T. J . Senden, and V . E . Shubin. J. Chem. Phys., 99:6098, 1993. [87] R. M . Pashley and J . N . Israelachvili. J. Colloid Interface Sci., 101:511, 1984. [88] P. M . McGuiggan and R. M . Pashley. J. Phys. Chem., 92:1235, 1988.  Bibliography  121  [89] R. Kjellander, S. Marcelja, R. M . Pashley, and J . P. Quirk. 92:4399, 1990.  J. Chem. Phys.,  [90 F. Lado. J. Comp. Phys., 8:417, 1971. [91 M . J . Gillan. Mol. Phys., 38:1781, 1979. [92 K . - C . Ng. J. Chem. Phys., 61:2680, 1974. [93 M . Kinoshita and D. R. Berard. J. Comp. Phys., 124:230, 1996. [94 W . H . Press, S. A . Teukolsky, W . T. Vetterling, and B . P. Flannery. Recipes in FORTRAN. Cambridge, 1992. [95 E. J-. W . Verwey and J . T. G . Overbeek. Colloids. Elsevier, 1948.  Numerical  Theory of the Stability of Lyophobic  [96 S. Engstrom and H . Wennerstrom. J. Phys. Chem., 82:2711, 1978. [97; J. C. Neu. Phys. Rev. Lett, 82:1072, 1999. [98 B. V . Derjaguin and L . Landau. Acta Phys. Chim. USSR, 14:633, 1941. [99; R. Kjellander and D. J . Mitchell. Chem. Phys. Lett, 200:76, 1992. [100 B. Berne, editor. Statistical Mechanics A. Plenum Press, New York, 1977.  Appendix A Numerical Algorithm  We first present the general method of numerical solution of the A H N C equations used in our program. Additional numerical complications for M M systems with effective ion-ion potentials based on hard-sphere ion models (i.e., results from Kusalik and Patey [52]) at finite bulk concentrations are discussed thereafter.  A.l  General Method  For many details, including the notation, we follow Kjellander [29]. For all species the space between the walls is partitioned into parallel layers, giving an overall number of N  layers with widths AZJ.  z  i,j = l , . . . , i V  2  are now indices for the individual layers,  representing the z-value in the middle of the layer, z and the particle type ojj. The wall i:  layer for each species is particularly thin ( A Z J « 0.005 A) and layers close to the wall are usually chosen to be slightly thinner than those near the midplane. Parallel to the wall we use N grid points, r^, with a cutoff r^ = r r  ~P i){y  r  u  r2  + (i ~ z  , beyond which we set hij(r) = 0 and Cy(r) =  the latter choice equating c^- with its asymptotic value. These  j) )i  z  m a x  2  long-range tails are handled analytically throughout the calculation. Typical parameter values used in the numerical calculations are N = 60 — 160, N « 300, and r z  r  m a x  « 18 A,  which are chosen to ensure that our results are independent of these values. Due to the symmetry of the model (i.e., two identical walls) and the inherent symmetry of the correlation functions [g (r, z^, Zj)= g {r, Zj, zA] one can reduce the required memory aiaj  ai0li  considerably (by up to 75%), but a complicated system of indices must be introduced. 122  Appendix  A.  Numerical  123  Algorithm  The OZ equation is conveniently solved in Fourier space. The Fourier transformation in two dimensions, called Hankel transformation functions [f(r)  (HT), for spherically symmetric  = /(|r|)] becomes poo  f{k) = 2?r / drr f(r)J {kr)  (A.l.l)  0  Jo  f(r) = ±- rdkkf(k)J (kr), 0  Z7T JO  where Jo(x) is the zeroth order Bessel function. transformed OZ-equation reads h (k,z z ) on  u  =  2  c (k,zi,z ) + a7  dz  2  <5  3  J  Applied to Eq. (2.6.1) the Hankel  c {k, aS  zi,  z )p (z )h {k, z , 3  s  3  Sl  3  z) , 2  (A.l.2)  which involves now only a one-dimensional integration, instead of the three-dimensional integral in Eq. (2.6.1). Thus, we have to perform two Hankel transforms in every iteration loop for which no fast "FFT-like" algorithm applicable to our problem is known. We employ the orthogonality-preserving Hankel transform quadrature proposed by Lado [90], which requires the grid points to be placed according to  = Ci^max/Civ, where Q is the r  i-th positive root of Jo(x). The corresponding grid points in Fourier space are ki = C /r . li  with a cutoff A ;  max  m&y  = CAT,./Vmax- The formulas for the discrete transformations are given by =  TJ-  E  "m ' ax =  fciW*  ^d  (A.1.3)  j=l  Z^-' LHki)W 1  ii  with  (A.1.4)  'max j=l W  _  JQ(CJ(J/(N ) T  where the matrix Wij has to be calculated only once at the beginning of the program. Using the trapezoidal rule to perform the integration, Eq. (A.1.2) can be written in the matrix form  H(k)  = C(k) +  H{k)RC(k) ,  (A.1  Appendix A. Numerical Algorithm  124  where H(k) = {hij(k)}, C(k) = {dij(k)}, and K= {piAzi 5ij} (<% is the Kronecker delta). Solving Eq. (A.1.5) for Ff yields H = C ( l - R C ) " = (R - R C R ) - - R 1  1  1  ,  (A.1.6)  where l = {5ij}. At hard-core contact r = a^ =  \J\{d ^) + d ^)  2  a  a  —  (zi  — Zj)  2  the functions Cij(r) and  hij(r) including their derivatives c^-(r) = dcij/'dr and ti^r) are discontinuous. Therefore the corresponding functions in Fourier space have long-range tails. This can be avoided by adding a second order polynomial of the form Qy(r) = Ac  tj  + ^ ^ ^ j  (A.1.7)  to the function in real space for r < a (inside the core) and subtracting the Hankel transformed part Qij(k) = 2TTa j[AcijJ (a jk)/k i  1  i  - Ac J {a jk)/k ) 2  ij  2  i  (A.1.8)  after performing the numerical H T , where J (x) is the n-th order Bessel function. Acij = n  Ahij and Ac'^ = Ah'^ can be determined by using hij(r < a) = — 1 and interpolating the continuous (even at hard-core contact) function yij( ) = hij{r) - c^r) r  ,  (A. 1.9)  as usually no grid point lies at r = a. With the H N C closure, hij(r) = exp{—ftuij(r)  +  Uij(r)} — 1, one finds Ac  i:j  = c°f(aij)  - <%-(aij) = expi-Buijtdij)  + yy(ay)} ,  and Ac'^ is given by a simple numerical derivative of Eq. (A. 1.10).  (A.l.10)  Appendix  A. Numerical  125  Algorithm  The actual iteration between the OZ-equation and the H N C closure is performed with the short-ranged functions  4 = v + Pv-ij ' c  a  n  Vij = VU ~ Pufj •  d  (A. 1.11)  ufj and its Hankel transform are given by  fig(r)  , and  =  4ire e^r + ^  {  K  )  6ij  where &  m m  =  =  M  M  (A.1.12)  bj  2  2  0  ^M-Kk}  2  max (b ,^(zi  w  .  t  h  - Zj) ^ , 2  min  > 0 has been introduced to avoid divergences for small r.  The terms in the density equation, Eq. (2.6.3), involving  and the integral over the  asymptotic part of c^-, —Qufp have to be evaluated together to yield the finite result [34] - / ^ ( * I ) - 2 T T £  jdz (z ) 2Pl  2  Jdrrpui^  + iz.-z,) )  ,  2  Cfdwall + £  2ee  <?  7  dz \z 2  1  (A.1.14) Z \p {z ) 2  7  2  0  When using the trapezoidal rule to evaluate the integral in Eq. (2.6.3) at least one more grid point at hard-core contact must be introduced to deal with the discontinuity of the functions Cy and hij. The required value C{j(r — 0) is determined by a third-order polynomial extrapolation. The general iteration scheme is depicted in Fig. A . l . If no good initial guesses from a previous calculation are available we use Poisson-Boltzmann density profiles (see Appendix C) and Cij(r) = —Bufj(r) outside the hard-core. For the low-density Coulombic systems considered in this thesis the initial guess for the correlation functions is not an important issue. Instead of using the output of the n-th iteration for a variable x, x ° , u t  Appendix  A. Numerical  126  Algorithm  if  As.< 8 HT  OZ-equation A  A  c(r) — y(r) initial  y  HT  X  density equation electroneutr. condition  guess c(r),p  HNC closure y(r)  c(r)  newp  Figure A . l : Iteration scheme for the algorithm.  if Ap <&2  Appendix  A.  Numerical  127  Algorithm  as input for the n + l-th iteration,  we employ the Picard mixing scheme [91] that  gives the new estimate as < i  = £<  u t  + (1 - O C - i ,  (A.1.15)  with the mixing parameter £, 0 < £ < 1- £ has to be adjusted to yield optimal convergence. For low-density systems £ « 0 . 5 — 0.8 turns out to be an adequate choice for both density profiles and correlation functions, whereas the high-density hard-sphere systems required a significantly lower £ ( « 0.1). Implementation of the straightforward procedure suggested by Ng [92] might improve the convergence drastically in the latter case.  The more  sophisticated Newton-Raphson technique [91,93] cannot be applied at the inhomogeneous level, as the calculation of the Jacobian becomes very time consuming. Defining Ax — :r°  —x^i,  ut  10  - 3  the system of equations is iterated until for successive iterations  Ay/y<ei^  and A p / p < £ ~ 1 0 . - 4  2  Many subroutines of the numerical program have been adapted from routines given by Press et al. [94], for example: calculation of an inverse matrix via L U decomposition, polynomial interpolation, Bessel functions ( J , J\, J ) , root finding with Ridders' method, 0  2  and solution of quadratic/cubic equations. In developing the algorithm one often has to trade off time efficiency with memory requirements, for instance, the values J i (<%&:;) and Jz(aijki)  (i,j = l,...,N  z  and I = 1,N ), r  needed in Eq. (A.1.8), are not stored in  three-dimensional matrices but are rather repeatedly calculated via linear interpolation from one-dimensional arrays. We tested the accuracy of our program by recalculating results for a 2:1 P M electrolyte between charged plates published by Kjellander et al. [42], with contact densities and midplane densities given to three significant figures. The two sets of results deviated by less than 1%.  Appendix A. Numerical Algorithm  A.2  128  M M Systems w i t h S t r o n g A n i o n - C a t i o n A t t r a c t i o n  The systems with anions as well as cations between the walls interacting via the effective ion-ion potentials from Kusalik and Patey [52] (discussed in Section 4.4.2) require extra care. Due to the strong attraction between oppositely charged ions, h+- and c+_ are very narrow and "spiked" close to ion-ion contact and the integrals in Eq. (2.6.3) must be evaluated using a very fine grid near contact, as this region gives the biggest contribution to p *(z)- Moreover, the spiked nature of these functions in real space expends the "tails" e  in Fourier space and a much finer grid in each layer becomes essential. This is because in the orthogonality-preserving H T a smaller grid width in real space corresponds to points at longer wavelengths in Fourier space. In the present calculations N = 1600 was found r  to be necessary for r  m  a  x  = 19.6A.  The like-ion functions y* {k, z\, z ) decay as C (zi, z )/k . This can easily been seen z  aa  2  a  2  by expanding a slightly modified Eq. (A. 1.6) using 1/(1 — x)zxl + x + x + ... 2  Y = C ( 1 - R C ) - - C « C R C + ... ,  (A.2.1)  1  where Y(k) — {yij(k)}. As a simple example we consider a calculation with only one layer for each ionic species, then c  V— y-+ . y+-  c-+  y++ .  ' P-  C++ . .  o" P+ .  0  p-c __ + p+c_+c+_ 2  _ p_c _c__ + p c +  +  + +  c _ +  c  c-+  .+a  p-c- c—+ +  (A.2.2)  C++ _  p c_ c +  +  p_c+_c_+ + p+c  ++  2 ++  For large k the only surviving contribution to c _ = c _ is the term proportional to +  Ji(a -k)/k +  zero (as A c  +  in Eq. (A. 1.8), which decays as fc / , whereas C++ and c -3 2  + +  and A c  long-range tails of y  t  are essentially  are very small). From this and Eq. (A.2.2) it follows that the and y  decay as k~ and that the corresponding functions for 3  ++  Appendix A. Numerical Algorithm  129  unlike species are of much shorter range. These tails are significant only if z\ ~ z , i.e., 2  for the functions yu(k). The proportionality constants Ci are determined for each layer, and corrections that take account of the contribution from points beyond the cutoff in Fourier space (A.2.3) are included in the H T . The required integrals  dkJ (kri)/k are evaluated once for 2  0  all rj. Finally, it is essential to use an identical set of layers for all ionic species even though coions are usually present at much lower concentration than counterions. This ensures that the density profiles are largely independent of the layering, although details of the correlation functions might still exhibit some artifacts due to the numerical method. In the present calculations a maximum of 101 layers per species (N = 202) was used, z  but these calculations can take several days to converge on modern workstations. We note, that most of the numerical problems can be avoided if "soft" rather than "hard" short-range interactions are used in the models.  Appendix B Formulas for  Three different routes to evaluate  P i s  i t  P H S  t  have been used in our numerical calculations:  the contact theorem, the midplane formula, and the derivative of the grand potential with respect to the wall separation. The latter has merely been employed for comparison and here we only state the formulas without proof. The midplane formula is an important tool in the interpretation of the pressure results and therefore we give a brief derivation, as Kjellander and Marcelja have not explicitly published these steps [35,39].  B.l  Midplane Formulation  In order to simplify the notation, we assume here that only one species is present between the walls, i.e., we derive the formulas for the case with only counterions (bulk at infinite dilution). We begin with the contact theorem a --  2  Psiit = k Tp(-z ) B  m&x  r^max  /  ZCCQ  Using the relationship J^_ dz(dp/dz)  J-z  sh>I  dzp(zf  1  ' .  (B.1.1)  OZ  m & x  = p(0) — p(—2 )  ZmAx  dv (z)  and a simplified version of the  max  B G Y equation, Eq. (2.3.4), k s T ^  =- P W ^  "  / * r  9  , (B.1.2)  ( r , * , * ) * 4 ^  we obtain  zeeo n  +2TT  f° /  ,  J-z  , . /"  m & x  2max  dz p(zi) l  oz J-z , . , r°° , , dz p(z ) drrg(r,z ,z )  m s u x  2  130  2  l  2  Oz .du(r,zi,z ) 2  -z  .  Appendix B. Formulas for P H S  131  t  The first term on the RHS is called Fkin, whereas terms three and four give P aii- As W  (dv /dz) = (dv el  BS  /dz) = 0 (for the integration range) one finds  Pwaii =  dlv^' (z) dzp(z)-+ 1  QyShJf^ Q dzp(z)—7^+/  /-.max  f  -/  /•»..,  dif^'iz)  r«  .  +  V^'" (z)}  -»  K  . .dv' -"(z) h  *™-ar Lj**' —5r -  -  )  L+  '  L  (B  Both terms in the expression for P u give an identical contribution if v ' (z) =v ' sh  T  L4)  (—z)  sh n  wa  (as in our model). The last term on the RHS in Eq. (B.1.3) can be simplified by splitting the integral over z into two terms, A and B, yielding 2  A +B  =  2TT f  dz ( )  +  c f ° 2?r /  lP  •/-Zmax  f ™dz p{z )  rdrrg(r,z z ) '  Z  Zl  2  JO  J-Zm&x  J - Z  m  u  2  JO  i i Nf ° dzip{zi)  , / \ f°° , / dz p{z ) drrg(r,z z ) 2  &  2  u  <OZ\  sdu(r,zi,z ) 2  2  JO  x  (B.1.5)  du{r Zl,Z2)  2  .  OZ\  B is zero for symmetry reasons, as can easily be shown: The assignment of the labels "1" and "2" in the expressions A and B is arbitrary and an exchange of coordinates, z\ < 4 z , 2  should not alter the value of the terms. For the same reason the pair distribution functions are also symmetric, g(r, z\, z ) = g(r, z , z\) . The particle pair potential depends only on 1  2  2  the separation R = \Jr + (z\ — z ) , which means that 2  2  2  du du dR d~z[ ~ d~R ~dz[  du z\ — z d~R R  2  =  =  du ~d~z ~ '  .  2  . (••)  After exchanging z\ and z in (B.1.5) the term B can be rearranged to the original form, 2  but multiplied by an additional factor (—1) stemming from (B.1.6), i.e., the term is overall antisymmetric and has to be zero. In term A, however, the integration boundaries are different for z\ and z and the additional symmetry operation Z\ —> —zi and z —»• — z 2  2  must be performed to bring the term back to its original form (due to the fact that *If more than one species is present the indices have to be exchanged, as well: ga-y{r,zi,z2) = 9ya{r, Z ,Z ). 2  1  2  Appendix B. Formulas for  P i s  132  i t  z = 0 is a symmetry plane for our model one has p(z) = p(—z) and g(r, —zi, —z ) = 2  g(r,zi,z )). 2  Performing only the latter transformation on term A and splitting up u into  its components, Eq. (B.1.3) becomes CT  — -Pkin + -Pvall ~  -Psiit  2  (B.1.7)  z e Co fZmai rv rzmax . . r° -1-E \ dz p(zi) / 1  -Zmax J— Z  JO  m  a  dz p(z ) 2  2  x  roo roo  /  d[u  RS  drrg(r,zi,z )  +u + u ] el  2  •'U •'O  sh  8Z\  P h is given directly by the part of the integral involving u . Combining the electrostatic sh  S  terms defines P i and yields e  r2  P„i  =  "^ f° ^ ^ Jo — — J-z^'^^Jo  f ° ° „  p( ) f° dz2p(z2) J-z .* Jo  rdrrh{r,z z )  rZm  2ee  0  r d &X  -2TT  Zl  The  eicl  ,. f°° r°°  2  ^  7 o  2 e  °  2  .  (B.1.8)  rdr  - 2 2 ) 7  e  2  OZ\  ma  last identity follows from h = g — 1, roo , f°° du Cw a cf ,. drr—— = - - — ( 2 1 2TT /  '  dueX{r ZuZ2)  u  Zl  JO  ,du*\r, ,z ) dz, Zl  ~"  7  o  ^ / r + (21 - 2 ) 2  =-3  2  2  q zi - z 2e e | 2 i - 2 | 2  (B.1.9)  2  0  2  and f0  /"Zmax  q  dzp(z) = q  dz p{z) = o ,  (B.1.10)  J—z .x  JO  mB  where the symmetry of the system and the electroneutrality condition have been used. The  remaining term in Eq. (B.1.7) is the pressure contribution due to hard-core  interactions across the midplane, P re- To simplify the expression for P CO  core  we first  define f = \Jd — (z\ — z ) (r-value for which the particles in layers 21 and 2 are in 2  2  2  2  contact, R(f) = d) and one further finds S(R-d)  = ^  e  - ^  OR dr  S  = -(3e-^ -^  (B.l.ll)  S( - r) = p(r - r) .  (B.1.12)  sd  -1 r  Appendix  B. Formulas for P n s  133  t  The equality (B.l.ll) follows from the fact that the Dirac delta function, S(x — XQ), is the derivative of the Heaviside (or step) function and Eq. (B.1.12) is given by a transformation property of 5 (f(x) — f(x )).  Using the last two equations we obtain  0  / drrg{r,z ,z )—— = rg{r,z ,z )^—^(-k T)-, Jo oz\ a 1  2  1  2  B  B.1.13  r  and finally P  coie  = 2Trk T  dzip(zi)  B  ^0  dz p(z )(zi-z )g(f,z ,z ) 2  2  2  l  .  2  (B.1.14)  J—Zm&x  Summing up all contributions yields Eq. (2.7.3) .  B.2  Calculation of £l  A  The grand potential per unit surface area and the corresponding Helmholtz free energy, F , for an electrolyte solution are related via A  Q = F -p±jdz A  A  \p (z) + p-(z)} .  (B.2.1)  +  The last term on the RHS can easily be written in the form p± fdz[p+(z)+  p-{z)}=  Pa fdzp (z) a  + k T—lnX,  (B.2.2)  B  where the definitions from Eq. (2.6.6) and the electroneutrality condition have been used, and a constant independent of  o?n has been dropped. X w a  depends on  d n and is wa  determined by solving Eq. (2.6.7). F has been given, for example, by Kjellander and Marcelja [34], exploiting again the A  HNC approximation and thereby avoiding a numerical coupling parameter integration . 2  The excess part of F is obtained via integration of n *(z) over the whole slit system, and then it is added to the appropriate ideal part. 2  A  e  Appendix B. Formulas for P n s  134  t  After simplification of their formulas one finds F  A  (B.2.3)  2ee k T 1 ^ Q,  kT B  0  B  dz dz p (z )p (z )\z x  dz dz x  -— Z  2  2  a  x  pa{zi)p (z ) 1  2  1  2  -  x  z  2  2irrdr  J dk kln[det{l + H R ) ] a  a  J  J  with H , R , and 1 defined following Eq. (A.1.5). Note that the terms J2  Q  Eqs.  B.3  (B.2.2)  and  (B.2.3)  cancel and therefore the values p  +  fdzp (z) a  in  and p_ are not needed.  R e m a r k s o n the N u m e r i c s  The use of the contact theorem requires an accurate determination of the contact densities. Taking the density of the layers closest to the wall for this is not sufficient, as the wall layers represent an average over a thin but finite region and the profiles close to contact are usually very steep. We find that a third-order polynomial extrapolation (using the four closest layers) improves the accuracy significantly and makes the contact densities somewhat independent of the layering. The contact theorem and midplane formulations should yield comparable values for Fsiit,  although small deviations must be expected as only approximate (HNC) correlation  functions are inserted in Eq.  (2.7.3).  This is a helpful check of whether the midplane  algorithm is implemented in a correct manner, or, as in the case of the M M calculations for finite bulk densities (Section  4.4.2),  whether the correlation functions obtained are  accurate enough to employ the midplane formula at all. From calculations with pure HS  Appendix B. Formulas for P nt  135  s  fluids, where the only terms for the midplane formula are P^ and P re, we deduced that n  CO  the integrals in jPcore must be evaluated with extra care. Overall, fluctuations in the P i;t s  values due to the layering and the numerical procedure are of the order of a few percent, without significant differences between the two methods of calculating the pressure. The accuracy of fi is comparable to the accuracy of P u from the contact theorem A  s  t  and midplane formula routes, but the differentiation with respect to d n has to be perwa  formed via a numerical finite difference operation, which increases the error considerably. This renders the third method of evaluating P for comparison.  slit  not practical and here we used it merely  Appendix C P B and D L V O Theory  Poisson-Boltzmann (PB) theory is the mean-field treatment for point ions in a dielectric continuum for an inhomogeneous system, which was developed independently by Gouy and Chapman around 1910. Compared with the exactly solved P M , at the P B level both ion size effects and correlated charge fluctuations are neglected. However, P B theory might still be the most widely employed theoretical level to study electrical double layers and to analyze the interactions between colloidal particles. Here, we merely give a brief outline of the ideas involved. The starting points are the Poisson equation, which describes the charge density as —*  source for the electrical potential, <j)(R),  and a Boltzmann ansatz for the particle density  p (R) = p exp{-Pq <p(R)},  *  0  a  a  a  (C.2)  where p° is a constant (usually the bulk density of species a). Eqs. (C.l) and (C.2) must a  be solved self-consistently with the appropriate boundary conditions for <j)(R) [1,95]. For the slit geometry with constant surface charges, with counterions as the only species present, the problem is reduced to solving (see for example [96])  2e ek Ts  qcrz  2  0  '  W  =  B  ««cos*(Ww)  max  W  '  136  t h  S t a n ( s )  =  - ^ T '  ( C  '  3 )  0  Appendix C. PB and DLVO Theory  137  Alternatively, one can obtain the P B p (z)  through iteration of Eq. (2.6.3) with only  a  the mean field electrostatic contribution, Eq. (A.1.14), included, i.e., by setting c ( r ) = Q7  - / 3 « Q ( r ) and h (r) = 0. 7  ai  The net pressure between like-charged walls obtained from P B theory  Q  a  B  (C.4)  lim p (0)  p {0) -  kT  •Pnet — Pkin ~ Pbulk —  «wall->°0  is always positive, which means that the double layer interaction is repulsive . For 1  ^waii -> co one finds [3]  = f° }  2 Za p  P  net  ~ exp{-/«Z  wall  }  with  K  .  (C.5)  K ~ is called the Debye length and arises also as the screening length from Debye-Hiickel 1  theory. A few examples of P [ P B ] in comparison with P M and M M results are shown in net  Section 4.4.1 . Note, that for the P M two additional contributions to P  net  arise: P  c o r e  >0  (due to ion-ion hard-core interactions) and P \ < 0 (due to correlated charge density e  fluctuations). As long as core-core interactions between ions are not important the former contribution can be ignored and therefore P [ P B ] > P [ P M ] . net  net  The "classical" theory for the interactions between colloids, the Derjaguin-LandauVerwey-Overbeek (DLVO) theory [95,98], includes also an attractive van der Waals interaction between the macroparticles. The corresponding pressure contribution for the slit system is (as a first approximation) [1]  A  (C.6)  vdW 1  27T<4all  where A is called the Hamaker constant. A depends on the fluid medium and on the nature of the walls. At small d \\ (a few Angstroms) P dw + P*net < Oand one expects an wa  v  This still holds for other macroparticle geometries and, in addition, for two macroparticles confined by two walls with constant <j) boundary conditions, as has recently been shown by Neu [97]. x  Appendix C. PB and DLVO Theory  138  effective wall-wall attraction, whereas at larger separations the repulsion outweighs the van der Waals term. For many systems studied experimentally D L V O theory has been found to be satisfactory down to wall separations as low as « 20 A [2], although this is often due to the procedure of fitting the results to an "effective" surface charge, which probably differs from the true value [99].  Appendix D McMillan-Mayer Theory  In the following we derive some results from the McMillan-Mayer solution theory with the inhomogeneous system in mind. The original paper was published in 1945 [68]. Excellent summaries including applications have been given by Hill [19] and by Friedman and Dale in Ref. [100]. Most of the derivation is easily done for a general multi-component mixture, but we restrict the discussion from the beginning to two species: a solvent component, c, and a solute species, s, which is for the moment assumed'to be uncharged (to avoid problems related to the electroneutrality of the system). For mixtures of overall N^ + Ng — N particles the definitions from Eqs. (2.2.7) - (2.2.9) can be extended to the n-particle density functions p ({n};a„a )  =  {n)  s  p^l^ {R ...,R \a^a ) n  = ^TVT E  u  n  T n T f \  s  f d { N - n }  exp[-8(V +U )] N  N  , (D.l)  where the dependence of p^ and E on the activities has been indicated (as it will be exploited later) and the abbreviation {n} for the indices and positions of n particles has been introduced. VV and U are the A^-particle potential energies due to the external N  potential and particle interactions, respectively. Assuming that n and n of the n indices ?  in  s  refer to solvent and solute particles, respectively, with n = n +n , (;  substitution  m  a  —  N  a  —  n  a  ,  t)  t)  and making the  Eq. (D.l) can be written in the form  // )({n};o o,)H(o a,) n  s  - " <" * E - ~ri / m  m  £  m>0  139  ^  [-Wm+n + Wm+n)] •  (D.2)  Appendix D. McMillan-Mayer  140  Theory  Integration over d{n} on both sides yields  A^r-  d{n}pM({n};a<,a )=  £  s  d{m+n} exp [ - / ? ( V  m + n  + U )}  .  M+N  (D.3) We now carry out a two-dimensional Taylor expansion of E ( a , a ) about the "reference f  s  acitivities" a* and a* for S ( a , a ) to obtain ?  S  K.««)  =  s  E S / W ^ N ^ + W J V ) ] v  (a.-a^K-q;)"-  ^  njn !  (D.4)  /^H(a„a )\ s  9 a ^ a a ? « J a, = a*  v  s  {  1  J  a = a* s  =  ^ K-o;)^(a.-a;)"« ^ n>o n .n . > q  s  m  0 0 /" m .m . J  0  q  d { m + n  } -/»Cv ^-w. .) c  B  +  i  ( D  .  6 )  s  where for the last step the partial derivatives in Eq. (D.5) have been evaluated using Eq. (D.4) and the substitution m = N —n a  a  has been made. Comparing Eqs. (D.3) and  a  (D.6) one finds a relationship between the partition functions for the activity sets (o , a ) ?  s  and (a*, a*). One has  SK, a.) = ~(a*, a:) £ < " u^X^ /<> (a  ^ "'  dn  (D  Analogous to Eq. (2.2.10), the n-particle distribution functions are defined via  g  <*(H;a ,a,) =f«">•«»->  .  t  (D.8)  nr=iPai(^;ac,a ) s  We further introduce the n-particle potential of mean force by setting wW({n};a„a ) = - A : T l n ^ ( { n } ; a , a ) . s  B  c  s  (D.9)  The latter definition is motivated by investigating Eq. (D.2) in the dilute gas limit a ,a —»0. In this limit only the term for m = 0 survives and we find c  s  p ({n};0,0) (n)  =  a»xp[-/?(V +Z4)] n  (D.10)  7)  Appendix D. McMillan-Mayer Theory  = =  141  n{p ( R ;0,0)exp[^ (P )]}exp[-/?(V +W )] Qi  J  i  Qi  i  n  (D.ll)  n  f[ {R ;0,0)exp[-pU ] , (D.12) t=l where S(0,0) = 1 [see Eq. (D.4)], lim ^ p (R) = a exp[0v (R)], and V„ = E?=i have been used. Equation (D.12) yields the interpretation for the n-particle potential of Pai  i  n  aa  0  a  a  a  mean force in the dilute gas limit, w^({n};0,0)=U ,  (D.13)  n  which is just the total potential energy due to particle interactions. Inserting Eqs. (D.9) and (D.8) into Eq. (D.7) and further expressing p (R; a , a ) a  c  s  via the (1-particle) potential of mean force W (R; a , a ) [see Eq. (2.3.7)] one finds the a  ?  s  McMillan-Mayer expression  E(a*,a* ) s  ^  Tj..in.i  n>0  x Jd{n}  exp  -P E ^ ( ^ ; <, <) + (")({n}; a*, a* ) W  \i=l  s  /  For discussion purposes, we consider the case where the solvent activities are the same for both systems, a^ = a*, and the solute activity is zero in the reference system, a* = 0 (i.e., the infinite dilution limit for the solute concentration in system "*"). This important situation arises for the slit system when the bulk is infinitely dilute with respect to the electrolyte. Then, the reference system "*" represents, for example, the uncharged wall case with a*± — a* = a*_ = 0 and with pure solvent between the walls. The system of +  interest (no superscript "*") has a finite negative wall charge density, which is cancelled by an equivalent number of counterions. If one assumes a very small but finite a±, then the surface potential due to the double layers decreases the coion activity, a_, even further (see Eq. (2.6.6) with p>0 and Fig. 2.3). Therefore a_ = 0 remains valid for the charged  Appendix D. McMillan-Mayer  Theory  142  slit in equilibrium with a bulk in the limit that the salt is infinitely diluted. In contrast, a  +  increases because of the surface potential and one has a+^>a_, i.e., the counterions  are the only solute species to be considered (a = a s  +  in the following). From Eq. (D.14)  we obtains for the case with a = a* and a* = 0 c  E(a*a )  ^  s  <  f .  s  -P £ W (R\; a*, 0) + w^\{n }; S  n >0  "«  s  s  a*, 0)  , (D.15)  \i=l  where the summation is now only over solute particles. Formally, the expression on the RHS is identical to the partition function of a one component system at finite concentration [set one of the a to zero in Eq. (D.4)], with the difference that the (vacuum) a  particle-particle interaction energies are replaced by "effective" interaction energies, given by the potentials of mean force of the infinitely dilute solute particles in the solvent at activity a*. Remarkably, the solvent does not enter Eq. (D.15) explicitly, but only via W  and u>( '). Of course, all changes in the solvent structure due to the presence of the n  s  solute are implicitly included in the solute potentials of mean force. Equation (D.15) can be further simplified as the number of counterions (per unit surface area) for the case described above is fixed by the wall charges to n = — 2o/q . s  +  A l l terms with n ^ n are zero because of the infinite electrostatic contribution in Yl W s  s  s  for non-electroneutral slit systems. Hence, the grand canonical expression (fixed chemical potential) reduces to the corresponding canonical expression (fixed number of particles) and one has 2 f l . ° «  .  s  s(o*,0) n  j/ d{n d{n\ } ^n" n\  S  tt  =  1 , 7 1 - 1 s s  s  exp  -(3 £ W (Ri; a*, 0) + w^({n }; 8  s  a*, 0)  (D.16)  \i=l  Equation (D.16) is exact if the contributions of practice, one introduces the pairwise approximation  i=l j>i  to all orders are included. In  Appendix  D. McMillan-Mayer  143  Theory  where w (Ri, Rj) is the pair potential of mean force. Unlike for the potential energy U , ai  N  the pairwise approximation for  is severe, as the effects due to changes in the solvent  structure around the ions cannot be expected to be additive. Still, at electrolyte concentrations up to several moles/litre the pairwise M M level for bulk electrolyte solutions gives reasonable results [71]. For the inhomogeneous system considered in this thesis another approximation is made, namely replacing w (Ri, Rj) by its bulk counterpart ay  w (Ri, Rj) « w™\\Ri ai  - Rj\) ,  (D.18)  and thereby ignoring effects due to solvent ordering close to the walls. The implications of these approximations for the appropriate choice of W (Rj,;a*,0) in Eq. (D.15) are s  discussed in Section 5.2 .  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0061490/manifest

Comment

Related Items