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 this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. 1 further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada DE-6 (2/88) 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 cho-sen 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 contin-uum 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 mono-valent 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 attrac-tion 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. ' i i i Table of Contents Abstract i i List of Figures v i i List of Tables xv List of Important Symbols and Abbreviations xv i Acknowledgement xx 1 Introduction 1 2 Equil ibrium 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.3.1 Particles of Equal Size 29 3.3.2 Larger Ions 39 3.3.3 Monovalent Counterions 44 3.4 Summary 49 4 Molecular Solvent Effects at the McMillan-Mayer Level 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 4.5 Summary 83 5 Further Investigations at the McMillan-Mayer 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 Ion-Wall Potentials of Mean Force 97 5.2.1 The Model 98 5.2.2 Results 102 5.2.3 Summary 104 5.3 Remarks on Some Experimental Results 106 6 Summary and Conclusions 111 v Bibliography 117 Appendices < 122 A Numerical Algorithm 122 A. l General Method 122 A. 2 MM Systems with Strong Anion-Cation Attraction 128 B Formulas for P s i i t 130 B. l Midplane Formulation 130 B.2 Calculation of ftA 133 B.3 Remarks on the Numerics 134 C P B and D L V O Theory 136 D McMil lan-Mayer Theory 139 vi List of Figures 2.1 The frame of reference used in the calculations. Note that this a two-dimensional representation of the three-dimensional system. For clarity, only the case where all particles have the same diameter da = d is considered. 15 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 dou-ble 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 vel = qa/A/(2ee0) • The corresponding case for a spherical macroparticle in b) yields vel = 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) 21 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 dion = 4s and Zion = +2. In (a) d w a n = 2.6dhs and in (b) d w a l l = 3.1 d h s - • • • 30 vii 3.2 Contact densities for (a) hard-sphere and (b) ion components for different models. The model labeled UCPM+HS is identical to PM+HS but the charges are switched off as described in the text. The remaining models are as in Fig. 3.1. In all cases d[0n = 4s and zT ion = +2 32 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 di0Ti = 4s and zT ion = +2 34 3-4 #ionion(7") parallel to the walls in the midplane and contact plane for the P M and PM+HS (rfion = 4s) with d w a l l = 2.05 dhs and Zion = +2 35 3.5 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. . . . 37 3.6 The electrostatic component of the pressure for the P M and PM+HS with ^ion = +2. Results for 4)n = 4s and for 4,n = 1.52 4s are shown 38 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 4m = 1-52 4s and Z\oxl = +2. In (a) 4aii = 2.5 4s and in (b) 4aii = 2.9 4s 40 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 d-l0n = 1.52 4s and Zion = +2 42 3.9 The net pressure acting between the walls obtained for the different models and with different superposition approximations. In all cases d10n —1.52 4s and zTj0n = +2. The labels are as in Fig. 3.2 and as discussed in the text. 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 Zion =+l. In (a) 4aii = 2.6 4s and in (b) 4aii = 3.0 4s- • • • 45 viii 3.11 Contour plots of p i o n i o n ^ , 0, z) for the P M and the PM+HS model with ^ion = and Z\on = +1 for d w an = 2.6 d n s . r and z are in units of d 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, —zma,x, z) for the P M and the PM+HS model with G?ion = ^ hs and zT;on = + l for d w a n = 2.6dhs- r and z are in units of dns- 47 3-13 ^ ionion(^) parallel to the walls in the midplane and contact plane for the P M and PM+HS (dion = dhs) with dVfan = 2.0dhs and Z 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 56 4.2 P n e t and its components for systems including only monovalent counterions with a = -0.1307 e/d2 = -0.267 C / m 2 . The Poisson-Boltzmann (PB) curve is shown for comparison (see Appendix C), where dw an[PB] is redefined to give the same accessible space (i.e., the same z m a x ) as the finite size ion models. P c o re[MM] is less than 2 x 10 - 5 kBT/d3 throughout and is therefore not shown. For comparison: 0.020fcBT/d3 = 1.5 M P T = 3.7 x 10 6 N / m 2 , where R is the gas constant 59 4.3 P n et and its components for systems including only monovalent counterions with a= -0.0653 e/d 2 = -0.134 C / m 2 . The Poisson-Boltzmann (PB) curve is shown for comparison (see Appendix C). P c o r e [ M M ] is less than 2 x 10 - 6 A;sT/d 3 throughout and is therefore not shown 61 ix 4.4 P n e t and its components for systems including only monovalent coun-terions with a = -0.1634 e/d2 = -0.334 C / m 2 . P c o r e [ M M ] is less than 4 x 1 0 - 5 kBT/d? throughout and is therefore not shown 62 4.5 The a dependence of P n e t and Peff for systems including only monovalent counterions at fixed d w an = 4 d 63 4.6 P n et and its components for systems including only divalent counterions with cr = -0.1307e/d2 = -0.267 C / m 2 . The Poisson-Boltzmann (PB) curve is shown for comparison (see Appendix C). P c o r e [ M M ] is less than 10~ 1 7 kBT/d3 and P c o r e [PM] is less than 4x 10~4 kBT/d3 throughout, there-fore these curves are not shown 64 4.7 Comparison of P s u t calculated via the contact theorem and via the mid-plane formula at a = —0.1307 e/d2. 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 w an is also shown 66 4.8 Density profiles for 1:1 electrolyte solutions between uncharged walls, p is LOOM = 0.0132/d3 and 1.03M = 0.0136/d3 for the P M and the M M model, respectively, p is indicated by horizontal lines. Results,are shown for d w an = 1 d, 2 d, and 8 d. The profiles for cations and anions are identical. 69 4.9 Density profiles for 1:1 electrolyte solutions between charged walls with a = —0.1307 e/d2 = —0.267 C / m 2 for (a) counterions and (b) coions. p is LOOM = 0.0132/d3 and 1.03M = 0.0136/d3 for the P M and the M M model, respectively, p is indicated by horizontal lines. Results are shown for d w an = 2.5 d, 4.0 d, and 6.0 d 70 x 4.10 Midplane densities of monovalent counterions and coions at <7nigh = —0.1307e/d2 and C T I o w = a h i g h / 2 for (a) the M M system and (b) the P M . The bulk den-sities p are as in Fig. 4.9 and are indicated by horizontal dotted lines. . . 7 4.11 P n e t for 1:1 electrolyte solutions at p = 1.0M for cr hj g h = —0.1307e/d2, ciow = Ch igh / 2 , and for uncharged walls 7 4.12 P n e t for M M systems in equilibrium with 1:1 electrolyte solutions at various bulk concentrations and o — —0.1307 e/d2 = —0.267 C / m 2 . The curves for p = 0.5 M and 1.0 M are indistinguishable within numerical error and therefore only the latter case is shown 7 4.13 The dependence of P n e t on a at fixed dw an — Ad. Results are shown for 1:1 electrolyte solutions at p = 1.0 M and for the the infinite dilution (counterions only) case 7 4.14 $ 0 7 ( r ) a t the midplane parallel to the walls (z\ — z^ — 0) for 1:1 electrolyte a = —0.1307e/d2. = is shown in (a) and g++ and g ("4-" denotes counterion and "—" the coion) are shown in (b) 7 4.15 Contour plots of g++(r, — z m a x , z) for systems where only monovalent coun-terions are present with d w au = 3.6d and cr= —0.1307 e/d2. 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 7 4.16 Contour plots of g++(r, —zmSiX, z) for systems where only divalent counte-rions are present with rfwall = 2.8rf and a = —0.1307 e/d2 8 4.17 Contour plots of p+(z)g^.+ (r, 0, z)dz, the average counterion density around a midplane coion. The ions are monovalent and p=1 .0M, d w an = 4d, and o — —0.1307 e/d2 8 xi 4.18 Contour plots of p-.(z)g+-(r, — zmax, z)d3, the average coion density around a counterion in contact with one wall. The system parameters are as in Fig. 4.17 82 5.1 u e f f (P) components of the ion-ion potentials of mean force at infinite dilu-tion. The curves labeled [KP], [GRP] and [PR] are from Refs [52], [74,75] and [76], respectively 86 5.2 The dependence of Pnet and its components on the wall separation for N a + [GRP] and N a + [PR] counterions. The wall charge density cr = —0.334 C / m 2 = — le/48A 2 . M represents moles/litre and R is the gas constant 89 5.3 The dependence of Pnet 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 = +le/48A 2 . The notation is as in Fig. 5.2 90 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 d10n = 2.94A. The notation is as in Fig. 5.2 92 5.5 Ion-ion pair distribution functions at the wall plane (z\ = z2 = 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 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 93 5.6 P n e t for the NaCl [GRP] model [74,75] at different salt concentrations, p, and cr = —0.334C/m 2 = - l e / 4 8 A 2 . The inset shows the tails on an expanded scale. The notation is as in Fig. 5.2 95 xii 5.7 Comparison of P n e t for four systems with monovalent counterions, includ-ing the P M with dion = 2.94.A, and a = - 0 .334C/m 2 = - l e / 4 8 A 2 . The 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 96 5.8 Two model ion-wall potentials for counterions, V+> a s functions of the distance to the planes of closest approach to the walls, £. For C= (zmax ± z)>ld and for coions at all distances tpa = 0. For x = C/d<l the curves are given by the polynomials 2.75a;3 — 7.75a;2 + 7.25a; — 2.25 (-01) and 71.5a;5 - 145.9a;4 + 91.0a;3 - 9.3a;2 - 7.7a; + 0.4 (^2) 101 5.9 P n e t for monovalent counterions and a wall charge o = —0.1307e/d2, with 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 (-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/d2. The lines and symbols indicate the results for effective ion-wall potential ipl and for the case with no additional short-range ion-wall potential (•0 = 0), respectively 105 xiii 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 separa-tion 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 cor-responding to a point between "F" and " G " and the local maximum at "E" cannot be detected experimentally 110 A . l Iteration scheme for the algorithm 126 xiv List of Tables 5.1 Fit parameters for the Pettitt-Rossky [76] ion-ion potentials of mean force. The units are chosen to yield ueS in fcBT if R is given in A 87 xv 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 w an 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 7Yr number of grid points parallel to the walls Nz total number of layers (grid points perpendicular to the walls) P pressure Pbuik bulk pressure Pnet pressure difference between the slit and the bulk xvi Psiit pressure in the slit q charge of a particle R position vector R particle separation r cylindrical coordinate parallel to the walls 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 ueS effective potential component of w due to molecular solvent uel Coulomb pair potential M H S hard-sphere pair potential ush 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) u e l Coulomb interaction with the walls vHS hard-sphere interaction with the walls vsh 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 xvii z coordinate perpendicular to the walls i-^ max z-coordinates of planes of closest approach to the walls 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 m a x p one-particle density p{n) n-particle density function a surface charge density grand partition function c distance from plane of closest approach to a wall —* V gradient operator / Hankel transform of f A H N C anisotropic hypernetted-chain B G Y Born-Green-Yvon HNC hypernetted-chain HS hard-sphere M moles/litre xviii M D molecular dynamics M M McMillan-Mayer M S A mean-spherical approximation P B Poisson-Boltzmann P M primitive model R H N C 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]. An 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 fun-damental 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 sur-face 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. Some-times 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 sep-aration independent dielectric constant, have received considerable attention and their influence on wall-wall forces is now reasonably well understood [4]. The classical Poisson-Boltzmann (PB) theory treats point ions at a mean-field level and always yields a re-pulsion between like-charged macroparticles, which decays exponentially with increasing wall separation. By 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 PB theory. For example, it is now well established [4] that like-charged walls immersed in a P M solution with divalent counte-rions 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 mech-anism 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 repul-sive 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 fi-nite macroparticle concentrations [7], and additional macroparticle attractions for models that allow for a variation of the wall charge via ion absorption potentials [8]. Recent ex-perimental 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 fluctua-tions) 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 separa-tion, 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 (^ 1M) . 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 ex-plicitly 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 prob-lem 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 McMillan-Mayer 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 the-ory 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,fluids1 however, 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 1 Above the critical temperature no distinction between gas and liquid states can be made and we refer to the state of matter as "fluid". 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 meth-ods 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 con-text 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 func-tional 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 Basic Definit ions 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 orien-tation of the particles. The definitions are straightforward adaptations from standard Chapter 2. Equilibrium Statistical Mechanics 8 ^ ^TJ^JdRN eM~PVN({P^}) - PUA{^})\ , (2-2.2) formulas given, for example, in [17,19]. The Hamiltonian, riN, of a system consisting of —* ./V particles, with coordinates Ri, momenta pi, particle types a,, and masses mai, is given by HN({Ri},{pi}) = £„({&})+ VN({&})+UN({Ri}) = £ ^ + £ ^ ) + £ £ W ^ ) > (2-2-1) i=l Z'' iai i=l i=l j>i where /CAT, VAT, and UN are the contributions due to the kinetic energy, an external potential (va(R)), and the pair particle interactions (uai(Ri,Rj)), respectively. The grand partition function, E(pA, pB,V,T), is defined by exp(NAppA + NB(3pB) r,^f,SN r anj (tB\ r - i v i - = £ N m iH3N dP d R em-^N({Ri},{pi})\ NA+Nb=0 ^A-^B-n j aA aB where NA + NB = N, j3 = (kBT)~l, and the activity for species a is given by aa = e^/Al. (2.2.3) V, T, kB, h, and pa 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 sepa-rations 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(pA,pB,V,T) = -kBT\nE(pA,pB,V,T) . (2.2.5) Chapter 2. Equilibrium Statistical Mechanics 9 With this relationship, the system's pressure, P, for example, is given by p==-{dv) T=kBTM T- (2-2-6) 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). A different approach employs the n-particle density functions, p^[„an(Ri, Rn), which are obtained by integration over (N — n) 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 oo oo NA NNB PA(RI) 4 E E ( N A ZN , dR2...dRN e-W"W . (2.2.7) The AA and AB 2-particle densities are given by i oo oo NA NB p!&(£i,&) = , E E * a , dR3...dRN e-flv*-"*) a n d ( 2 . 2 . 8 ) - Na=2nb=0 ("a ~ 2 j ! / V B ! j i oo oo NA NNB R fikRiX) 4 E E ( n tufr. 1 V dh\...dRN e~W»W .(2.2.9) - NA=I NB=I -iy.{NB -iy. J The pair distribution functions, gay, are defined via n (P P \ - J^S^ll^L (oo-in\ gai(Kl,H2) - —-=r- • (2.2.10) pa{Ri)p1{R2) 9a-y(Ri, ^2)^7(^2) is the particle density of species 7 at R2 under the condition that a particle of species a is located at Ri. If \Ri — R2\ is large compared to the typical length over which particle interactions are significant (and the fluid is not close to a Chapter 2. Equilibrium Statistical Mechanics 10 critical point) then ^ a 7 ( P i , P 2 ) ~ 1 and the particle densities at R\ and R2 are not correlated. Pair distribution functions can often be measured experimentally via X-ray or neutron scattering, at least for bulk fluids, where gay{Ri, R2) = ^ ( l - R i — -R2I) is called 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) calcu-lation of the functions pa(R) and ga^(Ri, R2), whereas higher order particle correlations 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 a i , hai(R\, R2) = gay{Ri, R2) — 1 j (2.3.1) and the direct correlation function, caj, via the integral relationship2 ha-y(Ri, R2) = ca~f(Ri, R2) + ^2 dR3cas(R\, R3)ps(R3)h$1(R3, R2) • (2.3.2) 5 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 (va(R) = 0) pa is translationally invariant and 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. 2Note, that ca^(Ri,R2) can be alternatively defined via the functional derivative of \n[pa(Ri)/aa] + 0va(Ri) with respect to p 7 (^)-Chapter 2. Equilibrium Statistical Mechanics 11 In order to calculate hai for a given set of pa, another equation relating hai and c a 7 is needed. Several authors studied these relationships for homogeneous systems in the late 1950's and 1960's, using diagramatic3 and functional differentiation techniques (for 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 /* Q 7 (P i , i ? 2 ) + l = e x p [ - / ? ^ ( ^ ^ , (2.3.3) where bai is called the bridge function. Although certain features of the bridge function are known (it is short-ranged for hard-sphere (HS) systems [23]), the exact calculation of ba^(Ri, R2) is not possible, as an infinite sum of "bridge diagrams" is involved. Hence, all 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 pa{R) is not a constant, requires an additional equa-tion for the average particle density. One option is the first member of the Born-Green-Yvon (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 [dR2py(R2)gaj(Ru^V.u^R,,R2) , (2.3.4) 7 J 3 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). Chapter 2. Equilibrium Statistical Mechanics 12 although the Wertheim-Lovett-Mou-Buff (WLMB) 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, / ^ X ( P L ) , defined via pa = p}^l(R) + pe*(R) = kBT[3\nAa + \nPa(R)} + va(R) + p™(R) . (2.3.5) pea{R) represents the change in £1(PA, PB,V,T) due to particle interactions when a parti-—* 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, uai(R\, i?2|£ = 0) = 0, and £ = 1 stands for the fully "switched on" particle, uay(Ri, R-2\l) = uai(R\, R2). We further define gai(Ri,R,21£) 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 R2 by d^duai(Ri, R2\^)/d^ and integration over the whole system for the total coupling process, £ = 0 —> £ = 1 , yields pe*(Ri) such that = E / d R 2 Pl(R2) J\ d u « ^ m ) g a , ( R u R2\0 • (2.3.6) Eq. (2.3.5) gives for pa(R) pa(R) = aaexp{-(3[va{R) + p,f(R)]} = aaexp{-0Wa(R)} , (2.3.7) where Wa(R) is the potential of mean force acting on particle a located at R. Introduction of the activity coefficient, ja(R), via aa = pa(R)ja(R) yields therefore 1*(R) = exp{/3Wa(R)} . (2.3.8) 2.4 Closure Approx imat ions 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 ba-y(Ri,R.2) = 0 in Eq. Chapter 2. Equilibrium Statistical Mechanics 13 (2.3.3) gives the hypernetted-chain (HNC) closure hay(Ri, R2) + l = exp[-puay(Ru R2) + hai{Ri, R2) - cai(Ri, R2)} . (2.4.1) This closure yields good results when the system is dominated by long-range forces [27,28] such as the 1/R2 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 / d R 2 p y { R 2 ) [ h l 7 ( R u R 2 ) - caj(RuR2)] - ^[c^RuRJ + 1] . (2.4.2) The HNC 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 bai(Ri,R2) 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^(RuR2) = ( l - vcp{(3uai{Ru R2)}) gai(Ru R2) (2.4.3) is preferable over the HNC closure [31] (although fewer diagrams are taken into account) and yields excellent results for inhomogeneous systems [32,33]. Chapter 2. Equilibrium Statistical Mechanics 14 When particles interact via a hard-core potential (u Q 7 ( i? i , R2) = oo, if \Ri — R2\ <dai), 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 cai(Ri, R2) = — (3uai(Ri, R2) outside the hard core and exploiting <?Q7(i?i, R2) = 0 inside the hard core. 2.5 General Model 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 da and the charge qa = Zae, where e is the elementary charge. Particles of species a and 7 at a distance R interact via the pair potential uai(R), which can be expressed in the form uai(R) = u™(R) + u%(R) + u%(R) (2.5.1) 00 , if R < \(da + d 7) , = < « 9 L ^ + U s h ( i 2 ) > otherwise, I AirtocR ' where e is the dielectric constant of the solvent and eo is the vacuum permittivity. The terms u^(R) and u^(R) are the hard-sphere and Coulombic interactions, respectively, while u^(R) is an additional short-range contribution. An electrolyte model comprising Chapter 2. Equilibrium Statistical Mechanics 15 d/2 t T Z max z = 0 D r Z2 -z t max 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 da = d is considered. 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 dWaii- Hard-core particle-wall in-teractions are assumed and therefore the space accessible to particle a has the width Da = dwa\i — da. Each wall is charged with a surface charge density a and taken to have 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. lC/m 2 ) . In our calculations we employ the reference frame shown in Fig. 2.1 with coordinates z and r = \/x2 + y2, where z = 0 and z = ±zmaiXiCt (zma,x,a = Da/2) 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 Mechanics 16 © 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 / ( — z ) is an additional short-range 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 function ipa, with ^Q(0) = ^ (-^max.a) = <h'/7(^max,a), vf'I(z)=ipa(\z + z max, a|), and vSa,n (z)=if)a(\z — zmaXiQ\). The assumptions implicitly made by writing this contribution in the above way is that the effects from both (identical) walls are superimposable and that ipa does not depend on dw an itself. 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, pa, for each species, Chapter 2. Equilibrium StatisiicaUMechanics 17 or, in the case of an electrolyte solution, p±. Specifying the pa is equivalent to fixing the bulk densities, pa. The walls are assumed to be wide and long enough, so that neither 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, pa(z), depend only on the z-coordinate, whereas the correlation functions, hay(r,Zi,Z2) and c Q 7 (r, 'z\, z2), are functions of three variables. The OZ equation, Eq. (2.3.2), takes the form haj(r,z1,z2) = cai(r,zuz2) + (2.6.1) with r3 = \Jx3 + y2 and r2S = ^J(x3 — r)2 + y\. The HNC closure, Eq. (2.4.1), becomes haj(r, zu z2) + 1 = exp[-/?uQ 7(i2) + hai(r, zu z2) - c Q 7 (r , Zi,z2)] , (2.6.2) with R = \jr2 + (z2 — Zi)2, and the equation for the density profiles is given by combining Eqs. (2.3.7) and (2.4.2), yielding pa(zi) = aaexp{ -8va(zi) - ^[cQ Q(0, zu zx) + 1] (2.6.3) " 2 7 r H Jdz2p1(z2) jdrr\^h2ai{r, zu z2) - cQ 7(r, zu z2)] } 7 We prefer this approach over the use of the B G Y or W L M B equations, as once the aa 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 w a i i to keep the fluid between the walls in equilibrium with the bulk. In addition, we apply the Chapter 2. Equilibrium Statistical Mechanics 18 electroneutrality condition to the slit system, which means that 2o + ^Qa fdzpa(z) = 0 . (2.6.4) a J Electroneutrality for the slit at small dwan (i.e., a few ion diameters) can only be expected 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) 4 . Once the 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, aa follows directly from Eq. (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\Z+_\XZ~, i.e., one unit of the salt dissociates into cations of valency Z+, and Z+ anions of valency Z _ . We start from the definition for the mean chemical potential, p±, t*± = Z + + 1 | Z _ | ( I Z - I ^ + + ' ( 2 - 6 - 5 ) 4 A study for equivalent systems using density functional methods has been recently published by Henderson et al. [37]. 5 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. Chapter 2. Equilibrium Statistical Mechanics 19 which is satisfied by setting A 3 A 3 p+= p± + Z+p + kBTln -j± and p,_ = p± + Z_~p + kBTln , (2.6.6) where A = + 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 Xz+a+ [-*$- + Xz-a„ 1° = 0 , (2.6.7) J 7+(z) J 7_(z) exp{/ fy ± }/A 3 with 7a(-2) — aa/pa(z) 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± ""(z) = xz"^)-w (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 (va = 0) applied to Eq. (2.6.5) gives - ^ ^ = 7 ± p + ( j ^ i J = 7 ± " - U r J • ( 2- 6 ' 9 ) where the mean activity coefficient 7± is defined via ( z + + | z _ | ) / | Z _ | Z + / o c i n \ 7± = V7+ 7- ; (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 w an the midplane densities must equal the desired bulk concentrations, P + and p_. The advantage of writing the constant exp{/3p±}/A3 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 Mechanics 20 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, Da, remains unchanged. For example, relocating the charged planes further apart (while keeping Da constant) changes the constant d w an —» d w a l i in the expression for u®1 [see Eq. (A. 1.14)], 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 walls6. The effect disappears when 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 HNC 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. 6The chemical potential of the counterions close to a charged wall for a PM electrolyte solution has been studied by Svensson and Woodward [40]. 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 vel = 0crA/(2eeo). The corresponding case for a spherical macroparticle in b) yields vel = 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). Chapter 2. Equilibrium Statistical Mechanics 22 2.7 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 component7. For the former component two contributions must be considered. One is due to the attractive short-ranged van der Waals interaction which decays as d~ln and will not be included in our discussions. The other, here called P sn t, arises from the interaction 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 / ° ° dzpa(z)^^- , (2.7.1) i.e., the derivative of the average interaction energy of the fluid with the wall. Inserting vIa(z) = vaS,I(z) + vtJ{z) + vSa''{z) i n t o E c l - (2.7.1) and using Eq. (2.6.4) one obtains ^ (J 2 - . /-2max,Q dVSh,I(z) Pm = ^ T ] > > a ( - z m a x , Q ) ~-0 £ / dzpa{z) 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 7The normal pressure component, Px, is constant throughout the system, whereas the parallel com-ponent, pH (z), is not uniquely determined and depends on the position between the walls. For limdw a l l_^0 O one has, of course, PH(0) = P x . and the quantity J02maxcte[P-L - P"(z)] is the surface tension [44]. Chapter 2. Equilibrium Statistical Mechanics 23 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 sij t is split into a term depending on the midplane density, Pkjn, and contributions representing the various particle interactions [see Eq. (2.5.1)] across the midplane. One finds (see Appendix B for details) -fslit = -Pkin + Pcore + Pel + P s h + P w a l l , with (2.7.3) Pki„ = A:BT5>Q(0), a k fZma.x,a fO -Pcore = /cBT2_,27r / dzipa(zi) / dz2p1(z2) x (zi — z2) <7Q7(rQ7, zi, z2) , 22 dzipa(zi) dz2p7(z2) a,7 /•oo 5 < ' ( r , 2 i , Z 2 ) v x / 2-ivrdr—2XL hal(r,zx, z2) , Jo azi Fsh = ~E / " ' " " " ' ^ l P o ^ i ) / dz2Pl{z2) x / 2nrdr a i \ gai(r,zuz2) , and wall The integrations over Z ! and z2 for P c o r e are to be carried out only if r2tJ = ^ (da + d 7 ) 2 — (zi — z2)2 is positive, i.e., if two particles at z\ and z2 can be in contact. Note that the calculation of P e i involves hay, and not gaj, as the second term on the RHS of Eq. (2.7.2) is absorbed into P e i ; only deviations from the mean field case (i.e., correlated electrostatic fluctuations) contribute and in general one has P e i < 0 8. The use of hai 8 The pressure component Pe\ 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 Chapter 2. Equilibrium Statistical Mechanics 24 instead of gai in an equivalent expression for Psh.given by Eq. (3) of Ref. [50] is due to a 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 ap-propriate total thermodynamic potential of the system between the walls, here the grand potential per unit surface area, flA, and obtains the pressure via ^=-(IR • <™> \ C a w a l l / { / 1 } ) T The HNC approximation yields an explicit formula for QA, which is given in Appendix B. The net pressure, Pnet, is given by the difference between the internal (slit) pressure 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 b u l k = lim P s l i t ) (2.7.6) "wall-* OO which is exploited in our calculations to obtain Pbuik- Wall separations yielding P n e t = 0 represent a fluid-wall (macroparticle) system in equilibrium (when no additional outside force is applied), although the equilibrium position is stable only if dPnet/<9dwaii < 0. of the ions in a simple model [49]. 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. An 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 Model 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), uai(R) = u™(R) + ue^(R) and (3.2.1) va(z) = v**(z) + v*(z), (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. = kBT [phs(z max,hs ) + Pion( 2max,ion)] ~~ ~ > a n o - (3.2.3) •Pslit = -Pkin + -Pel + -Pcore , (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 HNC approximation does not perform compara-tively 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 appro-priate 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 HNC approximations with exact results [31]. Surprisingly, they found the worst performance for the A H N C theory, with even the singlet HNC 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 D uik from the exact value is « 9 % , where we use wall separations of ~ 9 dhS or larger to determine Pbuik with sufficient accuracy. Over-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 Resul ts 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 PM+HS models except that all charges are switched off while keeping the number of "ions" between the walls fixed1. We label systems of this type UCPM+HS (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 n s = 2.8 A, e = 78.7 and a temperature T = 298 K. Thus, the hard-sphere 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; o n = +2) with diameters dj 0 n = dhs and d; o n = 4.25 A —1.52 d h s 2 . 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 1This can be achieved by scaling down all charges (a and gi0n) by a factor 1000, while still applying the electroneutrality condition. 2The reduced temperatures, T* = 4ne0edionT/q2, for the systems investigated are 0.1 {Z\on = 2, 4>n = 2.8 A), 0.15 (Zion = 2, d i o n=4.25 A), and 0.39 (Zion = l, dion = 2.8 A), but note that the significance of T* for the charged slit system is changed in comparison to the corresponding bulk. Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 29 basically due to the more interesting behaviour of P s n t for the corresponding P M (a range of attractive wall-wall interaction). A few results concerning monovalent counterions are also presented. The results reported are for phs = hm<2wal_>oo Phs(O) =0.492/d n s. To achieve this value, we used a bulk activity coefficient 7 n s = 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 HNC closure. We note that our bulk density is a little lower than that of water under ambient conditions (i.e., ~ 0.7/o!ns). However, the numerical solutions are 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 pa(z) 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 w au = 3.1 d h s plotted in Fig. 3.1(b) shows Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 30 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 z / d 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 i o n = d h s and Zion — +2. In (a) dWaii = 2.6dhs and in (b) dw an = 3.1 d 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 phs = 0.7/d3hs [55,56]. A n interesting point is that the average total particle density between the walls, ptot, (calculated by integrating the density profiles of all particles for a given model) is higher for PM+HS 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 ptot[PM + HS] and p tot[UCPM + HS] are always within 0.5% of each other, whereas ptot[HS] is up to 5% lower than the values for the binary mixtures. This indicates that the UCPM+HS model is not equivalent to HS even though all the species have identical interactions. The difference in ptot 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 w au for the different models. The contact densities (Fig. 3.2) for all components in the dense systems (all except the PM) show to some degree maxima (around d w aji = 2.05 dhS, 3.1 dhS) and minima (around dwaii = 1.6dhSj2.6dhs)- 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 dns neighboring particle layers interact 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 cox: T5 c a N Q 2.0 1.0 h 0.0 0.8 0.4 0.0 (a) PM+HS HS - -UCPM+HS •••• J I L J I I L J I I L J L PM+HS PM UCPM+HS - (b) J I I i J I I L J L J L 1 2 3 dwall / ^hs Figure 3.2: Contact densities for (a) hard-sphere and (b) ion components for different models. The model labeled UCPM+HS is identical to PM+HS 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 w an « 2.15 dns due to correlated ion-ion density fluctuations. The midplane densities (Fig. 3.3) for the dense systems exhibit layering features similar to those of the contact values. For d w a i i ^ 2.2 dns) the densities at the midplane and contact plane for PM+HS and 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. Fea-tures such as higher ion contact densities and a faster decrease of the ion midplane density with increasing d w an for the PM+HS (compared to the PM) come about because of the hard-core interactions. For larger wall separations the curves for the charged and un-charged 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 w an increases, whereas a constant nonzero contact density is reached for the PM+HS [see Fig. 3.2(b)]. The solution of the integral equations also yields the particle distribution functions 9a-y(r, zi, z2)- 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 #i0nion(r) parallel to the wall at the contact plane (z\ = z2 = zmajX) and midplane (zi = z2 = 0) for d w an = 2.05 d n s . The differences between the c/j 0 n i o n(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.7dns> which constitutes the global maximum as well (as seen in a full 2D contour plot), is split into two maxima in the mixture. The oscillations in <7i0nion(?") in the contact plane are smaller and no maximum 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-wn = d„s and Z i o n = +2. Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 35 Figure 3.4: #i O nion(0 parallel to the walls in the midplane and contact plane for the P M and PM+HS (o?j0n = d ns) with d w a n = 2.05 dns and Z ; o n = +2. 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 e i , which is always attractive as it arises because of correlated ion-ion fluctuations across the midplane. The net pressure for the system of interest, PM+HS, 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 PM+HS with simple addition schemes involving more basic components, different possibilities are available. The simplest suggestion is to add the values of Pnet for the pure systems (i.e. P n et[HS] + P n e t [PM]) . From Fig. 3.5(b), we see that for the equal-sized case this gives good agreement down to dWaii~l-8dhs-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 ei[PM]. The idea leading to this choice is that for small wall separations packing effects are crucial for the fluid structure and should be separated from the electrostatic contribution. P e i on the other hand shows 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 preced-ing paragraphs are of rather academic interest as we do not suggest that any procedure Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 37 Figure 3.5: The net pressure acting between the walls obtained for the different models and with different superposition approximations. In all cases d-lon = d n s and Zion = +2. The labels are as in Fig. 3.2 and as discussed in the text. Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 38 I i i i i I i i i i I i i i i L_i i i i I 0 1 2 3 4 ( d w a i i — d i o n ) / dhs Figure 3.6: The electrostatic component of the pressure for the P M and PM+HS with ^ion = +2. Results for di 0 n = <4s and for d i o n = 1.52dhs are shown. Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 39 of adding pressure components of different models will yield exactly P n e t [PM+HS] . Fur-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 n e t [HS] and Pnet [UCPM + HS] are not the same although the "ions" and hard spheres are identical in UCPM+HS. 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 sys-tems with identical diameters for the neutral HS component, but with ions which are 50% larger, or, more precisely, dj 0 n = 4.25 A = 1.52 d h s . With this choice, the ions are considerably larger than the solvent particles and deviations from the equal-sized case are expected to be significant. Moreover, this value of d; o n 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 dw an); the divalent counterions stay sufficiently apart from each other that ion-ion hard-core interactions are not important. See, for example, the P e i plots given in Fig. 3.6. For larger values of dj o n 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 o n = 1.52 d n s and Zion — +2. In (a) d w an = 2.5 <4S and in (b) d w a i i = 2.9tZhs-Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 41 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 w a n~2.5 dhS = dhs + ^ion, which means that ions are pushed 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 w a l l = 2 dhs [Fig. 3.8(a)]. Since hard spheres in PM+HS and UCPM+HS 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 Pnet for PM+HS 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 n e t [PM] is not good at any wall separation if the particle diameters are sufficiently different. On the other hand, the more complicated scheme, P n e t [ P M + HS] ^ P n e t [ U C P M + HS] + P ei[PM], works very well for smaller values of c/waii • Deviations do occur at larger wall-wall separations where the fluid structure in U C P M + H S differs significantly from that of PM+HS. Nevertheless, these results show that to a good approximation the net pressure in the PM+HS can be viewed as a super-position 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 o n = 1.52d h s and zT i o n = +2. 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; o n = 1.52 dhs and Z i o n = +2. The labels are as in Fig. 3.2 and as discussed in the text. Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 44 3.3.3 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 PM+HS model at d w an = 3.0d n s [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\ = — zm ax,ion) f ° r a wall separation of d w a i i = 2.6d n s- The white region around the particle is forbidden to other ions because of the hard-core interactions. The P M contour plots show that beyond R=^r2 + (z — Zi)2 = d\on the likelihood of encountering another ion is reduced through the electrostatic repulsion, c/i 0 n; o n(i?) < 1 is monotonic in all directions, and only a small degree of anisotropy is introduced by the walls. In contrast, the corresponding c/ i o n i o n(i?) for the PM+HS 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 w a i i = 2.0dhS shown in Fig. 3.13 demonstrates the same effective attraction for the midplane, whereas along the walls <7iOnion(0 remains <1 (similar to the P M curves). Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 45 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 z / d 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 i o n = d h s and Zion = + l . I n ( a) d w a i i = 2.6dhs and in (b) d w a l l = 3.0 d h s -Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 46 g(r,z) for counterion/counterion PM -1.0--6.0 -4.0 -2.0 0.0 +2.0 +4.0 +6.0 r g(r,z) for counterion/counterion PM+HS 1.2000 • 1.1000 • 1.0100 • 0.9900 I 1 0.9000 0.8000 • 0.7000 • 0.6000 • 0.0000 • Figure 3.11: Contour plots of g\on\on{r, 0, z) for the P M and the PM+HS model with dion = dhs and Z j o n = +1 for d w a n = 2.6 dns- r and z are in units of d^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. Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 47 r Figure 3.12: Contour plots of # 0nion ( r , -zmax, z) for the P M and the PM+HS model with 4>n = 4s and Zion = +l for dwaa=2.6dha. r and z are in units of 4s-Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 48 G o 1.6 \-1.2 5 0.8 0.4 \-0.0 dw«u=20 d hs s 1 midplane contact plane PM PM+HS PM PM+HS 2 3 r / d h s Figure 3.13: #i0nion(?") parallel to the walls in the midplane and contact plane for the P M and PM+HS (dj o n = 4s) with d w a i i = 2.0dhs and Z\on = +l. Chapter 3. Mixtures of Neutral and Charged Hard-Spheres - Packing Effects 49 3.4 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 calcula-tions [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 to-wards 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 oscil-lations 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 e i , of the P M . This scheme 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 e i component in the mixture differs little from that of the P M . Chapter 4 Molecu la r Solvent Effects at the M c M i l l a n - M a y e r Leve l 4.1 In t roduct ion Incorporating sophisticated models of solvent molecules into calculations of the interac-tion 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 simplifica-tion 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 HNC equation used to calculate the macroion-macroion potential of mean force. An important step was taken by Kinoshita et al. [66] who examined potentials of mean force between a pair of spherical macropar-ticles 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 homoge-neous 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. At 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)/S(a*, a*), where the indices <; and s stand for the 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 as^0 for the counterions because of 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)] , (4.2.1) where n is the number of counterions fixed by the wall charges, and Wn({R}; a f , as = 0) represents the "effective" potential energy of n counterions (including "effective" ion-wall potentials) at infinite dilution in the solvent at activity aq [compare with Eq. (D.16)]. If one substitutes the n-particle potential energy (V„ + Un) for W n , the RHS of Eq. (4.2.1) 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 n is set to Wn{RlfR\) = it W(Ki) + E E w(Ri> Rj) > . (4-2.2) i—l i = l 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) wa7(R) = -kBT\ngay(R) , (4.2.3) where ga-y{R) is the ion-ion radial distribution function obtained in the molecular solvent. 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. Na + ) 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 wai. 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 wai(R), while ignoring the effects of a molecular solvent on Wa(z). Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 55 4.3 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 solvent-averaged pair potential of mean force which can be conveniently expressed in the form wa,(R) = u™(R) + < ( P ) + v*(R) , (4.3.1) where ue^(R) represents the contribution from the molecular solvent and is not present in the P M 1 . This "effective potential" will depend on the solvent model employed. Here, we 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. At 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 wa^(R) for monovalent ions with equal diameters d+ — 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 1We 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). Chapter 4.. Molecular Solvent Effects at the McMillan-Mayer Level 56 10.0 5.0 ? 0.0 h — 1 MM +- -- \ MM ++&-- -— — \ PM + - -- - \ . PM + + & - - • - 1 1 1 i 1 • 1 , 1 i 1.0 1.5 2.0 2.5 R / d 3.0 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 Level 57 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 way) and the wai 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 ue^y(R) to be consistent with the correct dielectric constant. An important feature (discussed below) is the minimum in waj at RtzlAd. We also note that the solvent-induced contribution u^a(R) has a short-range attractive tail and hence waa 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^s(z) and v^(z), respectively. 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 • (4.3.2) Peff is calculated by using ueS(R) as ush(R) in the formula for Psh [see Eq. (2.7.3)]. 4.4 Resul ts 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/d3, and with ions and solvent of equal diameter d. 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 Level 58 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 w a l i and the higher the wall charge density the better the fluid structure (and therefore Psiit) is approximated by the infinite dilution result. However, it must be noted that increasing the bulk density increases Pbuik- Therefore, because Psi;t is less 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 PM case, except when the ions are very close to contact (P<1.2d) . 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 n e t [MM] to be less than Pnet[PM] (note that for the infinite dilution case Pbuik = 0 and -Pnet — -fslit) • We first discuss a system of monovalent counterions with the wall charge density cr = -0.1307e/d 2 = -0.267 C / m 2 = - l e / 6 0 A 2 . The density profiles are not shown because 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 PM result, but that the wall-wall interaction Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 59 E-CO PH 0.08 0.06 0.04 0.02 0/00 0.10 0.05 0.00 -0.05 -0.10 MM, PNET PM, PNEL &eroe*%% 0 / o MM, PEFF A PM, fcore JL 6 8 d w a l l / d Figure 4.2: P n e t and its components for systems including only monovalent counte-rions with o - -0.1307 e/d2 = -0.267 C / m 2 . The Poisson-Boltzmann (PB) curve is shown for comparison (see Appendix C), where dw an[PB] is redefined to give the same accessible space (i.e., the same z m a x ) as the finite size ion models. P c o r e [ M M ] is less than 2 x 10~5 ksT/d3 throughout and is therefore not shown. For comparison: 0.020 kBT/d3 = 1.5 MRT = 3.7 x 106 N / m 2 , where R is the gas constant. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 60 is attractive (P n et < 0) for 3.4 d < dw an < 5 d. A comparison of the different pressure components from Eq. (4.3.2) demonstrates that this behaviour is solely due to the new pressure component Peff with P k i n and P e i being basically unchanged from the P M case. One has Peff < 0 as long as ions interact across the midplane mainly through the attractive part of uf+{R). The minimum of Peff at d w an « 3.6 d corresponds to a configuration in 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 Peff[MM] RS P c o r e [PM] and moreover P n e t [ M M ] « P n e t [ P M ] . The corresponding picture for a system with only half the previous wall charge (a = —0.0653 e/d2) is given in Fig. 4.3. We note that Pes is still negative for most wall separations, but much smaller in magnitude such that now P n e t >0 for all <iwan. For the rather high wall charge o — —0.1634e/d2 = —0.334C/m 2 (achievable on mica sheets), P n e t [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 w an = Ad. The results are shown in Fig. 4.5 and reveal that P n e t [MM] becomes more and more negative with increasing \o\, again due to the pressure component Peff. Interestingly, P n e t [PM] remains comparatively constant over a wide range of a. These observations 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 Level 61 E-0.08 0.06 0.04 0.02 0.00 CO 0.10 0.05 0.00 -0.05 -0.10 \ \ \ • \ \ MM, Pnel PM, Pnet PB, Pnet v. MM, r^ in x PM, Pkin MM, Pel o PM, Pel * — * — * - H MM, Peff A PM p — — i i»i, i core 4 6 d w a l l / d 8 Figure 4.3: P n e t and its components for systems including only monovalent counterions with o= -0.0653 e/d2 = -0.134 C / m 2 . The Poisson-Boltzmann (PB) curve is shown for comparison (see Appendix C). P c o r e [ M M ] is less than 2 x 1 0 - 6 ksT/d3 throughout and is therefore not shown. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 62 03 CO CU 0.06 0.04 0.02 0.00 -0.02 0.10 0.05 0.00 -0.05 -0.10 MM, Pnet PM, Pnet MM, riin- x PM, PK M MM, PEL o PM, PEL - -/CPAAAAA MM, PEFF A > PM P i xu, i core 4 6 d w a l l / d 8 Figure 4.4: P n e t and its components for systems including only monovalent counterions with o = -0.1634 e/d? = -0.334 C / m 2 . P c o r e [ M M ] is less than 4 x 10~5 kBT/d3 throughout and is therefore not shown. 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\ds/e Figure 4.5: The o dependence of P n e t and Peff for systems including only monovalent counterions at fixed d w an = 4 d. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 64 0.10 0.00 -CO DH « • • 9-^* MM, It,, x PM, Pkin MM, Pel o PM, Pel MM, Peff A 4 6 dwai i / d 8 Figure 4.6: P n e t and its components for systems including only divalent counterions with o = -0.1307 e/d2 = -0.267 C / m 2 . The Poisson-Boltzmann (PB) curve is shown for comparison (see Appendix C). P c o r e [ M M ] is less than 10~ 1 7 kBT/d3 and P c o r e [ P M ] is less than 4 x 10~4 kBT/d3 throughout, therefore these curves are not shown. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 65 reasonable waj(R) for divalent ions by recognizing that in a linearized version of the theory (the L H N C approximation) wai(R) scales to a very good approximation as \ZaZ7\ [73]; this can be interpreted as the attainment of structural saturation. Then, since for monovalent ions the potentials of mean force given by the HNC and L H N C theories are very similar, we take W2+2+{R) = 4w++(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/d2 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 n e t [ M M ] is even more negative and has an extended attractive range due to the attractive component Peff. Also, we note that unlike the monovalent case, ion-ion 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 Peff is now at dwa\\Xi2.8d, where the favorable interaction via ue±+ 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), QA, see Appendix B.2 . One then obtains the pressure by taking its (negative) derivative with respect to d w an. This method has the conceptual advantage that the derivation of the expression for QA Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 66 E-PQ CO CO 0.10 0.05 0.00 h contact plane p=1.0M midplane n -aQ A /ad w a „ -0.05 \--0.10 only contact plane counterions midplane x A ^ — a — a — a • • d • _ n D ° O • Ijl D a 4 6 dwal l / d 8 Figure 4.7: Comparison of P s i i t calculated via the contact theorem and via the midplane formula at o = —0.1307 e/d2. Results are shown for a 1:1 electrolyte solution in equi-librium 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 w an is also shown. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 67 exploits the HNC approximation, and therefore the level of approximation involved is consistent with that of the correlation functions. Moreover, f i A depends "symmetrically" 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 —dftA/ddwa,\\. 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 c o r e and _Peff depend strongly on gai at, and close to ion-ion contact, respectively. From P M calculations using anisotropic integral equations it is known that the HNC 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 Ps\\t calculated via the contact 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 dw an, the midplane concentrations reach the bulk value p faster than for t r^O. For the 1:1 electrolyte, we find j± = 0.38, p = 0.0136/d3 = 1.03 M for the M M case, and j± = 0.58, p = 0.0132/d3 = 1.00 M for the P M . Some sample density profiles for these systems are shown in Fig. 4.8 (Note that due to the symmetry of the system pa(z) = 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 w an we observe a higher degree of depletion of ions from the slit in the M M model than in the P M . Corresponding density profiles for a wall charge density of —0.1307 e/d2 (in the fol-lowing discussion we shall refer to this value as (Thigh) 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 ion-ion 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 w an }t 4 d but always remains lower than the corresponding counterion density. For dwa\\> 10d we find p+(0) « P-(0) and this value 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 pQ(0) curves are very similar to the 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 dWaii> 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 S-O.OOS -MM PM 0.000 -3.0 I i -2.0 z / d -1.0 0.0 Figure 4.8: Density profiles for 1:1 electrolyte solutions between uncharged walls, p is 1.00 M = 0.0132/d3 and 1.03 M = 0.0136/d3 for the P M and the M M model, respectively. p is indicated by horizontal lines. Results are shown for d w an = Id, 2d, and 8d. The profiles for cations and anions are identical. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 70 CO T l N o Q. 0.3 0.2 -4 0.1 0.0 0.015 h 0.010 \-0.005 r-0.000 -2.5 -2.0 1.5 -1.0 z / d 0.5 0.0 Figure 4.9: Density profiles for 1:1 electrolyte solutions between charged walls with cr = —0.1307e/d2 = —0.267C/m 2 for (a) counterions and (b) coions. p is LOOM = 0.0132/d3 and 1.03M = 0.0136/d3 for the P M and the M M model, respectively. p is indicated by horizontal lines. Results are shown for dw&\\ = 2.5 d, 4.0 d, and 6.0 d. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 71 0 . 0 2 5 0 . 0 2 0 0 . 0 1 5 ^ 0 . 0 1 0 -(a) 0 . 0 2 0 0 . 0 1 5 0 . 0 1 0 A • . ^high A counterions A crlow * A A <7high D A A A coions • A • • A a1( a D n A A 6 . ow a _PJ_ 1 J I I I I I L (b) \\ counterions a h i g h \ * ^low 2 4 6 8 10 12 14 d w a i i / d Figure 4.10: Midplane densities of monovalent counterions and coions at fhigh = -0.1307 e/d2 and c r l o w = ahigh/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. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 72 - 0 . 0 2 1 1 ' « ' ' ' 1 ' 1 2 4 6 8 10 d w a n / d Figure 4.11: P n e t for 1:1 electrolyte solutions at p = 1.0M for a h i g h = -0.1307e/d2, o"iow = Chigh/2, and for uncharged walls. 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 our CTIOW and 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 n e t for crnigh is shown in Fig. 4.12. Curiously, 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/d2) this trend is reversed, as shown i n Fig. 4.13 where we compare the variation of P n e t 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/d2. 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 counterion-counterion 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 w an. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 74 Figure 4.12: P n e t for M M systems in equilibrium with 1:1 electrolyte solutions at various bulk concentrations and cr = -0.1307 e/d2 = -0.267 C / m 2 . The curves for p = 0.5 M and 1.0 M are indistinguishable within numerical error and therefore only the latter case is shown. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 75 0 . 0 2 0 . 0 0 0 . 0 5 0 . 1 0 0 . 1 5 0 . 2 0 0 . 2 5 k|d2/e Figure 4.13: The dependence of Pnet on o at fixed dwan = Ad. Results are shown for 1:1 electrolyte solutions at p = 1.0 M and for the the infinite dilution (counterions only) case. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 76 4.4.3 Pair Distribution Functions The anisotropic pair distribution functions, #Q7(r, Zi, z 2 ), contain information about the 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 2) is 27rrp7(z2)<jrQ7(r, z\, z 2 ) . As gai depends on three variables it is difficult to depict in full and here we select a few examples to illustrate some general features. <?a7(r, 0,0) (i.e., at the midplane, parallel to the walls) for a 1.0 M , 1:1 electrolyte in a slit with d w an = 4d is shown in Fig. 4.14. Whereas in the P M case all functions (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 3 around 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 w an as fewer coions are surrounded by more counterions. The coion-coion 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: gai{r) at the midplane parallel to the walls (zi = z2 = 0) for 1:1 electrolyte cr = —0.1307 e/d2. g+_=g_+ is shown in (a) and g++ and g ("+" denotes counterion and "—" the coion) are shown in (b). Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 78 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 gaj parallel and perpendicular to the walls are surprisingly similar and this is especially true for cation-anion pairs. Figures 4.15 and 4.16 show g++ for monovalent and divalent counterions, respectively, at the wall, with c?wan corresponding to the minima of Peff in Figs. 4.2 and 4.6. The white region around the particle is forbidden to other ions because of the hard-core 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 In order to interpret the contour plots in terms of average densities, gai must be multiplied by the density profile of ion 7, py{z). 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 w an « Ad. Finally, Fig. 4.18 shows the ordering of coions at preferred 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). separation R—Jr2 + (z + z m a x ) 2 ~ l . Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 79 g(r,z) for counterion/counterion MM -6.0 I •4.0 I -2.0 I 0.0 r +2.0 +4.0 +6.0 g(r,z) for counterion/counterion PM -6.0 I -4.0 -2.0 l 0.0 r ^ ^ ^ ^ +2.0 +4.0 +6.0 1.1500 m 1.0500 • 1.0100 • 0.9900 r ~ | 0.9000 • 0.8000 • • 04000 • 0.0000 • Figure 4.15: Contour plots of g++(r, -zmax,z) for systems where only monovalent coun-terions are present with d w a u = 3.6d and o = -0.1307e/d2. 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. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 80 g(r,z) for counterion/counterion MM • 1.1500 ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ ^ P ^ r ^ ^ ^ ^ ^ ^ ^ ^ ^ T ^ 0.0000 -6.0 -4.0 -2.0 0.0 +2.0 +4.0 +6.0 D Figure 4.16: Contour plots of g++(r, -zmax, z) for systems where only divalent counterions are present with d w a n = 2.8d and o = -0.1307 e/d2. 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)d3, the average counterion density around a midplane coion. The ions are monovalent and p=1 .0M, dwa,n = 4d, and a = -0.1307e/d2. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 82 [g(r,z) dens(z)J for counterion/coion MM [g(r,z) dens(z)J for counterion/coion PM 0.0500 • 0.0200 • 0.0100 • 0.0080 0.0060 • 0.0040 • 0.0020 • 0.0010 • 0.0000 • Figure 4.18: Contour plots of p_(z)p+_(r, -zmax, z)d?, the average coion density around a counterion in contact with one wall. The system parameters are as in Fig. 4.17. Chapter 4. Molecular Solvent Effects at the McMillan-Mayer Level 83 4.5 Summary The force between two like-charged walls immersed in electrolyte solution can be cal-culated 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 mag-nitude 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 w a i i . Nevertheless, the fluid structure around the coions that are present is drastically altered compared to the P M . Chapter 5 Further Investigations at the McMil lan-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. An 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 McMil lan-Mayer 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 NaCl 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 Level 85 5.1.1 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 ion-wall 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; o n and hard walls at a separation d w an (as in Chapter 4) Z) = d w a n — dj o n . 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 SPC (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 are set to zero for R>7.$A. The label P R denotes results obtained by Pettitt and Rossky [76] using RISM theory. M D simulations have shown that the RISM 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 RISM calcula-tions 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: ueff(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. Chapter 5. Further Investigations at the McMillan-Mayer Level 87 a / (kBTA) b / A~l c / {kBTA) d / 1 - 1 e / A " 1 f Na+/Na+ 6.95 • 1012 -8.685 7.96 -0.3172 3.01 -2.71 c i - / c r 770. -1.247 657. -0.9618 1.96 -3.94 Table 5.1: Fit parameters for the Pettitt-Rossky [76] ion-ion potentials of mean force. The units are chosen to yield u e f f in kBT if R is given in A. results) by the function uef!(R) = a exp(bR)/R + c exp(dR) cos(eR + f)/R , (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 e f f is only slightly deeper than the corresponding one for the N a + / N a + [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 way(R) =oo if waj(R) ^ 10 kBT. We stress that although the exact physical location of the walls, and 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 88 5.1.2 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 2 = le /48A 2 . The walls may be negatively or positively charged depending on 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.8kBT at R = 3.6 A (see Fig. 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 Peff. At small D, Pnet becomes . positive partly because the ion layers at the planes of closest approach to the walls interact via the repulsive part of w e f f. Also, the component Pki n grows with decreasing slit width as the average ion density between the walls increases. Results obtained for N a + [PR] counterions are also shown in Fig. 5.2. For the P R model ue++ has a shallow first minimum of ~ —0.6 kBT followed by a slowly decaying oscillatory tail. The resulting pressure component P eff 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 n e t 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 2 = - l e / 4 8 A 2 . M represents moles/litre and-R is the gas constant. Chapter 5. Further Investigations at the McMillan-Mayer Level 90 Figure 5.3: The dependence of P n e t 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 2 = +le /48A 2 . The notation is as in Fig. 5.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 P eff (and hence a repulsive Pnei) as the walls are brought together. In contrast, the midplane density and P^n remain rather constant as D decreases. In the C l ~ [PR] case, «!_ff_ has a strong attractive region with a mini-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 ueS 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 (Na + [PR]) or no attractive well at all (PM 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 ueS. Strong oscillations in u^a produce more pronounced oscillations and higher peaks in 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 gaa(r) is less oscillatory and lies closer to the P M result. 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 sepa-ration £> = 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 o n = 2.94A. The notation is as in Fig. 5.2. 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\ = z2 = 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 0 n = 2.94A (PM) at the same wall charge are also shown. The C l - [PR] curve reaches a maximum of 8.3. Chapter 5. Further Investigations at the McMillan-Mayer Level 94 the present discussion is restricted to one model. P n e t curves for the NaCl [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><13A at 1.0 M), and the repulsive tail found for £>>13A 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 e f f 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 e f f . These in turn can give rise to different ion "structure" between the walls and to 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 ueS does not necessarily reduce P n e t significantly. For example, E q + [KP] and N a + [PR] have comparable attractive wells, but E q + [KP] leads to an attractive wall-wall interaction whereas N a i [PR] does not. Chapter 5. Further Investigations at the McMillan-Mayer Level 95 0 5 10 15 20 D / A Figure 5.6: P n e t for the NaClJGRP] model [74,75] at different salt concentrations, p, and o = -0.334 C/m 2 = - l e / 4 8 A 2 . The inset shows the tails on an expanded scale. The notation is as in Fig. 5.2. Chapter 5. Further Investigations at the McMillan-Mayer Level 96 Figure 5.7: Comparison of P n e t for four systems with monovalent counterions, including the P M with d ion = 2.94.A, and o = -0.334 C / m 2 = - l e / 4 8 l 2 . 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. Chapter 5. Further Investigations at the McMillan-Mayer Level 97 Even calculations with a toy model that has 1.8 times the N a + [PR] w e f f (the depth of the first well is —1.0 ksT) still yield a positive P n e t 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 e f f 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, uef[ does display a strongly attractive well and/or an attractive tail [52,74,30]. Thus, solvent-mediated counterion-counterion interactions could lead to attractive wall-wall forces, and this possibility should be considered seriously in the interpretation of experimental observations for negatively charged mica sheets in aqueous NaCl 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 Level 98 here. 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, wai(R), calculated by 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, Wa(z). This "effective" ion-wall interaction can be written in the form [see Eq. (2.5.2)] Wa(z) = VUS(z) + V*(z) + V Q ( < W + Z) + MZmzx ~ z) , (5.2.1) where tpa(C) is the solvent-induced short-range interaction of an ion a with one wall as a function of the distance to the plane of closest approach, £ > 0. By using the ansatz in Eq. (5.2.1) we implicitly assume that both walls give identical superimposable contribu-tions and that ipa does not depend on the wall separation itself. These assumptions are reasonable for large wall separations but are clearly an approximation at short d w an. A good model for ipa (() is difficult to choose. 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]. At 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 im-mersed 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 solvent-averaged 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 sol-vent 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 calcu-lated 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 dwaii- In addition, the superposition ansatz made in 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 prob-lems 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 concentra-tion (0.1 M) [83] to the same macroion surface charge density as used in our calculations (le/60A 2). Then, only the part close to contact with the wall has been retained, rep-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 signif-icantly. 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 w a n [see Eq. (2.7.3)], Chapter 5. Further Investigations at the McMillan-Mayer Level 101 E-m + 0.5 0.0 0.5 1.0 1.5 2.0 2.5 \ . \ - \ f l ^ ' / / / / 1 \ 1 1 < 1 < i , 0.0 0.2 0.4 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 ( = (zmax ± z) > 1 d and for coions at all distances ijja = 0. For x = (,/d < 1 the curves are given by the polynomials 2.75a;3 - 7.75a;2 + 7.25x - 2.25 (ipl) and 71.bx5 - 145.9a;4 + 91.0a;3 - 9.3a;2 - 7.7x + 0.4 Chapter 5. Further Investigations at the McMillan-Mayer Level 102 and one has •Pslit = -fkin + -Pcore + Pel + -Peff + Pwall • (5.2.2) Note that Pwai i^O if 0+(£) does not act across the midplane. This is true for our model potentials if d w an > 3 d. 5.2.2 Results We present results for systems with monovalent ions at T = 298 K, e = 78.8, and a = —0.1307 e/d2 = —0.267 C / m 2 in Fig. 5.9. 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 w an, the potential "02" yields more attractive (less repulsive) values. Differences between the "^1" and uip2" cases come about because ions experiencing 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 .25 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 Pnei is more attractive than for the infinite dilution case (compare with Fig. 4.12), while at reduced d w an the pressure oscillations are slightly increased in magnitude. The extremely 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 E-DQ CO 0) e 0.00 H -1 6 7 d w a n / d Figure 5.9: Pnet for monovalent counterions and a wall charge a = —0.1307 e/d2, with 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. Chapter 5. Further Investigations at the McMillan-Mayer Level 104 is « —0.5d3/kBT) 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 addi-tional 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 Pe\[ipl]> Pe\[tp = Qi\ 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 Peff['0 = O] but the minimum is shifted to smaller rfwau, and there is an additional oscillatory tail. As long as P w a i i ~ 0 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 at the McMillan-Mayer Level 105 1.0 2.0 3.0 4.0 d w a n / 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/d2. 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. 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/d3. 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 experi-mentally 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 macroparti-cles). 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 Level 107 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 2 ) are achieved. With increasing bulk concentration 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 1. Only for 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 pre-dictions from DLVO theory (see Appendix C) for large surface separations and/or low (<10~ 5 M) salt concentrations. For surface separations ^$20 A and medium to high salt 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 1Schoen 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 approxi-mation, that this quantity is related to the net pressure between a corresponding system of two planar walls at a separation h via -27rPnet(/i) = (d/dh)[F(h)/R}. 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, prob-ably 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 homoge-neous HNC 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 concentrations2. Furthermore, their ob-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 fluctu-ations 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 2In 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. 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 fluctua-tions 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 kBT). 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 be-yond 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 s i j t , can be obtained either from the contact theorem or via an equivalent "midplane formula". The latter route divides Psnt into several pressure components, representing particle interac-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 Conclusions 112 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 n e t , is given by the difference between P s l i t and the outside bulk pressure. P n e t is a quantity of major interest, as it indicates whether the walls (macroions) show a net attraction ( P n e t < 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 PM) . 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 e i , 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), wa7(R), 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 Conclusions 114 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" ion-wall 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 waj(R) were calculated before by Kusa-lik and Patey. Initially the solvent-mediated wall potentials are ignored in order to in-vestigate 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, Pes, arises, which is due to solvent effects on the effective 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 Peg can be very important, leading to an attractive net wall-wall force at small separations (;$ 15 A), even for monovalent counterions if the wall charge den-sity is sufficiently high (see Fig. 4.5). With rising salt concentration in the bulk the Chapter 6. Summary and Conclusions 115 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 em-ployed 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 Conclusions 116 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. Hi l l . Statistical Mechanics. McGraw-Hill, 1956. [20] L. S. Ornstein and F. Zernicke. Proc. Akad. Sci. (Amsterdam), 17:793, 1914. 117 Bibliography 118 [21] H . L. Frisch and J. L. Lebowitz, editors. The Equilibrium Theory of Classical Fluids. W. A. Benjamin, Inc., 1964. [22 [23; [24; [25 [26; [27 [28 [29 [30 [31 [32 [33 [34 [35; [36 [37; [38 [39 [40 [41 [42; T. Morita and K . Hiroike. Prog. Theor. Phys., 25:537, 1961. P. Attard and G. N . Patey. J. Chem. Phys., 92:4970, 1990. J. S. Rowlinson and B. Widom. Molecular Theory of Capillarity. Clarendon Press, 1989. H. Greberg and R. Kjellander. Mol. Phys., 83:789, 1994. R. Kjellander and S. Sarman. J. Chem. Phys., 90:2768, 1989. J. C. Rasaiah, D. N . Card, and J. P. Valleau. J. Chem. Phys., 56:248, 1972. M . Plischke and D. Henderson. J. Chem. Phys., 88:2712, 1988. R. Kjellander. J. Chem. Phys., 88:7129, 1988. Erratum - 89:7649. M . Vossen and F. Forstmann. Mol. Phys., 86:1493, 1995. M . Plischke and D. Henderson. Proc. R. Soc. Lond. A, 404:323, 1986. R. Kjellander and S. Sarman. Mol. Phys., 70:215, 1990. R. Kjellander and S. Sarman. Chem. Phys. Lett, 149:102, 1988. R. Kjellander and S. Marcelja. J. Chem. Phys., 82:2122, 1985. R. Kjellander and S. Marcelja. Chem. Phys. Lett., 112:49, 1984. Erratum -114:124. M . Lozada-Cassou, W. Olivares, and B. Sulbaran. Phys. Rev. E, 53:522, 1996. D. Henderson, P. Bryk, S. Sokolowski, and D. T. Wasan. Phys. Rev. E, 61:3896, 2000. M . Lee and K . - Y . Chan. Chem. Phys. Lett:, 275:56, 1997. R. Kjellander and S. Marcelja. Chem. Phys. Lett, 127:402, 1986. B. R. Svensson and C. E. Woodward. Mol. Phys., 64:247, 1988. J . P. Valleau, R. Ivkov, and G. M . Torrie. J. Chem. Phys., 95:520, 1991. R. Kjellander, T. Akesson, B. Jonsson, and S. Marcelja. J. Chem. Phys., 97:1424, 1992. Bibliography 119 [43] H. Greberg, R. Kjellander, 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. J. Phys. Chem., 92:6489, 1988. [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. J. Chem. Phys., 107:3611, 1997. [64] R. M . Shroll and D. E. Smith. J. Chem. Phys., 111:9025, 1999. [65] A . Trokhymchuk, D. Henderson, and D. T. Wasan. J. Colloid Interface Sci., 210:320, 1999. Bibliography 120 [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. J. Chem. Phys., 92:4399, 1990. [90 [91 [92 [93 [94 [95 [96 [97; [98 [99; [100 F. Lado. J. Comp. Phys., 8:417, 1971. M . J . Gillan. Mol. Phys., 38:1781, 1979. K . - C . Ng. J. Chem. Phys., 61:2680, 1974. M . Kinoshita and D. R. Berard. J. Comp. Phys., 124:230, 1996. W. H . Press, S. A . Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in FORTRAN. Cambridge, 1992. E. J-. W. Verwey and J. T. G. Overbeek. Theory of the Stability of Lyophobic Colloids. Elsevier, 1948. S. Engstrom and H. Wennerstrom. J. Phys. Chem., 82:2711, 1978. J. C. Neu. Phys. Rev. Lett, 82:1072, 1999. B. V . Derjaguin and L. Landau. Acta Phys. Chim. USSR, 14:633, 1941. R. Kjellander and D. J. Mitchell. Chem. Phys. Lett, 200:76, 1992. 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 Nz layers with widths AZJ. i,j = l , . . . , i V 2 are now indices for the individual layers, representing the z-value in the middle of the layer, zi: and the particle type ojj. The wall 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 Nr grid points, r^, with a cutoff r^r = r m a x , beyond which we set hij(r) = 0 and Cy(r) = ~Pui){yr2 + (zi ~ zj)2)i the latter choice equating c^ - with its asymptotic value. These long-range tails are handled analytically throughout the calculation. Typical parameter values used in the numerical calculations are Nz = 60 — 160, Nr « 300, and 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 [gaiaj(r, z^, Zj)= gai0li{r, Zj, zA] one can reduce the required memory considerably (by up to 75%), but a complicated system of indices must be introduced. 122 Appendix A. Numerical Algorithm 123 The OZ equation is conveniently solved in Fourier space. The Fourier transforma-tion in two dimensions, called Hankel transformation (HT), for spherically symmetric functions [f(r) = /(|r|)] becomes poo f{k) = 2?r / drr f(r)J0{kr) ( A . l . l ) Jo f(r) = ±- rdkkf(k)J0(kr), Z7T JO where Jo(x) is the zeroth order Bessel function. Applied to Eq. (2.6.1) the Hankel transformed OZ-equation reads hon(k,zuz2) = ca7(k,zi,z2) + dz3 caS{k, zi, z3)ps(z3)hSl{k, z3, z2) , (A.l.2) <5 J 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/Civr, where Q is the i-th positive root of Jo(x). The corresponding grid points in Fourier space are ki = Cli/rm&y. with a cutoff A ; m a x = CAT,./Vmax- The formulas for the discrete transformations are given by = TJ- E fciW* ^ d (A.1.3) "'max j=l = Z^-'1LHki)Wii with (A.1.4) 'max j=l W _ JQ(CJ(J/(NT) 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.5) 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 ) " 1 = (R - R C R ) - 1 - R 1 , (A.1.6) where l = {5ij}. At hard-core contact r = a^ = \J\{da^) + da^)2 — (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) = Actj + ^ ^ ^ j (A.1.7) to the function in real space for r < a (inside the core) and subtracting the Hankel transformed part Qij(k) = 2TTaij[AcijJ1(aijk)/k - AcijJ2{aijk)/k2) (A.1.8) after performing the numerical HT, where Jn(x) is the n-th order Bessel function. Acij = Ahij and Ac'^ = Ah'^ can be determined by using hij(r < a) = — 1 and interpolating the continuous (even at hard-core contact) function yij(r) = hij{r) - c^r) , (A. 1.9) as usually no grid point lies at r = a. With the HNC closure, hij(r) = exp{—ftuij(r) + Uij(r)} — 1, one finds Aci:j = c°f(aij) - <%-(aij) = expi-Buijtdij) + yy(ay)} , (A.l.10) and Ac'^ is given by a simple numerical derivative of Eq. (A. 1.10). Appendix A. Numerical Algorithm 125 The actual iteration between the OZ-equation and the HNC closure is performed with the short-ranged functions 4 = cv + Pv-ij ' a n d Vij = VU ~ Pufj • (A. 1.11) ufj and its Hankel transform are given by fig(r) = , and (A.1.12) 4ire0e^r2+ b2j ^ { K ) = M M 2^M-Kk} w . t h 6ij = max (bmin,^(zi - Zj)2^ , where & m m > 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 £ jdz2Pl(z2) Jdrrpui^ + iz.-z,)2) , (A.1.14) Cfdwall + £ <?7 dz2\z1- Z2\p7{z2) 2ee0 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 C y 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 Algorithm 126 HT OZ-equation A A c(r) — y(r) initial guess c(r),p if As. y HT <8 X HNC closure y(r) c(r) density equation electroneutr. condition newp if Ap <&2 Figure A . l : Iteration scheme for the algorithm. Appendix A. Numerical Algorithm 127 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 ° u t —x^i, the system of equations is iterated until for successive iterations Ay/y<ei^ 1 0 - 3 and A p / p < £ 2 ~ 1 0 - 4 . 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 0 , J\, J 2 ) , root finding with Ridders' method, 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,...,Nz and I = 1,Nr), 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 128 A.2 M M Systems w i t h Strong 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 pe*(z)- Moreover, the spiked nature of these functions in real space expends the "tails" in Fourier space and a much finer grid in each layer becomes essential. This is because in the orthogonality-preserving HT a smaller grid width in real space corresponds to points at longer wavelengths in Fourier space. In the present calculations Nr = 1600 was found to be necessary for r m a x = 19.6A. The like-ion functions y*aa{k, z\, z2) decay as Ca(zi, z2)/kz. This can easily been seen by expanding a slightly modified Eq. (A. 1.6) using 1/(1 — x)zxl + x + x2 + ... Y = C ( 1 - R C ) - 1 - C « C R C + ... , (A.2.1) where Y(k) — {yij(k)}. As a simple example we consider a calculation with only one layer for each ionic species, then V— y-+ c c-+ ' P- o " c c-+ . y+- y++ . C++ . . 0 P+ . . a+- C++ _ (A.2.2) p-c2__ + p+c_+c+_ p-c-+c—+ p+c_+c++ _ p_c +_c__ + p + c + + c + _ p_c+_c_+ + p+c2++ For large k the only surviving contribution to c + _ = c_ + is the term t proportional to Ji(a+-k)/k in Eq. (A. 1.8), which decays as fc-3/2, whereas C++ and c are essentially zero (as A c + + and A c are very small). From this and Eq. (A.2.2) it follows that the long-range tails of y and y++ decay as k~3 and that the corresponding functions for Appendix A. Numerical Algorithm 129 unlike species are of much shorter range. These tails are significant only if z\ ~ z2, i.e., 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 are included in the HT. The required integrals dkJ0(kri)/k2 are evaluated once for 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 (Nz = 202) was used, 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. (A.2.3) Appendix B Formulas for P S H t Three different routes to evaluate P s i i 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 a2 r^max dvsh>I(z) Psiit = kBTp(-zm&x) - - / dzp(zf 1 ' . (B.1.1) ZCCQ J - z m & x OZ Using the relationship J^_ZmAxdz(dp/dz) = p(0) — p(—2max) and a simplified version of the B G Y equation, Eq. (2.3.4), k s T ^ = - P W ^ " / * r 9 ( r , * , * ) * 4 ^ , (B.1.2) we obtain zeeo J - z m & x oz J - z m s u x Oz n f° , , . /"2max , . , r°° , , .du(r,zi,z2) +2TT / dzlp(zi) dz2p(z2) drrg(r,zl,z2) -z . 130 Appendix B. Formulas for P S H t 131 The first term on the RHS is called Fkin, whereas terms three and four give PWaii- As (dvel/dz) = (dvBS /dz) = 0 (for the integration range) one finds /-.max QyShJf^ fQ dlv^'1 (z) + V^'" (z)} Pwaii = - / d z p ( z ) — 7 ^ + / dzp(z)-+ K-» /•».., dif^'iz) r« . . .dv'h-"(z) - *™-arL+Lj**')—5rL- (B'L4) Both terms in the expression for P w a u give an identical contribution if vsh'T (z) =vsh'n (—z) (as in our model). The last term on the RHS in Eq. (B.1.3) can be simplified by splitting the integral over z2 into two terms, A and B, yielding A + B = 2TT f dzlP(Zl) fZ™dz2p{z2) rdrrg(r,zuz2)du{r'Zl,Z2) (B.1.5) •/-Zmax JO JO <OZ\ c f ° i i N f ° , / \ f ° ° , / sdu(r,zi,z2) + 2?r / dzip{zi) dz2p{z2) drrg(r,zuz2) . J-Zm&x J - Z m & x JO 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\, z2) = g(r, z2, z\)1. The particle pair potential depends only on the separation R = \Jr2 + (z\ — z2)2, which means that du du dR du z\ — z2 du . . d~z[ ~ d~R ~dz[ = d~R R = ~d~z2~ ' (••) After exchanging z\ and z2 in (B.1.5) the term B can be rearranged to the original form, 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 z2 and the additional symmetry operation Z\ —> —zi and z2 —»• — z2 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, Z2,Z1). Appendix B. Formulas for P s i i t 132 z = 0 is a symmetry plane for our model one has p(z) = p(—z) and g(r, —zi, —z2) = g(r,zi,z2)). Performing only the latter transformation on term A and splitting up u into its components, Eq. (B.1.3) becomes C T 2 -Psiit — -Pkin + -Pvall ~ (B.1.7) z e Co r z m a x . . r° roo d[uRS + uel + ush] f Z m a i rv roo -1-E \ dz1p(zi) / dz2p(z2) / drrg(r,zi,z2) JO J— Z m a x •'O -Zmax •'U 8Z\ PSh is given directly by the part of the integral involving ush. Combining the electrostatic terms defines P ei and yields r 2 rZm"^ f° ^ ^ f ° ° „ ,du*\r,Zl,z2) P„i = 2ee 0 Jo — — J-z^'^^Jo ~" dz, - 2 T T r&XdZl p(Zl) f° dz2p(z2) rdrrh{r,zuz2)dueX{r'ZuZ2) . (B.1.8) JO J-zma.* Jo OZ\ The last identity follows from h = g — 1, duei a2 , , r°° rdr roo 2TT f°° , Cwcl cf . . f/ drr—— = - - — ( 2 1 - 2 2 ) 7 = 7 o ^ 2 e ° e 7 o ^ / r 2 + (21 - 2 2 ) -3 2 q2 zi - z2 2e0e | 2 i - 2 2 | (B.1.9) and /"Zmax f0 q dzp(z) = q dz p{z) = o , (B.1.10) JO J—zmB.x 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, PCOre- To simplify the expression for P c o r e we first define f = \Jd2 — (z\ — z2)2 (r-value for which the particles in layers 21 and 2 2 are in contact, R(f) = d) and one further finds S(R-d) = ^ e - ^ S = -(3e-^sd-^ ( B . l . l l ) OR dr - 1 S(r - r) = p(r - r) . (B.1.12) Appendix B. Formulas for Psnt 133 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(x0)). Using the last two equations we obtain / drrg{r,z1,z2)—— = rg{r,z1,z2)^—^(-kBT)-, B.1.13 Jo oz\ a r and finally Pcoie = 2TrkBT dzip(zi) dz2p(z2)(zi-z2)g(f,zl,z2) . (B.1.14) 0^ J—Zm&x Summing up all contributions yields Eq. (2.7.3) . B.2 Calculation of £lA The grand potential per unit surface area and the corresponding Helmholtz free energy, F A , for an electrolyte solution are related via QA = FA-p±jdz \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 fdzpa(z) + kBT—lnX, (B.2.2) where the definitions from Eq. (2.6.6) and the electroneutrality condition have been used, and a constant independent of o?wan has been dropped. X depends on dwan and is determined by solving Eq. (2.6.7). FA has been given, for example, by Kjellander and Marcelja [34], exploiting again the HNC approximation and thereby avoiding a numerical coupling parameter integration2. 2The excess part of FA is obtained via integration of ne*(z) over the whole slit system, and then it is added to the appropriate ideal part. Appendix B. Formulas for Psnt 134 After simplification of their formulas one finds FA (B.2.3) kBT 2ee0kBT 1 ^ Q, dzxdz2 pa{zi)p1(z2) 2irrdr dzxdz2pa(zx)p1(z2)\zx - z2 -— J dk kln[det{l + HR)] Z a J a J with H , R , and 1 defined following Eq. (A.1.5). Note that the terms J2Q fdzpa(z) in Eqs. (B.2.2) and (B.2.3) cancel and therefore the values p+ and p_ are not needed. B . 3 Remarks on the Numer ics The use of the contact theorem requires an accurate determination of the contact den-sities. 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 Psnt 135 fluids, where the only terms for the midplane formula are P^n and PCOre, we deduced that the integrals in jPcore must be evaluated with extra care. Overall, fluctuations in the Psi;t 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 fiA is comparable to the accuracy of P su t from the contact theorem and midplane formula routes, but the differentiation with respect to d w an has to be per-formed via a numerical finite difference operation, which increases the error considerably. This renders the third method of evaluating P s l i t not practical and here we used it merely for comparison. Appendix C P B and D L V O Theory Poisson-Boltzmann (PB) theory is the mean-field treatment for point ions in a di-electric continuum for an inhomogeneous system, which was developed independently by Gouy and Chapman around 1910. Compared with the exactly solved P M , at the PB level both ion size effects and correlated charge fluctuations are neglected. However, PB 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 pa(R) = p0aexp{-Pqa<p(R)}, * (C.2) where p°a is a constant (usually the bulk density of species a). Eqs. (C.l) and (C.2) must 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]) 2e0ekBTs2 qcrzmax ' W = « « c o s * ( W w ) W ' t h S t a n ( s ) = - ^ T ' ( C ' 3 ) 136 0 Appendix C. PB and DLVO Theory 137 Alternatively, one can obtain the P B pa(z) through iteration of Eq. (2.6.3) with only the mean field electrostatic contribution, Eq. (A.1.14), included, i.e., by setting c Q 7 (r) = - / 3 « Q 7 ( r ) and hai(r) = 0. The net pressure between like-charged walls obtained from P B theory •Pnet — Pkin ~ Pbulk — kBT pa{0) - lim pQ(0) « w a l l - > ° 0 (C.4) is always positive, which means that the double layer interaction is repulsive1. For w^aii -> co one finds [3] P n e t ~ exp{- /«Z w a l l } with K2=Zaf°p} . (C.5) K ~ 1 is called the Debye length and arises also as the screening length from Debye-Hiickel theory. A few examples of P n e t [PB] in comparison with P M and M M results are shown in Section 4.4.1 . Note, that for the P M two additional contributions to P n e t arise: P c o r e > 0 (due to ion-ion hard-core interactions) and Pe\ < 0 (due to correlated charge density fluctuations). As long as core-core interactions between ions are not important the former contribution can be ignored and therefore P n e t [PB] > P n e t [ P M ] . The "classical" theory for the interactions between colloids, the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory [95,98], includes also an attractive van der Waals in-teraction between the macroparticles. The corresponding pressure contribution for the slit system is (as a first approximation) [1] A vdW 127T<4all (C.6) where A is called the Hamaker constant. A depends on the fluid medium and on the nature of the walls. At small dwa\\ (a few Angstroms) P vdw + P*net < Oand one expects an xThis 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]. 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 DLVO 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 Hi l l [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)({n};a„as) = p^l^n{Ru...,Rn\a^as) = ^ T V T E T n T f \ f d { N - n } exp[-8(VN+UN)] , (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 UN are the A^-particle potential energies due to the external potential and particle interactions, respectively. Assuming that n ? and n s of the n indices in refer to solvent and solute particles, respectively, with n = n(;+ns, and making the substitution m a — N a — n a , Eq. (D.l) can be written in the form //n)({n};o t )o,)H(o t )a,) - "m<"m* E -£~ri / ^ [-Wm+n + Wm+n)] • (D.2) m>0 139 Appendix D. McMillan-Mayer Theory 140 Integration over d{n} on both sides yields A^r- d{n}pM({n};a<,as)= £ d{m+n} exp [ - / ? ( V m + n + UM+N)} . (D.3) We now carry out a two-dimensional Taylor expansion of E(a f , a s ) about the "reference acitivities" a* and a* for S(a ? ,a s ) to obtain S K . « « ) = E S / W ^ N ^ + W J V ) ] (D.4) v ( a . - a ^ K - q ; ) " - / ^ H ( a „ a s ) \ ^ n j n s ! v 9 a ^ a a ? « J a, = a* { 1 J a s = a* = ^ K - o ; ) ^ ( a . - a ; ) " « ^ 0 0 / " d { m + n } c - / » C v B ^ - w . + . ) i ( D . 6 ) n>o nq.ns. m>0 mq.ms. J where for the last step the partial derivatives in Eq. (D.5) have been evaluated using Eq. (D.4) and the substitution ma = Na—na has been made. Comparing Eqs. (D.3) and (D.6) one finds a relationship between the partition functions for the activity sets (o ?, as) and (a*, a*). One has SK, a.) = ~(a*, a:) £ ( a< " u^X^ /d<n> ^ " (D'7) Analogous to Eq. (2.2.10), the n-particle distribution functions are defined via g<*(H;at,a,) =f«">•«»-> . (D.8) nr=iPa i (^ ;ac ,a s ) We further introduce the n-particle potential of mean force by setting wW({n} ;a„a s ) = - A : B T l n ^ ( { n } ; a c, a s) . (D.9) The latter definition is motivated by investigating Eq. (D.2) in the dilute gas limit a c,a s—»0. In this limit only the term for m = 0 survives and we find p ( n ) ({n};0,0) = a » x p [ - / ? ( V n + Z 4 ) ] (D.10) Appendix D. McMillan-Mayer Theory 141 = n{p Q i ( J R i ;0 ,0)exp[^ Q i (P i )]}exp[- /?(V n +W n )] ( D . l l ) = f[Pai{Ri;0,0)exp[-pUn] , (D.12) t=l where S(0,0) = 1 [see Eq. (D.4)], limaa^0 pa(R) = aaexp[0va(R)], and V„ = E?=i have been used. Equation (D.12) yields the interpretation for the n-particle potential of mean force in the dilute gas limit, w^({n};0,0)=Un, (D.13) 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 pa(R; a c, as) via the (1-particle) potential of mean force Wa(R; a ? , as) [see Eq. (2.3.7)] one finds the McMillan-Mayer expression E(a*,a*s) ^ T j . . i n . i n>0 x Jd{n} exp -P E ^ ( ^ ; <, <) + W(")({n}; a*, a*s) \i=l / 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 (as = a+ in the following). From Eq. (D.14) we obtains for the case with a c = a* and a* = 0 E(a*as) ^ < s f . ns>0 "« -P £ WS(R\; a*, 0) + w^\{ns}; 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 concen-tration [set one of the aa to zero in Eq. (D.4)], with the difference that the (vacuum) 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 Ws and u>(n'). Of course, all changes in the solvent structure due to the presence of the 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 ns = — 2o/q+. A l l terms with ns ^ ns are zero because of the infinite electrostatic contribution in Yl Ws 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 . ° « ttS 1 , 7 1 - 1 . n s = n^" / d{ns\ exp s(o*,0) ns\ j d{ns} -(3 £ W8(Ri; a*, 0) + w^({ns}; a*, 0) (D.16) \i=l Equation (D.16) is exact if the contributions of to all orders are included. In practice, one introduces the pairwise approximation i=l j>i Appendix D. McMillan-Mayer Theory 143 where wai(Ri, Rj) is the pair potential of mean force. Unlike for the potential energy UN, 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 con-centrations 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 way(Ri, Rj) by its bulk counterpart wai(Ri, Rj) « w™\\Ri - 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 Ws(Rj,;a*,0) in Eq. (D.15) are discussed in Section 5.2 .
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Forces between like-charged walls immersed in electrolyte...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Forces between like-charged walls immersed in electrolyte solution Otto, Frank 2000
pdf
Page Metadata
Item Metadata
Title | Forces between like-charged walls immersed in electrolyte solution |
Creator |
Otto, Frank |
Date Issued | 2000 |
Description | 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 MM 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 MM level, the force between like-charged 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 Na+ and CI- 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. |
Extent | 8653171 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0061490 |
URI | http://hdl.handle.net/2429/11305 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2000-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
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
Cite
Citation Scheme:
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>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0061490/manifest