Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Structure and properties of fluids near inert and metallic surfaces Bérard, Daniel Robert 1994

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

Item Metadata

Download

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

Full Text

STRUCTURE AND PROPERTIES OF FLUIDS NEAR INERT AND METALLIC SURFACES By Daniel Robert Bérard B. Sc., Hon. (Chemistry) McGill University, 1989  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  December 1994  ©  Daniel Robert Bérard, 1994  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. I 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  of  The University of British Columbia Vancouver, Canada Date  DE-6 (2/88)  Dec q  Abstract  General reductions of the Ornstein-Zernike equation are given for multicomponent molecular fluids near an isolated planar surface and between planar surfaces. This al lows integral equation approximations such as the hypernetted-chain (HNC) and refer ence hypernetted-chain (RHNC) theories to be solved numerically. Theoretical methods for treating multipolar particles near inert, dielectric and metallic surfaces are consid ered. The metal is represented using a jellium model together with density-functional (DF) theory and the two phases interact electrostatically. Mean-field theories which re duce the many-body electrostatic wall-particle interactions to effective pair potentials are described for dielectric and metallic surfaces. The two interfacial phases are solved self-consistently for the wall-particle distribution function and, in the case of metal sur faces, for the electron density distributions. Explicit results are given for dipolar hard sphere fluids and for electrolyte solutions. The wall-induced fluid “structure”, electro static potential drop across the interface and electron density distribution of the metal are discussed in detail. Close to metallic surfaces, it is found that a highly ordered region exists. The dipoles are strongly ordered normal to the surface with the positive ends out. This is because the solvent structure effectively dictates the ion distributions near the surface. Hard sphere, Lennard-Jones and dipolar hard sphere fluids are investigated between inert walls. Lennard-Jones fluids are also considered between attractive walls. RHNC and HNC results for the fluid structure and the force acting between surfaces are compared with computer simulations. With the exception of Lennard-Jones fluids confined between inert walls, it is found that the integral equation theories show good agreement with  11  simulation results. Integral equation and computer simulation results for Lennard-Jones fluids show better agreement as the distance of the state point from liquid-vapour bulk coexistence is increased. A detailed investigation of a Lennard-Jones fluid confined between inert walls using the grand canonical Monte Carlo method is performed. Capillary evaporation is found for liquid subcritical bulk states. General methods are given for simulating a metastable fluid. The force acting between the walls is found to be attractive and increases rapidly as the spinodal separation is approached. On the equilibrium liquid branch, the net pressure appears significantly more attractive than the van der Waals attraction at small separations. This might be the origin of the experimentally measured attraction between hydrophobic surfaces in water.  111  Table of Contents  Abstract  ii  Table of Contents  vi  List of Tables  vii  List of Figures  ix  Table of Symbols  xvii  Acknowledgement  xxii  1  Introduction  1  2  Statistical Mechanics  5  2.1  Two-body functions and their properties  5  2.2  Integral equation theories  8  3  2.2.1  The Ornstein-Zernike equation  2.2.2  Blum’s reduction of the bulk Ornstein-Zernike equations  12  2.2.3  The reference hypernetted-chain theory  14  8  2.3  Monte Carlo methods  17  2.4  Models and pair potentials  20  Fluids near planar surfaces  23  3.1  23  Introduction  iv  3.2  Reduction of the planar wall-fluid Ornstein-Zernike equation  25  3.3  Mean-field theories for polarizable surfaces  33  3.3.1  Dielectric surfaces  33  3.3.2  Jellium surfaces  41  3.4  4  Dipolar fluids and electrolyte solutions near inert, dielectric and metallic surfaces  47  3.4.1  Fluid model  47  3.4.2  Asymptotic behaviour of h(01) near an inert wall  52  3.4.3  Results for a pure dipolar fluid  56  3.4.4  Results for electrolyte solutions  74  Fluids between planar surfaces  94  4.1  Introduction  94  4.2  Theory  96  4.2.1  4.2.2  5  The interaction free energy between two walls separated by a molec ular fluid  96  The wall-solvent-wall correlation functions  99  4.3  Grand canonical Monte Carlo method for Lennard-Jones particles  4.4  Results for pure fluids between planar walls  106  4.5  Cavitation of a Lennard-Jones fluid between hard walls  116  Summary and Conclusion  .  .  102  126  Bibliography  133  Appendices  143  V  A Properties of the Wigner symbols, rotation matrices and rotational invariants.  143  B The electrostatic interaction potential between a molecule and a uni form surface charge distribution.  V  147  C Numerical considerations  149  vi  List of Tables  3.1  The first few coefficients of the Legendre polynomial expansion, Pi(x)  =  ax 3.2  31  Reduced ion diameters d±/d together with the examples and labels used in this work. The solvent diameter d was taken to be 2.8A  3.3  Systems considered in the present work. In all cases the surface is in con tact with a solution of total reduced density p dipole moment and ionic charge are  =  0 where a 2.65a 0  =  =  8, respectively. For  0.529177A is the Bohr radius  52  The self-consistent electrostatic potential drop across the jellium-vacuum =  consists of dipolar hard spheres with p T  =  32A. When present, the liquid =  0.7,  =  d  2.8A and  298K  74  The self-consistent electrostatic potential drop across the jellium-solvent interface for r 3  =  0 and z 2.65a  dipolar hard spheres with 3.6  0.7. The reduced solvent  298K and, unless otherwise stated,  and jellium-solvent interfaces for z  3.5  =  and q*  =  the metal slab, we assume that T  3.4  51  =  =  32A. In all cases the liquid consists of at T  =  298K  74  Contact values of the wall-solvent pair distribution function near an inert surface  3.7  78  The potential drop in volts across inert surface-electrolyte solution inter faces.  The dielectric constant is also included. In all cases p  4=/andq*=8  =  0.7, 80  vii  3.8  The potential drop in volts across jellium surface-electrolyte solution in terfaces. The dielectric constant is also included. In all cases T = 298K, 3 = 2.65a 0 d = 2.8 and r  4.1  92  Details of the grand canonical Monte Carlo simulations for a Lennard Jones fluid confined between hard walls. The last column labeled “opera tions”, denotes the total number of particle moves, insertions and deletions included in the simulation. (The grand canonical Monte Carlo simulations for the bulk were performed in a central cell 9u x 9u x 9g.)  4.2  116  Details of the grand canonical Monte Carlo simulations for a Lennard Jones fluid with T* = 1.0 and res/kBT = —3.29 confined between planar hard walls. (The grand canonical Monte Carlo simulation for the bulk, performed in a central cell 9u x 9a x 9o-, has (ex)/kBT =  3.00, Ps0 3 =  0.75 + 0.01 and 3 Pbu / kBT = 0.44 ± 0.05.) 4.3  119  Details of the grand canonical Monte Carlo simulations for the second state point (T* = 0.9 and res/kBT = —3.84). (The grand canonical Monte Carlo simulation for the bulk, performed in a central cell 7u x 7cr x 7cr, has  ex  3.60, Ps /kBT = 0.4 ± 0.1.) 3 3 = 0.79 ± 0.01 and Pbu  viii  120  List of Figures  3.1  The coordinate system used in the reduction of the wall-particle Ornstein Zernike equation  3.2  27  The coordinate system used in Sec. 3.3.1 for a dipolar fluid near a dielectric surface showing two particle-image pairs. The images are shown to the left of the interface with a prime included in the labels  3.3  35  The coordinate system used in Sec. 3.3.2 for a dipolar fluid near a jellium metal slab  3.4  42  RHNC results for the projections occuring in the wall-solvent pair distri bution function for a dipolar hard sphere fluid with u.  =  and p  0.7.  In (a)-(d), the solid curves represent the present planar wall calculations and the dotted curves are results obtained previously (Ref. 10) for a macrosphere of diameter 30d 8 3.5  57  The wall-solvent pair distribution function  (O, 5 g  z) for the pure dipolar  fluid near an inert wall. The solvent parameters are as in Fig. 3.4 3.6  .  .  .  59  The asymptotic behaviour of (a) h° and (b) h . The dipole moment 2 and density are as in Fig 3.4. The solid curves are the full RHNC re sults obtained numerically. The dotted and dashed curves are the RHNC [Eq. (3.4.16)] and classical continuum asymptotes. The dash-dot curve in cluded in (a) represents the exact asymptotic behavior given by Eq. (3.4.19) and evaluated as discussed in the text  ix  60  3.7  RHNC results for the projections occuring in the wall-solvent pair corre lation function for a dipolar hard sphere fluid with The solid curve represents results for an inert wall  =  (=  /‘  and p  =  0.7.  1). The long- and  short-dashed curves were obtained in the cavity approximation for  =  cc  and 5, respectively. The long-dash-dot and short-dash-dot curves are the results obtained with the superposition approximation for  =  cc and 5,  respectively 3.8  62  The wall-solvent pair distribution function near a dielectric surface with  ,  =  gws(&,  z) for a pure dipolar fluid  5. The results in (a) and (b) are for  the cavity and full superposition approximations, respectively. The model parameters for the solvent are given in Fig. 3.7 3.9  64  The wall-solvent pair distribution function g(O, z) for a pure dipolar fluid near a conducting surface  (  =  cc). The results in (a) and (b) are for  the cavity and full superposition approximations, respectively. The model parameters for the solvent are given in Fig. 3.7  65  3.10 Oscillations in the net potential drop across the jellium-vaciium interface for r  0 2.65a  66  3.11 The self-consistent electron density (a) and effective potential (b) of a jellium slab, 64A wide, immersed in the pure dipolar hard sphere fluid. The dotted, solid and dashed curves are for r 8  =  2.0, 2.65 and 3.0a , 0  respectively. In (b), the dash-dot curve shows the self-consistent effective potential when no fluid is present and r 3  =  . The distribution of 0 2.65a  electrons in (a) is very similar to the distribution when the fluid is absent. The inset in (a) shows the change in the distribution due to the presence of the fluid (r  =  ) 0 2.65a  68  x  3.12 RHNC results for the projections of the wall-solvent pair correlation func tion for a dipolar hard sphere fluid in contact with a jellium metal. The dotted, solid and dashed curves are as in Fig. 3.11. The dash-dot curve is for an inert surface (r 3  oo)  69  3.13 The probability density p(8, z) for solvent particles at distances (a) 0.5db, (b) 0.8db, (c) 1.2d and (d) 1.5db from the jellium edge. The curves are as in Fig. 3.11  70  3.14 The wall-solvent pair distribution function wall with r  =  g(O,  z) near a jellium metal  . The solvent parameters are as in Fig. 3.7 0 2.65a  71  3.15 The reduced electrostatic potential across the interface when r 3  . 0 2.65a  (a) The net potential for the jellium slab in a dipolar hard sphere fluid (pa rameters as in Fig. 3.11) is represented by the solid curve and is the sum of contributions from the jellium (dotted curve) and the orientationally po larized fluid (dashed curve). The net potential across the jellium-vacuum interface is represented by the dash-dot curve. The scale on the right shows the electrostatic potential çb(z) in volts at 298K. (b) The electro static potential for z” > 0 on an expanded scale; curves are the same as in (a)  73  3.16 Comparison of the projections occuring in the wall-solvent pair correlation function for a dipolar hard sphere fluid near a jellium metal (r 3  =  ) 0 2.65a  surface and a dielectric surface. The short-dash curve is the jellium metal result. The long-dash and dash-dot curves are the results near a dielectric surface using the cavity and superposition approximations, respectively. The solid curve is result near an inert surface  xi  75  3.17 Projections of the wall-solvent pair distribution function for inert walls. The parameters are as in Tables 3.2 and 3.3. The solid curve is the pure solvent result and the dashed curve is for 0.945M CsF  77  3.18 Wall-ion pair distribution functions for 0.945M electrolyte solutions in con tact with an inert surface. The parameters are as in Tables 3.2 and 3.3. In  (----) Li in LiT; (---) Eq in EqEq; (_ —) Cs in CsT. In (b) the curves are: (—) Cl— in NaC1; (———) Eq in EqEq; (.-•—.) F— in CsF; (— —) 1 in CsI (a) the curves are:  (—) Na  in NaC1;  79  3.19 Projections of the wall-solvent pair distribution function for NaC1 solu tions in contact with a jellium slab. The solid curve is the pure solvent result. The dotted and dashed curves are for 0.945M and 7.56M solutions, 81  respectively. The parameters are as in Tables 3.2 and 3.3 3.20 Projections of the wall-solvent pair distribution function for RbF solu tions in contact with a jellium slab. The solid curve is the pure solvent result. The dotted and dashed curves are for 0.945M and 7.56M solutions,  82  respectively. The parameters are as in Tables 3.2 and 3.3 3.21 The wall-solvent pair distribution function  g(O,  z) for the 0.945M NaC1  solution in a dipolar solvent and near a jellium metal wall  83  3.22 The wall-solvent pair distribution function gws(O, z) for the 7.56M NaC1 84  solution in a dipolar solvent and near a jellium metal wall 3.23 The wall-cation  (—)  and wall-anion  (----)  pair distribution functions for  several 0.945M electrolyte solutions in contact with a jellium slab. The parameters are as in Tables 3.2 and 3.3. In (h), the wall-cation distribution functions for RbF, CsI, CsF and ILi are shown on an expanded scale.  .  87  3.24 An illustration of the distortions of the solvation shell due to the jellium induced interfacial polarization  88 xii  3.25 Wall-ion pair distribution functions for NaC1 solutions at 0.945M. The dotted, solid and dashed curves are for r 8  =  , 2.65a 0 4.0a 0 and 2.0a , 0  respectively  89  3.26 The wall-ion pair distribution function for NaC1 solutions near a jellium surface. The curves are labeled as in Fig. 3.23  90  3.27 The wall-ion pair distribution function for RbF solutions near a jellium surface. The curves are labeled as in Fig. 3.23  91  3.28 The self-consistent mean-field electrostatic potential for 0.945M NaC1 in contact with a metal slab (z  =  32A). The long-dash-short-dash, dash-  dot and long-dash curves are the jellium, dipole and ion contributions, respectively. The short-dash curve is the sum of the dipole and ion con tributions and the dotted curve is the dipole contribution obtained for the pure solvent. The solid line in (a) is the total electrostatic potential. The parameters are as in Tables 3.2 and 3.3 4.1  93  Net pressure between hard walls in a hard sphere fluid at a density of p  =  0.681. (a) The solid curve is the HNC result, the almost coincident dotted curve is the HNCP result and the dashed curve is the result obtained using the Verlet-Weis [27] and Henderson-Plischke [31] bridge functions.  (b)  Comparison on an expanded scale with the simulation data of Karlström [123] (squares); the curves are the same as in (a) 4.2  The wall-hard sphere distribution function  gws(z;  separations 1. The parameters are as in Fig. 4.1  xlii  107  1) for several wall-wall 109  4.3  The (a) HNC and (b) RHNC contact density minus the bulk density as a function of wall-wall separation 1. pg(d/2; 1)  —  The crosses denotes values of  Ps• The curves are the net HNC and RHNC pressures  calculated from Eq. (4.2.9) 4.4  110  HNCP pressure between hard walls in dipolar hard sphere fluids at a den sity of p’ = 0.7. In (a), the solid, dashed and dotted curves are for reduced dipoles moments of  j  = 0 (hard spheres),  and  =  respec  tively. In (b), the HNC asymptotes calculated from Eq. (4.4.5) (smooth curves) are compared with numerical HNCP results on an expanded scale 111 4.5  RHNC results for the angle-averaged wall-particle pair distribution func tion g°(z; 1) for a dipolar hard sphere fluid confined between hard walls. Results are shown for several wall-wall separations 1. The parameters are as in Fig. 4.4  4.6  112  RHNC results for the projection g 2 of the wall-particle pair distribution function for a dipolar hard sphere fluid confined between hard walls. Re sults are shown for several wall-wall separations 1. The parameters are as in Fig. 4.4  4.7  113  Wall-solvent distribution functions for a Lennard-Jones fluid confined be tween Lennard-Jones 10—4 walls, Eq. (4.3.3) (T* = 1.21, p = 0.5925 and 2 = 0.8). The wall-wall separations are (a) l = 2cr and (b) 1* = 3cr. Psurf The solid and dotted curves are the HNC results using the Lennard-Jones, Eq. (2.4.10), and truncated Lennard-Jones, Eq. (4.4.6), potentials, respec tively. The HNC results were obtained using the bulk direct correlation function of perturbation theory [124,125]. The circles are results from molecular dynamics simulations [108]  xiv  114  4.8  Wall-solvent distribution functions for a Lennard-Jones fluid confined be tween hard walls for three state points and two wall-wall separations. The solid curves are results from grand canonical Monte Carlo simulations. The dashed and dotted curves are the HNC and HNCP results. In (a) T*  2.0 and p”  =  in (c) and (d) T* 4.9  0.710 ± 0.003; (b) T*  =  1.0 and p  =  =  1.35 and p  =  0.70 ± 0.01; and  0.75 ± 0.01  115  The grand canonical Monte Carlo average density  , between 2 (N)/1L  hard walls for a Lennard-Jones fluid in equilibrium with the bulk. Results for the first state point (res/kBT = —3.29 and T* (a) and for the second state point (res/kBT  =  =  1.0) are shown in  —3.84 and T* = 0.9)  in (b). Crosses represent points that are thermodynamically stable and open squares represent those that are metastable. Open circles are used when the equilibrium and metastable fluids could not be distinguished. Except where illustrated, the standard deviations in 5 are smaller than the symbols. The inset shows an expanded view of the gas branch  118  4.10 Grand canonical Monte Carlo density distributions for the equilibrium and metastable liquid branches for wall-wall separations l = 5, 6, 7, 8 and 10. The first (a) and second (b) state points are defined as in Fig. 4.9. Here, the midplane between the walls is at z  0  121  4.11 The net pressure between hard walls. The symbols are the same as in Fig. 4.9. The solid curve is the conventional van der Waals attraction, Eq. (4.5.3), and the dashed curve is the fitted mean-field result, Eq. (4.5.2). The fitting parameters for the first state point (a) are A0/kBT = 0.15, =  0.02, and 10 = 3.Ou. For the second state point (b), Au/kBT = 0.13,  = 0.02, and 10 = 3.0cr. For the gas branch, the standard deviations in the net pressure are smaller than the symbols xv  122  4.12 The interaction free energy per unit area measured between hydrophobic surfaces in water [128]. The solid curve is the conventional van der Waals attraction, with a Hamaker constant of 2 x 10 J, and the dashed curve is 20 the fitted mean-field result, Eq. (4.5.1) (A and a  =  4 x 2 nm 4 10 )  =  1.2 x 1 Jm lo 2 10’ ,  =  l2nm, 125  xvi  Table of Symbols  Symbol  Description  0 a  Bohr radius  b(12)  Bridge function between particles of a and /3.  b(12)  Reference bridge function.  =  O.529177A.  In this thesis, either the  Henderson-Plischke or HNCP bridge functions are used.  B(l)  Wall-wall bridge function per unit surface area.  c(l2)  Direct correlation function between particles of a and  c(r)  Projection of the direct correlation function.  C(r)  x  e  Elementary charge.  dc,  Diameter of species a.  ffll  Nonzero constant. Here,  ) 2 g(l  Rotationally and translationally invariant pair distribution  transform of c. (r). Also C,(r) when 1  fmnl  =  =  0.  (2m + 1)(2n + 1).  function between particles of a and g(1, 2)  x  /3.  /3.  General pair distribution function between particles of a and  ‘3. g(l; 1)  Wall-particle pair distribution function for a particle of species  /3 confined between hard parallel walls separated by a distance g(1; 2)  Three-body distribution function defined after Eq. (3.3.5).  g(l2)  Projection of the pair distribution function.  xvii  h  Planck constant divided by 2ir.  ha(12)  Pair correlation function between particles of a and /3.  h(12)  Projection of the pair correlation function.  mnx Hpj,;(r)  x  jn(r)  Spherical Bessel function of order ri.  kB  Boltzmann constant.  1  Separation between walls and length of simulation cell along  mnl mn transform of h,j,;(r). Also H,;a(r) when  x = 0.  z axis. L  Length of grand canonical Monte Carlo simulation box along the x and y axes.  m  Mass. Average number density of elementary charges in metal from atomic cores.  N  Number of particles.  N(l)  iiww(l) per unit area.  , z) 0 pw(  Probability density that a particle of species /3 is at a distance z from a wall and at an angle 0 with respect to the surface  P  normal vector. Pressure. Specifically the net pressure P(l) = P(l)  Pt  Internal pressure. Pressure at a surface due to fluid confined between two surfaces.  Pb  Bulk fluid pressure.  Pj(cos 0)  Legendre polynomial of order 1.  qa  Charge on a particle of species a.  xviii  —  Pb.  Electric multipole moments of order  mq.I  of species o in spher  ical tensor notation and measured in the reference frame of the molecule. Electric multipole moments of order mqi of species o in spher ical tensor notation and measured in the interparticle refer r  ence frame. General vector. r 1 is the vector denoting the position of par ticle labeled 1; r 12 Unit vector i  =  =  2 r  r/r.  —  . 1 r  When used as an argument for the  Wigner spherical harmonics, i  =  (q, 0, 0) denotes the Euler  angles. 3 r  Wigner-Seitz radius. r 3  =  (3/47rn)’/.  Wigner generalized spherical harmonic. S(x)  Polynomial. See page 30  T  Temperature.  ut(z)  Mean-field potential due to fluid.  usI(z)  Self-image potential. See Eq. (3.3.13).  u°’(z)  “Other-image” potential. See Eqs (3.3.22) and (3.3.23).  iel()  Potential due to jellium background. Exchange and correlation potential. Pair potential between species o and  .  u(12)  Projection of the a(12).  V  Volume.  wa(12)  Potential of mean force potential between species a and 8.  w(12)  Projection of the w(12).  xix  y  y  Y(O, q)  Spherical harmonic.  z  With reference to the partition functions and Monte Carlo  =  47rps/9kBT  simulations, z is the activity. Otherwise, z denotes a position along the z-axis, normal to the surface. z  Vector directed normal to the surface. Unit vector in the direction of z. Fermi energy. Eigenvalue of the vth eigenfunction &(z). Eigenvalue of the kth eigenfunction  lk(r).  Dielectric constant of species a. Series function.  =  h  —  c.  Isothermal compressibility. Chemical potential. ex  Excess chemical potential.  id  Ideal chemical potential.  res  Called the residual chemical potential. Defined here as 3 where z is the activity. kBT in zu  Dipole moment of species a. Grand canonical partition function.  E(, V, T)  Number of species composing system. Average bulk number density of species a.  p p(’)(1, 2,. Pth  ..  ,  n)  n-particle density. Threshold density defined on page 104. Lennard-Jones collision parameter.  xx  =  , 112) 2 Z  Translational invariant basis function. One of the azimuthal Euler angles.  cb(z)  Scalar potential.  x  One of the azimuthal Euler angles. 2,  12 r )  Rotational invariant basis function. One dimensional one-electron wavefunction for state v.  I1k(r)  Three dimensional one-electron wavefunction for state k. Euler angles  (çf,O,x).  Accents and other symbols Used to denote —m. h  Fourier or Hankel transform of h.  ()  Grand canonical ensemble average.  xxi  Acknowledgement  It is my pleasure to thank my academic supervisor, Prof. G. N. Patey, whose has given me guidance, encouragement and independence in my studies over the last five years. I would also like to thank the Natural Sciences and Engineering Research Council (NSERC) of Canada and the University of British Columbia Chemistry Department for financial support and Dr. Phil Attard, Dr. Masahiro Kinoshita, Dr. Glen Torrie, Dr. Dongqing Wei and Dr. Amalendu Chandra for numerous interesting discussion pertaining to my work.  xxii  Chapter 1  Introduction  The solid-liquid interface has received a great deal of interest in the past few decades from both theoreticians and experimentalists. The principle reason for this interest is that the properties of the individual phases at the interface can be significantly different from those of the bulk. Furthermore, the interfacial properties depend on the nature of both phases and the properties of a liquid, for instance, can vary significantly at different surfaces. A familiar example at the macroscopic level are the wetting properties of water at hydrophobic and hydrophilic surfaces. However, fundamental to an understanding of any physical system is knowledge of the underlying microscopic “structure”. Statistical mechanics provides a means of calculating the structure and other ob servable properties of a given system based only upon a description of the microscopic interactions and the state parameters. One of the principle tasks in any statistical me chanical theory is the description and development of microscopic models and two general approaches are often followed. One of these employs potentials that have been fitted to describe the particle interactions in specific physical systems. Such potentials are often restricted to a particular state of the system and may only model a small set of properties accurately. Alternatively, simple model potentials which represent certain intrinsic prop erties of molecular interactions can be employed. Such an approach is directed towards understanding the relationship between general features of intermolecular interactions and observable properties. The simplest model potential which can reproduce an observ able property often yields the most insight into the nature of a system.  1  Chapter 1. Introduction  2  Particle distribution functions, in particular the singlet and pair distributions, are of central importance to an understanding of the microscopic structure of fluids and solids. The singlet distribution function describes the probability density of a particle being at a given position and orientation. The pair distribution function describes the probability density of two particles simultaneously having given sets of positional and orientational coordinates. If we assume that the potential energy can be written as the sum of pairwise potentials, all thermodynamical properties of a system can be calculated from the singlet and pair distribution functions. Integral equation methods provide one approach for calculating pair distribution func tions. Methods based on the Ornstein-Zernike equation are commonly used in liquid state studies and will be employed here. The Ornstein-Zernike equation relates the pair distri bution g(12) to the direct correlation function c(12) and may be considered as a definition of c(12). In order to employ the Ornstein-Zernike equation, a second closure relation be tween g(12) and c(12) is necessary. An exact closure can be formally written but this introduces a third unknown function b(12). The Ornstein-Zernike equation together with an approximate closure form an integral equation theory. Approximate closures can be viewed as different approximations for the unknown function b(12). Computer simulations provide an alternative approach to the liquid state. In order to understand the influence of a surface on a fluid, it is particularly useful to compare the fluid properties at the interface with those of the bulk. Grand canonical Monte Carlo calculations are especially useful in this regard. Monte Carlo methods are based on sampling system configurations. Starting with an initial configuration of particles, additional configurations are generated by translating, rotating, inserting and removing particles according to their relative Boltzmann probabilities. The equilibrium distribu tion functions and other properties of the system are calculated as ensemble averages. The only inputs required by Monte Carlo and other computer simulation methods are  Chapter 1. Introduction  3  the interaction potentials and state parameters. Hence, for a given model, computer simulation methods yield essentially exact results. In this thesis, the liquid side of the interface is of primary concern and we will study the effects of surfaces and surface properties on the fluid structure. Experimentally, there is almost no structural data on which to compare theoretical results at interfaces. In fact, the liquid structure at the interface has only recently been probed using tunneljunction devices [1] and x-ray diffraction measurements [2]. A few limited results have been obtained for water at mercury and silver surfaces. There has been a considerable amount of theoretical work focussed on inhomogeneous systems. These include dense fluids near simple surfaces [3] and electron density distri butions at metallic surfaces [4]. However, these methods emphasize a detailed model and theory of one phase and simplify the other phases enormously. For example, metal surfaces have been studied primarily in vacuum or in dilute gases (emphasizing adsorp tion phenomena) [4]. When the liquid phase is of interest, detailed classical liquid state models and theories have been employed at inert and simple attractive walls. This thesis is divided into five chapters. Chapter 2 is devoted to the fundamentals of liquid state theory and begins with a discussion of the properties of pair functions in homogeneous and isotropic systems. The integral equation theories and grand canonical Monte Carlo method are presented and the chapter concludes with a discussion of pair potentials. Chapter 3 is devoted to molecular fluids near isolated planar surfaces and contains an introduction to the relevant literature. A general reduction of the Ornstein Zernike equation for fluids near planar surfaces is presented and theories for molecular fluids near inert, dielectric and metallic surfaces are derived. Simple dipolar hard sphere fluids and electrolyte solutions in contact with these surfaces are considered. The asymp totic behaviour of the pair distribution function for the pure dipolar fluid at long range is derived and compared with analytical continuum and exact asymptotic profiles. Chapter  Chapter 1. Introduction  4  4 is concerned with the structure of fluids confined between surfaces. General reductions of the Ornstein-Zernike equation describing the pair correlation functions in these sys tems are presented. The grand canonical Monte Carlo method used is described in detail and a general method of simulating metastable states with Monte Carlo is included. A variety of model fluids are studied between inert and attractive surfaces. The integral equation results are compared with simulations and a discussion concerning the appli cability of integral equation methods to different interfacial models is given. A detailed Monte Carlo study of a Lennard-Jones fluid confined between inert walls is also described. Equilibrium and metastable simulation results are presented including density profiles, average densities, and the force between walls as a function of surface separation. The possible relevance of the Monte Carlo results to the experimentally measured attraction between hydrophobic surfaces is discussed. Finally, a summary of conclusions is given in Chapter 5.  Chapter 2  Statistical Mechanics  In this chapter, the basic equilibrium statistical mechanical tools and theories used in this thesis will be outlined.  Two-body functions and their properties  2.1  In classical statistical mechanics, the grand canonical n-particle density is defined as N  p(1, 2,.  .  .  ,  n) =  (N  where the arguments 1, 2, N,  z  n)!  —  ...,  is the activity, U(1, 2,.  .  .  f  NTd(n 2 e_U(1  + 1) d(n +2).  d(N)  ,  (2.1.1)  N denote the coordinates of particles labeled 1, 2, ,  N) is the N-particle potential, kBT is the Boltzmann  constant times the temperature and  is the grand canonical partition function. The  pair distribution function will be of fundamental importance and is defined in terms of )(1, 2) as 2 the pair density p( (2)(1  g(1,2)  2) p(’)(1)p(’)(2)  (2.1.2)  For a homogeneous and isotropic system of spherically symmetric particles, the pair dis 12 tribution is simply a function of the interparticle separation r  =  2 r  —  1 and describes r  12 from particle 1. Conversely, the probability density of finding a particle 2 at distance r for a molecular system, g(1, 2) can be a function of as many as nine independent co ordinates in a translationally invariant environment and twelve in an inhomogeneous environment. In order to work with these functions, it is necessary to first introduce the mathematical tools needed to reduce them to a tractable form. 5  Chapter 2. Statistical Mechanics  6  The translationally invariant pair function for a molecular system can be written as g(1,2)  =  (qj,Oj,xj)  1 r where 1 12 , 1 g(1Z , 2 ),  ing the orientation of molecule  j  and r 12  =  2 r  represents the set of Euler angles denot  —  . In this form, the two-body function 1 r  can be expanded in an orthogonal basis set which spans the anisotropic orientational space of the two particles. Specifically, g(1,2)  =  l r 12 j 1 g(Z , 2 )  = mnl  12 where r  =  12 r 12 ,i 2 Iri 12 = /r  >  g,)1(ri2) T  ,I(11,12,i12)  ,  (2.1.3a)  and the basis functions are  , cz 1 T),(1Z , 112) 2  = ffl 1  (i 0 R, ) ) R,(fZ 1 R( ) 12 2  Here, R(fZ) is a Wigner generalized spherical harmonic [5] and  .  fmnl  (2.1.3b) is an arbitrary  nonzero constant chosen for mathematical convenience. When used as an argument of the Wigner spherical harmonics, i denotes the Euler angles of r, (q, Or, 0). The eight indices in Eq. (2.1.3) take on all integer values for which the Wigner generalized harmonics are defined. That is,  m,n,l  =  (2.l.4a)  0,1,2,...  —m<<m  and  —m’<m  n  and  —n  —n  v  —l  )‘  <1  i/  .  n  (2.1.4b) (2.1.4c) (2.1.4d)  Provided the expansion is rapidly convergent and the sum can be truncated for relatively small values of m, n and 1, the one-dimensional coefficients  are relatively easy  to calculate. In this work, only translationally and rotationally invariant pair functions will be con sidered and, in general, these can be described by six (or fewer) independent coordinates.  Chapter 2. Statistical Mechanics  7  In this case, an appropriate basis set expansion can be constructed from Eq. (2.1.3) when expressed as a contracted set of the coefficients,  gIl(r12),  and basis functions,  2). The rotational invariance of g(1, 2) implies that an arbitrary rotation of the molecular pair (while keep the relative orientation and interparticle spacing fixed) leaves g(1, 2) unchanged. Mathematically, the concerted rotation is achieved by apply ing a rotational operator R(f) to each particle and to the interparticle vector 112 in Eq. (2.1.3) [6]. This gives  R,,(1Z) R ) ( 1  ) 2 fmnlgnn,,(ri  g(12) mnl  z)’  x () 11 R  (i, 0 R ) RIA,l(f) 2  where Eq. (A.13) has been used’. Integrating over all orientations  (2.1.5) and using Eq. (A.1O)  yields  ) 2 g’(r,  g(12) =  ,f, mnl(f,, ) 2  ,  (2.1.6a)  mnl  where the so-called rotational invariants are defined as g(r,2)  =  gsl(r,2),  (2.1.6b)  /I,V,)  mnl( 2,  ) 2 r,  = ?I,v,),  )  (  R(,) R( )2 2 (, 0 R, )  I  (2.1.6c)  is the usual 3-j symbol [5] and, in order to emphasize that the pair functions  depend only on the interparticle coordinates, the commas have been removed from the arguments of g(12). The contracted basis functions have five indices (one for each degree of rotational freedom). These indices can take on the values given in Eq. (2.1.4); however, the 3-j symbols are nonzero only if  Im iI —  1  m  + n.  properties of the Wigner functions are summarized in Appendix A.  (2.1.7)  Chapter 2. Statistical Mechanics  8  The angle-averaged distribution function is simply g°(r). The expansion (2.1.6) provides a general mechanism for describing two-body func tions which only imposes the constraints that the functions are invariant to rotations and translations. The pair functions of interest for molecular fluids must satisfy a few additional constraints [6]. Specifically, the two-body functions of interest here are also real, invariant to the interchange of particles and invariant under the symmetry oper ations of molecules 1 and 2. As a result of these properties, many of the coefficients become degenerate or vanish, and the number of unique coefficients is reduced. The first conditions that the functions are real implies rnnl*( 1 r g  where i  = —ii,  —  —  (  ‘m+n+1+#+v mnl(r’ \  2 1  ,  and is obtained by requiring that the product of two real rotationally and  translationally invariant functions is real. The second conditions can be expressed as g(12)  =  1) 2 g(  (2.1.9)  .  (i) 0 Using this, the parity of the generalized spherical harmonics R  = (_)l  (—i) and 0 R  the symmetry of the 3-j symbols (see Appendix A) gives the identity g(r)  = (_)m+ng1l(r)  .  (2.1.10)  The molecular symmetry properties further reduce the number of unique coefficients and identify those which must be zero (see Ref. [6] for a discussion). This reduction will be discussed in Chapter 3 for dipolar particles.  2.2  2.2.1  Integral equation theories The Ornstein-Zernike equation  The Ornstein-Zernike (OZ) equation is the basis for many of the integral equation theories of fluids and was first introduced in 1914 [7] to study density fluctuations in fluids near  Chapter 2. Statistical Mechanics  9  the critical point. The OZ equation relates the pair and direct correlation functions, h and c, respectively, through a convolution and is given in its most general form for a one-component molecular fluid as h(1,2)  c(1,2)+--Jp(’)(3)h(1,3)c(3,2)d(3),  where p(’)(3) is the singlet density, h(1, 2) = g(1, 2)  —  (2.2.1)  1 and fd(3) = fdr 3d 3 denotes  the integral over all positions and orientations of particle labeled 3. The OZ equation is often viewed as the defining relation for the direct correlation function; however, c can also be defined as the functional derivative [3] —  i 2  6ln[p()(1)/z*(1)] 222  6p(’)(2)  —  where z*(1) = z e_()1IT and q(1) is an external potential specified by the problem at hand. Using this definition, and starting from the grand canonical ensemble of classical statistical mechanics, the OZ equation can be derived if we assume that all particles interact through a pairwise additive potential such that [3] N  U(1,2,. ..,N)  = u(i,j)  .  (2.2.3)  i<j  The assumption that the fluid behaves according to classical mechanics implies that A =  2rh  N rnkBT  <<a,  (2.2.4a)  where A is the thermal de Broglie wavelength, h is Planck’s constant divided by 27r, m is the particle mass and a represents the average nearest-neighbour separation between particles. For molecular fluids, the additional constraint that  Orot = is required, where  °rot  h 2 <<T 21kB  (2.2.4b)  is the characteristic rotation temperature and I is the molecular  moment of inertia. These conditions are satisfied for typical bulk fluids and, in particular,  Chapter 2. Statistical Mechanics  10  are satisfied for liquid water. In the process of deriving the OZ equation from Eq. (2.2.2) [3], the specification of the chemical potential is replaced by the average singlet density p(1).  Eq. (2.2.1) is usually referred to as the inhomogeneous OZ equation since the  correlation functions are solved in an external field. When the two-body correlation functions are invariant to translations and rotations and the functional derivative in Eq. (2.2.2) is evaluated at q(l)  =  0, the OZ equation  has the simpler form 2 h(12)  =  c(12) +  J  h(13) c(32) d(3)  ,  (2.2.5)  where the singlet density reduces to the average grand canonical density p and the com mas have again been removed from the function variables to emphasize that they depend only on the interparticle coordinates. Equation (2.2.5) is referred to as the homogeneous OZ equation and can not be solved for both h and c until a second closure relation is introduced. The inhomogeneous OZ equation requires two additional relationships. To date, many approximate closures have been proposed and some of these will be considered in Sec. 2.2.3. The designation “homogeneous” and “inhomogeneous” are somewhat inappropriate and have apparently led to some confusion in the literature concerning the applicability of Eq. (2.2.5) [8,9]. Hence, it is useful to define these terms in the present context. Since the OZ equation relates two-body functions, the distinction between Eqs (2.2.1) and (2.2.5) must be made at the pair level. Hence we define a system as homogeneous (or pairwise homogeneous) if the interactions and correlations between a pair of particles depends only on the interparticle coordinates (i.e., translationally and rotationally invariant). In this way, homogeneous also implies that the system is pairwise isotropic and is valid whenever q(l)  =  0; and  The homogeneous OZ Eq. (2.2.5) is also be obtained if particle 1 is assumed to be the source of q(l). 2  Chapter 2. Statistical Mechanics  11  inhomogeneous if the two-body correlations are translational and/or rotational vari ants. The external field in the inhomogeneous equation can often be defined in terms of the field of one or more particles. In this case, the distinction between a pairwise homogeneous and a pairwise inhomogeneous system is academic. In practice, the inhomogeneous equation not only provides more information about the microscopic system structure (effectively as the three-body correlation functions), but is typically more accurate for a given closure approximation. However, it is worthwhile to note that the inhomogeneous theories based on Eq. (2.2.1) are significantly more complicated and time-consuming to solve. For the systems that will be considered in this thesis, all species are modeled explicitly and the simpler homogeneous Eq. (2.2.5) is used. Before proceeding, we generalize the homogeneous OZ equation to a multicomponent system of  + 1 species to obtain  h(12)  =  ca(12) +  f  ha(13) c(32) d(3),  (a,  =  0,1,.  .  .  , ),  (2.2.6)  Chapter 2. Statistical Mechanics  12  where p is the average bulk number density of species 7. In this work, the species labeled w, for wall, is considered to be infinitely dilute (Pw  =  0) and will be used to represent a  surface species. As a result, the coupled set of Eqs (2.2.6) is reduced to four smaller sets  h(12) h(12) haw(12) h(12) where the labels (cr, /3  —  —  —  —  =  c(12) c(12)  =  c(12)  =  c(12)  =  1, 2,..  .  ,  J J f 2f f 2  (13) c(32) d(3), 7 h  (2.2.7a)  h(13) c(32) d(3),  (2.2.7b)  h(13) c(32) d(3),  (2.2.7c)  h(13) c(32) d(3)  (2.2.7d)  ,  now exclude the dilute species w. The first set of  Eqs (2.2.7a) describe the correlation functions for the bulk fluid in the absence of a wall and have been the topic of much earlier work [10—16]. The solutions to Eqs (2.2.7b) and (2.2.7c) are equivalent—the two-body correlations between w and the remaining species 1, 2,.  .  .  ,  c—and only one is required. The last equation describes the correlations  between two dilute w particles. The solution of Eq. (2.2.7a) is required as input to Eqs (2.2.7b) and (2.2.7c) which in turn are required as input to Eq. (2.2.7d). 2.2.2  Blum’s reduction of the bulk OZ equations  Although we are primarily interested in the solution to Eqs (2.2.7b)—(2.2.7d), we require the bulk pair correlation functions before any of these equations can be solved. Hence, we briefly discuss the solution of Eq. (2.2.7a) for nonspherical particles. As we have seen, the rotationally and translationally invariant two-body functions vary with six independent variables and it is obviously not realistic to solve Eq. (2.2.7a) in its present form. The reduction of the OZ equation for pairwise homogeneous and isotropic fluids to a set of  Chapter 2. Statistical Mechanics  13  coupled algebraic equations in one variable is due to the work of Blum and Torruella [6] and of Blum [17,18]. The reduction can also be found in Ref. [10] and only a brief review is given here. Since the OZ equation is in the form of a convolution, the spatial integral can be eliminated by taking the three-dimensional Fourier transform. This gives 7a(12)  —  J  cyp(12) =  , 3 cv() y(32) d  (2.2.8)  where the Fourier transforms are defined as (Zi, 2, 3 iai  J  k) =  e’’  ha(i, 12, r) dr  ,  r)dr feirc&(i,1z , 2  =  (2.2.9a) (2.2.9b)  ,  and k is a reciprocal space vector. The k-space correlation functions can be expanded in a rotational invariant basis set yielding Jia(12)  =  E  k)  , 2, 1(12) 1 ‘(fZ  k)  mnl()  ,  (2 2 lOa)  mnl  a(12)  =  (22 lOb)  mnl j.LV  The coefficients  are related to the real space coefficients by the Hankel  and  transforms defined as k)  =  1 47ri  j  jj(kr)h(r)dr 2 r  ,  (2.2.11)  where ji(kr) is the spherical Bessel function of order 1. The Hankel transforms can be evaluated using Fourier transform methods [10,6,17,19]. Substituting these expansions into Eq. (2.2.8) and using properties of the rotational invariants (see Appendix A), the bulk OZ equation reduces to relations of the form (2.2.12)  =  -y  12fl 1  Chapter 2. Statistical Mechanics  where  h  =  —  c,  f  ‘712111  Li mnn  —  —  —x and  1  2 fmnilifninl 12 11 l /1211 l’\ (i)m+n+ni 21 + 1 fmnl 1 +1 2n mn nif 0 o oJ fmnhllfnlnl2 2l + 1 1 f +1 1 2n rjl /mnl” 1 /n ’ ’\ 2 l 1 mn  (  x  =  14  (2.2.13a)  (2.2.13b)  { } is a 6-j symbol [5]. If we use the definition fl  n  fmnl  =  /(2m + 1)(2n + 1),  (2.2.14)  referred to henceforth as Blum’s notation, and define the so-called ymn;x (k’ iw;a/3  /  k)  ,nn 1  -mnl c,(k)  =  x transforms  as  ,  (2.2.15a)  ,  (2.2.15b)  l) =  o (x /mnl)  mnl  then Eq. (2.2.12) simplifies considerably to give  fmn;x(j ILv;a/3\  Blum  ()x+iii  /  -y  where  Blum =  Ffmn1;xfk+mn1;x(k] wi;ay\ j I. 1wl;a7 I  nin;x  1) k  i’;8 1 i  (2.2.16)  1’1  denotes the use of Eq. (2.2.14). These equations can also be conveniently  assembled into matrix form as shown in Refs [17], [18] and [10]. 2.2.3  The reference hypernetted-chain theory  The OZ equation (2.2.5) relates the unknown pair and direct correlation functions. In order to make use of the relationship, a second closure condition relating h and c is required. Such a relation was simultaneously developed by several authors [20—24] around 1960 and can be derived as a functional Taylor series expansion of Eq. (2.2.2) [3,25] to give h(12)  —  c(12)  =  in [h(12) +  +  u(12) kBT  —  b  ‘12),  (2.2.17)  15  Chapter 2. Statistical Mechanics  where u is the pair potential between particles of species c and 3. The function bafi represents the bridge (or elementary) diagrams and consist of an infinite series of multicentered integrals. With the exception of the first two- and three-centered integrals (to order p 2 and p 3 in density) [26], the bridge functions are too complicated to evaluate. . However, 2 They are usually considered to be short ranged; decaying at long range as h for fluids near macroscopic surfaces it is expected that b may become very important in certain cases (for example, as the equilibrium bulk phase approaches the coexistence curve). In order to obtain a useful integral equation theory, an approximation of the ex act relation (2.2.17) must be constructed. The approximation bai = 0 is called the hypernetted-chain (HNC) theory since the resulting expression for c consists of diagrams called simple chains, netted chains and bundles. Another approximation is to replace b.o with those of a simpler reference system b, such as a hard sphere fluid. The resulting closure ha(12)  —  c(12) = ln [ha(12) +  +  u(12) —  b(12),  (2.2.18)  is called the reference hypernetted-chain (RHNC) approximation. Since the bridge dia grams for hard spheres are also difficult to evaluate, several indirect methods of estimating the bulk hard sphere-hard sphere [26—30] and the hard wall-hard sphere [31,26,32] func tions have been developed. Most of these methods are based on analytically fitting the pair distribution function. The bridge function is then calculated from Eq. (2.2.18) using the OZ definition for c(12). In the method of Attard and Patey [26], the hard sphere bridge diagrams are calculated explicitly to second, b2(r), and third, b(r), order in density. The bridge function is then approximated by the Padé form b(r)  Ps  1  —  WS\  /  b (r)/bws (r) 5 p  (2.2.19)  16  Chapter 2. Statistical Mechanics  This equation is referred to as the HNCP approximation. Although the resulting bridge functions from the various methods can differ significantly, most of these methods give very good estimates of the hard sphere-hard sphere and hard wall-hard sphere pair distri bution functions. Unless otherwise stated, we use the Verlet-Weis [27] method generalized to multicomponent fluids by Lee and Levesque [28] for the bulk bridge functions  and  either the Henderson-Plischke fit [31] or the HNCP approximation for the wall-particle functions b. In order to make use of a closure for nonspherical particles, Eq. (2.2.17) must be expanded in a rotationally invariant basis set. One method of achieving this is to rear range the closure equation and take the partial derivative with respect to the interparticle separation r 12 [10]. After rearranging the result, this gives Oca(12) 12 Or  18a,(12) +3(12) 12 Or 12 kBT Or  hafl(12)Owa, ( 3 12) 12 Or kBT  =  ,  (2.2.20)  where wa(12)  =  (12) 3 —kBTlng,  =  —kBTi(12) + ua(12)  —  (2.2.21)  kBTba(12)  is the potential of mean force. Integrating from an arbitrary radius r to oc and assuming that ca(12)  *  ua(12)/kBT  —*  1  c(12) =  kBT  12 0 as r ha(12)  oo (as is usually the case), one finds that  —*  Ow  112  (12) dr ‘r  (12)  + ba(12)  .  (2.2.22)  —  Expanding the two-body functions in rotational invariants using Eq. (2.1.6) and reducing yields pmnl  jj  c(r) =  OWm2n2t2  oc  h’(r’)  l 22 f m . t 2 L  mnl  ( up,;af3r1 —  1  ‘  m  i  jmn1  Y v;a/3T)  ,  (r’)  dr’  17  Chapter 2. Statistical Mechanics  where fmn l fm2n212  (2m + 1)(2n + 1)(2l +  =  X  1 ni  1 n  ‘1  Tfl2  2  12  1)(_)m+n+1+1+1+2+M2  m m (ni n 2 n (l 12 1 (m 2 1 1 - 11 - 11 1 2 VI J 1 V 1 I2 /1) \i’ \IL 000 \  ,  (2.2.23b)  mnl and 131  J  33  34 35 36 18  39  is a 9-j symbol [5]. Since the reference fluid consists of spherical particles in all current applications of the RHNC theory, the only nonzero projection of b(12) is b(r). Similar methods in which partial derivatives are taken with respect to angular coordi nates have also been derived [33,34]. In this approach, the assumption that ca(12) as  12 —* i’  2.3  _+  0  oo is not necessary.  Monte Carlo methods  The grand canonical Monte Carlo method provides an alternative approach for liquids in bulk or confined systems. The method basically consists of sampling configurational phase space points occurring in the partition function Z  E(, V, T) =  J  -U(1,2  N)/kuTd(l)  d(2)... d(N).  (2.3.1)  d(1)d(2)... d(N),  (2.3.2)  The distribution functions and observables, co  (A) =  N 2  N!f  A(1 2 E(t,’çT)  together with the associated statistical error, are then calculated as ensemble averages. In principle, the Monte Carlo method is exact for a given model system and better estimates of the observables can be obtained by sampling more phase space points.  Chapter 2. Statistical Mechanics  18  For nearly all systems of interest (N)  c’-’  Avogadro’s number NA and a method of  reducing the number of independent particles N is necessary. The standard method of doing this while still considering a large system size is to introduce periodic boundary conditions. In this approach, a given configuration of N particles in a finite volume V is replicated periodically to form an infinite lattice. The N particles in the central cell are free to translate over all space provided their replicate “images” follow identical trajectories. Even when (N) is relatively small, direct sampling of phase space is inefficient since most phase points make little or no contribution to the integral. In an alternative tech nique called importance sampling, phase points are chosen from a weighted distribution. The Metropolis method implements this technique by starting with an initial configura tion and follows a Markov chain with a limiting distribution of states given by PiVT  =  z V N!  N e_U(1,2,N)PBT  E(t,V,T)  ‘  for spherical particles. The Metropolis method in the grand canonical ensemble consists of a series of particle translations, insertions and deletions. In a translation from coilfiguration m to n, a particle is chosen at random and a new set of test coordinates are selected from within a volume dxdydz centered at the old position. The move is accepted if the new configuration is lower in energy Un  —  Urn <0; otherwise, it is accepted with a  probability N  rm  =  z N e _Ufl/kBT/N!’ / .  zNe_Um/IcBT/N!’’ /  =  e_(_Jm)T  < 1  .  (2.3.4)  —  .  In either case, the resulting configuration contributes to the ensemble average.  In a  deletion step, a particle is again chosen at random. In an insertion step, a new set of coordinates within the central cell are selected at random for a test particle N + 1. In either case, if the system energy is lowered by a deletion or insertion, the new  19  Chapter 2. Statistical Mechanics  configuration of N  —  1 or N + 1 molecules is accepted; otherwise, it is accepted with a  probability N  Pn Pm  N C(UnUm)/kBT zV  =  (2.3.5)  ,  for a deletion and N  Z  —  V  (2 3 6)  e__Um)T  N+1  —  ‘  for an addition. The principle task in a Monte Carlo simulation is the continual calculation of the potential energy U in the Markov chain of configurations. When periodic boundary con ditions are used, the potential is usually partitioned into short and long ranged contri butions. The short-range contribution involves direct summation of the pair interactions using the minimum image convention where a particle i interacts with the nearest image of particle  j.  In addition, the short-range contributions to U is typically truncated using  a spherical cut-off. The longer ranged contributions from particles beyond this truncation sphere can be included using various mean-field techniques [35]. The grand canonical Monte Carlo calculation is usually carried out at constant re 3 duced activity z  =  /A 3 ehh/Tcr  (  has units of length and is an intrinsic property  of the molecular model) instead of constant chemical potential holding  it  .  This is equivalent to  fixed but has the advantage that the deBroglie wavelength A does not have to  be defined. Using z, the excess chemical potential is given by itex  =  =  where it’ itid =  =  ititid itres  —  3 kBTln  (2.3.7)  ,  3 is often called the residual or configurational chemical potential, kBT ln za  kBT ln+ 3 in A is the chemical potential for an ideal gas and  =  (N)/V. Another  Chapter 2. Statistical Mechanics  20  method, due to Adams [36,37,35], can also be used. Here the simulation is performed for constant B, V and T where B=lnzV ex/kT+l_V  =  (2.3.8)  This method also has the advantage that the deBroglie wavelength does not have to be defined. However, even when the chemical potential and activity are constant, B is a function of V, the dimensions of the simulation cell. The Monte Carlo method has the advantage that it can provide an exact solution for a given model provided that a sufficient number of statistical sample configurations have been averaged. However, Monte Carlo simulations can require several orders of magnitude more computer time than the integral equation theories.  2.4  Models and pair potentials  In order to use the integral equation theories, Monte Carlo method or any other statisti cal mechanical theory, it is necessary at some point to consider the interaction potential. Fundamentally, the potential is governed by the dynamic electronic and nuclear struc ture of the constituent molecules. At the pair level, the potential can be partitioned artificially into two basic components including a short-range repulsive term,  US,  arising  from electron and nuclear repulsion as the molecular charge distributions overlap, and a long-range electrostatic term,  el,  to get  up(l2)  =  u(12) + uj(l2)  .  (2.4.1)  The short-range potential is the least-well characterized contribution to the pair poten tial and either fitted potentials or simple generic models must be used to represent it.  21  Chapter 2. Statistical Mechanics  Here, we are not interested in modeling specific systems and will employ some simplified approximations for u including the hard-sphere potential  I  HS u(r)  r<(d+d)/2,  oo,  (0,  (2.4.2)  r>(da+d)/2,  and the algebraic inverse power-law potential zt(r)  cx  -—  (2.4.3)  ,  where n is typically chosen between 9 and 12. Although both forms are chosen in part for mathematical convenience, they provide a reasonable approximation of us’ for “classical fluids” and give good results for many properties [3,11—16,38]. The electrostatic pair interaction between particles of species a and 9 in electrostatic units (esu) is given by u(12) = e2JJ n(r1)n(r2)dd  (2.4.4)  where e is the elementary charge and na(r) is the molecular charge distribution of species a. In a point charge model, na(r) is simply n(r) where  j  (2.4.5)  =  is the number of point charges q at r of species a. Equation (2.4.4) can be  conveniently arranged in the form of a rotational invariant expansion [17], u(12)  =  ) 2 u(rl  i 12 , 2 M 1 ( )  ,  (2.4.6a)  mnl in,  mnl;el  (_)m6m÷ni  12 —  \j  (21 + 1)! (2m)!(2n)!  13 Q;c Qi;  appropriate for use in the integral equation theories. Here  1 r  2.4.  represents the electric  multipole moments of species a in spherical tensor notation and in the reference frame  22  Chapter 2. Statistical Mechanics  are constant and given  of the molecule. For a rigid molecular model, the moments by eJrma(r)R(I)dr  =  ,  (2.4.7)  where ia(r’) is the charge distribution measured in the molecular reference frame. The multipole moments in the space-fixed,  Q,,  and molecule-fixed,  reference frames  are related by a rotational transformation where =  ) 1 R(f  (2.4.8)  Q,,  and  Q(a) =  eJrmn (r)R m (I)dr  .  (2.4.9)  The average contribution to the pair potential from molecular polarizability can be in cluded in a mean-field dipole-dipole contribution to the electrostatic potential. When the molecules have no permanent multipole moments, the Lennard-Jones potential, [(a)12  u’(r)  =  4e  —  j 6 (cT) ,  (2.4.10)  where e is the depth of the potential well and a, the collision parameter, is the parti cle separation when interaction.  uLJ(r) =  0, provides a reasonably good approximation of the pair  Chapter 3  Fluids near planar surfaces  3.1  Introduction  As discussed in Chapter 1, the solid-liquid interface is of importance in many physi cal systems and has received a considerable amount of theoretical attention in recent years [39—79]. The objective of much of this work is to develop a completely molecular description of the interface and, towards this goal, two fundamental approaches arise naturally. One of these approaches uses quantum density-functional theory to study the metal-vacuum interface [80]. The second is based on classical statistical mechanical meth ods including integral equation theories [39—53] and computer simulations [54—61]. The statistical mechanical methods have been used to study liquids and electrolyte solutions near simple non-metallic surfaces. During the last decade, additional advances have been made towards combining the density-functional and statistical mechanical approaches, and developing more realistic descriptions of the metal-liquid interface. However, with the exception of work based on the formalism described in Sec. 3.2 [52,81], recent theoretical methods of the liquid phase structure near planar surfaces have been limited to dipolar solvent models. Furthermore, these calculations have been confined to the mean spheri cal approximation (MSA) or linearized hypernetted-chain (LHNC) theory and have been solved only near inert surfaces. The more robust RHNC theory has never been solved for molecular fluids near planar surfaces. Hence, before considering model potentials for the dielectric-fluid and metal-fluid interfaces, the reduction of the OZ theory will be given  23  Chapter 3. Fluids near planar surfaces  24  for molecular fluids near planar surfaces. Among work dedicated to the study of interfacial models, there has been some con tribution treating the surface as a dielectric continuum [60,61,9,64—66]. However, these studies have been limited to computer simulations of water-like fluids [60,61] and to fluids of spherical ions [9]. A formal analysis of the long-range asymptotic profile of dipolar fluids near dielectric surfaces has also been given [62,63]. In addition, there has been a considerable amount of work devoted to the development of theories for the metal-liquid interface [64—79]. In these studies, particular attention was given to the metal-electrolyte solution interface and the importance of explicitly including the metal side of the in terface was clearly demonstrated. In these calculations, the metal was modeled as an electron gas in a uniform background of positive charge. This so-called “jellium” model was then solved at varying levels of approximation using the general approach of the Hohenberg-Kohn-Sham density-functional (DF) formalism [82,83].  However, we note  that the self-consistent method developed for the jellium-vacuum interface [80] was only applied to the metal side of the interface by Halley, Price and coworkers [73,75]. Hence, although the metal-liquid calculations demonstrated the importance of treating the liq uid side at a molecular level, they did not properly couple the solution structure at the interface to the electronic structure of the metal. That is, while the solvent and ion fields were allowed to influence the electron density distribution at the metal surface, the liquid structure itself could not react to the field generated by the electron distribution in the metal. For polarizable interfaces, the principal problem hindering the development of fully self-consistent theoretical methods arises from the fact that the forces involved are funda mentally many-bodied in nature. For example, the electrostatic interaction of a particle in solution with a polarizable surface depends in principle upon the coordinates of all  Chapter 3. Fluids near planar surfaces  25  other particles in both the solid and liquid phases. Unfortunately, at present most sta tistical mechanical theories applicable to interfaces are only useful if formulated at the level of the translationally and rotationally invariant pair potential. The purpose of this chapter is to propose a nontrivial practical solution to this prob lem. We develop self-consistent mean-field theories for the dielectric-liquid and metalliquid interfaces which reduce the many-body wall-particle interaction to an effective pair potential suitable for use in integral equation theories such as the RHNC approximation. The formalism is given for general multipolar models and explicit results are obtained for a pure dipolar hard sphere fluid and electrolyte solutions near inert and metallic surface. The dipolar hard sphere fluid is also studied near dielectric surfaces characterized by different dielectric constants. The remainder of this chapter is divided into four sections. In Sec. 3.2, the re duction of the OZ equation is given for multicomponent molecular fluids near a planar wall. Sec. 3.3 is devoted to the wall-particle potentials and we begin by developing a general self-consistent mean-field theory for the dielectric wall-particle interaction. A self-consistent mean-field theory is then presented for the metal-particle interaction. In Sec 3.4, results are presented for pure dipolar fluids near inert, dielectric and metallic surfaces.  Results for model electrolyte solutions near inert and metallic surfaces are  discussed in Sec. 3.5.  3.2  Reduction of the planar wall-fluid OZ equation  The Ornstein-Zernike equation for a multicomponent fluid was introduced in Sec. 2.2 and the reduction of the bulk equations was briefly outlined. For molecular fluids in contact with a planar surface, the reduction of either Eqs (2.2.7b) or (2.2.7c) is required.  Chapter 3. Fluids near planar surfaces  26  If Eq. (2.2.7c) is written as h(12) —c(12) =Jcw(13)h(32)d(3), where we have used the property h(12) = ha(21) [see Eq. (2.1.10)], it should be clear that the reduction of either equation is in principle equivalent. However, in practice from a numerical solution perspective, there are small differences in the reduction as discussed by Kinoshita and Harada [51] using the present formulation. In this section the reduction of Eq. (2.2.7b) is presented. Differences in the reduction of Eq. (2.2.7c) will also be noted. A discussion of the numerical solution is given in Appendix C. The equations for a planar surface in contact with a fluid of spherical particles were first considered by Perram and White [84] and by Henderson et al. [85,86]. In their pioneering work, the planar surface is obtained in the limit that a spherical particle becomes infinitely large. Subsequently, Sullivan and Stell [87] provided a fast and efficient Fourier transform method for solving the resulting wall-particle equation.  Since the  surface is uniform, it has no orientational dependence and the rotationally invariant expansion (2.1.6) reduces to g(01)  = = gw(fZi,roi)  (3.2.la)  Chapter 3. Fluids near planar surfaces  27  2  21 r  ,4 ) , i (r i i = 2  1 1 z  Figure 3.1: The coordinate system used in the reduction of the wall-particle Orn stein-Zernike equation. where we have used 0 to denotes the wall coordinates. The wall-particle correlation func tions will be considered exclusively in a reference frame directed along the interparticle 01 vector where r  =  the Euler angles ioi  1 and z z 1 =  =  0 at the surface (see Fig. 3.1). Since  (O) 0 (0,0,0), R,  = ‘n  ,  0) =  corresponds to  &‘o and the basis set reduces considerably to Onn  +  (3.2.lb)  R(q, 0, x)  (_)nfOnn (__4.__1/2 =  2n+1 2n+1}  Y(x,O).  For the present choice of wall-particle reference frame, the basis set reduces to the spher ical harmonics Y(x, 0) [5].  Chapter 3. Fluids near planar surfaces  28  Starting with Eq. (2.2.7b) and expanding the correlation functions in rotational invariants yields  ( 77W/3(01)  )nfOnn Onn  /  (3.2.2)  =  (  )fli  (1Z 2 R ) 1  = =1  fOni ni fm2 2 12  /(2ni + 1)  i”i mnl 12 ‘2  22 Im n 12) v 2 2 2  j R,  I h°  x  (12)  2 (1 d11 R ) 2  1 n 2 m  (3.2.3)  .  j  Since the coefficients of the correlation functions depend only on interparticle spacing, the definition R, (i’) 0  =  eiA’d,o(O) is used to reduce the spatial integral to  J—  m2n212 / ) C;.yi 2 h.(z (i dr 0 R ) 2 iT21) 21 3  J  —  where  921  1 I n 2 m \ 112 / \ 0 h T2l)ao(O2l)’2 7 ; v 2 ovi;w2)c  ,  (3.2.4)  21 (see Fig. 3.1). Subsequent application of the orthogonal is the polar angle of r  properties of the Wigner functions (A.9) yields relations of the form Onn / o;wl)  Blum  v’  x where Pi(cos 0)  (_)fl+fla+i  =  f  (2n + 1) 1 In n 1 + 1) o o o) (2n  , 2 (z C(r21) Pi(cos 0) dr 7 h ) 2  (3.2.5)  (0) is the Legendre polynomial of order 1 and we have used Blum’s 0 d  notation [17] for the  [See Eq. (2.2.14)].  We now take the Fourier Transform of Eq. (3.2.5). However, since not vanish in the limit z  —*  (or c) does  —oo, its transform diverges. This problem is easily solved  for a hard wall interaction where h,(Z,z)  —1,  z < d/2  ,  (3.2.6)  Chapter 3. Fluids near planar surfaces  29  and d is the diameter of species 3. For a soft wall interaction, the same condition will hold at a smaller separation z and is nearly always satisfied for z < 0; hence, for simplicity and generality we use h(f,z)  =  —1,  z <0  (3.2.7)  .  In terms of the expansion coefficients, Eq. (3.2.7) implies that h(z)  z <0  =  (3.2.8)  .  Taking advantage of this condition, Eq. (3.2.5) can be further simplified. Towards this end, it is convenient to define h(z)  h(z)H(±z)  =  (3.2.9)  ,  where 11(x) is the Heaviside step function 11(x)  (1,  x>0,  10,  x<0.  (3.2.10)  =  Equation (3.2.5) can now be expressed as the sum z)  =  5z) +  (3.2.lla)  ,  where p(_)fl+f1+M1) (nini) xJ h Here, the Fourier transform of z  —*  +oo and the remaining integral  below.  4  (z2)c  2 i)Pi(cos O)dr 2 (r  (z) is now well defined since  (3.2.llb)  .  (z)  —*  0 as  (z) can be evaluated in real space as shown  Chapter 3. Fluids near planar surfaces  30  If we started with Eq. (2.2.7c) in the reduction, Eq. (3.2.5) would be unchanged except that the function labels h and c inside the integral need to be switched, and a division of the functions similar to those yielding Eq. (3.2.llb) is necessary. However, in this case the division must be made in terms of the direct correlation function cp(01) but, unlike h(01), c,3(O1) does not satisfy an exact condition inside the surface. Nev ertheless the direct correlation function does approach a limit and we can make use of the approximation c(fZ,z)  where ITd/kBT  =  Z  <  Zt  <  0  (3.2.12)  ,  —(1/V)(OV/OP)T is the isothermal compressibility of the bulk fluid.  This approximation is exact in the limit zt  —* —QO.  Using this approximation, Eq. (2.2.7c)  can be arranged in the form of Eq. (3.2.llb) but again with the function labels h and c  inside the integral switched. In this case, we emphasize that the equal sign is strictly  valid only in the limit  Zt  —  —°°.  5,(z). Expanding the Legendre polynomials in Eq. (3.2.llb)  Let us now consider explicitly as  Pi(cosO)=acosiO=a(Z1Z2)  21 r  i=O  (3.2.13)  ,  where for convenience the coefficients a are given in Table 3.1, and using cylindrical coordmates,  Onn(—) 7 ? Ov;w  (z) reduces to  Onn(—) 7OW/3  1 00 Onn 21rp,j r 2 Sn(zi/r)cou;.y(r)dr Izil -y Izi I —4ir6oH(—zi) c?(r) dr 2 r p  (zi)  j  ,  (3.2.14)  where x—1, a  S(x)  (xi+1  —1)  =  2 (x  —  1)/2,  n=0, n  =  1,  =  n=2,4,6,...  (3.2.15)  Chapter 3. Fluids near planar surfaces  Table 3.1: Pi(x) = Pi(x) Po 1 P 2 P 3 P 4 P 5 P 6 P 7 P  31  The first few coefficients of the Legendre polynomial expansion,  1 1 -1/2  3/2 -3/2  5/2 -30/8  3/8  35/8  15/8  -70/8  63/8  105/16  5/16  -315/16  -35/16  315/16  231/16 -693/16  429/16  and the particle labels have been dropped from the coordinates. The Fourier transform of  Onn(+) iOv.w  (z) is given by  1 Onn(+) 100 e ikz )dzi 1 lOii;w (z 7  _onn(+)(k)  0v;w/3 71  J —00  nñ 1 (2n+1) /n (2ni+1) o oo)  (_)fl+fl1+  =  niai  x  If  1 hofhfh(+) ikz 0vi;wy  ) 2 (z  nl 1 n  j,;y,(r21) 1 c  Pi(cosO) dr 2 dz 1  Since in these equations k spans the reciprocal space of z, we define k e1  etl’21.  =  .  (3.2.16) k to get  Using the Rayleigh expansion [6,88] 00  eik.r =  (21+1) ji(kr) Pi(cosO)  (3.2.17)  1=0  and the identity  dr 2 dr , ff1  dr 1 f f dr  together with the orthonormal properties of  the Legendre polynomials we obtain  OL’;w13 7  )  (_)fl+fl1+v1  =  vil 1 n  Y  X  (0  n l 0  oJ  jOnini+  (2n+1) 1 + 1) (2n (k)  ,  (3.2.18)  Chapter 3. Fluids near planar surfaces  32  where  I  =  Ov,w/  00  Onn;(+)  0 e’h  Onn;(+) is the Fourier transform of 0 h (z) and ;  (z) dz  (3.2.19)  is the Hankel transform defined in  Sec. 2.2.2. Taking the inverse Fourier transform of Eq. (3.2.18) and using Eqs (3.2.11)  1 > 0 can be expressed in the form and (3.2.14), the zi) coefficients for z P00  Onn  /  ‘  7 2irp  =  1ov;w*Z1)  /  (r) dr 3 r) 1 S 2 r / (z cg  1 Jz  -y  (_)f11  x  J  ’ 1 e  Onini(+) Ovi ;w  (2n+1) (ninl\ 1 + 1) 0 0 (2n  (k) k)dk  (3.2.20)  .  z < 0 will also be useful in Chapter 4 and are given by Eq. (3.2.20) The equations for 1 after replacing the first term by Eq. (3.2.14). When all indices are zero, it is easy to verify that these equations reduce to the known result for spherical particles. It is useful to write out the final form we would have obtained if we had started with Eq. (2.2.7c). The approximation (3.2.12) gives (  Onn ( \ CowZ)  Onn(+) C , 0 ;y (z),  {  —  —  z > Zt  (3.2.21)  Z <Zt  where Onn(+) COV;W  (z)  Onn = Cov;w(Z)H(Z — Zt)  (3.2.22)  .  Following the above analysis, we obtain P00  z)  —*  S(z/r) h(r) dr 2 r  2irpJ -y  (_)fl+fli+Vi  navil  l /n n I(2n+1) 1 1 + 1) 0 0 (2n  P00  X  for z> 0.  /  e_i0fh71 CoL,;.q3  ( k ) h4y,y(k)dk  ,  Zt  ,  (3.2.23)  Chapter 3. Fluids near planar surfaces  33  Closure relations for the OZ equation were considered in Sec. 2.2.3 including a general reduction in rotational invariants to yield Eq. (2.2.23) [10]. When c is a planar wall, only coefficients in the rotational invariant expansion of the closure relation Eq. (2.2.17) with rn  =  1 m  =  c(z)  2 m  =  Blum  0 contribute and Eq. (2.2.23a) reduces considerably to give (71171271) (nln2n)  (—)“(2n+1)  V1 ‘2 V 1 V V 2  °° 1 x_jj h(z’)  8w°fl2fl2 ‘‘  ( ) dz’  —  tt  Onn w/3Z  ( )  + bD(z)  (3.2.24)  where the property of the 9-j symbol [5] 1 li•i (On 0  (0  2 71  12  1  =  J  1 11 0 In 2 12 0 (71 1 0  6fl16fl1116fl212  ,  .  J  (3.2.25)  2 + 1) 1 + 1)(2n /(2n + 1)(2n  has been used to reduce p to Onn  Poz  —  —  n 2 fOnin, fOn J J “ fOnn  fl-i-  Z’l+’2 n i 2 ) ij  6n11,6n212  ) (nin n 2  1)(2n + 1)(2n + 1) +2 (_y.(2n+1)(nhn2)  3.3  2 V \“i V  ’ (nin n 2 \ 0 0 01  (71171271)  (3.2.26) (3.2.27)  .  Mean-field theories for polarizable surfaces  In this section, two mean-field approximations for the wall-particle potentials are consid ered. In one, the wall is modeled as a dielectric continuum and in the other as a jellium metal.  3.3.1  Dielectric surfaces  In the present model, the particles of the fluid (immersed in a vacuum) are near a semiinfinite dielectric continuum, characterized by its dielectric constant  .  The wall-particle  Chapter 3. Fluids near planar surfaces  34  potential is given by sr  ttwj3  I  \  UwZ)  diel  +  U  ‘  1 represent the short-range repulsive and electrostatic contributions, where u(z) and u respectively. The functions  and u 1 depend on the orientational and positional  coordinates of all fluid particles. In the dielectric model, when a charge q resides at distance z from the dielectric (at z), an electric field results from both the charge and the polarization of the surface. The field can be calculated using Maxwell’s equation of electrodynamics [89,90]; however, the boundary conditions (and hence the electric field) are equivalent to those if the surface is replaced by a second charge q’ located at —z inside the wall. The so-called image charge q’ which satisfies this condition is given by [89,90] =  U  (3.3.2)  q.  In the presence of the fluid, the polarization of the dielectric surface, and hence depends upon the positions and orientations of all particles in the fluid. Here, we will develop a mean-field approximation which reduces the many-body potential to an effec tive wall-particle pair interaction suitable for use with the homogeneous integral equation theories. However, it is worthwhile to note that an exact formulation of this many-body problem could be obtained using the inhomogeneous equations when the dielectric surface is viewed as an external potential [9]. In modeling the wall as a dielectric continuum, the mean potential generated by the surface is calculated after considering the image-charge distribution of all particles in the fluid (see Fig. 3.2). The electrostatic interaction between a solvent particle and the wall consists of a self-image (SI) term plus contributions from the images of all remaining particles in the fluid (denoted 01 for other images). Therefore, one can write diel_  3 1 uw  —  SIir  \  i  01  13 1,z1) 1- tL 3 tw/  ,  .  Chapter 3. Fluids near planar surfaces  Wall,  35  Fluid  E  2 z  2’ 12 r  —  2  2 2z  12 r 1  Figure 3.2: The coordinate system used in Sec. 3.3.1 for a dipolar fluid near a dielectric surface showing two particle-image pairs. The images are shown to the left of the interface with a prime included in the labels. 1 is the separation of the solvent particle from the surface and fi where z  = (i,  Oi, Xi)  represents the Euler angles. The self-image term is a function of only the coordinates of a particle. If we denote the image to be a second particle of species /3’, then u(O1) can be expressed in the form u(1i,zi)  =  ,—2zi)  ,  12 is the electrostatic pair potential [Eq. (2.4.6)] and r where u(fZi, f2, )  (3.3.4) denotes  the orientation of the image 1’. Here, the self-image potential is reduced by a factor of one half from the electrostatic pair potential to account in part for the work required to polarize the surface. In Eq. (3.3.3), the term u represents the instantaneous interaction of a molecule in  Chapter 3. Fluids near planar surfaces  36  the fluid with the polarized surface for a given configuration of the remaining particles. We propose a mean-field approximation in which the instantaneous interaction between a molecule and the surface is replaced by the result obtained after averaging over the positions and orientations of the “other images”. This yields the effective 01 potential ) 1 u(i, z  J  =  1, r 1  —  ) g(l; 2) d(2) 2 2z  ,  (3.3.5)  where g,(1; 2) is the probability density that particle 2 of species y has coordinates (12,  2 ), and hence 2’ has (fZ, r 2 r  —  ), given that particle 1 has coordinates (1Z, zi). 2 2z  The 01 potential is also reduced by a factor of one half to account for the work required to polarize the surface. Note that we distinguish the effective potential u(Zi, z ) from 1 the instantaneous result u by specifying that it depends only on the coordinates of a single particle. In order to proceed with the reduction of the wall-particle potential to a rotational invariant expansion u(O1)  =  we need to define the multipole moments of the images. To do this, consider the charge distribution n(r) of a molecule of species  ci  in the space-fixed reference frame. The  distribution of charge of the “image species” cV (denoted na’) is obtained after reflecting n,. through the surface at z  =  0. Since the image of a vector r  =  (r, 0, so) is (r, r  —  , 9 t  so),  we have —  =  U ::)  na(r,,so),  (3.3.6)  Chapter 3. Fluids near planar surfaces  37  where the prefactor arises because the charge-image relationship given in Eq. (3.3.2). Using Eqs (3.3.6), (2.4.9) and (A.6d), the multipole moments of the image are given by (a’) =  frna’(r)Rso()dr  —  /1  ( +) frna(r)Rio()dr (1_— +  (_)fl+V’  w\  (a)  (3.3.7)  ‘.  However, in the pair potential we require the multipole moments in the molecular-fixed reference frame. The moments in the space-fixed and molecule-fixed reference frames are related by a rotational transformation, R7(’)  =  2 Q  (3.3.8)  ,  1 Li  where 1’ denotes the orientation of the image. If a has orientation Z =  (q, ir  —  0, r  — x)  =  (, 0, x) then  and we can use Eq. (3.3.7) and the property [5] (_)‘“  ,  R()  (3.3.9)  to obtain  \1 +  —  Qp;a  (3.3.10)  •  The coefficients of the particle-image pair potential can now be written as mnl;el l 3 UILi;/  (r)  ( —) =  =  m’  Om+n,l  fmnl  Qv;’ 13 j (21 + 1)! Qi; (2rn)!(2n)! 11 r  oh (1_— \1+ Ew)  mnl;el ujp;(r).  (3.3.11)  and the self-image potential is given by u(01)  1 i’l  —  =  +  mnl;el (_)Vup; (zi) 1313 mnl  pa,  m .  (3.3.12)  Chapter 3. Fluids near planar surfaces  38  With the help of several properties of the Wigner spherical harmonics, including Eqs (A.8b), (A.11) and (3.3.9) the self-image potential can be arranged in the form of a wallsolvent rotational invariant expansion appropriate for use with the wall-solvent RHNC theory. Doing this gives uj((01)  =  u  =  1 (1  ‘(zi)8(O,1i,O) —  :w)  U’(zi)  (—)  n+L.’  (3.3.13a)  ,  (2fl+1) 3/2  A  u11(2z)  ,  (3.3.13b)  m Ii n 1 l.1 ‘1  where A  (::‘;) [(-)‘ (  =  )(  g)].  (3.3.13c)  Turning now to u(01), the distribution function g(1; 2) used in Eq. (3.3.5) is difficult to obtain exactly but can be approximated by the superposition r ,zi;fZ 1 g( , 2 )  (3.3.14)  .  If only the excluded volume contribution to the particle-particle distribution function is retained, one has the simpler approximation z , g-( ) 2 H(ri  ,r 2 g(i,zi;1Z )  —  d)  (3.3.15)  ,  where H(x) is the Heaviside step function defined in Eq. (3.2.10) and d  =  (da + d)/2.  In the remainder of this chapter, the approximations defined by Eqs (3.3.14) and (3.3.15) will be referred to as the full superposition and cavity approximations, respectively. Starting with the simpler cavity approximation (3.3.15), and expanding the particle particle and wall-particle functions as in Eqs (2.1.6) and (3.3.13a), relationships of the form ug°’(zi)  =  !i  ( ::)  )+ 2 x [h(z  (_)n+vFmnl  J  u(\/r +  )dr O 12 — d,)P,(cos 2 6] H(r o  4z1z2)  (3.3.16)  Chapter 3. Fluids near planar surfaces  39  Here Pj(x) is the Legendre polynomial of order 1, cos O2  are obtained.  =  1 + —(z  2 and z 1 )//r + 4z 2 z fOnnfmnl/  m+n  F  mm f° Blum  m+m  12 Noting that one can write H(r Omm;OI (z ) UO;w/3 1  =  —  -y Ti  —mnl;el  =  =  /1  6w”  1  —  H(d  , a trivial rearrangement yields ri ) 2  —  (_)n+um nI —mnl;el  2 ) + o] dz 2 (z  1+z z )’[h 2  rz +d  where  d) —  d Jz 1  (1  1—lr  +  2)  mnl;el (r), x r 1-H fl,;-y  =  (3.3.17)  V  ,00  xL-J 2 0 —  /2m+1 “mni” 2n + 1 0 0  ( + ) .ç-  7 irp  2m+1 (mnl\ (2n+1) o oo)  Onn  =  1  j(x) dx x P 2  J —1  ]dz 0 +6 [07(2) 2  cos °2’ and A  /  x)dx f’xl_2p ( 1 JA  (3.3.18)  ,  ]  . Using the 2 z 1 /d + 4z z ) 1 +2 (z  orthogonality property of the Legendre polynomials,  I  1  P(x)Fj(x)dx  —1  =  (2j+  )6i3  (3.3.19)  ,  (3.3.20)  together with the expansion [88] 1—2  2 x 1  =  (x) b P 2 i=0  one can demonstrate that the first integral over Legendre polynomials is zero when 1 (Note when 1 > 2, i  1  —  i(x)dx x P 2 < 1 and the integral j 2  =  2  0). The remaining  integrals are evaluated using the related expansion Pi(x)  =  ax  (3.3.21)  ,  i=0  where the coefficients a are listed in Table 3.1. Combining Eqs (3.3.18)—(3.3.21) gives Omm;OI  3 1 Uoj;  ) 1 (z  =  7rYZpy  /1  (+  jzi +d,  —  z1—ci 1 -y 3  mn;cav  Onn )dz 2 hoz,;w(z2)T,;fry (zi,z  ,  (3.3.22a)  40  Chapter 3. Fluids near planar surfaces  where  T;c(x, y)  a (1+ i  mç  = Sm+ni()m  (x+y) —  (3.3 .22b)  1) ,/(d + 4xy) ’ 1  This is the final form of the mean-field interaction using the cavity approximation (3.3.15). Eq. (3.3.22) is valid for a system of uncharged particles. Consider now u(01) using the full superposition approximation of Eq. (3.3.14). Again expanding the particle-particle and wall-particle potential and correlation func tions and reducing with the use of Eqs (3.3.11), one obtains the general result Omm;OI UO;W  (1  (zi) =  P00  I  —  mmn;fu  Onn  0 j  )dz ,(3.3.23a) 2 (zi,z  where mn;ful  F00  r  (zi,z2)  n m 1 m n212 2  2 dr 1 r 12  r  nili;el 1 m  Livi;  ‘  2 (Vr +4z1z2) 12  ’ z’ 2 ? A  Vj 1 fL  X  /  n212 2 m 2M2/3Y()  (3.3.23b)  ] O 2 ( 0 d, O) 12 ( 0 d, )  and =  £Onn (2m  l j 2 ()m+n+fll+LIl+ILfmlnhllfm2n  +  fOmm,  1)3/2  +1  1 n\ 2n 1 m) (n (m m 2 2 l  2 1  \ /rn n 2 m n\ 1 n 2 m’ ,‘n 1 m 2 ’ 1 l n 1 /m ,  )  Here, d() = R(0,O,0) and cos 012  .  o) = (z 2  —  . 12 zi)/r  Eqs (3.3.23) is time consuming. However, the functions  The numerical evaluation of  mmflfUl  lated only once for a given model of the bulk fluid, and when  (3.3.23c)  ) 2 , z 1 (z  need to be calcu  is known,  Omm;OIi UO;W  ) 1 z  can be obtained rapidly from Eq. (3.3.23a). We will discuss results for dipolar fluids in Sec. 3.4.  Chapter 3. Fluids near planar surfaces  3.3.2  41  Jellium surfaces  In this subsection, we consider a simple quantum mechanical model for a metallic surface and assume that the surface and fluid interact only through an electrostatic potential. As in earlier sections, this preempts any mass or charge transfer between the two phases and hence represents a model for what is known as an “ideally polarizable electrode” in electrochemistry. Model, potential and general theoretical approach We represent the metal surface by the so-called “jellium model”, in which the atomic cores (including the nucleus and core electrons) are replaced by a uniform background of positive charge [80,91]. This model is characterized by specifying the average number density of the atomic cores n+, or equivalently, the Wigner-Seitz radius 3 r  = (3/4.)1/3 .  (3.3.24)  Using density-functional theory, the solution of the jellium model gives the ground state (T  =  0) electron density distribution. In this section, the surface is modeled as a metal  slab of finite width 2z (see Fig. 3.3) with the origin chosen at the center. This is the model used by Gies and Gerhardts [91] to study the metal-vacuum interface and from a numerical perspective is a more convenient and accurate model than the semi-infinite wall of Lang and Kohn [80]. In particular, the boundary conditions are well defined in the finite slab (either the wavefunction or its derivative is zero at the center of the slab and the function decays to zero at long range). In the Lang and Kohn method for semi infinite metals, the wavefunction also decays to zero at long range outside of the metal. However in the metal, the self-consistent wavefunction near the surface must match the bulk (plane wave) wavefunction. This matching condition makes it difficult to satisfy charge neutrality with sufficient accuracy.  Chapter 3. Fluids near planar surfaces  42  2  T21  =  ,6 ) i,ci (r i 2 i 1  1 z  r  I  I  0  z  L  Figure 3.3: The coordinate system used in Sec. 3.3.2 for a dipolar fluid near a jellium metal slab. It is well-known from early work on metallic surfaces [80,91] that the electrons of the metal spill several angstroms out from the surface and provide an electric field Ejel strong enough to polarize the liquid near the interface. The liquid in turn polarizes the metal through a potential which depends in principle upon the position and orientation 1 also depends on the instantaneous of all particles in the fluid. Hence, in general, E configuration of the liquid. Here, as in the dielectric continuum case, we will use a mean-field approximation which reduces these inherently many-body terms to an effective one-dimensional interaction.  Chapter 3. Fluids near planar surfaces  43  The interfacial structure (including the electron and molecular distributions) is ob tained by the self-consistent solution of two coupled theories. The liquid phase is de scribed by the RHNC approximation [92,51], and the metal side by the Hohenberg Kohn-Sham implementation of density-functional theory [82,83], as further developed by Gies and Gerhardts [91]. The RHNC theory has already been described in detail. The DF theory consists of the effective one-electron Schrödinger equation _V2wk  +  [e  +  xc  —  k] k  0,  (3.3.25)  where k and T ’k are the one-electron eigenvalue and normalized eigenfunction for the  kth state, m is the mass of an electron, and  Ue  is the exchange and correlation potential  represents the instantaneous interaction energy of an electron in the field of the  jellium and surrounding liquid for a given configuration of the particles in the fluid. We use the local density approximation to solve Eq. (3.3.25) with Wigner’s expression for the correlation energy’ [93] such that [0.611 + (r) 5 r where r (r) 3  =  147 ( 8 4r r) + 23.4ao] (r) + 7.8a 5 [r ] 0  (3.3.26)  ,  (4n(r)/3)/3, e is the magnitude of the charge of an electron,  n(r)  =  2  ‘T!k(r) 12  is the electron density distribution and F is the Fermi energy.  (3.3.27) The factor of 2 in  Eq. (3.3.27) is to account for the electron spin degeneracy. In addition, we propose a mean-field approximation in which the instantaneous in teraction of a particular electron with the external potential is replaced by the result obtained after averaging over all positions and orientations of the fluid particles. This ‘Note that there is a typographical error in Eq. (2.15) of Gies and Gerhardts [91]. In the numerator 7.8 should read 23.4.  44  Chapter 3. Fluids near planar surfaces  yields the mean-field electrostatic potential Ue(Z) = el(Z)  +  distinguished from the instantaneous potential the linear coordinate u(z)  z.  In Eq. (3.3.28),  (Z) 1 UJe  ect(Z)  Ue  and  (3.3.28)  ,  b specifying that it depends only on u(z)  denote the contributions to the  due to the jellium background and fluid, respectively. The Schrödinger equation is  then simplified considerably and reduced to the one-dimensional form [91] 2 d h 2 ——-----i/,,(z) + 2mdz  [Ueff(Z)  —  €,] i/,(z)  =  0 ,  (3.3.29)  ,  (3.3.30)  where Ueff(Z)  + u(z)  =  Ue(Z)  =  iel()  + ut(z) +  —  (3.3.31)  —(k + kg), ‘‘IJk(r) ,  1 2’re  =  u(z)  (3.3.32)  and hk, hk are the x and y components of the electron momentum. The electron density can be reduced to n(z)  =  (F  -  -  e)  I (z)  12.  (3.3.33)  Using normalized wavefunctions, the Fermi energy can be determined from the charge neutrality condition n(z)dz  =  (3.3.34)  2zn÷,  which together with Eq. (3.3.33) gives = CF  fl+Zw 2 (2n  +  /,  (3.3.35)  45  Chapter 3. Fluids near planar surfaces  where  F 11  is the number of eigenstates with energy e, <  F•  The potential contribution from the metal is obtained from the solution of Poisson’s equation (multiplied by the charge of an electron) 2 [nW(z) 4e  =  where W(z)  —  n(z)]  (3.3.36)  ,  11, HI<z,  (3.3.37)  =  10, HI>z. Integrating once and dividing by e, we get the electric field Eiel(z) =  j  4e  [nW(z’)  —  n(z’)] dz’.  (3.3.38)  A second integration gives the interaction potential —  uiel(O)  =  2 j(z 4re  =  e  4e  j  [ZEe1(Z)  —  —  z’)  [nW(z’)  D(z)]  —  n(z’)]  dz’ (3.3.39)  ,  where D(z)  =  z’ [nW(z’)  —  (3.3.40)  n(z’)] dz’.  In order to complete the theory, we must first consider the effective metal-particle, u(01), and electron-fluid, ut(z), interaction potentials needed in the RHNC and DF theories, respectively. Consider first the mean scalar potential =  qt(z)  = _uext(z)/e  Jgwa(02) q(2,1)d(2)  where q( , 1) is the scalar potential at a position 2  1 z  at ,  z  defined by (3.3.41)  from the center of the slab due to  a single molecule of species a with coordinates 2. Since Eq. (3.3.41) is in the form of a convolution, Fourier transform methods can be used to calculate qt(z). However, since  Chapter 3. Fluids near planar surfaces  46  a(2, 1) is known analytically, Eq. (3.3.41) can be reduced significantly in real space. Towards this end, we expand q(2, 1) in rotational invariants, (2, 1)  (_)m2m + 1  =  2,  i0  m+1 r 2 1  fmom  mi  -m Q;amom,  O) 21 r  (3.3.42)  (this is equivalent to the pair potential expansion [Eq. (2.4.6)] when  to  =  obtain (_)m+a  g(z Pm(cos 021) ) 2 . 2 dr 2m+1 r 1  f  eXt() =  2m + 1  (3.3.43)  Using the Legendre polynomial expansion Eq. (3.2.13), the integral reduces to  J  Ommf  Pm(cosO2l) 2 dr 21 r  go;wz2,  m+1  2ir  =  am (m + i  ( (m +ai  27r  +(_)m  I  1 Jz  where we have used the property that  Table 3.1). When m m  1) Jo  —  m  /  =  r°°  z 2 1 0mm 2 iIm+i_1Y0;wa(Z2) dz 2 z 1  ){jzi —  1)  0mm  2 g 1 ) (z dz ] 2  i 2 Iz  (_)i = (_)m  j  0  ,  (3.3.44)  for nonzero coefficients ar (see  2, we have a  1 m 2 X Pm(X)  (m +  =  dx  =  0,  m  2,  (3.3.45)  That is, a uniform planar distribution of quadrupole and higher order moments make no  contribution to qSt(z) and a test charge interacts only with wall-ion and wall-dipole pair distributions functions. Hence,  ext()  is given in general as  r ) 1 ext(  =  —47r  j  (zi  —  z’)  L  4 (YPc;a 3 c4t  paqa2wa(z’)]  dz’  ) dz 2 , 2 j’ g.(z  (3.3.46)  where we have added the constant P00  —27r  I Jo  (zi  —  z’)  [  Pa92(z’)j dz’  —  (YPa;a  j  g,(z2) dz . 2  47  Chapter 3. Fluids near planar surfaces  The interaction potential of a particle with the metal wall is derived in Appendix B and can be written as ua(O1)  =  u  a(zi)(O,i,O)  ,  (3.3.47a)  where —n  Onnel uov;wa(Z)  Qv;a  /2n + 11  0ieI(z)  (3.3.47b)  and iel() =  _uiel(z)/e  (3.3.48)  is the scalar potential due to the charge density of the jellium metal.  3.4  Dipolar fluids and electrolyte solutions near inert, dielectric and metallic surfaces  In this section, we consider three model surfaces in contact with a fluid of dipolar hard spheres and in contact with solutions of charged hard spheres in a dipolar hard sphere solvent. The effects of salt concentration and ion diameter on the interfacial structure are examined. The long-range asymptotic behaviour of a pure dipolar fluid near an inert surface as implied by the RHNC theory is also derived and compared with classical expressions and exact results [62]. 3.4.1  Fluid model  ; 5 The solvent (s) is characterized by the dipole moment p and hard sphere diameter d the cation  (+)  and anion  (—)  are characterized by charges qj, and diameters  4.  If we  choose a molecular frame in which the z-axis is directed along the solvent dipole moment, then a rotation  x  about the principle axis does not change the particle configuration or  Chapter 3. Fluids near planar surfaces  48  the correlation functions and, by symmetry, only coefficients in the rotational invariant expansion with (3.4.1) are nonzero. For this choice of molecular reference frame  =  =  ,u, is the only  multipole moment needed to describe the solvent, and the solvent-solvent, ion,  solvent  and ion-ion, u±±(r), interactions are given by the pair potentials =  + u(r) 2(12)  u(r) + u±(r)  u±(12) u±(12)  uHS(r)  =  iO1(12)  u(r) + ug?(r)  (3.4.2a)  ,  (3.4.2b)  ,  (3.4.2c)  ,  respectively, where the hard sphere potential u and multipole moment interactions are given by Eqs (2.4.2) and (2.4.6). The condition (3.4.1) also implies a significant reduction of the wall-solvent pair func tions and the expansion (3.2.1) can now be written simply as gw(O1) where the subscripts (v angle between  ji  =  =  (_)?2gTh(zi)Pn(cos0i)  (3.4.3)  ,  0) have been dropped from the coefficients and O is the  and the surface normal vector  .  Since the wall-particle interactions  also include a hard repulsive contribution, the wall-solvent,  and wall-ion,  pair  potentials can be written as ‘u(01) u±(01)  =  (Do,  IzH(zw+ds/2,  1. u(z), 1°°  I  z  Izl<z+d±/2,  =  I  z  (3.4.4a)  z + d/2  z  (3.4.4b)  + d/2  where u(z) denotes the electrostatic contribution to the wall-particle potential defined for the dielectric and metallic surface models and we have used z to denote the z  Chapter 3. Fluids near planar surfaces  49  coordinate of the surface. For an inert wall, u(z) = 0. In the case of a jellium metal wall, u(01) is given by Eq. (3.3.47). For dielectric walls, we consider only the pure dipolar hard sphere solvent and u(O1)  u(01) + u(01)  (3.4.5)  ,  where the self-image term is given by u(O1)  (  =  w)  (2)  [2Po(cosOi) + P (cosOi)]. 2  (3.4.6)  For dipolar fluids, one finds that the cavity approximation Eq. (3.3.22) implies that u(01)  =  0  (3.4.7)  .  This is only true for simple dipolar fluids. In general, there are contributions from the images of other particles, even when the cavity approximation is employed. Using the symmetry properties of the dipole-dipole interaction (i.e., u  =  v  =  = ill =  /12 =  =  0), the full superposition result given by Eq. (3.3.23) simplifies somewhat to yield HmuOrnm;OI(z) Pm(cos  u(O1)  O)  (3.4.8)  ,  where ) 1 um;OI(z and  =  Ps 2  (  ) 7;(zi, z 2 f°°g(z ) dz 2 , 2  (3.4.9a)  ) has only the dipole-dipole contribution 2 z I’’(zi, z ) 2  r  = m2n212  x  J  2 zil 1z  ) 2 g2fl2l2(ri  12 r  12 dr  2 [u2;e1(\/r?  (O (O 0 d, ) 12 10 d ) 2 j  .  ) 2 z 1 + 4z (3.4.9b)  Chapter 3. Fluids near planar surfaces  50  In general, the properties of the bulk solution are completely determined by specifying the reduced parameters 1LS/\/kBTd  =  q*  pd  (3.4.lOb)  ,  ±q±/kBTd pd = p_d  4  (3.4.lOa)  ,  =  8 d/d  (3.4.lOc)  ,  (3.4.lOd)  ,  (3.4.lOe)  ,  where Pa denotes the number density of species a. When we consider fluids near metallic surfaces, it will also be necessary to specify the fluid temperature and the solvent diameter since the reduced elementary charge e*  =  e/V”kBTd is needed to define the electrostatic  metal-particle and electron-fluid potentials. (* For all systems considered, the total reduced density of the bulk fluid  =  p + 24),  the solvent dipole moment and the ion charges are constant and have the values p’ = 0.7, and q* fixed, the only and q* 8. For monovalent electrolytes with pa”, = remaining solution variables are the concentration and the ion diameters. The reduced ion diameters considered are given in Table 3.2 with some examples of real ions of similar “size” (assuming d  =  2.8A, the “diameter” of water). These values are based on x-ray  density measurements of ionic crystals [94] and were chosen to be consistent with earlier work [11—13,47—49,95]. The ions labeled “Eq” have the same diameter as the solvent particles. Below we will also refer to a salt labeled “ILi” where the ion diameters are the same as model Lii, but the charges are switched. The characteristic parameters for all systems considered are summarized in Table 3.3. It should also be noted that the model ion pair Na+ and F- have the same diameter; this is also true of Rb+ and C1. Results for the pure dipolar hard sphere fluid and for the salt EqEq in a dipolar hard sphere solvent have previously been obtained using the RHNC theory [10,45] and  Chapter 3. Fluids near planar surfaces  51  Table 3.2: Reduced ion diameters d±/d together with the examples and labels used in this work. The solvent diameter d was taken to be 2.8A.  d±/d 0.68 0.84 1.00 1.16 1.28 1.44  cation Li Na Eq Rb Cs  anion F Eq C1 Br  the reader is referred to these articles for details. A basis set study was also included in Ref. [10] and it was found that truncating the basis at n, m  2 was sufficient to  calculate the pair distribution functions. In the present work, the basis set is truncated with n, rn  <  4 or 6 to ensure accurate results. Results for n, m  from those obtained with  ii,  4 are indistinguishable  m <6.  8 extending at least All RHNC calculations were performed on a grid of 0.01d  212  points from the surface. The bulk particle-particle reference bridge functions were ob tained from the parameterization of Lee and Levesque for hard sphere mixtures [27,28]; the one-component wall-solvent reference bridge functions from the parameterization of Henderson and Plischke [31]. Unfortunately, wall-particle bridge functions are not avail able for hard sphere mixtures. Hence for electrolyte solutions near a wall, the necessary bridge functions were constructed by “shifting” the Henderson-Plischke results [31]. That is, we have used the Henderson-Plischke function (for a reduced density of 0.7) with the contact position shifted to correspond to the appropriate ion diameter. We would expect this to be a reasonable approximation for the present systems where the reduced ion densities and size differences are not very large. The density functional calculations for the jellium metal used a grid width of 0.02A  Chapter 3. Fluids near planar surfaces  52  Table 3.3: Systems considered in the present work. In all cases the surface is in contact with a solution of total reduced density p” = 0.7. The reduced solvent dipole moment = and ionic charge are and q* = 8, respectively. For the metal slab, we assume 3 0 where a 2.65a that T = 298K and, unless otherwise stated, r 0 = 0.529177A is the Bohr radius. salt none LiT NaC1  EqEq  RbF  CsF CsI ILi  p 0.7 0.675 0.675 0.6 0.5 0.675 0.6 0.5 0.688 0.675 0.6 0.5 0.675 0.675 0.675  p 0 0.0125 0.0125 0.05 0.1 0.0125 0.05 0.1 0.006 0.0125 0.05 0.1 0.0125 0.0125 0.0125  molarity 0 0.945 0.945 3.78 7.56 0.945 3.78 7.56 0.453 0.945 3.78 7.56 0.945 0.945 0.945  and an infinite potential barrier was placed 8A out from the jellium surface. At this distance, the electron density distribution is essentially zero in the absence of the potential barrier and is unaltered by its presence. 3.4.2  Asymptotic behaviour of h(01) near an inert wall  It is possible to determine analytically the long-range asymptotic behaviour of the wall particle pair correlation function within the RHNC (or HNC) approximation when the bridge functions rapidly decay at long range. Near an inert wall, the RHNC equations for simple dipolar fluids require that only even-n projections occur in the pair expansion  53  Chapter 3. Fluids near planar surfaces  Eq. (3.2.1) such that Eq. (3.2.20) leads to the simplified relationship s 2 P  =  (zi/r)c(r)dr r S 2  j  (2n+1) (ninl +PsN(2+1) 0 00  c(r )Pi(cosO)dr ) xJhmnh(z i 2 , Since c(01)  —*  —u(O1)/kBT as z  than h (01) and h(01) 5  —*  —*  >0.  1 z  (3.4.11)  oo, we would expect c(01) to be shorter ranged  (01) in the same limit. Hence the required asymptotic  limit can be obtained from Eq. (3.4.11). Explicitly, for those projections which dominate asymptotically, it is not difficult to show that as z h(z)  2rp°  j  —*  oc (3.4.12)  (z/r)c(r)dr, r S 2  where =  {i  0)//2n p2 ( 0 —  +  [i + pssTh0(0)//2Tj + 1] The OZ equation at k  =  (3.4.13)  .  0 has been used to obtain the second equality in Eq. (3.4.13).  From Eq. (3.4.12), the limiting behaviour of h(z) as z  —  oc is determined by the  limiting behaviour of cz(r). Within the HNC/RHNC approximation, the solvent-solvent direct correlation func tion has the asymptotic limit 12) c ( 8  —  h(12)  u(12) ,  r  —*  (3.4.14a)  oo.  r) and h ( For the dipolar model the leading terms are determined by the projections 2 , and Eq. (3.4.14a) yields 3 r) which decay as r u’ ( 2  c(12)  1  [h2(r)j2(12)]  112(  2 —  kBT  ) 2(12)  ,  r  —*  oc  .  (3.4.14b)  Chapter 3. Fluids near planar surfaces  54  112 we obtain  For mril  15  c’(r)  pmnl  1( — )  L  2 (47Tps)  -‘-‘oo  “2I2  r  6 r  Ys 3  —+  oo,  (3.4.15a)  where 2 mnl*(ç , 1  I (ii2(f 2, —00 , 2, ri2) 1 1 U f ‘(1Z  pmnl — .LJ00 —  fz 2 f 12 1 I d 1 )dc d 2 fz i’ 2 1 d 1 )dfz d 2 (ni, 2, 12 2,  (3.4. 15b)  Also, we have used the exact limiting expression [96]  V3  (r) 2 h  —  1) 1  7r,osEsY 4  r  where e is the dielectric constant of the solvent and y we find that all the B’ vanish except for B°°° 00 and (3.4.15b) we find that as z  —*  oo  —*  (3.4.15c)  —  4Kp/9kBT. From Eq. (3.4.15b)  =  p022 ‘‘OO — —  —  1/5. Using Eqs (3.2.15), (3.4.12)  co (_—_1)212 000 1 I z 3 5 327rp 3yc  /  220  F(_—_1)212 I 64rp YEs 3 1  z 22 h° “  ]  L  /  (3.4.16a)  j  L  WS  (3.4.16b)  .  Equations (2.1.6) and (3.4.16) imply that within the HNC/RHNC approximation h(Ofl as  z  F(  1 —  Ps 32  [  —  1)212  Ys 3  1  (  000  +  1 220P(cosO))  (3.4.17)  ,  oo.  —*  It is interesting to compare the HNC/RHNC expression with the exact asymptote derived by Badiali [62] and with the result implied by the “image-dipole” potential of classical continuum theory.  For comparison purposes it is convenient to express the  results as follows hHNC(Ol) WS \ lexact( It1v  \  01)  hdl(01) WS \  —  i 3 kBTz  (S_1)2(S_12  1  4 FjXss  9  000+  3y  \.  (cosO)] 2 X220P  I—1 (—1) 1 ) 16 Z7(n)x P(cosO) 1 3y c(c + kBTz 3 2 ( — 1) 1 4 + P (cos 0)] 2 — kBTz 3 + 1)  )  .  ,  ,  (3.4.18a) (3.4.18b) (3.4.18c)  Chapter 3. Fluids near planar surfaces  55  In Eq. (3.4.18b), n must be even and 7(n)  + a(n, 0),  +  = n,0  (3.4.18d)  where a(n, 0) is related to the three-body direct correlation function given in Ref. [62]. We note that the HNC/RHNC and exact results both depend upon  terms but that  the cx(n, 0) contributions do not contribute in the HNC/RHNC case. Also, it is obvious from Eqs (3.4.18a) and (3.4.18b) that the HNC/RHNC “prefactor” differs from the exact  x°  result. The classical continuum expression does not include  contributions which is  anticipated in the absence of a structured solvent. It is possible to show [62] that when averaged over orientations Eq. (3.4.18b) yields 1  1 WS  000  i  16K c(c + 1) op z 3  (3 4 19) ‘  which can be compared with the HNC/RHNC result given by Eq. (3.4.16a). Furthermore, using the expansion 1 + 3y + 3y (1 + b) + Q(p) 2  (3.4.20)  ,  where b is related to the second dielectric virial coefficient, together with the result XSS —  where  T 1  =  1 (T)p+ 2 1+2B  3421 )  ‘  .  .  (T) is the second virial coeffi 2 is the isothermal compressibility and B  cient, it is possible to derive low density expansions for the h°(z) asymptotes. Explicitly, in the low density large z limit we obtain b 2 {aps + [4a 8T —  where a  =  in density.  8T  —  (T)a 2 2B  (T)a 2 b —2B 2 {aps + [3a b 2 {aPs + [a  —  —  —  ]p + Q(p)} 2 2a a2]p + O(p)}  ]p + Q(p)} 2 a  ,  (3.4.22a)  ,  ,  (3.4.22b) (3.4.22c)  4irz/9kBT. We note that all three results are identical only to linear order  Chapter 3. Fluids near planar surfaces  3.4.3  56  Results for a pure dipolar fluid  The structure of dipolar fluids near inert planar and spherical surfaces has been the topic of several publications [54,45,50,51]. However, these systems have never been solved near planar surfaces using the full RHNC theory and it is useful to understand the structure of the fluid in this case before considering fluids near dielectric and metallic surfaces. The long-range asymptotic behaviour implied by the RHNC theory will also be considered. Pure solvent near inert surfaces The four nonzero projections of the wall-solvent pair distribution function g(z) in cluded in the numerical calculations are shown in Fig. 3.4 (note that z  =  (z  —  z)/d  where z denotes the z coordinate of the wall). The projections for the dipolar fluid near a large inert sphere of diameter 30d [50] and with a hard sphere fluid near a planar inert wall are included. Only even-n projections occur in the expansion because symmetry dictates that  gws(0,  z)  = gws(lr  —  0, z) and because there is equal probability for a dipole  to be directed towards or away from the surface. The oscillation in the angle-averaged distribution function g° indicates ordering of the fluid perpendicular to the surface and arises from the inhomogeneous pressure exerted on a particle near the interface. In the case of dipolar hard sphere fluids, the contact value drops from the pure hard sphere value as a result of electrostatic interactions among the solvent particles and corresponds to a reduction in the bulk pressure Pb. This result is anticipated from the contact value theorem which states for a hard wall interaction that the contact density is proportional to the bulk pressure pg°(d/2)  =  Pb/kBT  (3.4.23)  .  At contact, dipoles favour orientations parallel to the surface (i.e., cos 0  i—’  0 and g 2 <0),  since this configuration lowers the interaction energy among particles in the dense surface  Chapter 3. Fluids near planar surfaces  57  5 0 4 3  —1  2 —2  1 0 0  1  2  3z*  0  1  2  3z*  0  1  2  3z*  0 0.4 —0.02 —0.04  0.2  —0.06 0  —0.08 0  1  2  3z*  Figure 3.4: RHNC results for the projections occuring in the wall-solvent pair distribution and p = 0.7. In (a)-(d), the solid function for a dipolar hard sphere fluid with j = curves represent the present planar wall calculations and the dotted curves are results obtained previously (Ref. 10) for a macrosphere of diameter 30d.  Chapter 3. Fluids near planar surfaces  58  layer. Further out, as indicated by g , the orientations oscillate between parallel and 2 perpendicular configurations and the surface-induced structure decays rapidly to zero. As a result, the full wall-solvent pair distribution function gws(O, z) (shown in Fig. 3.5) is symmetric about cos 0 = 0. As noted in Sec. 3.2, the planar surface can be viewed as an infinitely large sphere. From the projections near the planar and macrosphere surfaces, we see that the distribu tion functions are almost identical. Thus, for a simple dipolar fluid, the 30d macrosphere results already approach the flat plate limit at short range. However at long range, this (r) (r is 2 cannot be true since the macrosphere-solvent (ms) projections h(r) and h , whereas the corresponding 6 the center-center separation) decay asymptotically as r 3 as discussed in Sec. 3.4.2. wall-solvent projections decay as z (z) is shown in Fig. 3.6. The RHNC 2 The long-range behaviour of h°(z) and h asymptotes were obtained from Eq. (3.4.16) using the values and  x°  = 49.24,  x°  0.0738  = 1.349 determined from the solution of the RHNC equations for the bulk  dipolar hard sphere fluid. The small discrepancy discernible in Fig. 3.6b between the analytical and numerical RHNC asymptotes gives an indication of the accuracy of the numerical solutions (i.e.,  i0). The classical continuum limits are apparent from  Eq. (3.4.18c) and were evaluated using the RHNC dielectric constant. Also included in Fig. 3.6a is the curve given by the formally exact Eq. (3.4.19) using RHNC values for the bulk fluid properties. It should be noted that this result is “exact” only in the sense that the formally exact expression is used. The true exact asymptote could only be found if exact values of  Es,  Es/öP, and 9  x°  were known and this is not the case for the present  model. It can be seen from Fig. 3.6 that in the full RHNC result the oscillatory structure characteristic of molecular solvents persists until z  . At larger separations the 5 13d  projections obey the analytically determined asymptotic limits. For h°(z), we observe that both the RHNC and “exact” curves lie below the continuum result, but that the  Chapter 3. Fluids near planar surfaces  59  (O,z) 5 g  z/d __%.,_• •5_-  cos9  Figure 3.5: The wall-solvent pair distribution function g(9, z) for the pure dipolar fluid near an inert wall. The solvent parameters are as in Fig. 3.4  Chapter 3. Fluids near planar surfaces  60  0.0 0.0  —0.2 —0.4  —0.1  —0.6 —0.8 —1.0  —0.2 5  10  15  5  z’  10  15  z  . The dipole moment and 2 Figure 3.6: The asymptotic behaviour of (a) h° and (b) h density are as in Fig 3.4. The solid curves are the full RHNC results obtained numerically. The dotted and dashed curves are the RHNC [Eq. (3.4.16)] and classical continuum asymptotes. The dash-dot curve included in (a) represents the exact asymptotic behavior given by Eq. (3.4.19) and evaluated as discussed in the text. RHNC approximation overestimates the magnitude of the coefficient occurring in the z) h ( asymptotic expression. From Fig. 3.6b, it is obvious that the RHNC result for 2 also lies well below the continuum curve. Unfortunately, unlike the case for hO;ct(z), ;et(z) in the absence of explicit there does not appear to be any way of obtaining h 2 results for the three-body direct correlation function. Pure solvent near dielectric surfaces We now consider the dipolar fluid near hard walls having dielectric constants of 5 and cc using both the cavity approximation for wall-solvent pair potential, u(01)  =  u(1Zi,  zi),  z >0  (3.4.24)  Chapter 3. Fluids near planar surfaces  61  which amounts to retaining only the self-image contribution in the present model, and the full superposition result u(O1)  =  ) + u(fZi,zi), 1 u(fi,z  z > 0  (3.4.25)  ,  where interactions with images of other particles are included. The four projections of the wall-solvent pair distribution function retained in the cal culations are shown in Fig. 3.7 for both approximations and dielectric constants c. The results near an inert surface (c  1) are also included. In general, results for the cavity  =  and superposition approximations are significantly different whereas those obtained for a given approximation with  =  5 and  €  =  co are similar. Also, it is apparent from  Fig. 3.7a that g°(z) does not have a strong dependence upon the value of  .  except  very near contact where the density is nearly double that of an inert surface. The large increase in g°(z) at contact is a direct consequence of the self-image potential which is attractive for all molecular orientations. For the g (z) projection, the superposition 2 and cavity approximation curves are significantly different even several diameters from the surface. Also note that the full superposition results at short range are similar to those obtained for an inert wall. This is important because this projection makes a large contribution to g(01). The curves obtained for g (z) and g 4 (z) with the full super 6 position and cavity approximations are also plotted in Fig. 3.7 and essentially bracket the inert wall results. The physical significance of the results obtained for different models and approxi mations is most easily understood by considering the full wall-solvent pair distribution function g(01) plotted in Figs 3.8 and 3.9 for  =  5 and  w =  respectively. In  the cavity approximation (Figs 3.8a and 3.9a), particles near the surface show relatively little orientational preference when  =  5 and prefer orientations parallel to the surface  Chapter 3. Fluids near planar surfaces  10  000  62  (a) I  2  I  8 8.0  _\  6  0  4.0 I  I  0  I  —2  I  0  2 —4 0 0  1  2  3z*  0  1  2  3z*  0  1  2  3z*  0.0 1.0 —0.1 0.5 —0.2 0.0 0  1  2  3z*  Figure 3.7: RHNC results for the projections occuring in the wall-solvent pair correlation = and p = 0.7. The solid curve function for a dipolar hard sphere fluid with 1). The long- and short-dashed curves were represents results for an inert wall ( = cc and 5, respectively. The long-dash-dot obtained in the cavity approximation for curves are the results obtained with the superposition approximation and short-dash-dot for = cc and 5, respectively.  Chapter 3. Fluids near planar surfaces  normal  when  =  63  oc. Further from contact, the particles shows an oscillatory struc  ture similar to the case for an inert wall but with orientations parallel to  more strongly  favoured. In contrast, with the full superposition approximation (Figs 3.8b and 3.9b), the solvent structure near contact is very similar to the inert case shown in Fig. 3.5. It is clear (for orientations parallel to  )  that the images of other particles act to oppose  the influence of the self-image interaction and again orientations perpendicular to  are  preferred. Further from the surface, the structure is also similar to that found for the inert wall but now dipolar orientations parallel to  are not favoured at any separation.  Comparing the two approximations, it is evident that for dense fluids the self-image term alone is clearly insufficient to describe the interaction with a dielectric continuum surface. Pure solvent near jellium metal surfaces Here we consider the molecular and electronic structure of the metal-solvent interface and compare results for several jellium densities. However, before proceeding, it is necessary to consider the finite slab model used in the jellium calculations since many properties of the metal (such as the Fermi energy, the potential drop and the electron density distribution) vary as function of the width 2z. In Fig. 3.10, the potential drop LS across the jellium-vacuum interface is plotted as a function of z, when r 3 Wigner-Seitz radius of 2.65a , where a 0 0  =  =  . (A 0 2.65a  0.529177A is the Bohr radius, corresponds  approximately to the bulk electron density of Hg.) The potential oscillates with a period of  2A with each period corresponding to the inclusion or exclusion of one eigenstate  (z). For narrow slab widths (z  =  8—hA), the magnitude of the oscillations are  r’.I  0.1  volts (approximately the order of magnitude of weak van der Waals interactions). For larger widths (z  =  29—32A), the oscillations are much smaller (i.e.,  “.‘  0.01 volts) and,  as z, increases, the potential drop approaches the value obtained for the semi-infinite  Chapter 3. Fluids near planar surfaces  64  gws(O,z)  (a)  z /d .-  cos 0  —--.  gws(6,z) 0 f1  (b)  z /d  •  cos 0  Figure 3.8: The wall-solvent pair distribution function gws(O, z) for a pure dipolar fluid near a dielectric surface with E = 5. The results in (a) and (b) are for the cavity and full superposition approximations, respectively. The model parameters for the solvent are given in Fig. 3.7.  Chapter 3. Fluids near planar surfaces  65  gws(6,z)  1.0 10.0  (a)  z Id 8  cos 0  .  gws(O,  z)  100  (b)  z .  cos 0  Figure 3.9: The wall-solvent pair distribution function g(O, z) for a pure dipolar fluid near a conducting surface ( = co). The results in (a) and (b) are for the cavity and full superposition approximations, respectively. The model parameters for the solvent are given in Fig. 3.7.  Chapter 3. Fluids near planar surfaces  66  3.32  3.32 :  3.28  3.28  3.24  3.24  3.2  3.2 8  9  10z(A)  I  -Lçb  II  I  liii  —  (volts)  :11  29  30  31z(A)  Figure 3.10: Oscillations in the net potential drop across the jellium-vacuum interface for r 8 = 2.65a . 0 . 2 j ellium-vacuum interface The Fermi energy and potential drop across the jellium-liquid interface show similar behaviour. In particular the magnitude of the oscillations in the Fermi energy is volts at narrow slab widths and  i-.’  0.016  0.002 volts at the larger widths. In comparison to the  solvent effects considered below, the dependence on slab width is very small, particularly at the larger widths, and provides a good approximation for a semi-infinite metal. All results considered for the jellium-liquid interfaces are for a slab 64A wide. Also in these calculations, the solvent particles have a diameter of 2.8A and the temperature is taken to be 298K. The solvent diameter and the temperature are required in order to specify the jellium-solvent interaction. The self-consistent electron density distributions for a jellium slab immersed in the dipolar fluid are shown in Fig. 3.lla for three values of r. The electron density shows the characteristic Friedel oscillations inside the metal and quickly decays to the bulk density The potential drop across the semi-infinite jellium-vacuum interface was calculated using N.D. Lang’s 2 original self-consistent program [801.  Chapter 3. Fluids near planar surfaces  67  further in the slab. The magnitude and frequency of the oscillations depend strongly . Just outside the jellium surface, the density decays exponentially and drops to 5 on r 0.1% of the bulk value at  3A.  It is interesting to compare the electron distribution at the jellium-solvent interface with the jellium-vacuum case. The inset in Fig. 3.lla shows the change in the electron density distribution resulting from the presence of the dipolar fluid. It can be seen that the solvent only slightly perturbs the distribution and that the maximum deviation is approximately three orders of magnitude smaller than the bulk electron density. Also from the inset, we see that the electron density spills further out of the metal due to the presence of the solvent. The effective potentials obtained for the jellium-solvent systems are plotted in Fig. 3.llb. The result for the jellium-vacuum interface with r  =  0 2.65a  is also shown. The oscillations which occur just outside the jellium edge arise from the fluid structure. The effective potential is dominated by the exchange-correlation energy. The electrostatic contribution to  Ueff  will be considered in detail below (see discussion  related to Fig. 3.28). We now consider the self-consistent fluid structure at the interface. The first few nonzero projections of the pair correlation functions are shown in Fig. 3.12 and the results are compared with those for an inert surface. (The inert surface can be thought of as a semi-infinite jellium metal with r  =  oo). A direct consequence of the interaction  of the two phases is that the probabilities of a dipole pointing into and out of the surface are no longer equal. Hence, odd-ri projections now occur in the Legendre polynomial -dependent trends and rapidly 3 expansion. Near the surface, the curves exhibit general r approach the inert wall results at larger separations. From Fig. 3J2a, it is apparent that near the surface g°(z), and hence the local solvent density, increases with deceasing . We also note that the second peak in the function moves closer to the surface as 8 r . The 3 r decreases, implying that the solvent packs more tightly for smaller values of r  Chapter 3. Fluids near planar surfaces  68  0  1 0.8  —0.5  0.6 0.4  —1  0.2  —1.5  0 —15 —10 —5  0  5 z  —15 —10 —5  0  5 z  Figure 3.11: The self-consistent electron density (a) and effective potential (b) of ajellium slab, 64A wide, immersed in the pure dipolar hard sphere fluid. The dotted, solid and dashed curves are for r 5 = 2.0, 2.65 and 3.0a , respectively. In (b), the dash-dot curve 0 shows the self-consistent effective potential when no fluid is present and r 3 = 2.65a . The 0 distribution of electrons in (a) is very similar to the distribution when the fluid is absent. The inset in (a) shows the change in the distribution due to the presence of the fluid 8 2.65a (r ). 0 projection g’(z) [Fig.3.12b] is a measure of the local orientational polarization of the solvent. For inert surfaces, g (z) = 0 since symmetry requires that the solvent has no 1 (z) [Fig. 3.12c], the most striking feature is the change of sign 2 net polarization. For g near the surface as the solvent becomes polarized through interactions with the jellium. It is interesting to introduce the orientational probability density p(O, z) = pws(O, z). For the present dipolar model,  p(O, z)  p(O,z) =  is given simply by  1 (_)  Onn  O()Pn(cos)  (3.4.26)  .  Orientational probability densities for four wall-solvent separations are compared in 3 Fig. 3.13. The probability density near an inert surface is included. For r the complete pair distribution function, in the interval 0.5d  z  —  z  =  , 0 2.65a  3.5db, is plotted  69  Chapter 3. Fluids near planar surfaces  10  0  8  (b)  011  6  —5  4  —10  2  —15  0 10  111111111  g  —10  0.5  3z*  2  1  0 __L  III  (c)  10  0.6  3z*  2  1 1 0  5  —1 0  —2 1111111  I  3z*  2  1  11111  0  1  2  3z*  0.6 0.4 0.2 0 —0.2 —0.4 0  1  2  3z*  Figure 3.12: RHNC results for the projections of the wall-solvent pair correlation function for a dipolar hard sphere fluid in contact with a jellium metal. The dotted, solid and 5 = oo). dashed curves are as in Fig. 3.11. The dash-dot curve is for an inert surface (r  Chapter 3. Fluids near planar surfaces  70  I  2.0  2.0  1.5  1.5  1.0  1.0  0.5  0.5  0.0  0.0 0  —1 I  2.0  111111  111111  1.5  1.0  1.0  0.5  0.5 I  ii  —1  I  11111 I I  0  iii  cosO  ii  I  1  Ill  III  I  I  I  0 I  I  I  I  I  I  I  0.0 —l —1  I  I  (b)  I  2.0 Ep(,1.5d ) 8  1.5  0.0  II  p(, O.8d)  I  III  (c) Z  -p(6, 1.2d ) 5  II  —1  1  cosO  III  I  cosO I  I  I  I  I  I  1 I  I  (d)  I  I  0  I  I I I  cos9  1  Figure 3.13: The probability density p(O, z) for solvent particles at distances (a) 0.5d, (b) 0.8db, (c) 1.2db and (d) 1.5db from the jellium edge. The curves are as in Fig. 3.11.  Chapter 3. Fluids near planar surfaces  71  gws(,z)  10.0  0  z /d  e5  cos6  Figure 3.14: The wall-solvent pair distribution function g(O, z) near a jellium metal wall with i’ . The solvent parameters are as in Fig. 3.7. 0 3 = 2.65a  Chapter 3. Fluids near planar surfaces  72  in Fig. 3.14. The structure of the dipolar fluid near an inert surface was discussed ear lier (see page 56 and Fig. 3.5). In contrast, near the jellium-liquid interface the solvent orientations are highly asymmetric. At contact, the dipoles orient preferentially along the surface normal vector (i.e., the positive end of the dipole points away from the metal surface) giving rise to significant surface induced polarization of the solvent. This is qual itatively consistent with calculations for dipolar monolayers [67,68]. As we move further from the surface the most probable dipole orientation oscillates. For example, at O.8d and 1.2d 5 [Figs. 3.13b and 3.13c] the dipoles have rotated more towards the surface, but at 1.5d [Fig. 3.13d] they again prefer to align along the surface normal. We note that the solvent structural features again decay rapidly with distance from the surface although the correlations are a little longer ranged than in the inert case (see Fig. 3.14). 3 Electrostatic potentials obtained for r  =  0 are plotted in Fig. 3.15. The total 2.65a  potential across the jellium-liquid interface (solid curve) is the sum of the contributions from the electron density (dotted curve) and from the solvent dipoles (dashed curve). We see that as expected the net potential drop across the interface is reduced by the presence of the solvent. Structural features related to solvent granularity and orientational order are also evident. The potential drop, /q  =  4(z  oo)  —  q!(z  =  0), across the jellium-dipolar fluid  interface is given by = =  471-  S  Sf  3 +4ej  g’(z’)dz’  0  z’[n(z’)  —  nW(z’)]dz’  (3.4.27)  where /qfd and zSS are the dipolar and jellium metal contributions. Results for  LSd,  Lq!j and Z.4 at selected r 3 values are given in Table 3.4. The potential drop across the jellium-vacuum interface is also included. We note that the results obtained for dipolar  Chapter 3. Fluids near planar surfaces  73  3 100 2 50  1  0  0 —5  0  5z*  20 10 0 —10 —20 0  2  4  6  z  Figure 3.15: The reduced electrostatic potential across the interface when r 3 = 2.65a . 0 (a) The net potential for the jellium slab in a dipolar hard sphere fluid (parameters as in Fig. 3.11) is represented by the solid curve and is the sum of contributions from the jellium (dotted curve) and the orientationally polarized fluid (dashed curve). The net potential across the jellium-vacuum interface is represented by the dash-dot curve. The scale on the right shows the electrostatic potential çS(z) in volts at 298K. (b) The electrostatic potential for z’ > 0 on an expanded scale; curves are the same as in (a). monolayers are of comparable magnitude and exhibit similar trends [68]. The effect of varying the solvent diameter and/or density on the total potential drop is shown in Table 3.5. If the solvent diameter is reduced while holding p’ constant, the dipoles can approach the jellium edge more closely and the potential drop decreases. Also, as we would expect, the potential drop decreases if p’ is increased with the diameter held fixed. It is interesting to compare the solvent structure at the jellium and dielectric surfaces. The wall-solvent pair distribution functions at a jellium surface with r 3  =  0 and at 2.65a  a continuum conducting surface using the cavity and superposition approximation are plotted again in Fig. 3.16 for comparison. The dielectric wall can only influence the solvent’s preference to point either parallel to the surface or in and out from the surface but can not induce a net polarization. Hence, only even-n projections are shown in  Chapter 3. Fluids near planar surfaces  74  Table 3.4: The self-consistent electrostatic potential drop across the jellium-vacuum and jellium-solvent interfaces for z = 32A. When present, the liquid consists of dipolar hard = spheres with p = 0.7, d = 2.8A and T = 298K. L4 (volts)  j ellium-vacuum j ellium-solvent 2.00 2.65 3.00  —6.301 —2.950 —2.045  —6.787 —3.275 —2.318  Table 3.5: The self-consistent electrostatic potential drop across the jellium-solvent in terface for r 3 = 2.65a 0 and z = 32A. In all cases the liquid consists of dipolar hard spheres with t = at T = 298K. p 0.663 0.7 0.7 0.7 0.75  8 d  (A) 2.8 2.6 2.8 3.0 2.8  Liq!’  39.9 50.3 50.3 50.3 77.7  (volts) —2.969 —2.898 —2.950 —3.007 —2.940  Fig. 3.16. In both approximations the contact density increases from the inert surface value as a result of the attractive wall-solvent interaction. Comparing the (022), (044) and (066) projections, we note that the solvent structure is significantly different in all cases. In addition, we note that the jellium curves for even-n projections (short-dash curves) resemble the inert wall results (solid curves) except very near contact where they are most similar to the cavity approximation (long-dash curves).  3.4.4  Results for electrolyte solutions  In this section we will discuss electrolytes near inert and metallic surfaces.  Chapter 3. Fluids near planar surfaces  10  000  (a)  75  2  8 8.0  -\\  0  6 4.0  —2  4  0  2 —4 0 0  1  2  3z*  0  1  2  0  1  2  3z*  0.1 1.0  0.0  0.5  —0.1  0.0  —0.2 0  1  2  Figure 3.16: Comparison of the projections occuring in the wall-solvent pair correlation function for a dipolar hard sphere fluid near a jellium metal (r 8 = 2.65a ) surface and a 0 dielectric surface. The short-dash curve is the jellium metal result. The long-dash and dash-dot curves are the results near a dielectric surface using the cavity and superposition approximations, respectively. The solid curve is result near an inert surface.  Chapter 3. Fluids near planar surfaces  76  Electrolyte solutions in contact with an inert surface The first three wall-solvent projections obtained for a O.945M CsF solution near an inert wall are shown in Fig. 3.17. The results for pure dipolar hard spheres are included. In general the solvent structure is only slightly modified by the electrolyte. As the anion and cation diameters become equal, the influence of the electrolyte is further reduced. As for pure dipolar fluids, ion-dipole mixtures with ions of equal size cannot be po larized by the presence of an uncharged inert wall and the potential drop across the interface is zero. However, if the ions are of unequal size the symmetry is broken and a small surface-induced solvent polarization results. For such systems, dipoles will have a preference to point either into (g 1 > 0) or out of (g’ < 0) the surface. For CsF at 0.945M (Fig. 3.17) we see that g’ is positive near contact indicating that the dipoles prefer to have their positive ends directed towards the surface. Further from the surface, g’ oscillates and decays to zero from below indicating a small but long-ranged polar ization of the solvent. At higher concentrations, the decay to zero becomes more rapid as ion screening grows in importance. If the salt concentration is not too high it is possible to predict the sign of g’ (d/2) from the relative ion sizes. For example, from Table 3.6 we note that at 0.945M g (d/2) 1 is negative for Lii and NaC1 and positive for CsF. This is as one would deduce from the sign of the smaller ions; i.e., small positive ions located between the surface and the solvent particles will direct the dipoles outward and small negative ions will do the reverse. For CsI, where both ions are larger than the solvent but Cs+ is smaller than 1  (see Table 3.2), g’(d/2) is positive. Again this is predictable if one considers a  “layer” of Cs+ ions lying just outside the solvent “layer” in contact with the surface. However, from Table 3.6 we also note that the sign of g’(d /2) is not so easily predicted 8 at high concentrations. For example, for NaC1 solutions g’(d/2) becomes positive at  Chapter 3. Fluids near planar surfaces  77  2.0 0.02 1.5 0.00 1.0 —0.02 0.5 0  1  23  4z*  0  1  2  4  0  1  2  3  0.0  —1.0  3  z  Figure 3.17: Projections of the wall-solvent pair distribution function for inert walls. The parameters are as in Tables 3.2 and 3.3. The solid curve is the pure solvent result and the dashed curve is for O.945M CsF.  Chapter 3. Fluids near planar surfaces  78  Table 3.6: Contact values of the wall-solvent pair distribution function near an inert surface. model pure solvent LiT NaC1  EqEq CsF CsI  molarity 0 0.945 0.453 0.945 3.78 7.56 0.945 0.945 0.945  g° 4.39 4.89 4.50 4.58 4.87 5.23 4.43 4.73 5.46  g°” 0 —0.073 —0.017 —0.024 0.007 0.114 0 0.025 0.005  22 g° —2.58 —2.41 —2.39 —2.27 —1.92 —1.78 —2.18 —2.37 —3.06  higher concentrations indicating that correlation effects become more important as the concentration is increased. Apart from these relatively small polarization effects, the addition of salt does not strongly influence the interfacial solvent structure. We see from Fig. 3.17 that the even-n projections are only slightly influenced at 0.945M and in fact salt effects are not very significant even at 7.56M. The value of g°(d/2) (Table 3.6) is increased for all systems considered, most significantly for the 7.56M NaC1 and 0.945M CsI (where both ion di (d/2) decreases in all cases except for 2 ameters are larger than d) solutions. Also, g (d/2)/g°(d/2), which is 2 0.945M CsI. The more relevant quantity for comparison g proportional to (P 2 (cos 0)) per particle, decreases in all cases indicating that in the pres ence of salt the dipoles have a reduced tendency to align parallel with the surface. This is obviously consistent with the polarization effects discussed above. The wall-cation and wall-anion distribution functions for a few salt solutions at 0.945M are shown in Fig. 3.18. At 0.945M the anion and cation distributions are nearly independent of the coion diameter; for example, the distributions of Cs in CsF (not shown) and CsI are similar. Also note, that in the simple dipolar solvent the distribution  Chapter 3. Fluids near planar surfaces  79  5  5  4  4  3  3  2  2  1  1  0  0  0  1  3z*  2  0  1  3z*  2  Figure 3.18: Wall-ion pair distribution functions for O.945M electrolyte solutions in con tact with an inert surface. The parameters are as in Tables 3.2 and 3.3. In (a) the curves are: (—) Na in NaC1; (----) Li in Lii; (----) Eq in EqEq; (— —) Cs in CsI. In (b) the curves are: (—) Cl— in NaC1; (---) Eq in EqEq; (-.-.) F— in CsF; (— —) 1 in CsI. functions are identical for ions of equal size and opposite charge. The most significant feature evident in Fig. 3.18 is a sharp increase in the contact values as the ion diameter is increased. This comes about simply because the smaller ions (which can approach the solvent particles more closely) are better solvated than the larger ones in the dipolar solvent. The potential drops across the inert wall-electrolyte solution interface  (/-4  =  /qd  +  including the ion and dipole contributions, are given in Table 3.7. The dielectric constants of solution are also included. If the ion pairs are of equal diameter but opposite charge (e.g., RbF and NaC1), that although  IL4  has the same magnitude but opposite sign. We note  iI and ILci range between 0 - 146mV, the ion and dipole contributions  are always opposite in sign and cancel such that qS is always less than we would expect,  ILq  5OmV. As  generally increases with the difference in ion size. Further, Lq  passes through a minimum at  ‘S-’  0.945M. It is also evident that  ILq  is not simply related  Chapter 3. Fluids near planar surfaces  80  Table 3.7: The potential drop in volts across inert surface-electrolyte solution interfaces. q* = 8. = /ä and The dielectric constant is also included. In all cases p = 0.7, model pure solvent LiT NaCl  EqEq RbF  CsF CsI ILi  molarity 0 0.945 0.453 0.945 3.78 7.56 all 0.453 0.945 3.78 7.56 0.945 0.945 0.945  c 50.3 30.23 37.22 30.64 17.23 11.47 37.22 30.64 17.23 11.47 30.94 35.36 30.23  Lq 0 —0.060 —0.034 —0.032 —0.068 —0.146 0 0.034 0.032 0.068 0.146 0.042 —0.019 —0.060  0 0.048 0.036 0.029 0.042 0.099 0 —0.036 —0.029 —0.042 —0.099 —0.036 0.014 0.048  to the dielectric constant of solution. For example, although  f  0 —0.012 0.002 —0.003 —0.026 —0.047 0 —0.002 0.003 0.026 0.047 0.006 —0.005 —0.012  decreases with increasing  salt concentration, Lq does not. Electrolyte solutions near a metal surface We begin by considering the first three projections of the metal wall-solvent pair distri bution function. Results are shown in Figs 3.19 and 3.20 for NaCl and RbF solutions, respectively, at concentrations of 0.945M and 7.56M. Pure solvent curves are also in cluded. The full pair distribution function for 0.945M and 7.56M NaCl are shown in Figs 3.21 and 3.22, respectively. We note that in all cases the general structure is qual itatively similar to that of the pure solvent. However, the salt does induce significant perturbations which depend on both the cation and anion diameters. The g’ and g 2 projections are of most interest here.  Chapter 3. Fluids near planar surfaces  81  0.5  2.0  0.0  1.5  —0.5 1.0 —1.0 0.5 0  1  23  0  1  2  4z*  0  1  2  3  0.0  —1.0  3  Figure 3.19: Projections of the wall-solvent pair distribution function for NaC1 solutions in contact with a jellium slab. The solid curve is the pure solvent result. The dotted and dashed curves are for O.945M and 7.56M solutions, respectively. The parameters are as in Tables 3.2 and 3.3.  Chapter 3. Fluids near planar surfaces  82  0.5  2.0  0.0  1.5  —0.5 1.0 —1.0 0.5 0  1  2  3  4  z*  0  1  2  3  4  z  0  1  2  3  0.0  —1.0  Figure 3.20: Projections of the wall-solvent pair distribution function for RbF solutions in contact with a jellium slab. The solid curve is the pure solvent result. The dotted and dashed curves are for 0.945M and 7.56M solutions, respectively. The parameters are as in Tables 3.2 and 3.3.  Chapter 3. Fluids near planar surfaces  83  gws(O,z)  10.0  z /d 5  cos -j.  Figure 3.21: The wall-solvent pair distribution function g (O, z) for the O.945M NaC1 5 solution in a dipolar solvent and near a jellium metal wall.  Chapter 3. Fluids near planar surfaces  o.o  gws(O,  84  z)  1 1  z Id 5  cos 6 -.-  Figure 3.22: The wall-solvent pair distribution function g(O, z) for the 7.56M NaC1 solution in a dipolar solvent and near a jellium metal wall.  Chapter 3. Fluids near planar surfaces  85  The g’ projections show the strong dipolar polarization induced by the metal sur face. The large negative values at contact indicate that the dipoles strongly prefer to orient perpendicular to the surface. From Figs 3.19 and 3.20 we see that this effect is magnified as the salt concentration is increased. Further from the surface g’ oscillates and approaches zero. For the pure solvent these oscillations are centered about zero and quickly decay. At 0.945M, g’ is longer ranged and oscillates about negative values for both NaC1 and RbF solutions. At 7.56M the solvent is more strongly ordered near 2 projections ex contact, but g’ is highly screened and decays very rapidly. The g hibit qualitatively similar trends and become more positive at small separations as the concentration is increased. These wall-solvent results will be discussed further after we consider the wall-ion distributions. The short-range potential arising from the jellium is attractive to anions and repulsive to cations. Thus if the solvent were not explicitly present (e.g., a primitive model) we would expect a high anion and low cation density near the surface. However, from Fig. 3.23 we see that this is clearly not the case for the present system. At 0.945M the most obvious features, common to all solutions considered, are the large cation peak near the metal surface and a major anion peak located a full solvent diameter from contact (near z  =  1 + (d_/d )/2). This structure comes about because of the strong ordering of the 3  solvent particles near the metal surface and can be understood by considering the ion dipole interactions. In bulk solution (see Fig. 3.24), the dipoles in the first solvation shell prefer to be directed towards an anion and away from a cation. As an ion is moved towards the surface, an increasingly strong force tends to orient the dipoles normal to the surface. For anions the orienting fields of the surface and of the ion act more or less constructively on dipoles located between them for separations  8 + d_ /2. (Note that d  the surface field is very short ranged and also partially screened on the far side of the anion.) Closer to the surface the fields act destructively and destabilize the solvation  Chapter 3. Fluids near planar surfaces  86  shell. These observations explain all but the very short-range behaviour of the anion distributions shown in Fig. 3.23. At small separations the direct metal-anion potential is strongly attractive and consequently the anion density increases rather sharply at contact. For cations essentially opposite arguments apply. At intermediate separations, the surface and cation fields act to destabilize the solvation shell. Nearer the surface (i.e., at  ‘‘  ) the polarization of the dipoles creates a new region of stability (see Fig. 3 O.6d  3.24). However, at still smaller separations the direct metal-cation interaction begins to dominate and the cations are repelled from contact. The full cation distributions near contact for RbF, CsF, CsI and ILl are shown in Fig. 3.23h. We see in general that the contact peak heights for cations near the metal follow the same trend as in Fig. 3.18 for an inert surface (i.e., the contact density increases sharply for larger ions as a result of the weakened ion-solvent interaction). We also note that cation “dewetting” due to the repulsive interaction with the metal occurs only for smaller ions (d peak at  “-‘  1.2d), an indication of the range of the jellium potential. The anion  3 + d_/2 is relatively insensitive to the ion diameter, but increases in height d  as the cation density near the surface increases. Cation and anion distribution functions obtained for different values of r 3 are shown in Fig. 3.25. We observe that the qualitative features are similar for the different values 3 is decreased from 4.0a of r , but that the structure sharpens considerably as r 3 0 to 2.0a . 0 This reflects the increasing degree of orientational order in the solvent contact layer as the metal-dipole interactions increase in strength. The results discussed above show that at O.945M the ion-solvent interactions largely determine the ion-wall distribution functions. The metal orders the solvent which in turn directs the ions. The direct ion-metal and ion-ion interactions tend to be dominated by the ion-solvent forces. The only important exception to this occurs for small ions near contact where the direct ion-metal interaction sharply repels cations and attracts anions.  Chapter 3. Fluids near planar surfaces  4  11111  liii 11111  87  111111 I I  I  I  I  I  I  I  i  i i  i  j  I  I  I I  I  I  I  Ii  I  II I  (a)LiI.  2  2  0  0  I  I •iI  I  i  i  i  I  i  i i  I  i i  i  III\,IIIIIIii;IiiiII  2  2  0  0  2  2  0  0  I I I Iii I I I  I  i  i  I  i  i  I  i  :L’>’ IIIIIIIIIIIIIIIIIIIIIIII  0 20  1  I  2  3  1111111  I  11111  (h)  2  10  —CsF  —  RbF  0 0  1  2  3  4z*  0  I  0  I  I  I  —I—1  1  I  2  I  I  z  Figure 3.23: The wall-cation (—) and wall-anion (---.) pair distribution functions for several O.945M electrolyte solutions in contact with a jellium slab. The parameters are as in Tables 3.2 and 3.3. In (h), the wall-cation distribution functions for RbF, CsI, CsF and ILi are shown on an expanded scale.  Chapter 3. Fluids near planar surfaces  88  e  G 0 0  eØ Bulk Figure 3.24: An illustration of the distortions of the solvation shell due to the jel hum-induced interfacial polarization. Results for higher concentrations (up to 7.56M) of NaC1 and RbF are shown in Figs 3.26 and 3.27. Although screening and other concentration dependent effects can be seen in the plots, we note that the basic structural features persist at the higher concentrations. It is also worth noting that the ion-wall distribution functions are in full accord with the solvent structure illustrated in Figs. 3.19 and 3.20. The excess anion concentration near z  —  1.5 induces an additional polarization of the solvent and hence g 1 is negative  in this region and approaches zero from below. Further, for 0.945M RbF (d_ g’ is more negative than in the case of 0.945M NaC1 (d_  =  =  ) 3 0.84d  1.16db). For 7.56M RbF  Chapter 3. Fluids near planar surfaces  5H 4-  g  (a)  Il  89  5 4 3  2z  2  1-  1  0_I  0  1  2  3z*  0 0  1  2  3z*  Figure 3.25: Wall-ion pair distribution functions for NaC1 solutions at O.945M. The dotted, solid and dashed curves are for r 5 = 4.0a , 2.65a 0 0 and 2.0a , respectively. 0 the region near z’  =  2 is also of interest. Here g’ has a positive peak indicating that  dipoles prefer to orient towards the surface. This peak is presumably associated with the strongly bound first solvation sphere of the small F— ion. The potential drops across the metal-solution interface are given in Table 3.8. The contributions from the dipoles,  ions,  and jellium,  j,  are included. The  net contribution from the electrolyte solution is also included and denoted by LSd. We observe that the presence of salt increases zg,j by about a factor of three above the pure solvent value. However, this increase is largely canceled by zj such that the total solution contribution to that  dj 4  remains near that obtained for the pure solvent. Also, we note  is roughly an order of magnitude larger than corresponding values for an inert  wall (see Table 3.7). The partial cancellation of /qf’d and qj serves to emphasize that all effects must be taken into account in a self-consistent manner if accurate theoretical results are to be obtained. The ion and solvent contributions are strongly coupled and one cannot focus separately upon one or the other as is frequently done in double layer theories. Several other “trends” are evident from Table 3.8. For example, we see that  Chapter 3. Fluids near planar surfaces  90  4  4  3  3  2  2  1  1  o  0-j0  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  (c) 7.56M  r’ii;iiiiiiIi,ii: 1  2  3  4z*  3 2 1 0 0  1  2  3  4z*  Figure 3.26: The wall-ion pair distribution function for NaC1 solutions near a jellium surface. The curves are labeled as in Fig. 3.23. Lf is not strongly concentration dependent. Also, the change in Lq with respect to the jellium-vacuum interface,  tends to decrease with increasing concentration,  although for RbF a minimum occurs at  r’.i  3.78M. Furthermore, metal induced charge  asymmetry effects are evident in that the potential drops for NaC1 and RbF, Lii and ILi, are significantly different. We recall that for an inert wall the magnitude of the potential drop is the same for salts consisting of ions of equal size but opposite charge. The mean electrostatic potential and its various contributions for a 0.945M NaCl solution are plotted in Fig. 3.28. Similar results were obtained for the other systems considered. As in the pure solvent case (see Sec. 3.4), the jellium contribution decays very rapidly on the solution side of the double layer and is more or less negligible for z  0.6. However, it is the rapid decay of this potential which gives rise to the strong  electric field orienting the dipole in the contact layer. The separate contributions of the  Chapter 3. Fluids near planar surfaces  91  8  8  6  6  4  4  2  2  0  0 0  1  2  3  4z*  6 4 2 0 0  1  2  3  4z*  Figure 3.27: The wall-ion pair distribution function for RbF solutions near a jellium surface. The curves are labeled as in Fig. 3.23. ions and dipoles are large but cancel to such an extent that the net result does not differ greatly from the pure solvent curve also shown in the plot. This demonstrates again the strong coupling of the ionic and dipolar potentials emphasized above. osta  Chapter 3. Fluids near planar surfaces  92  Table 3.8: The potential drop in volts across jellium surface-electrolyte solution interfaces. The dielectric constant is also included. In all cases T = 298K, d = 2.8 and r 5 = 2.65a . 0 model vacuum pure solvent LiT NaC1  EqEq  RbF  CsF CsI ILi  molarity 0 0.945 0.945 3.78 7.56 0.453 0.945 3.78 0.453 0.945 3.78 7.56 0.945 0.945 0.945  1.0 50.3 30.23 30.64 17.23 11.47 37.57 31.08 18.17 37.22 30.64 17.23 11.47 30.94 35.36 30.23  0 0.412 0.955 1.047 1.191 1.013 1.003 1.187 1.310 1.199 1.376 1.438 1.299 1.524 2.011 1.520  LS4 0 0 —0.550 —0.635 —0.836 —0.716 —0.574 —0.774 —0.925 —0.768 —0.956 —1.008 —0.852 —1.099 —1.592 —1.084  /q’d1  0 0.412 0.405 0.412 0.355 0.297 0.429 0.413 0.385 0.431 0.420 0.430 0.447 0.425 0.419 0.436  —3.275 —3.362 —3.375 —3.375 —3.384 —3.379 —3.373 —3.382 —3.387 —3.374 —3.382 —3.412 —3.421 —3.387 —3.396 —3.393  —3.275 —2.950 —2.969 —2.964 —3.028 —3.082 —2.944 —2.969 —3.002 —2.943 —2.962 —2.982 —2.974 —2.962 —2.976 —2.957  LLqf 0 0.325 0.306 0.311 0.247 0.193 0.331 0.306 0.273 0.332 0.313 0.293 0.301 0.313 0.299 0.318  Chapter 3. Fluids near planar surfaces  93  40 3 100 2 50  1  0  .IIiIIIIIIIIIIIIIIIIII.!  (b) 20  —  -  0 -  —20  0 —40  —2  —1  0  1  2 z  0  1  2  3  4z*  Figure 3.28: The self-consistent mean-field electrostatic potential for O.945M NaCl in contact with a metal slab (z = 32A). The long-dash-short-dash, dash-dot and long-dash curves are the jellium, dipole and ion contributions, respectively. The short-dash curve is the sum of the dipole and ion contributions and the dotted curve is the dipole contribution obtained for the pure solvent. The solid line in (a) is the total electrostatic potential. The parameters are as in Tables 3.2 and 3.3.  Chapter 4  Fluids between planar surfaces  4.1  Introduction  In chapter 3, the effect of surface properties on the interfacial fluid structure was ex amined. This represents one approach towards understanding fluids in inhomogeneous environments. Another approach, taken in experiments [97—101], computer simulations [102—108] and theories [109—118] over the last couple of decades, is to confine the fluid, for example between two surfaces, and study the separation-induced structure, phase behaviour and pressure of the fluid. In these confined systems, integral equation theories have been used primarily to study primitive model electrolytes (spherical ions in a dielec tric medium) [109—112]. Lennard-Jones fluids have also been considered near an isolated inert surface [119] using homogeneous integral equation theories and between soft planar walls using the inhomogeneous theories [113,114]. In addition, there has been an exten sive amount of work devoted to wetting (or drying) transitions at isolated surfaces [120]. To a lesser extent, capillary condensation [103—105,115,116] and evaporation [115] have also been investigated. In this chapter we consider fluids confined between planar walls and study the sep aration-induced fluid structure and fluid-mediated force between surfaces. General re ductions of the wall-wall and wall-solvent OZ equations are obtained for molecular fluids confined between planar surfaces. In addition, a general method is given for simulat ing metastable fluids beyond coexistence. Results are discussed in two parts. In the  94  Chapter 4. Fluids between planar surfaces  95  first, a survey of various integral equation methods for different model fluids between surfaces is given. These results are compared with computer simulations to verify the accuracy of the HNC and RHNC theories for different systems. In the second part, grand canonical Monte Carlo is used to simulate subcritical Lennard-Jones fluids between inert walls as a function of the wall-wall separation and to study capillary evaporation. The fluid, which is liquid in the bulk, evaporates at small separations. These results shed light on separation-induced phase changes. The possible relevance of these results to the experimentally measured attraction between hydrophobic surfaces is discussed. The remainder of this chapter is organized as follows. The reduction of the planar wall-wall OZ equation for molecular fluids is presented in Sec. 4.2.1. In Sec. 4.2.2, a method of calculating the wall-particle pair distribution function for a fluid confined between planar walls (using the homogeneous OZ equation) is presented. In Sec. 4.3, we elaborate on the grand canonical Monte Carlo method and apply it to a fluid of spherically symmetric particles between planar walls. Here, general methods of including long-range interactions for Lennard-Jones particles and of simulating metastable fluids are included. Results for hard sphere, Lennard-Jones and dipolar hard sphere fluids using the HNC and RHNC theories are presented in Sec. 4.4 and compared with simulation results. A general discussion of the applicability of the integral equation theories near surfaces and between surfaces is presented. Results for equilibrium and metastable simulations of Lennard-Jones fluids near coexistence are given in Sec. 4.5. These results are compared with a simple theory and experimental data.  Chapter 4. Fluids between planar surfaces  4.2  96  Theory  4.2.1  The interaction free energy between two walls separated by a molecular fluid  In this subsection, we consider correlations between planar parallel walls immersed in a molecular fluid. The relevant OZ equation is given by Eq. (2.2.7d) where the two walls, considered to be infinitely large solutes at infinite dilution, have coordinates denoted by 1 and 2 respectively. For smooth surfaces, the wall-wall pair functions vary only with 1 the separation 1 and we use a coordinate system in which z surface and 1  =  =  0 (see Fig. 3.1) at the left  0 at wall-wall contact. The reduction of the OZ relation follows from the  analysis of the wall-solvent OZ equation presented in Sec. 3.2. Using the identity c(32)  =  (23) [see Eq. (2.1.9)], the OZ equation is first written 7 c  in the equivalent form  h(l)  —  J  c(l) =  h(O,  )c 3 z (O, 7  3 (z  —  l)) d(3).  (4.2.1)  Expanding the wall-solvent pair correlation function as in Eq. (3.2.1) and the wall-solvent direct correlation functions as (23) 7 c  =  =  c(l  — Z3)  (_)nf  ()  —) c(l  -  Z3)  , R( ) 3  (4.2.2)  Chapter 4. Fluids between planar surfaces  97  where Eq. (A.8b) has been used, we get’ fOmm\ /Onn\  0mrn0nn  2 8ir  )  o  mj.. 7W  p00  i00  /  x  /  dx  J—oo  [00  0mm  dy J—oo  J—oo  /  o o Onn  ‘  ) dz 3 z 3  /  hOi;w.yIZ3) Co,,;w..y(l  —  ) dfz 3 *(ç . 3 -oi o, i3  (4.2.3)  x(—)JRm(c  Applying the orthogonality property of the Wigner spherical harmonics, Eq. (A.9), this reduces immediately to =  (2n+1)  ,  x/  ico  a0  Here we notice that  =  /  dx  J—oo  p00  dy  J—oo  h  (  °mn ji  ,  ui(1)  —  c,,  1—00  On  2 \1  o)i  iOnn I \ Onn ‘7 7 0 1 3 ;w ) COz,;wyt  —  ) dz 3 z 3  diverges with the surface area.  .  (4.2.4) Hence, using  Eq. (2.2.21), we introduce the potential of mean force per unit area as F(l) =  where  w(l) f dxfdy  =  U(l)  —  kBTB(l)  —  kBTN(l)  is the direct wall-wall pair potential per unit area,  (4.2.5)  is the wall-wall bridge  function per unit area and -  N(l) —  l) fdxfdy (_)n_v  =  (2n+1)  =  xJ  00  hOni  ‘  flN 0 Ffonn( ooo]2  L  Onn  /  op;wy(Z3) COv;wl  -00  For hard walls, the bridge function B(l)  0 when 1  —  z3)  3 dz  .  (4.2.6)  1 (i.e., when the walls are  separated by one or more solvent diameters) [3]. In order to solve Eq. (4.2.6), we need ‘The orientational dependence of c.(23) is incorrectly defined in Ref. [121]. This led to an incorrect sign on the lower index of homm O ;wy in Eq. (4.15) of that reference but has no effect on the results presented for dipolar fluids.  Chapter 4. Fluids between planar surfaces  (z) and 7 hg  c(z)  for positive and negative values of z. For (z) 7 h  For  Onn  we use  98  c =  h  —  =  Z  nO&,O,  0  Onn(+)  Onn  where  <  0v;wy 77  we have (4.2.7)  .  (z) +  Onn(—) iOu;w 7  (z).  Onn(+) 77 0v;wy  (z) is  obtained from the inverse Fourier transform of Eq. (3.2.18) for both positive and negative values of z and  Onn(—)  lOuw7 1  (z) is obtained from Eq. (3.2.14).  Since the OZ equation is formulated in the grand canonical ensemble, the walls, although infinitely large, are “immersed” in a molecular fluid with bulk surrounding them. When the walls are at infinite separation, no net force acts upon them; at small separations, a net force results from the pressure of bulk fluid pushing the walls into contact and from fluid confined between the walls pushing them apart. A positive value of F corresponds to a repulsion between surfaces. The excess potential of mean force per unit area is given by F(l)  =  F(l)  —  U(l)  (4.2.8)  ,  and represents the fluid-mediated contribution to F(l). Notice that, since h,(z) and c,(z) —*  0 as z  —*  oc, we have F(l) and F(l)  0 as 1  —*  —  oc, and neither the bulk  free energy nor the surface energies contribute to F(l) or F(l). The net pressure between the walls, mediated by the fluid, is given as P(l)  dF ( = —  (4.2.9)  .  Specifically, the net pressure is the normal component of the internal pressure at one wall due to fluid P(l) minus the bulk pressure Pb, P(l)  =  (l) 1 P  —  Pb  .  (4.2.10)  Chapter 4. Fluids between planar surfaces  4.2.2  99  The wall-solvent-wall correlation functions  The wall-wall OZ equation only describes correlations between surfaces. The three-body wall-solvent-wall correlation functions are also of interest and can be defined from the inhomogeneous OZ equation h(l,3)  —  c(l,3)  ) h(l,4) 4 ps(  =  Css(4,3)  d(4)  (4.2.11)  if one wall [we use the left wall in Eq. (4.2.11)] is treated as an external field and defines a space-fixed coordinate frame. In this approach, the inhomogeneous direct correlation function 6 c ( 4, 3) and the singlet density distribution ps( ) near an isolated wall are re 4 quired as input, and a new solution is obtained for all wall-wall separations 1. The homogeneous wall-particle OZ equation can be modified to obtain similar results with relatively little effort. A similar procedure has been applied by Lozada-Cassou [109] to study primitive model electrolyte solutions (spherical ions in a dielectric medium) between hard walls. In the present homogeneous approach, the isolated wall-particle potential  uW,(C3,  ) 3 z  is replaced with the wall-particle-wall potential, uw,3(3; 1)  =  u(O,  ) + u( 3 z , 0, 1 3  —  ) 3 z  ,  (4.2.12)  and the resulting equations are solved for a given separation of the walls 1. The relevant OZ equation is now written as hp(3; 1)  —  c(3; 1)  =  2  J  (4; 1) 43) d(4), 7 h  (4.2.13)  where the two walls are treated as a single confining particle and the coordinate system . The only modification in the reduction of Eq. (2.2.7b) 2 defined in Sec. 4.2.1 has been used Comparing Eqs (4.2.11) and (4.2.13), we see that h(3; 1) 2 h(l, 3), c(3; 1) c(1, 3) and the inhomogeneous singlet density p(3) psgw(13) in an exact theory. Hence, it is interesting to point out that substitution of the results from the homogeneous calculations into Eq. (4.2.11) provides a relatively simple means of calculating the inhomogeneous correlation functions.  Chapter 4. Fluids between planar surfaces  (see Sec. 3.2) arises because h, 3 hg(z; 1)  100  —1 inside both surfaces,  =  z < 0  =  and in place of the definitions for  and  or  z > 1,  (4.2.14)  we use  h(3; 1)  =  3 h,(3; 1) W(z  —  h(3; 1)  =  3 W(z  1,  —  1/2)  —  1/2)  (4.2.15)  ,  (4.2.16)  where W(z)  (1,  z<l/2,  l0,  z>l/2.  =  (4.2.17)  .  For z > 0, relations of the form z; 1)  =  2irpyf  S[(z + z)/r] c(r) dr 2 r  Zyq-f-Z  jr2Sn[(zw 7 +2rp  2r  x  —  z)/r] c(r) dr  (2n + 1) (n 1n 1 [ (2n 1 + 1) 0 0 0  ‘  nl  j°° e(k; 1)  k) dk  (4.2.18)  ,  are obtained. Before proceeding, we must consider closure relations for Eq. (4.2.18). For spherical particles, the best method is direct application of Eq. (2.2.17) for some approx imation of the bridge functions. For molecular fluids between two walls, the approach of Fries and Patey [see Eq. (3.2.24)] is no longer convenient since u(3; 1) and the assumption that c(3; 1)  —*  1 0 as z  —+  —*  oc as z 3  —÷  +oo  oc is no longer valid. In this case, we  follow an analysis similar to Caillol [33,34] and take the derivative of Eq. (2.2.17) with respect to the Euler angles of the fluid particle h(12) 2 J  =  In general, this gives w(12), 2 J  (4.2.19)  Chapter 4.  101  Fluids between planar surfaces  or =  h(12)  c(12) 2 J  w(12) 2 J  (4.2.20)  u(l2) + J 2 -J b(12), 2  —  kBT  where  2 J  represents a linear combination of partial derivative operators with respect to  . Choosing J 2 f 2 as the step operator (see Eq. (3.45) of Ref. [5]) 1 \ sin O  8  cot  = ie  =  +  X2 8  and applying  J  +  81 ij  (4.2.21)  to the wall-particle basis functions yields J± R, ) 2 ( 1  /n(n  =  + 1)  —  v(v + 1)  ) 2 R(i  (4.2.22)  .  Expanding the closure Eq. (4.2.19) in rotational invariants for the wall-particle pair functions and applying the step operator yields J±h(12)  =  =  0ThTh f  ) 2 (°)Vn(n+1)_v(v+1)h(z)R±i(  00  —/n + ( ± 2 1)—v n v 1) 2 fl  xf°’’  (Oni  fOfl2fl2  2n (On 2  01  Onini 1 ;(z) w(z) R , 0 xg 1 (2)  21 (fZ R ) 2  (4.2.23)  .  Using either Eq. (A.9) or Eq. (A.11) then gives relationships of the form +V(2n + 2 (_)fl+Thl+fl  Blum v;w(z) 0 h  =  —  1)  (n + n 2  kBT 2 fl  x  12 In In n n n\ 1 2 n 1 ± 2 1v+1 o  1)  /n(n + 1) ‘  —  —  (v + v 2 v(v  1000  ± 1)  Onn o 1  A similar relation can be obtained starting with Eq. (4.2.20). When 3 Oow  1)  z>O. (4.2.24) n =  0, the projection  is independent of orientation and Eq. (4.2.24) is not valid. For this special case,  we use the definition of the potential of mean force z)  = e_()hhJT  —  1  (4.2.25)  Chapter 4. Fluids between planar surfaces  for the choice  102  (0, 0, 0). This gives  h(z) = exp  _w(z)) -1-  h(z).  (4.2.26)  (-  Together Eqs (4.2.24) and (4.2.26) form a complete closure for the wall-particle OZ equation.  4.3  Grand canonical Monte Carlo method for Lennard-Jones particles  The grand canonical Monte Carlo method was introduced in Sec. 2.4. In this section, it is considered in detail for a fluid of Lennard-Jones particles confined between two planar walls. The computer simulations are performed in a central cell of dimensions L x L x 1 in which periodic boundary conditions are imposed along the x and y axes and two planar walls are located at z = 0 and z = 1. The solvent-solvent interactions are given by the Lennard-Jones potential [(u)12  —  (7)6]  u(r) = 4€  ,  (4.3.1)  where e is the depth of the potential well and o is a measure of the solvent diameter. The solvent interaction with the two walls is calculated assuming a pairwise additive potential [see Eq. (4.2.12)]. The wall-solvent interaction is modeled predominantly by the hard wall potential u(z) =  I  co,  z</2,  (0,  (4.3.2)  z>a/2.  However, we also consider the Lennard-Jones 10—4 potential u(z) = 2 Psur  2  g  I  \1O  + u/2)  ,  — (z + u/2)  (4.3.3)  where Psin is the surface density. Equation (4.3.3) represents the interaction energy of a Lennard-Jones particle in the fluid (at a separation z from the surface) with a uniform  Chapter 4. Fluids between planar surfaces  103  distribution of Lennard-Jones particles in a plane centered at z  =  —u/2 [122]. To get the  particular form of Eq. (4.3.3), the Lennard-Jones wall-solvent parameters  (  and u)  are equal to the solvent-solvent values (€ and u). The planar confinement of the fluid suggests a cylindrical truncation of u (r) instead of the usual spherical truncation used 5 for homogeneous bulk fluids. In this way, the potential  u(r)  I =  ss(v”5  ) 2 +z  S  ,  10, is defined where  <  St  (4.3.4)  u(s,z) = S>St,  L/2 is the cylindrical cut-off radius. The cylindrical truncation  st  scheme has the advantage that the mechanical properties of the fluid can be corrected for long-ranged structural effects in a straightforward manner. The average configurational energy of the central cell, u(s, z) +  U =  where zj  zj  —  z,  sj =  i  u(z; 1)) + U,  (4.3.5)  j.<i —  s and the angular brackets denote a grand canonical  average, is calculated using the minimum image convention. The long-range correction U is given by U  =  2Jv  1J2  ;l) 2 p)(ri,r  +z? 2 {uss(Vs? ) 2 z _u(si ) ] i  1 is restricted to the volume V where the integral over r  =  ,  (4.3.6)  1 of the central cell 2 L  2 (f r  ,r 1 ; 1) is the inhomogeneous solvent-solvent pair distribution 2 in unrestricted) and p)(r function between two planar surfaces separated by a distance 1. truncation cylinder is centered on particle 1 such that s  =  In Eq. (4.3.6), the  s 2 > .s 0 and only values of 1  contributes. For the inhomogeneous fluid, the pair distribution function is required in order to evaluate the correction U. If the truncation radius s is chosen sufficiently large, spatial  Chapter 4. Fluids between planar surfaces  104  correlations can be ignored such that ,r 1 p(r ; 1) 2  =  ; 1) p(’)(z 1 p(’)(z ; 1) g)(r 2 ,r 1 ; 1) 2 p(’)(zi; 1) 2 )(z 1), 1 p( ;  12 > s s  where p(’)(z; 1) is the singlet distribution function. Further, if  p(1)(z;  (4.3.7)  ,  1) is assumed to be  uniform then the assumption used by Scheon and coworkers [107] is obtained. A slightly better approximation is p’>(z;  1)  po(z; 1)  (4.3.8)  where, for example, po(z; 1) is the grand canonical Monte Carlo singlet distribution func tion when U  =  0. This is the approximation used in the present calculations. The  correction to the energy can then be evaluated either in integral form, , 1 j po(zi; 1)1(zi) dr  U  (4.3.9a)  or as an ensemble average, (4.3.9b) where 1(zi) =  f  ; 1) [u(s + z? 2 po(z ) 2  —  (s )] 8 u , 2 12 dr z . 2  (4.3.10)  Explicitly, using the asymptotic limit of the Lennard-Jones potential, u(r)  ‘-.‘  , 6 —4(u/r)  one obtains U  _2  (j  po(z2;  (  _  +  Z)2)  ) 2 dz  .  (4.3.11)  The integral 1(z) is easily evaluated as a function of z for a given approximation po(z; 1). Here, U is calculated from Eq. (4.3.11) and po(z; 1) is obtained from at least 10 million configurations on a grid of /z  =  particle insertions or deletions.  0.Olu. This tail correction is included when attempting  Chapter 4. Fluids between planar surfaces  105  The contact value theorem [3] is used to calculate the internal pressure and relates the internal pressure on one wall to the contact density of the fluid, =  lim  p(1)(z;  1)  =  lim  p(1)(z;  1)  .  (4.3.12)  In the simulations, P (l) was obtained from the average contact value at the two walls, 1 each contact value is linearly extrapolated using values of  p(1)(z;  1) within 0.lu of contact.  The grand canonical Monte Carlo calculations were carried out at constant  =  , as discussed in Sec. 2.3. The equilibrium state of the system is obtained after 3 kBT in zo  attempting a large number of particle operations (translations, insertions and deletions) at constant  res,  V, T and 1.  Grand canonical Monte Carlo has several advantages over canonical Monte Carlo cal culations. It is particularly useful for studying inhomogeneous systems and for studying phase behaviour. In this work, we are interested in the pressure between walls as a func tion of separation at constant chemical potential, and further, in this pressure near phase transitions. Conventional grand canonical Monte Carlo is suited for finding the equilib rium coexistence point of two phases, but it can be adapted to study the metastable fluid beyond the coexistence point, and to find the spinodal point itself. Consider a fluid close to the liquid vapour coexistence curve where a phase transition can be induced by small variations in  V or T. In this region, there exists two minima in the free  energy as a function of density at constant  res,  V, and T. The global minimum corre  sponds to the equilibrium state and a second minimum corresponds to the metastable state. In general, for metastability to exist, the two minima are separated by a free energy barrier. The disappearance of this barrier corresponds to the spinodal point. The metastable branch can be simulated by preventing the system from crossing the barrier. In this way, configurations near the local metastable minimum are selectively sampled.  Chapter 4. Fluids between planar surfaces  106  Since the two branches (the equilibrium state, and the metastable state) correspond to different densities, one can confine the fluid to the metastable branch by imposing a threshold volume-averaged density,  öth, 1  and by rejecting any attempted particle inser  tions or deletions that would cross this threshold. By monitoring fluctuations in the number of particles during the simulation, one can ascertain whether the average density corresponds to the metastable state, at the local minimum, or whether it corresponds to the threshold density. The spinodal point in which cavitation occurs is approximated by the minimum separation at which a local metastable minimum is observed; the spin odal point in which condensation occurs is approximated by the maximum separation in which a local metastable minimum is observed.  4.4  Results for pure fluids between planar walls  This section is concerned with the structure of simple fluids and dipolar hard sphere fluids confined between walls and the solvent-mediated wall-wall pressure. The pressure between hard walls immersed in a hard sphere fluid [calculated from Eq. shown in Fig. 4.1 as a function of the reduced wall-wall separation 1*  =  (4.2.9)] is  l/d > 1. Here,  results obtained from the HNC theory, RHNC theory and grand canonical Monte Carlo simulations [123] are compared for a reduced hard sphere density of p  =  0.681. The  RHNC calculations were carried out using the Verlet-Weis [27] and Henderson-Plischke [31] parameterization for the hard particle-particle and hard wall-particle bridge func tions, respectively. The direct hard wall-wall bridge functions are zero when 1* > 1 and are not required. The integral equation results generally show little dependence on the approximation used. At narrow separations, the pressure is oscillatory with a period of about one hard sphere diameter. The oscillations arise from particle correlations between the solvent layers at the surfaces.  Chapter 4. Fluids between planar surfaces  10  III  II  1111111  I  I  107  II  0.5  (a)  P(l)d/kBT  5  0.0 —0.5  0 I  1  I  I  I  2  I  I  I  3  I  I  41*5  —1.0 2  3  4  5  6  71*  Figure 4.1: Net pressure between hard walls in a hard sphere fluid at a density of (a) The solid curve is the HNC result, the almost coincident dotted curve = 0.681. is the HNCP result and the dashed curve is the result obtained using the Verlet-Weis [27] and Henderson-Plischke [31] bridge functions. (b) Comparison on an expanded scale with the simulation data of Karlström [123] (squares); the curves are the same as in (a). Figure 4.lb compares the integral equation results with grand canonical Monte Carlo simulation results of Karlström [123] on an expanded scale (note that the bulk pressure of  3 has been subtracted from the internal pressure data of Karlström in order 3.695kBTd that the net pressure approaches zero as z  —*  oc). Although noise in the simulation  data makes it difficult to assess the effects of including bridge functions, the RHNC results generally show slightly better agreement with simulations. This system of hard spheres between hard walls has also been solved using the inhomogeneous OZ equation together with the Percus-Yevick (PY) approximation [113]. The homogeneous RHNC results are about as accurate as the inhomogeneous results. The largest disagreement between the RHNC and Monte Carlo results occurs near the troughs in the pressure. This disagreement is also seen in the inhomogeneous PY results [113]. The wall-solvent distribution functions  ; 2 gws(z  1) calculated from Eq.  (4.2.13) are  Chapter 4. Fluids between planar surfaces  108  shown in Fig. 4.2 for several wall-wall separations. The results for 1*  =  10 are com  pared with grand canonical Monte Carlo results [123]. Again the agreement with simu lation results is good except very close to the wall where the HNC theory overestimates (ds/2; 1). At 1* 8 gw  =  3.1 and 1*  =  10, the distribution function can also be compared  with the inhomogeneous PY results [113] (not shown). In both cases, the HNC results are about as accurate as the inhomogeneous PY theory except near contact. In an exact theory, the contact theorem [see Eq. (3.4.23)] tells us that pg(d/2; 1) is equal to the in ternal pressure when the walls are separated by a distance land P(l) In Fig. 4.3, we plot pg(d/2; 1)  —  =  Ps calculated from the HNC and RHNC theories and  compare the results with the net pressure results from Eq. (4.2.9). These results show that while the contact density in the HNC theory is too high, variations in the contact density due to confinement of the fluid are proportional to the net pressure and the con tact density for hard spheres near a hard wall can also be used to obtain the pressure between the surfaces. The net pressure between hard walls mediated by a dipolar hard sphere fluid was also calculated using Eq. (4.2.9) for the RHNC theory. In this case, the RHNC calculation was carried out using the Verlet-Weis parameterization [27] and HNCP approximation [26] for the solvent-solvent and wall-solvent bridge functions, respectively. Figure 4.4 shows that the results are qualitatively similar to the hard sphere fluid (Fig. 4.1) and that the magnitude of the oscillations are reduced as the dipole moment is increased. Using the analysis presented in Sec. 3.4.2, the long-range asymptotic behaviour of N(12) between planar walls immersed in a dipolar fluid can be determined. For a dipolar fluid near an inert wall, h (12) decays asymptotically as [see Eq. (3.4.17)] 5 [53_2]2  (12) 5 h where  —*  =  32p  (x + X  (cosO)) 2 P 5  z  —+  (4.4.1)  ‘r is the isothermal compressibility. In the HNC/RHNC approximation,  Chapter 4. Fluids between planar surfaces  109  6  6  4  4  2  2 0  0 0  1234  1  0  z*5  2  3  4z*5  10  10  10  10  8  8  8  6  6  6  6  4  4  4  4  2  2  2  0  0  0  0  1  2z*3  I  :‘  0  Ljd  :  :  Figure 4.2: The wall-hard sphere distribution function arations 1. The parameters are as in Fig. 4.1.  gws(z;  :‘  I  8  1.5  0  2 lz  :  2 :  0  lz  Ozl  1) for several wall-wall sep  the wall-solvent direct correlation function has the limit c(21)  h(21)  —  6 far from an inert wall. In addition, since c(21) and, hence, decays as 1/z z  —*  (4.4.2)  (21)/kBT, 8 u —*  as  cc, we make use of the approximation c(21)  1  Z  to obtain F(l)  kBT F(_l)2 3YEs  —L  2  j  I  .11  —dz  (4.4.3)  Chapter 4. Fluids between planar surfaces  110  III  0.5  0.5  0  (b)Z  —  I  I I  -  z>j7  —0.5  —0.5  1(11111  -1 _i  0  1111111111  -  I  2  3  4  5  6z*  2  I  i  i  i  I  4  3  i  i  i  I  I  5  I  I  i  6z’  Figure 4.3: The (a) HNC and (b) RHNC contact density minus the bulk density as a function of wall-wall separation 1. The crosses denotes values of pg(d/2; 1) Ps• The curves are the net HNC and RHNC pressures calculated from Eq. (4.2.9). —  2 32l  [  2 i)2] ,  1  -  oo,  (4.4.4)  where y = 47rpS/9kBT. The asymptotic limit of the net pressure, plotted in Fig. 4.4b, is given by P(l)  =  dF(l)  (4.4.5)  The first two projections of the wall-solvent distribution functions  gws(z2;  1) are shown  in Figs 4.5 and 4.6 for several wall-wall separations. Simulation results are available for dipolar hard spheres near isolated surfaces [54] and show good agreement with RHNC theories. The distribution functions for a fluid of Lennard-Jones particles confined between two Lennard-Jones 10—4 walls at separations of 1* = l/o = 2 and l = 3 are shown in Fig. 4.7. Here results from the HNC theory and molecular dynamics calculations [108] are compared for a fluid with a reduced density of p  = 0.5925 and a reduced  Chapter 4. Fluids between planar surfaces  111  0.5  111111111111111111111111  0.4 -  0.2  P(l)d/kBT x  io  (b)  -  0.0  0.0  -  1  -  —  It -t  —0.2  —  — .—.  St.  -  —0.5  LiI\J-  —1.0  ‘IIIIIIIIIIIIIIIIIIIIIIIII  —  —0.4 —0.6 2  345  61*  8  9  10  11  12 l  Figure 4.4: HNCP pressure between hard walls in dipolar hard sphere fluids at a density of p’ = 0.7. In (a), the solid, dashed and dotted curves are for reduced dipoles moments of = respectively. In (b), the HNC asymptotes and u = 0 (hard spheres), calculated from Eq. (4.4.5) (smooth curves) are compared with numerical HNCP results on an expanded scale. temperature of T*  =  kBT/e  =  3 were carried 1.21. The molecular dynamics simulations  out using a spherically truncated and shifted Lennard-Jones potential (u’(r)  —  nLJ(3.5a),  r <3.5a (4.4.6)  u(r) =  0,  r>3.5u.  In the HNC theory for Lennard-Jones particles, the bulk direct correlation function  c(r) was obtained employing the perturbation theory of Weeks, Chandler and Anderson [124,125]. In order to compare with the simulation results, the bulk direct correlation function was calculated using both the full Lennard-Jones, Eq. (2.4.10), and the trun cated Lennard-Jones, Eq. (4.4.6), potentials, respectively. Using Eq. (4.4.6), the HNC results show good agreement with simulations. Using the full Lennard-Jones potential, the results show a small but significant drop in the distribution function near contact. Using a cylindrical cut-off and correcting the interaction energy as described in Sec. 4.3, The state parameters in the molecular dynamics simulation [108] were chosen such that the confined 3 fluid is in equilibrium with the bulk having an average density 0.5925.  Chapter 4. Fluids between planar surfaces  112  6  6  4  4  2  2  o  0 4z*5  0  0  1  2  6  6  6  4  4  4  2  2  2  o  0  0  0  1  2z*3  0  lz*2  3  4z*5  0  lz*2  Figure 4.5: RHNC results for the angle-averaged wall-particle pair distribution function g°(z; 1) for a dipolar hard sphere fluid confined between hard walls. Results are shown for several wall-wall separations 1. The parameters are as in Fig. 4.4. grand canonical Monte Carlo can be used to simulate the full Lennard-Jones potential. Figure 4.8 shows the distribution function for a Lennard-Jones fluid confined be tween hard walls calculated from the HNC theory, RHNC theory with the HNCP hard sphere bridge function and grand canonical Monte Carlo simulations. The grand canon ical Monte Carlo simulation method is described in Sec. 4.3 including a correction for long-range Lennard-Jones interactions. Results for a wall-wall separation of 1* considered for three dense bulk state points at reduced temperatures of T*  =  =  5 are  1.0, 1.35  and 2.0. For the low temperature state point, the distribution function for a wall-wall  Chapter 4. Fluids between planar surfaces  113  :1 I I I  0  0  —1  —1  —2  —2  —3  —3 —l  lilt  II  3  2  1  0  4z*5 I  I I I I  I I I I  I  1*  I 1  f—I  0 —l  Ill  I  I  I  I  I  l*=2.5  =  i  2  —1  —1  —1  —2  —2  —2  0  I  11111  I  2z*3  1  —3  -f-I  0  lilt  III  I  i-I  2 lz  I  i  I  —I I I I  0  -f-I I I  I  3  0  —3  I I I  I 1 I  5  II  —  I  I-I  4z*5 I  I I I I—  0  —3  ÷1111111÷  o  2 i  Figure 4.6: RHNC results for the projection g 2 of the wall-particle pair distribution between hard walls. Results are shown confined sphere fluid function for a dipolar hard for several wall-wall separations 1. The parameters are as in Fig. 4.4. separation of 1*  =  10 is also included. Details of the computer simulations are shown in  Table 4.1. In general, the HNC and RHNC results show relatively little dependence on the state parameters, whereas the grand canonical Monte Carlo results clearly show that the Lennard-Jones fluid is in the liquid phase for the first two state points with 1* 4.8a and 4.8b) and in the gas phase for the third state with 1*  =  5 (Figs  5 (Fig. 4.8c). In addition,  results for the third state point (Figs 4.8c and 4.8d) indicates that the phase transition is a function, not only of the state parameters, but also of the wall-wall separation. The  Chapter 4. Fluids between planar surfaces  114  6  6  4  4  2  2  o  0.0  0.5  1.0  1.5z*2.0  0 0  1  2  z*3  Figure 4.7: Wall-solvent distribution functions for a Lennard-Jones fluid confined between 1.21, p = 0.5925 and Ps Lennard-Jones 10—4 walls, Eq. (4.3.3) (T* 2 = 0.8). The wall-wall separations are (a) 1* = 2 and (b) 1* = 3u. The solid and dotted curves are the HNC results using the Lennard-Jones, Eq. (2.4.10), and truncated Lennard-Jones, Eq. (4.4.6), potentials, respectively. The HNC results were obtained using the bulk direct correlation function of perturbation theory [124,125]. The circles are results from molecular dynamics simulations [108]. grand canonical Monte Carlo results are characteristic of a drying transition as a fluid with attractive interactions near an isolated “solvophobic” surface approaches the bulk coexistence curves [102,119]. The phase diagram for the Lennard-Jones fluid [126] shows that the density and temperature considered in Fig. 4.8c are for a state point close to the coexistence curve. The failure of the HNC theory to account for dewetting at an isolated solvophobic surface is well-known [119,8] and is a result of the approximation used in the closure. It is not surprising that the RHNC theory also fails in this case since the hard sphere bridge functions are used as an approximation of the state dependent hard wall—Lennard-Jones bridge function. At high temperatures, the strength of the Lennard-Jones attraction is  Chapter 4. Fluids between planar surfaces  4  _i11Fj11IIJIIIJIIIijiiii_  115  3  I  1 (a)  g(z;5u)  J  II  I  J  Illij  liii  (b)  g(z;5u)  3 2 2 1 1  0  I  I  I  o  3  I  I  I  1  2  3  liii  I  4z*5  (c)I  I I  —  —  2  I I I I I  —  —  —  I I I I I  —  1  f\  — -  /  ‘ /  /\ \/  \/ ‘  f‘  \  I  -  I  2  1  0 3  1111111111111111111  g(z;5a)  -  I  0  I  I  I  I  gws(z;  I  J  I  I  4z*5  3 I  I  I  I  I  lOu)  I  I  I  I  I  (d)  —  —  2  —  —  —  —  —  I  1  -  -  0  IIIIIIIIIIIIIIIIII11_II  o  1  2  3  4z*5  0  I I I  0  I  1  I I I  I  2  I  I  3  I 4z*5  Figure 4.8: Wall-solvent distribution functions for a Lennard-Jones fluid confined between hard walls for three state points and two wall-wall separations. The solid curves are results from grand canonical Monte Carlo simulations. The dashed and dotted curves are the HNC and HNCP results. In (a) T* = 2.0 and p = 0.710 + 0.003; (b) T* = 1.35 and p = 0.70 ± 0.01; and in (c) and (d) T* = 1.0 and p = 0.75 + 0.01.  Chapter 4. Fluids between planar surfaces  116  Table 4.1: Details of the grand canonical Monte Carlo simulations for a Lennard-Jones fluid confined between hard walls. The last column labeled “operations”, denotes the total number of particle moves, insertions and deletions included in the simulation. (The grand canonical Monte Carlo simulations for the bulk were performed in a central cell 9u x 9u x 9a.)  i/u 5 10 5 5  L/u 15 7 7 7  St/J  3.5 3.5 3.5 3.5  T* 1.0 1.0 1.35 2.0  preS/kT  -3.287 -3.287 -1.941 1.525  0.751 + 0.007 0.751 + 0.007 0.700+0.006 0.710+0.003  Operations 6 x10 10 10 10 10  attenuated and the bridge functions, defined as integrals over Mayer-f functions, f(12)  =  e_  12)/kT —  (4.4.7)  1,  approach the hard wall results and better results are obtained. For narrow wall-wall separations, even at relatively high temperatures, the HNC and RHNC approximation are less reliable. For the Lennard-Jones fluid near an attractive 10—4 wall, it is worthwhile to note that the state point considered in Fig. 4.7 (T* = 1.21 and p  =  —.5925) is also  very close to the coexistence curve yet the integral equation theories show relatively good agreement with simulations.  4.5  Cavitation of a Lennard-Jones fluid between hard walls  As pointed out in the last section, the HNC and RHNC theories, using hard wall-hard sphere bridge functions, are inadequate to describe Lennard-Jones fluids between hard walls for state points near the coexistence curve. In this section, grand canonical Monte Carlo is used to investigate the phase behaviour of this system for two bulk state points near the coexistence curve. In addition to studying the equilibrium states of the fluid,  Chapter 4. Fluids between planar surfaces  117  we also consider the metastable liquid and gaseous branches. Tables 4.2 and 4.3 contain details of the simulations. The threshold density th, defined in Sec. 4.3, is included for the metastable states. Table 4.2 is for the confined Lennard-Jones fluid with T* = 1.0 and ,t1e8/kBT eX/kBT  =  =  —3.29, which corresponds to a bulk having excess chemical potential,  —3.00; density,  3 = Ps  0.75 ± 0.01; and pressure, Pbu /kBT 3  Table 4.3 is for the second state point, T* = 0.9 and u/kBT to ex/kBT  =  =  3 a 8 —3.60, p  =  0.79 ± 0.01, PbJ /kBT 3  =  =  =  0.44 + 0.05.  —3.84, which corresponds  0.4 ± 0.1.  Again, the grand  canonical Monte Carlo simulations for the bulk liquids were performed in a central cubic cell using periodic boundary conditions. The Lennard-Jones pair potential was truncated spherically at half the cell length and the potential energy was corrected assuming a uniform distribution of particles beyond the truncation sphere. Figure 4.9a shows the average density between the walls as a function of wall separa tion for the first state point. The equilibrium phase points, denoted by crosses, show that a high density fluid is stable at large separations, and that a distinct low density fluid is stable at small separations. The two fluids coexist at a separation of  ‘-.  6.Ou. It was pos  sible to follow the metastable liquid branch, denoted by open squares, into separations as small as 4.Ou. Hence this separation approximates the spinodal point at which cavitation occurs. Although decreasing, the average density of the metastable liquid remains clearly greater than the density of the gas. Moreover, this density is larger than the threshold density (Table 4.2), which shows that the system is in a true local free energy minimum. As illustrated in Fig. 4.9a, it was also possible to simulate the metastable gas out to separations of  15o- and 20u. Although there is a significant L dependence in the larger  separation result, it is clear that the spinodal point for condensation is at very large separations Figure 4.9b shows the average density for the second, higher density, lower temper ature, state point. Here, the gas phase is also stable at narrow separations, 1  5  6.5a.  118  Chapter 4. Fluids between planar surfaces  I  0.8  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  (a)  -  0.6  X  X  X  X  0.4 0.2 0  I  I  :  I  I  I  I  I  I  I  I  I  I  I  I  15 I  I  I  I  0  0  I  20 I  I  *  I  I  I  II  11111 I  .  I  I  C  —  E  0.04-  -  0.02—  -  —  -  11111111111  0  0  I  I  (b)  0.6  0.2  I  10  —  0.4  I  I  5 0.8  C  L:J  XXX  x  20  40  60 0  XX I  I  I  I  5  I  I  10  I  15  I  I  I  I  20  i*  , between 2 Figure 4.9: The grand canonical Monte Carlo average density = (N)/1L Results for the first hard walls for a Lennard-Jones fluid in equilibrium with the bulk. state point (&/kBT = —3.29 and T* = 1.0) are shown in (a) and for the second state point (‘/kBT = —3.84 and T* = 0.9) in (b). Crosses represent points that are thermodynamically stable and open squares represent those that are metastable. Open circles are used when the equilibrium and metastable fluids could not be distinguished. Except where illustrated, the standard deviations in are smaller than the symbols. The inset shows an expanded view of the gas branch.  Chapter 4. Fluids between planar surfaces  119  Table 4.2: Details of the grand canonical Monte Carlo simulations for a Lennard-Jones fluid with T* = 1.0 and p7/kBT = —3.29 confined between planar hard walls. (The grand canonical Monte Carlo simulation for the bulk, performed in a central cell 9u x 9o x 9, has (pex)/kBT = —3.00, PsU /kBT 0.44 ± 0.05.) 3 3 = 0.75 ± 0.01 and Fbo  i/u 3.0 5.0 5.0 5.2 5.5 5.5 5.8 6.0 6.0 6.5 7.0 8.0 10.0 10.0 15.0 20.0  L/u 7 10 10 10 10 10 10 10 10 10 10 10 10 15 15 15  St/U  3.5 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0  3 iiu 0.0493 + 0.0002 0.0621 + 0.0005 0.47 ±0.01 0.50 +0.01 0.065 +0.001 0.551 ±0.007 0.572 ±0.006 0.068 ± 0.001 0.588 ±0.006 0.607 +0.006 0.624 + 0.007 0.643 ±0.007 0.668 ±0.006 0.084 +0.002 0.095 ±0.001 0.0988±0.0002  PIflt(l)cx / 3 kBT 0.048 + 0.002 0.051 ± 0.002 0.29 ±0.01 0.31 ±0.03 0.052+0.002 0.36 ±0.02 0.38 +0.02 0.052 + 0.002 0.38 ±0.02 0.40 ±0.03 0.405 + 0.02 0.41 ±0.03 0.42 ±0.03 0.052±0.002 0.053+0.03 0.041±0.03  0.40 0.40 0.45 0.48  0.10 0.10 0.10  Operations 6 x10 10 40 40 20 30 20 20 30 30 30 30 30 60 30 30 30  At large separations and in the bulk, the liquid phase is stable and the gas phase is metastable. For separations 6.5u  5  1  51u, the grand canonical Monte Carlo method  could not distinguish the equilibrium and metastable states after at least 20 million con figurations, and thus the coexistence separation was not determined. Compared to the state point shown in Fig 4.9a, the spinodal separation for cavitation is slightly smaller at 4.9u and although the coexistence separation is not known, it is at least slightly larger at 1  6.5cr. The density profiles for the equilibrium and metastable liquid branches are shown in  Fig. 4.10 using the midplane between the two walls as the origin. The sharp decrease in  Chapter 4. Fluids between planar surfaces  120  Table 4.3: Details of the grand canonical Monte Carlo simulations for the second state point (T* = 0.9 and res/kBT = —3.84). (The grand canonical Monte Carlo simulation for the bulk, performed in a central cell 7o’ x 7o x 7o-, has ex = 3.60, Ps°’ 3 = 0.79±0.01 0.1.) = 0.4 and Pb0’ ± /kBT 3  l/ 4.9 5.0 5.0 5.2 5.5 5.8 6.0 6.0 6.5 7.0 8.0 10.0 21.0 51.0  L/ 10 10 10 10 10 10 10 10 10 10 10 10 16 10  St/0’  5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 5.0 8.0 5.0  0.496 ±0.024 0.0287 ± 0.001 0.562 ±0.01 0.591 ±0.01 0.619 +0.01 0.634 +0.01 0.0294 ± 0.0001 0.644 ±0.007 0.662 ±0.006 0.674 ±0.005 0.695 ± 0.007 0.718 +0.006 0.0320 ± 0.0001 0.0323±0.0002  /kBT 3 Pt(l)cr 0.25 ±0.02 0.023 ± 0.002 0.31 ±0.02 0.34 ±0.02 0.36 ±0.02 0.37 ±0.02 0.023 ± 0.002 0.38 ±0.02 0.39 ±0.02 0.39 ±0.02 0.41 + 0.02 0.41 ±0.02 0.023 ± 0.002 0.023+0.002  3 Pth0’  0.40 0.40 0.40 0.42 0.42 0.45  Operations 6 x10 30 30 20 20 20 20 30 20 20 20 20 20 30 30  Chapter 4. Fluids between planar surfaces  1  0.8  I I I I (z;l) 5 pg I  I  I  I  I  I  I  I  I (a) I  0.6  121  1 0.8  0.6  0.4 0.2  0.2  0  0 —4 —2  0  2  LO—  —  0.4  678  l*=5  I 111111111 II pg(z;l) (b)  1111  4z*  78  liii I  —4 —2  I  0  I  liii 4z* 2  Figure 4.10: Grand canonical Monte Carlo density distributions for the equilibrium and metastable liquid branches for wall-wall separations l = 5, 6, 7, 8 and 10. The first (a) and second (b) state points are defined as in Fig. 4.9. Here, the midplane between the walls is at z = 0. the density adjacent to the walls is characteristic of a liquid close to drying. The density profiles show the oscillations that are characteristic of a liquid, and their amplitude is slightly larger for the higher bulk density state point (Fig. 4.lOb). The net pressure between the walls is displayed in Fig. 4.11. By definition, the net pressure must approach zero at large separations. In the simulations, the net pressure at large separations is equal to the difference in two nearly equal quantities (the total pressure and the bulk pressure) and is strongly influenced by errors in the bulk value. These can arise from statistical uncertainty, the finite size and shape of the simulation cell, and also from the different potential truncations used (spherical in the bulk and cylindrical between the walls). In order to ensure that the net pressures shown in Fig. 4.11 correctly approach zero, P(l) was calculated using the total pressure at 1  =  10 on  the liquid branch rather that the bulk pressure obtained by direct simulation. This only shifts the curve by a few percent. In general, Fig. 4.11 shows that the force between the walls is attractive, and on  Chapter 4. Fluids between planar surfaces  122  0  0  —0.2  —0.2  —0.4  —0.4 2  4  6  8  101*  2  4  6  8  101*  Figure 4.11: The net pressure between hard walls. The symbols are the same as in Fig. 4.9. The solid curve is the conventional van der Waals attraction, Eq. (4.5.3), and the dashed curve is the fitted mean-field result, Eq. (4.5.2). The fitting parameters for the first state point (a) are A0/kBT = 0.15, cw 2 = 0.02, and l = 3.Ou. For the second state point (b), AJ/kBT = 0.13, a 2 = 0.02, and 10 = 3.0g. For the gas branch, the standard deviations in the net pressure are smaller than the symbols. any one branch the magnitude decays monotonically with increasing separation. It is possible that the pressure on the liquid branch is oscillatory due to molecular packing between the walls [103,105], but the magnitude of any such oscillations appears small. Oscillations are not expected on the gas branch. At the coexistence separation, there is a discontinuous jump in the pressure as the liquid evaporates to the gas, in qualitative agreement with earlier results [104,116]. At small separations, the force is dominated by a large attraction. Along the gas branch, the total pressure in the slit is very small, and hence, the net pressure is essentially the negative of the bulk pressure. This represents a large attractive force, only weakly dependent upon the separation between the walls. It is interesting to compare the simulation data with a simple mean-field theory due to Attard [127]. In this theory, the fluid at the midplane between the surfaces is assigned an effective chemical potential l(l) equal to the chemical potential of the bulk liquid minus  Chapter 4. Fluids between planar surfaces  123  the excess mean-field contribution from fluid displaced by the presence of the walls. (The fluid between the walls and the walls are each assumed to be as polarizable as the bulk.) At some separation 1, the effective chemical potential will equal the coexistence chemical potential. At a smaller separation, it will equal the bulk spinodal value which defines the spinodal separation  10.  Close to the spinodal separation, the potential of mean force per  unit area, ‘(l), and pressure, P(l), between walls decay as w(l)  “  i1e__b0)1,  P(l)  1  .‘  1  10, ‘  lo,  (4.5.1)  (4.5.2)  where A and a are positive parameters of the theory. Strictly speaking, the pressure between the walls must always be finite. The simulation data on the liquid branch in Fig. 4.11 has been fitted to Eq. (4.5.2). These fits are also plotted in Fig. 4.11 as dashed curves. It can be seen that the simple theory gives a rather good representation of the data for the liquid branch. In both cases the values of lo are somewhat smaller than those estimated directly from the grand canonical Monte Carlo simulations. This is probably due to the simplicity of the meanfield theory which completely neglects the structure of the fluid adjacent to the walls. Asymptotically, the net pressure between hard walls mediated by a Lennard-Jones fluid, is given by [121] =  P(l)  6 2irpu  (4.5.3)  This is the pressure that arises from the van der Waals forces and is included in Fig. 4.11 as solid curves. It may be seen that the van der Waals contribution greatly underestimates the attractions for the metastable liquid branch. This is unambiguous at separations near the spinodal point and is expected to hold at larger separations. Figure 4.12 shows experimental measurements of the potential of mean force per unit area between two hydrophobic surfaces in water [128]. The measurements were  Chapter 4. Fluids between planar surfaces  124  obtained at standard temperature and pressure. In Fig. 4.12, the solid curve represents the conventional van der Waals attraction between crossed cylinders of radius R [129], = —  IIR 61  (4.5.4)  for a Hamaker constant of H = 2 x 10 J. This value is an upper limit for a hydrocarbon20 water-hydrocarbon or a mica-water-mica system. We see that the measured attraction is significantly larger than the van der Waals force and increases sharply as the surfaces are brought into contact. The dashed line represents the mean field result, with parameters A = 1.2 x 10’ J/m, 2  ‘0  = l2nm, and a = 4 x 2 nm The measured force is very 4 10 .  well fitted by this choice of parameters. At a narrow separation (1  10  ‘-.  l2nm), the  attraction between surfaces exceeds the counter force of experimental apparatus and the surfaces jump into contact. When the surfaces are subsequently separated, cavitation can persist in some cases for separations on the order of micrometers [99], suggesting that water is metastable between these surfaces at relatively large separations. Currently, the origin of this long-ranged attraction between hydrophobic surfaces is unknown. However, comparison of the experiments and simulations with the mean-field theory, suggest that the measured attraction might result from a separation-induced phase transition. The simulation results show that cavitation will occur for water at subcritical state points near the coexistence curve and between hydrophobic surfaces at sufficiently narrow separation. In the simulation results of a Lennard-Jones fluid between hard walls at a particular bulk state point showed a cavitation separation of about 4 diameters, whereas water at standard temperature and pressure between hydrophobic walls apparently cay itates at separations of about 40 water diameters.  Chapter 4. Fluids between planar surfaces  125  111111111  I  lii  w(l) (mJ/m ) 2  0  -  -  -  ctJ  0.05:  d  -  -  —0.1 —0.15  -  -  I  20  I  I  I  I  I  I  40  Figure 4.12: The interaction free energy surfaces in water [128]. The solid curve with a Hamaker constant of 2 x 10 J, 20 result, Eq. (4.5.1) (A = 1.2 x 1 Jm 2 10’ ,  I  60  I  I  I  80  I  I  I  I  100  I  I  120 1 (nm)  per unit area measured between hydrophobic is the conventional van der Waals attraction, and the dashed curve is the fitted mean-field nm 4 10 ) . o = l2nm, and = 4 x 2 t  Chapter 5  Summary and Conclusion  This thesis has been concerned with nonspherical particles near planar surfaces and with both spherical and nonspherical particles between planar parallel surfaces. General theoretical methods for obtaining the wall-particle pair distribution function for mul ticomponent fluids near isolated inert, dielectric and metallic surfaces were presented. These were used to study simple dipolar hard sphere fluids and simple model electrolyte solutions near such surfaces. General methods for finding the fluid-mediated net pres sure between planar walls and for calculating the wall-particle correlation functions for nonspherical particles between surfaces have also been derived. Fluids of hard spheres, dipolar hard spheres and Lennard-Jones particles were considered between inert and Lennard-Jones surfaces. The theories were used to investigate the fluid-mediated net pressure between surfaces and the effects of confinement on the fluid structure. Com parison of HNC and RHNC results with grand canonical Monte Carlo calculations and molecular dynamics simulations allowed us to establish the validity of these integral equation theories for certain systems. Finally, a grand canonical method for studying metastable fluids was given and applied to a Lennard-Jones fluid between hard inert walls. These later results allow us to suggest an explanation for the experimentally measured attraction between hydrophobic surfaces. General reductions of the OZ equation and the HNC/RHNC closure relation were given for multicomponent fluids of nonspherical particles near an isolated wall. The equations are presented in a formalism similar to that previously derived for bulk fluids 126  Chapter 5. Summary and Conclusion  127  [6,10—13,17,18], and provide a relatively efficient method of calculating the structure of molecular liquids near planar walls.  In addition, fully self-consistent theories are  derived for multipolar particles near dielectric and metallic surfaces. In these theories, the many-body interfacial potentials are reduced to an effective mean-field pair interaction appropriate for use in integral equation theories. For a pure dipolar solvent near an inert surface, analytical expressions for the asymp totic behaviour of the wall-solvent pair correlation function at large wall-solvent sepa rations were obtained. Within the RNC approximation, the resulting equations were compared with the classical continuum (image-dipole) expression and with an exact lim iting result obtained by Badiali [62]. It is shown that all three limiting expressions agree only to linear order in density. The HNC and continuum theories give different second order terms and both are inexact. Numerical solutions of the RHNC equations were obtained for a dense system and these indicate that solvent structural effects cease to be important and the z 3 asymptotic limiting behaviour is reached at about 11-13 solvent diameters. Comparison with the exact asymptotic result shows that the RHNC approx 3 in the limiting imation overestimates the magnitude of the coefficient multiplying z expression. At a semi-infinite dielectric slab, the interaction energy of a fluid particle with the polarizable material is approximated by a uniform mean-field pair potential. Formal expressions have been obtained for general multipolar models. For the metal-liquid in terface, a metal slab is modeled at the jellium level and the electron density distribution across the interface is obtained from density functional theory. The liquid and metal in teract electrostatically and the resulting electron density distribution and metal-particle distribution functions obtained are completely self-consistent. This self-consistency is a principle feature of the present theory for the liquid-metal interface and is a significant  Chapter 5. Summary and Conclusion  128  improvement upon earlier treatments based upon the mean spherical closure approxima tion where the liquid only sees the metal as an inert wall. Explicit RHNC results were given for a simple dipolar hard sphere fluid near inert, dielectric and metallic surfaces, and for simple model electrolyte solutions near inert and metallic surfaces. Symmetry requirements do not allow the simple dipolar solvent to be polarized by an inert surface. This remains true if the dissolved anions and cations are of the same size. If the ions are of unequal size a net solvent polarization arises from the asymmetric density distributions of anions and cations across the interface, but this effect is not very large for simple ions. For inert walls, the ion distributions are governed principally by the solvation energy and are directly related to the ion diameters. That is, ions which are small relative to the solvent have a low contact density and those which are large have a high contact density. This holds independent of the sign of the ion charge and the particular coion involved does not have a large influence on the distribution. The symmetry conditions at inert surfaces also apply at dielectric surfaces. For the pure dipolar case, it was shown that images of other solvent particles make a large contri bution to the interaction between a particular particle and the wall. If these contributions are ignored and one includes only the simple self-image potential then rather dramatic structural changes occur as the dielectric constant of the wall is increased. However, the image contribution from all remaining fluid particles acts to largely cancel the self-image term and, as a result, the qualitative solvent structure near the surface is not strongly influenced by the dielectric constant of the wall. Indeed, for dipolar hard sphere solvents, polarization effects near the surface act to enhance the structure already present for inert walls. At a metal surface the solvent and ion distributions differ fundamentally from the inert and dielectric wall cases. The metal induces a very strong short-range electric field and a net polarization of the solvent distribution results. This was most dramatic in the  Chapter 5. Summary and Conclusion  129  layer contacting the metal where the dipoles tend to align parallel to the surface normal (i.e., positive ends towards the fluid). This is in agreement with earlier calculations for dipolar monolayers [67,68] and contrasts with the inert wall case where the dipoles prefer to lie perpendicular to the surface normal [54,50]. The polarization of the solvent at the surface plays an important role in determin ing the ion distribution. This gives rise to a “layered” structure with anions located at  3 + d_ /2 (measured from the surface) and cations nearer the surface. These qualita d tive structural features were found for all solutions considered and were not particularly sensitive to ion size or salt concentration. It is worth noting that the direct ion-jellium interaction is attractive for anions and repulsive for cations and hence would imply an “in verse” order with more anions than cations near the surface. Thus the observed ion-wall distributions are to a large extent determined by the ion-solvent forces and are greatly influenced by the jellium-induced solvent order. These calculations emphasize that the solvent response to the metal potential is of crucial importance and must be included in double layer theories. The electron density distribution of the metal is also altered by the presence of the fluid and extends a small distance further from the jellium edge in the presence of the liquid. The electrostatic potential drop across the interface can be written as the sum of jellium, ion and dipole terms. In all cases, it was found that the jellium contribution increases due to the presence of the fluid while, as expected, the net electrostatic potential drop across the interface is reduced. Compared with the vacuum values the potential reductions for the pure solvent range from 6—10% (i.e., 0.42—0.24 volts) as 0 /a is varied 3 r from 2.0—3.0. It is interesting to observe that these values do not differ greatly from those given by the monolayer model [67,68]. In the presence of ions, the jellium contribution to the potential drop is not greatly altered from the pure solvent value. On the other hand, the ion and dipolar contributions are relatively large and are sensitive to both  Chapter 5. Summary and Conclusion  salt concentration and ion size.  130  However, these terms cancel to give a net solution  contribution comparable in magnitude with that of the pure solvent. This demonstrates that the ion and dipole contributions are very strongly coupled and should never be treated as independent quantities in theories of the double layer. Fluids between planar surfaces were also considered using the HNC and RHNC the ories and grand canonical Monte Carlo methods. A general reduction of the wall-wall OZ equation was presented for nonspherical particles between walls. This allows the wall-wall potential of mean force per unit area and the net pressure to be calculated. The associated OZ equation for a fluid confined between planar walls was also reduced for fluids of nonspherical particles. Previously this wall-solvent OZ equation had only been considered for spherical particles [109]. These equations are used to investigate the effects of confinement on the net pressure between surfaces, fluid structure and phase behaviour. For simple hard sphere fluids between hard walls, the net pressure between surfaces as a function of the wall-wall separation was obtained from the HNC and RHNC theo ries. It was found that the pressure had little dependence on the closure approximation. The HNC and, in particular, the RHNC results showed very good agreement with grand canonical Monte Carlo simulations [123] and with results obtained from the inhomoge neous PY theory [113]. The small discrepancies between the RHNC and grand canonical Monte Carlo results were comparable to those between the inhomogeneous PY and grand canonical Monte Carlo results. HNC results for the hard sphere distribution function be tween hard walls were also obtained for a few wall-wall separations in the range of 1.2—10 solvent diameters. The distribution function at a separation of 10 solvent diameters was in good agreement with grand canonical Monte Carlo results except very near contact. The same is true for a separation of 3.1 solvent diameters when compared with results from the inhomogeneous PY theory. While the HNC contact densities p g(d/2; 1) are 5  Chapter 5. Summary and Conclusion  131  incorrect, the contact theorem was found to still give good results for the net pressure between hard walls. Dipolar hard sphere fluids between hard walls were considered using the RHNC ap proximation with HNCP bridge functions. The net pressure was calculated for two values of the reduced solvent dipole moment. It was found that the solvent-mediated pressure is similar to the hard sphere-mediated pressure and the net pressure is smaller for larger dipole moments. The wall-solvent structure was also examined for several wall-wall sep arations. The distribution functions for fluids of Lennard-Jones particles confined between walls were also calculated from the HNC and RHNC theories. The bulk direct correlation function needed in the wall-solvent integral equation theories was obtained from pertur bation theory [124,125]. These results were compared with computer simulations. The Lennard-Jones parameters for the wall-particle pair potential were equal to the values used for the particle-particle pair interaction. For a Lennard-Jones fluid confined between Lennard-Jones 10—4 plates, the HNC distribution functions showed very good agreement with simulations. It is worthwhile noting that these results were obtained for a bulk Lennard-Jones state point near the coexistence curve. When the Lennard-Jones fluid was confined between hard walls, grand canonical Monte Carlo results indicate that the wall-solvent distribution function is highly sensitive to the bulk state parameters and wall-wall separation. In particular, it was found that the fluid, which is liquid in the bulk and at large wall-wall separations, evaporates at narrow separations. This separation-induced phase transition arises from the loss of attractive material (displaced by the inert walls) and corresponds to a decrease in the interaction energy per particle relative to the bulk. The HNC and RHNC results showed relatively good agreement with simulations for bulk state points far from the coexistence curve. As the state point approaches the coexistence curve, the IINC and RHNC theories fail  Chapter 5. Summary and Conclusion  132  to describe solvent dewetting at a hard solvophobic wall. Grand canonical Monte Carlo was used to study confined Lennard-Jones systems in detail. For one subcritical state point, (T* from grand canonical Monte Carlo was at l  =  =  1.0, p  0.75), the coexistence separation  6. The metastable branches on either side  of coexistence were also explored. Spinodal cavitation was found at a separation of 1* and spinodal condensation occurred near  =  =  5  15. The net pressure between walls was of  particular interest. As the wall-wall separation was decreased, the pressure was found to be attractive and increase monotonically on the liquid branch as the spinodal cavitation separation was approached. As previously reported [115], a discontinuous increase in the pressure was observed at the coexistence separation. On the metastable gaseous branch, the attraction decayed extremely slowly. In general, for systems in which a fluid with attractive forces is confined between solvophobic surfaces, separation-induced cavitation is certain to occur in the limit that the state parameters approach the bulk coexistence curve. The separations at which coexistence and spinodal points occur depend upon the specific system and the distance of the bulk state point from coexistence. Although experimental results for aqueous systems at standard temperature and pressure indicate a much longer ranged attraction (measurable to lOOnm), both experiments and simulations show qualitatively similar behaviour. Our simulation results suggest that a separation-induced phase transition might be the origin of the experimentally measured attraction between hydrophobic surfaces immersed in water.  Bibliography  [1] J. D. Porter and A. S. Zinn. Ordering of liquid water at metal surfaces in tunnel junction devices. J. Phys. Chem., 97:1190—1203, 1993. [2] Michael F. Toney, Jason N. Howard, Jocelyn Richer, Gary L. Borges, Joseph 0. Gordon, Owen R. Melroy, David G. Wiesler, Dennis Yee, and Larry B. Sorensen. Voltage-dependent ordering of water molecules at an electrode-electrolyte interface. Nature, 368:444, 1994. [3] J.P. Hansen and I.R. McDonald. Theory of Simple Liquids. Academic Press, Lon don, second edition, 1986. [4] J. A. Alonso and N. H. March. Electrons in Metals and Alloys. Academic Press, New York, 1989. [5] Richard N. Zare. Angular Momentum. John Wiley and Sons, New York, 1988. [6] L. Blum and A. J. Torruella. Invariant expansion for two-body correlations: Ther modynamic functions, scattering, and the Ornstein-Zernike equation. J. Chem. Phys., 56:303—310, 1972. [7] L. S. Ornstein and F. Zernike. Proc. Akad. Sci. (Amsterdam), 17:793, 1914. [8] R. Evans, P. Tarazona, and U. Marini Bettlol Marconi. On the failure of certain integral equation theories to account for complete wetting at solid-fluid interfaces. Mol. Phys., 50:993—1011, 1983. [9] Roland Kjellander and Stjepan Marelja. Inhomogeneous Coulomb fluids with image interactions between planar surfaces. I. J. Chem. Phys., 82:2122—2135, 1985. [10] P. H. Fries and G. N. Patey. The solution of the hypernetted-chain approximation for fluids of nonspherical particles, a general method with applications to dipolar hard spheres. J. Chem. Phys., 82:429—440, 1985. [11] P. G. Kusalik and 0. N. Patey. The solution of the reference hypernetted-chain approximation for water-like models. Mol. Phys., 65:1105—1119, 1988. [12] P. G. Kusalik and 0. N. Patey. On the molecular theory of aqueous electrolyte solutions. I. The solution of the RHNC approximation for models at finite concen tration. J. Chem. Phys., 88:7715—7738, 1988. 133  Bibliography  134  [13] P. G. Kusalik and 0. N. Patey. On the molecular theory of aqueous electrolyte so lutions. II. Structural and thermodynamic properties of different models at infinite dilution. J. Chem. Phys., 89:5843—5851, 1988. [14] P. 0. Kusalik and G. N. Patey. On the molecular theory of aqueous electrolyte solutions. IV. Effects of solvent polarizability. J. Chem. Phys., 92:1345—1358, 1990. [15] P. G. Kusalik. The structural, thermodynamic and dielectric properties of elec trolyte solutions: A theoretical study. PhD thesis, University of British Columbia, Vancouver, Canada, 1988. [16] J. Perkyns. A theoretical study of ammonia-salt mixtures. PhD thesis, University of British Columbia, Vancouver, Canada, 1990. [17] L. Blum. Invariant expansion. II. The Ornstein-Zernike equation for nonspherical molecules and an extended solution to the mean spherical model. J. Chem. Phys., 57:1862—1869, 1972. [18] L. Blum. Invariant expansion. III. The general solution of the mean spherical model for neutral spheres with electrostatic interactions. J. Chem. Phys., 58:3295—3303, 1973. [19] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes: The Art of Scientific Computing. Cambridge University Press, New York, 1986. [20] J. M. van Leeuwen, M. J. Groeneveld, and J. de Boer. New method for the calcu lation of the pair correlation function. Physica, 25:792—808, 1959. [21] E. Meeron. Nodal expansions. III. Exact integral equations for the particle corre lation functions. J. Math. Phys., 1:192—201, 1960. [22) T. Morita. Theory of classical fluids: Hypernetted-chain approximation. III. Frog. Theor. Phys., 23:829—845, 1960. [23] G. S. Rushbrooke. On the hypernetted-chain approximation in the theory of clas sical fluids. Physica, 26:259, 1960. [24] L. Verlet. Nuovo Cimento, 18:77, 1960. [25) Richard N. Stell. The Equilibrium Theory of Classical Fluids, pages 11—171. W. A. Benjamin, New York, 1964. [26] Phil Attard and 0. N. Patey. Hypernetted-chain closure with bridge diagrams. Asymmetric hard sphere mixture. J. Chem. Phys., 92:4970—4982, 1990.  Bibliography  135  [27] L. Verlet and J.-J. Weis. Equilibrium theory of simple fluids. Phys. Rev. A, 2:939— 952, 1972. [28] L. L. Lee and D. Levesque. Perturbation theory for mixtures of simple liquids. Mol. Phys., 26:1351—1370, 1973. [29] Anatol Malijevsky and Stanislav Labik. The bridge function for hard spheres. Mol. Phys., 60:663—669, 1987. [30] Stanislav Labik and Anatol Malijevsky. Bridge function for hard spheres in high density and overlap regions. Mol. Phys., 67:431—438, 1989. [31] D. Henderson and M. Plischke. Pair correlation function in a fluid with density inhomogeneities: parameterization for hard spheres near a hard wall. Proc. R. Soc. Lond., A400:163—174, 1985. [32] Anatol Malijevsky, R. Pospisil, W.R. Smith, and Stanislav Labik. The Ornstein Zernike equation for hard spheres near a hard wall. a rapid method of numerical solution and an accurate new RHNC theory. Mol. Phys., 72:199—213, 1991. [33] J. M. Caillol. A method for solving the hypernetted-chain approximation for molec ular fluids. Chem. Phys. Lett., 121:347—350, 1985. [34] J. M. Caillol. A method for solving the hypernetted-chain approximation for molec ular fluids. Chem. Phys. Lett., 156:357—362, 1989. [35] M.P. Allen and D. J. Tildesley. Computer Simulations of Liquids. Oxford Univer sity, Oxford, 1987. [36] D. J. Adams. Grand canonical ensemble Monte Carlo for Lennard-Jones fluids. Mol. Phys., 29:307—311, 1975. [37] D. J. Adams. Calculating the high-temperature vapour line by Monte Carlo. Mol. Phys., 37:211—221, 1979. [38] P. G. Kusalik and G. N. Patey. On the molecular theory of aqueous electrolyte so lutions. III. A comparison between Born-Oppenheimer and McMillan-Mayer levels of description. J. Chem. Phys., 89:7478—7484, 1988. [39] D. J. Isbister and B. C. Freasier. Adsorption of dipolar hard spheres onto a smooth hard wall in the presence of an electric field. J. Stat. Phys., 20:331—345, 1979. [40] J. C. Rasiah, D. J. Isbister, and J. Eggebrecht. Polarization density profiles for dipoles against an electrified wall in the MS and RLHNC approximations. J. Chem. Phys., 75:5497—5502, 1981.  Bibliography  136  [41] Steven L. Carnie and Derek Y.C. Chan. The structure of electrolytes at charged surfaces: Ion-dipole mixtures. J. Chem. Phys., 73:2949—2957, 1980. [42] L. Blum and D. Henderson. Mixtures of hard ions and dipoles against a charged wall: The Ornstein-Zernike equation, some exact results, and the mean spherical approximation. J. Chem. Phys., 74:1902—1910, 1981. [43] J. P. Badiali. Structure of the dipolar hard sphere fluid near a neutral hard wall. The orientational structure. Mol. Phys., 55:939—949, 1985. [44] C. E. Woodward and S. Nordholm. Dipoles at charged solid-liquid interfaces. Mol. Phys., 60:415—439, 1987. [45] W. Dong, M. L. Rosinberg, A. Perera, and G. N. Patey. A theoretical study of the solid-electrolyte solution interface. I. Structure of a hard sphere ion-dipole mixture near an uncharged hard wall. J. Chem. Phys., 89:4994—5009, 1988. [46] 0. N. Patey and 0. M. Torrie. Water and salt water near charged surfaces: A discussion of some recent theoretical results. Chemica Scripta, 29A:39—47, 1989. [47] 0. M. Torrie, P. G. Kusalik, and 0. N. Patey. Molecular solvent model for an electrical double layer: Reference hypernetted-chain (RHNC) results for solvent structure at a charged surface. J. Chem. Phys., 88:7826—7840, 1988. [48] G. M. Torrie, P. 0. Kusalik, and G. N. Patey. Molecular solvent model for an electrical double layer: Reference hypernetted-chain results for ion behaviour at infinite dilution. J. Chem. Phys., 89:3285—3294, 1988. [49] 0. M. Torrie, P. G. Kusalik, and G. N. Patey. Molecular solvent model for an electrical double layer: Reference hypernetted-chain results for potassium chloride solutions. J. Chem. Phys., 90:4513—4527, 1989. [50] 0. M. Torrie, A. Perera, and G. N. Patey. Reference hypernetted-chain theory for dipolar hard spheres at charged surfaces. Mol. Phys., 67:1337—1353, 1989. [51] M. Kinoshita and M. Harada. Numerical solution of the RHNC theory for fluids of non-spherical particles near a uniform planar wall: A study with dipolar hard spheres as an example system. Mol. Phys., 79:145—167, 1993. [52] M. Kinoshita and M. Harada. Numerical solution of the RHNC theory for waterlike fluids near a large macroparticle and a planar wall. Mol. Phys., 81:1473, 1994. [53] M. Moradi and G. Rickayzen. An approximate density functional for an inhomo geneous dipolar fluid. Mol. Phys., 68:903—915, 1989.  Bibliography  137  [54] D. Levesque and J.-J. Weis. Orientational structure of dipolar hard spheres near a hard neutral wall. J. Stat. Phys., 40:29—38, 1969. [55] 5. H. Lee, J. C. Rasaiah, and J. B. Hubbard. Molecular dynamics study of a dipolar fluid between charged plates. J. Chem. Phys., 85:5232—5237, 1986. [56] C. Y. Lee, J. A. McCammon, and P. J. Rossky. The structure of liquid water at an extended hydrophobic surface. J. Chem. Phys., 80:4448—4455, 1984. [57] J. P. Valleau and Anne A. Gardner. Water-like particles ar surfaces. I. The un charged, unpolarized surface. J. Chem. Phys., 86:4162—4170, 1987. [58] S. B. Zhu and G. W. Robinson. Structure and dynamics of liquid water between plates. J. Chem. Phys., 94:1403, 1991. [59] Zixiang Tang, L. E. Scriven, and H. T. Davis. Structure of a dipolar hard sphere fluid near a neutral hard wall. J. Chem. Phys., 96:4639—4645, 1992. [60] Anne A. Gardner and J. P. Valleau. Water-like particles ar surfaces. II. In a double layer and at a metallic surface. J. Chem. Phys., 86:4171—4176, 1987. [61] J. Hautman, J. W. Halley, and Y.-J. Rhee. Molecular dynamics simulation of water between two ideal classical metal walls. J. Chem. Phys., 91:467—472, 1989. [62] J. P. Badiali. Structure of a polar fluid near a wall. Exact asymptotic behavior of the profile, relation with the electrostriction phenomena and the Kerr effect. J. Chem. Phys., 90:4401—4412, 1989. [63]  Q. Zhang, J. P. Badiali, and W. H. Su. Polar molecules in planar interfaces. Role of the triplet direct correlation function in the asymptotic behavior of the profile. Chem. Phys. Lett., 92:4609—4619, 1990.  [64] J. P. Badiali, M. L. Rosinberg, and J. Goodisman. Effects of solvent on properties of the liquid metal surface. J. Electroanal. Chem., 130:31—45, 1981. [65] J. P. Badiali, M. L. Rosinberg, and J. Goodisman. Contribution of the metal to the differential capacity of an ideally polarizable electrode. J. Electroanal. Chem., 143:73—88, 1983. [66] J. P. Badiali, M. L. Rosinberg, and J. Goodisman. The metal in the polarizable interface coupling with the solvent phase. J. Electroanal. Chem., 150:25—31, 1983. [67] Wolfgang Schmickler. The potential of zero charge of jellium. Chem. Phys. Lett., 99:135—139, 1983.  Bibliography  138  [68] Wolfgang Schmickler. A jellium-dipole model for the double layer. J. Electroanal. Chem., 150:19—24, 1983. [69] J. P. Badiali, M. L. Rosinberg, F. Vericat, and L. Blum. A microscopic model for the liquid-ionic solution interface. J. Electroanal. Chem., 158:253—267, 1983. [70] Wolfgang Schmickler and Douglas Henderson. The interphase between jellium and a hard sphere electrolyte. a model for the electric double layer. J. Chem. Phys., 80:3381—3386, 1984. [71] Douglas Henderson and Wolfgang Schmickler. The capacitance for a single crystal face of an FCC metal/electrolyte interface. J. Chem. Phys., 82:2825—2829, 1985. [72] Wolfgang Schmickler and Douglas Henderson. The interphase between jellium and a hard sphere electrolyte: Capacity-charge characteristics and dipole potentials. J. Chem. Phys., 85:1650—1657, 1986. [73] D. L. Price and J. W. Halley. A new model of the differential capacitance of the double layer. J. Electroanal. Chem., 150:347—353, 1983. [74] J. W. Halley, B. Johnson, D. L. Price, and M. Schwalm. Quantum theory of the double layer: A model of the electrode-electrolyte interface. Phys. Rev. B, 3 1:7695— 7709, 1985. [75] J. W. Halley and D. L. Price. Quantum theory of the double layer: Model including solvent structure. Phys. Rev. B, 35:9095—9102, 1987. [76] D. L. Price and J. W. Halley. Electronic structure of the metal-electrolyte surfaces: Three-dimensional calculation. Phys. Rev. B, 38:9357—9367, 1988. [77] S. Amokrane, V. Russier, and J. P. Badiali. A model for the metal-solvent coupling at a charged interface: Effects on the differential capacitance. Surface Science, 217:425—444, 1989. [78] S. Amokrane and J. P. Badiali. A new analysis of the differential capacitance of an ideally polarized electrode. J. Electroanal. Chem., 266:21—35, 1989. [79] S. Amokrane and J. P. Badiali. A new analysis of the differential capacitance of an ideally polarized electrode. Part II. Non-aqueous solvents. J. Electroanal. Chem., 297:377—398, 1991. [80] N.D. Lang and W. Kohn. Theory of metal surfaces: Charge density and surface energy. Phys. Rev. B, 1:4555—4568, 1970.  Bibliography  139  [81] M. Kinoshita and D. R. Berard. Numerical methods for analyzing the bulk and surface-induced structure of electrolyte solutions using integral equation theories. submitted to JCP, 1994. [82] P. Hohenberg and W. Kohn. Inhomogeneous electron gas. Phys. Rev., 136B:864— 871, 1964. [83] W. Kohn and L. J. Sham. Self-consistent equations including exchange and corre lation effects. Phys. Rev., 140A:1133—1138, 1965. [84] J. W. Perram and L. R. White. Structure of the liquid/vapour and liquid/solid interfaces. Discuss. Faraday Soc., 59:29—37, 1975. [85] D. Henderson, Farid F. Abraham, and John A. Barker. The Ornstein-Zernike equation for a fluid in contact with a surface. Mol. Phys., 31:1291—1295, 1976. [86] D. Henderson, Farid F. Abraham, and John A. Barker. The Ornstein-Zernike equation for a fluid in contact with a surface (errata). Mol. Phys., 32:1785, 1976. [87] II. E. Sullivan and G. Stell. Structure of a simple fluid near a wall. I. Structure near a hard wall. J. Chem. Phys., 69:5450—5457, 1978. [88] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions. Dover, New York, 1970. [89] J. II. Jackson. Classical Electrodynamics. Wiley, New York, second edition, 1975. [90] David J. Griffiths. Introduction to Electrodynamics. Prentice Hall, New Jersey, second edition, 1989. [91] P. Gies and R. R. Gerhardts. Self-consistent calculation of electron-density profiles at strongly charged jellium surfaces. Phys. Rev. B, 33:982—989, 1986. [92] D. R. Berard and G. N. Patey. The application of integral equation theories to fluids of nonspherical particles near a uniform planar wall. J. Chem. Phys., 95:5281—5288, 1991. [93] F. Wigner. On the interaction of electrons in metals. Phys. Rev., 36:1002—1011, 1934. [94] D. F. C. Morris. Ionic radii and enthalpies of hydration of ions. Struct. Bonding, 4:63—82, 1968. [95] G. M. Torrie, P. G. Kusalik, and G. N. Patey. Theory of the electrical double layer: Ion size effects in a molecular solvent. J. Chem. Phys., 91:6367—6375, 1989.  Bibliography  140  [96) G. Stell, G. N. Patey, and J. S. Høye. Dielectric constants of fluid models: Statis tical mechanical theory and its quantitative implementation. Adv. Chem. Phys., 48:183—328, 1981. [97] T. D. Blake and J. A. Kitchener. J. Chem. Soc. Faraday Trans. 1, 68:1435, 1972. [98] J. N. Israelachvili and R. M. Pashley. J. Colloid Interface Sci., 98:500, 1984. [99] P. M. Claesson and H. K. Christenson. Very long range attractive forces be tween uncharged hydrocarbon and fluorocarbon surfaces in water. J. Phys. Chern., 92:1650—1655, 1988. [100] Y. I. Rabinovich and B. V. Derjaguin. Colloid Surf., 30:243, 1988. [101] J. L. Parker and P. M. Claesson. Direct measurements of the attraction between solvophobic surfaces in ethylene glycol and mixtures with water. Langmuir, 8:757— 759, 1992. [102] Farid. F. Abraham. The interfacial density profile of a Lennard-Jones fluid in contact with a (100) Lennard-Jones wall and its relationship to idealized fluid/wall systems: A monte carlo simulation. J. Chem. Phys., 68:3713—3716, 1978. [103] J. E. Lane and Thomas H. Spurling. Monte Carlo simulation of the effects of adsorption on interparticle forces. Aust. J. Chem., 33:231—239, 1980. [104] J. E. Lane and Thomas H. Spurling. The thermodynamics of a new type of surface transition. Aust. J. Chem., 34:1529—1533, 1981. [105] W. van Megen and I. K. Snook. Physical adsorption of gases at high pressure III. Adsorption in slit-like pores. Mol. Phys., 54:741—755, 1985. [106] I. K. Snook and W. van Megen. Solvation forces in simple dense fluids. I. J. Chem. Phys., 72:2907—2913, 1980. [107] M. Schoen, D. J. Diestler, and J. H. Cushman. Fluids in micropores. I. structure of a simple classical fluid in a slit pore. J. Chem. Phys., 87:5464—5476, 1987. [108] J. J. Magda, M. Tirrell, and H. T. Davis. Molecular dynamics of narrow, liquid filled pores. J. Chem. Phys., 83:1888—1901, 1985. [109] M. Lozada-Cassou. The force between two planar electrical double layers. J. Chem. Phys., 80:3344—3349, 1984. [110] M. Lozada-Cassou and E. Dcaz-Herrera. Three point extension for hypernetted chain and other integral equation theories: Numerical results. J. Chem. Phys., 92:1194—1210, 1990.  Bibliography  141  [111] M. Lozada-Cassou and E. Dcaz-Herrera. Three point extension hypernetted chain, conventional hypernetted chain, and superposition approximations: Numerical re sults for the force between plates. J. Chem. Phys., 93:1386—1398, 1990. [112] E. GonzJez-Tovar, M. Lozada-Cassou, and W. Olivares. Electrokinetic flow in ultrafine slits: Ionic size effects. J. Chem. Phys., 94:2219—2231, 1991. [113] Roland Kjellander and Sten Sarman. On the statistical mechanics of inhomoge neous fluids in narrow slits. An application to a hard sphere fluid between hard walls. Chem. Phys. Lett., 149:102—108, 1988. [114] Sten Sarman. The influence of the fluid-wall interaction potential on the structure of a simple fluid in a narrow slit. J. Chem. Phys., 92:4447—4455, 1990. [115] P. Tarazona, U. Marini Bettlol Marconi, and R. Evans. Phase equilibria of fluid interfaces and confined fluids. non-local versus local density functionals. Mol. Phys., 60:573—595, 1987. [116] R. Evans and U. Marini Bettlol Marconi. Phase equilibria and solvation forces for fluids confined between parallel walls. J. Chem. Phys., 86:7138—7148, 1988. [117] P. C. Ball and R. Evans. Structure and adsorption at gas-solid interfaces: Layering transitions from a continuum theory. J. Chem. Phys., 89:4412—4423, 1988. [118] Ziming Tan and Keith E. Gubbins. Selective adsorption of simple mixtures in slit pores: A model of methane-ethane mixtures in carbon. J. Phys. Chem., 96:845— 854, 1992. [119] D. E. Sullivan, D. Levesque, and J. J. Weis. Structure of a simple fluid near a wall. II. Comparison with Monte Carlo. J. Chem. Phys., 72:1170—1174, 1980. [120] R. Evans. Liquids at Interfaces. North-Holland, Amsterdam, 1990. [121] Phil Attard, D. R. Berard, C. P. Ursenbach, and G. N. Patey. Interaction free energy between planar walls in dense fluids: An Ornstein-Zernike approach with results for hard-sphere, Lennard-Jones and dipolar systems. Phys. Rev. A, 44:8224— 8234, 1991. [122] W. A. Steele. The Interaction of gases with solid surfaces. Pergamon Press, Oxford, 1974. [123] Gunnar Karlström. Solvation forces studied by grand canonical Monte Carlo sim ulations of hard spheres between hard walls. Chem. Scripta, 25:89—91, 1985.  Bibliography  142  [124] John D. Weeks and Davis Chandler. Role of repulsive forces in determining the equilibrium structure of simple liquids. J. Chern. Phys., 54:5237—5247, 1971. [125] Hans C. Andersen, John D. Weeks, and Davis Chandler. Relationship between the hard-sphere fluid and fluids with realistic repulsive forces. Phys. Rev. A, 4:1597— 1607, 1971. [126] J.-P. Hansen and L. Verlet. Phase transitions of the Lennard-Jones system. Phys. Rev., 184:151—161, 1969. [127] D. R. Berard, Phil Attard, and G. N. Patey. Cavitation of a Lennard-Jones fluid between hard walls, and the possible relevance to the attraction measured between hydrophobic surfaces. J. Chem. Phys., 98:7236—7244, 1993. [128] H. K. Christenson, J. Fang, B. W. Ninham, and J. L. Parker. Effects of divalent electrolyte on the hydrophobic attraction. J. Phys. Chem., 94:8004—8006, 1990. [129] Jacob Israelachvili. Intermolecular and surface forces. Academic Press, New York, second edition, 1992.  Appendix A  Properties of the Wigner symbols, rotation matrices and rotational invariants.  In deriving much of the theory presented in this thesis, the orthogonality and symme try properties of the Wigner 3-j symbols, rotation matrices and rotational invariants are used extensively. For convenience, many of the properties of these functions are listed here. All of these relationships can be found in Ref.  [5J or derived from relationships  therein. The matrix R associated with rotating a rank-rn spherical tensor has the elements  R(f) (known as the Wigner generalized spherical harmonics or the symmetric-top wavefunctions when properly normalized) where  =  (q, 0, x) denotes the Euler angles.  The elements R() form an orthogonal basis set that spans the angle space of an asymmetric particle and can be written as  ) R ( 4  =  e’d (0)&  (A.l)  ,  where the functions d(0) are real. The index m takes on all positive integer and half integer values;  and  t’  are given by =  —rn, —rn + 1,  ...,  rn  =  —m, —rn + 1,  ...,  m  Since the rotation matrices are unitary (RRt  =  —  —  1, m,  (A.2a)  1, m.  (A.2b)  I, where Rt is the complex transpose  of R and I is the identity matrix), we have =  [RtA,(x,o,)j* 143  =  [R(,o,x)j*  ,  (A.3)  Appendix A. Properties of the Wigner symbols,  144  ...  where R is the complex conjugate of R, and the elements R must satisfy the sum rules’ (A.4)  d(0) d,A,I(0) =  >,,(qf, 9, x) 1 R(q, 0, x) R  A’  A’  d,A(0)d,,A(0)  =  =  A’  .  (A.5)  A  A  The functions d(0) have several symmetry properties which will be of use. These are d(0)  =  d(r  —  0)  (A.6a)  = d,(0) (—)‘d(0)  (A.6b)  —0) = (—)‘d(—0) d ( 1  (A.6c)  = ()m+LLdm(0)  (A.6d)  where j = —p. Using Eq. (A.6b), we see that R() = (—)‘R(f)  (A.7)  .  In addition, (A.8a)  dL(O) =  d(+r) = Hm±IL’6,  (A.8b)  A few useful integrals, sums and products of the Wigner spherical harmonics follow:  f  R()  87r 2 6mn6,’a,”5,w , 2m+1 Imn i fmnl 2 = 87r  (A.9)  dfz =  f  () = (2l 1 R(f) R  --  1)  IA’A  RA(fZ)  =  (mn 1) Rl*()  (A.1O)  (A.n)  ,  (2l + 1) I-I’ L  m n 1 1  ,  ‘mn 1” frnnl\ ,  x R() R,(f)  ,  ‘The sum rules given in Eqs (3.77) and (3.78) of Zare [5] are incorrect.  (A.12)  Appendix A. Properties of the Wigner symbols,  where  145  ...  ) is a Wigner 3-j symbol.  ) then R() is equivalent 1 Applying two rotations in succession and in the order R(f to a single rotation R(f), =  R(fZ) 1 R,,L(Zl) () 1 R  =  .  (A.13)  IL”  Using the unitary property, the inverse expression is R(fZ)[ft(1i)]  =  ]. [R( ) 1  =  (A.14)  Ii”  These relations are useful when analyzing the pair functions for molecular fluids (see Appendix B). The 3-j symbols have a few useful symmetry properties. Exchanging two column is equivalent to multiplication by (m fl  —  m+n+l  (_)m+n+l.  (m 1  For example, (n m l  m+n+l  = m+n+l  (1 n m  (A 15  V)  Two such exchanges gives (mn l  —  (n 1 m  —  (1  A 16  V)VIL)ILV)  In addition,  (:)  = m+n+l  The orthogonality properties of the 3-j symbols are (inn l” (mn l’’\ jt  (2l +1)  (: ) (  —  —  ,  =  6ll’5)\’ (21+1)  ‘  ( A 18 ) (A.19)  Appendix A. Properties of the Wigner symbols,  146  ...  If, in Eq. (A.18), the sum over,\ is included, we get frnnl\ frnni’\  6.  (A.20)  =  A general expression for evaluating the Wigner symbols, as well as FORTRAN code, can be found in Appendix A of Ref. [5]. The orientational space of one asymmetric particle with respect to another is spanned by the rotational invariants mni  mnl(, 2, p21) /W  The complex conjugate of mnl*(, 2,  21 r )  m R i/  mnl  (A.21)  .  is (m  mni  =  21 ( 0 R, )) 2 () R,(  IL’1)’’  fl  /2 V  l Rrn* 21 ( 0 R )) 2 /t () R( 1 X) L  .  (A.22)  A few useful integrals over the rotational invariants are  J J  ‘(1, 2, ) 21 r  4minili*(f \  /L1111  2 512ir (fmnl) 6 (2n + 1) (2m + 1) (21 +  2,  =  6 5127r  \  )  (A.23)  1)8mm161hh111  Jrnfll(j2) mmfuh1(21\ m2z2l2(12) /L1L’l  21 df r ) 1 dfZ 2  /.42L’2  fmnlfmlnhllfm2n2l2  1 dfZ dfZ 2  (m m 1m 2 i 1L2)  (  1 2 2)  / i 11 i 2 2)  {  1 m 2 m m rj  1 j  11  12  }  .(A.24)  Since the rotational invariants are independent of the orientation of the pair, the in tegration over the 121 is usually omitted from the orthogonality relation (A.23) and, . 2 consequently, the normalization constant is reduced by a factor of 87r  Appendix B  The electrostatic interaction potential between a molecule and a uniform surface charge distribution.  In this appendix, the metal wall-particle pair potential is reduced to a rotational invariant expansion. Consider the static charge distribution n(r) of a molecule of species near a planar polarizable surface having a charge distribution n(z). The electrostatic interaction energy between the surface and molecule is given by e2JJfldrdt,  =  up(O1)  (B.1)  where 0 and 1 denote the coordinates of the wall and molecule 3, respectively, and Zt =  t  .  We assume that the charge distribution of the molecule is localized near the  center of mass and further that the surface and molecular charge distributions do not 1 and r overlap. Choosing an origin such that the molecule is at z  =  1 + r’, we can use z  the Taylor series expansion  —  where t’  =  t  —  z, t’  =  00  1  1  ri  =  1t9,  r’  t’ =  —  r’  (t)m+l  =  Ir’I and  m  Pm(cos ),  (B.2)  is the angle between the vectors t’ and r’.  Using the spherical harmonic addition theorem, ) 7 Fi(cos  =  i’) (’)R, 0 R ( “I  147  ,  (B.3)  Appendix B. The electrostatic pair potential  Eq. (B.2) and the identity drdt u(O1)  dr’dt’, Eq. (B.1) reduces to  J(r)mnp(r + ziê) R (’) dr’ 0  2 e  =  =  148  ...  mp’  f  2e  =  mi’  n(zi +  f  Zt’)  (t)m+l  + ft’ dt’ n(zi (j!)m+l  R(’) dt’  Zt’)  dz’  (B.4)  Pm(CO5O)  where Eq. (2.4.9) has been used to defines the molecular multipole moments cos 8  =  z’/t’  =  (zi  —  and  )/t’ denotes the polar angle of t’. Writing the multipole moments 1 z  in the molecular frame [using Eq. (2.4.7)], we get u,3(O1)  n(zt)  =  t_zim+1Pm050t  ,  (B.5)  mi  where  denotes the orientation of the molecule. When m  =  0, the integral in Eq. (B.5)  is simply the scalar potential qf(zi) a distance z from the surface and arising from the charge distribution n(z),  =  When m  =  f  z’  dt.  (B.6)  1, we get n(z) OZ  It—ziI  dt =  f  n(zt) 2 1 t—z  cosOdt  =  —E(zi)  (B.7)  where E(zi) is the electric field normal to the surface. In general,  I  1  n(zt)  m a  n(zt)  (B.8)  =  Thus, we have u(01)  1 =  Om(z)  t9zr  ) 1 R(1Z  (B.9)  or, using the rotational invariant expansion Eq. (3.2.1), 0mm ju 0 ;w(01)  —  —  (_)mI2m + 1 —m 3 L 1 Q ;/ mm m! f°  8m(zi)  Oz  (B.10)  Appendix C  Numerical considerations  In this appendix, the computer algorithms used to obtain the wall-particle pair and direct correlation functions are outlined. The computer programs used were written for a general multicomponent fluid composed of particles of C 2 or higher symmetry. The algorithm used to solve the fluid structure near an isolated surface and between two surfaces are similar irrespective of the surface model. Since the equations are non-linear, the correlation functions must be solved iteratively. A detailed method for treating longrange pair functions in the integral equation theories has been described elsewhere [15]. The basic algorithm used to solve the wall-particle correlation functions is summarized below. 1. Solve the bulk OZ equations using the computer program written by Kusalik [15]. 2. Calculate the wall-particle bridge functions using the methods of Henderson and Plischke [31] or Attard and Patey [26]. 3. Calculate  Onn(—) Ov;w  (z) using Eq. (3.2.14).  4. Calculate the Hankel transforms of the projections cj(r),  k) [18,10].  5. Obtain an initial guess for h(z). When no initial guess is available, the choice h(z)  (—1,  z<d/2,  10,  z>d/2,  (C.1)  =  is sufficient for the systems considered. Onn(+) Onn(+) 6. Calculate the Fourier transforms of the projections hv;w (z), 0 h ( ;  7. Calculate  (k) using Eq. (3.2.18). 149  150  Appendix C. Numerical considerations  8. Inverse Fourier transform  .-Onn(+) l7 Ov•w  (k) to get  Omm(+) iOu.w  (z).  z) and the projections of potential of mean force, w(z) [see 9. Calculate Eq. (2.2.21)]. 10. Using the RHNC closure Eq. (3.2.24), calculate c(z). Onn(+) 11. Calculate a new estimate for hv;w (z) from  Wth(z) + (1  —  (z) + c(z)] 3 W) [q,  ,  (C.2)  and c(z) are the most recent estimates and W is where h(z), interval 0 W < 1. a mixing parameter in the 12. Repeat steps 6—12 until converged. Convergence was obtained when the pair cor relation function was constant to five significant figures. Step 11 is characteristic of the Picard iteration strategy [52]. A faster convergence strat egy has recently been developed by Kinoshita and Harada [52]. When the wall-particle potential depends on the surface and fluid structure in a self-consistent manner, an ad ditional iteration cycle is required. At a dielectric surface, the algorithm is summarized below. 1. Calculate the wall-particle pair and direct correlation functions at an inert wall. 2. Calculate the dielectric wall-particle pair potential, Eqs (3.3.3) and (3.3.13), using either the cavity, Eq. (3.3.22), or superposition, Eq. (3.3.23), approximations. 3. Calculate the wall-particle pair and direct correlation functions using the new wall particle pair potential. 4. Repeat steps 2—4 until converged pair potential is self-consistent. At a metallic surface, the Schrödinger Eq. (3.3.29) must also be solved self-consistently. The basic method of Gies and Gerhardts for jellium-vacuum interfaces was followed using the iteration strategies of Kinoshita et al. [52]. When the jellium metal is immersed in fluid, the basic algorithm is summarized below.  Appendix C. Numerical considerations  151  1. Calculate the wall-particle pair and direct correlation functions at an inert wall as defined above. 2. Calculate the electron density distribution, electrostatic potential and effective po tential at a jellium-vacuum interface using Eqs (3.3.33), (3.3.39) and (3.3.30) with ext() = 0. 3. Calculate the wall-particle pair potential using Eq. (3.3.47). 4. Solve the OZ equation for the wall-particle direct and pair correlation functions using the above algorithm. 5. Calculate the external potential ut(z) from Eq. (3.3.46). 6. Calculate the electron density distribution, electrostatic potential and effective po tential for the jellium metal using ut(z). 7. Repeat steps 3—7 until the electron density distribution and fluid structure are self-consistent.  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

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

Comment

Related Items