STRUCTURE AND PROPERTIES OF FLUIDS NEAR INERT ANDMETALLIC SURFACESByDaniel Robert BérardB. Sc., Hon. (Chemistry) McGill University, 1989A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE STUDIES(Department of Chemistry)We accept this thesis as conformingto the required standardTHE UNIVERSITY OF BRITISH COLUMBIADecember 1994© Daniel Robert Bérard, 1994In presenting this thesis in partial fulfilment of the requirements for an advanceddegree at the University of British Columbia, I agree that the Library shall make itfreely available for reference and study. I further agree that permission for extensivecopying of this thesis for scholarly purposes may be granted by the head of mydepartment or by his or her representatives, It is understood that copying orpublication of this thesis for financial gain shall not be allowed without my writtenpermission.Department of________________The University of British ColumbiaVancouver, CanadaDate Dec qDE-6 (2/88)AbstractGeneral reductions of the Ornstein-Zernike equation are given for multicomponentmolecular fluids near an isolated planar surface and between planar surfaces. This allows integral equation approximations such as the hypernetted-chain (HNC) and reference hypernetted-chain (RHNC) theories to be solved numerically. Theoretical methodsfor treating multipolar particles near inert, dielectric and metallic surfaces are considered. The metal is represented using a jellium model together with density-functional(DF) theory and the two phases interact electrostatically. Mean-field theories which reduce the many-body electrostatic wall-particle interactions to effective pair potentialsare described for dielectric and metallic surfaces. The two interfacial phases are solvedself-consistently for the wall-particle distribution function and, in the case of metal surfaces, for the electron density distributions. Explicit results are given for dipolar hardsphere fluids and for electrolyte solutions. The wall-induced fluid “structure”, electrostatic potential drop across the interface and electron density distribution of the metalare discussed in detail. Close to metallic surfaces, it is found that a highly ordered regionexists. 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 thesurface.Hard sphere, Lennard-Jones and dipolar hard sphere fluids are investigated betweeninert walls. Lennard-Jones fluids are also considered between attractive walls. RHNC andHNC results for the fluid structure and the force acting between surfaces are comparedwith computer simulations. With the exception of Lennard-Jones fluids confined betweeninert walls, it is found that the integral equation theories show good agreement with11simulation results. Integral equation and computer simulation results for Lennard-Jonesfluids show better agreement as the distance of the state point from liquid-vapour bulkcoexistence is increased.A detailed investigation of a Lennard-Jones fluid confined between inert walls usingthe grand canonical Monte Carlo method is performed. Capillary evaporation is foundfor liquid subcritical bulk states. General methods are given for simulating a metastablefluid. The force acting between the walls is found to be attractive and increases rapidlyas the spinodal separation is approached. On the equilibrium liquid branch, the netpressure appears significantly more attractive than the van der Waals attraction at smallseparations. This might be the origin of the experimentally measured attraction betweenhydrophobic surfaces in water.111Table of ContentsAbstract iiTable of Contents viList of Tables viiList of Figures ixTable of Symbols xviiAcknowledgement xxii1 Introduction 12 Statistical Mechanics 52.1 Two-body functions and their properties 52.2 Integral equation theories 82.2.1 The Ornstein-Zernike equation 82.2.2 Blum’s reduction of the bulk Ornstein-Zernike equations 122.2.3 The reference hypernetted-chain theory 142.3 Monte Carlo methods 172.4 Models and pair potentials 203 Fluids near planar surfaces 233.1 Introduction 23iv3.2 Reduction of the planar wall-fluid Ornstein-Zernike equation3.3 Mean-field theories for polarizable surfaces3.3.1 Dielectric surfaces3.3.2 Jellium surfaces3.4 Dipolar fluids and electrolyte solutions near inert, dielectric and metallicsurfaces3.4.1 Fluid model3.4.2 Asymptotic behaviour of h(01) near an inert wall3.4.3 Results for a pure dipolar fluid3.4.4 Results for electrolyte solutions4 Fluids between planar surfaces4.1 Introduction4.2 Theory4.2.1 The interaction free energy between two walls separated by a molecular fluid4.2.2 The wall-solvent-wall correlation functionsGrand canonical Monte Carlo method for Lennard-Jones particles . .Results for pure fluids between planar wallsCavitation of a Lennard-Jones fluid between hard walls5 Summary and Conclusion 126Bibliography 133Appendices 14325333341474752567494949696991021061164.34.44.5VA Properties of the Wigner symbols, rotation matrices and rotationalinvariants. 143B The electrostatic interaction potential between a molecule and a uniform surface charge distribution. 147V C Numerical considerations 149viList of Tables3.1 The first few coefficients of the Legendre polynomial expansion, Pi(x) =ax 313.2 Reduced ion diameters d±/d together with the examples and labels usedin this work. The solvent diameter d was taken to be 2.8A 513.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 solventdipole moment and ionic charge are = and q* = 8, respectively. Forthe metal slab, we assume that T 298K and, unless otherwise stated,= 2.65a0 where a0 = 0.529177A is the Bohr radius 523.4 The self-consistent electrostatic potential drop across the jellium-vacuumand jellium-solvent interfaces for z = 32A. When present, the liquidconsists of dipolar hard spheres with p = 0.7, = d 2.8A andT = 298K 743.5 The self-consistent electrostatic potential drop across the jellium-solventinterface for r3 = 2.65a0 and z = 32A. In all cases the liquid consists ofdipolar hard spheres with = at T = 298K 743.6 Contact values of the wall-solvent pair distribution function near an inertsurface 783.7 The potential drop in volts across inert surface-electrolyte solution interfaces. The dielectric constant is also included. In all cases p = 0.7,4=/andq*=8 80vii3.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 r3 = 2.65a0 924.1 Details of the grand canonical Monte Carlo simulations for a LennardJones fluid confined between hard walls. The last column labeled “operations”, denotes the total number of particle moves, insertions and deletionsincluded in the simulation. (The grand canonical Monte Carlo simulationsfor the bulk were performed in a central cell 9u x 9u x 9g.) 1164.2 Details of the grand canonical Monte Carlo simulations for a LennardJones fluid with T* = 1.0 and res/kBT = —3.29 confined between planarhard 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, Ps03 =0.75 + 0.01 and Pbu3/kBT = 0.44 ± 0.05.) 1194.3 Details of the grand canonical Monte Carlo simulations for the second statepoint (T* = 0.9 and res/kBT = —3.84). (The grand canonical MonteCarlo simulation for the bulk, performed in a central cell 7u x 7cr x 7cr,has ex 3.60, Ps3 = 0.79 ± 0.01 and Pbu3/kBT = 0.4 ± 0.1.) 120viiiList of Figures3.1 The coordinate system used in the reduction of the wall-particle OrnsteinZernike equation 273.2 The coordinate system used in Sec. 3.3.1 for a dipolar fluid near a dielectricsurface showing two particle-image pairs. The images are shown to the leftof the interface with a prime included in the labels 353.3 The coordinate system used in Sec. 3.3.2 for a dipolar fluid near a jelliummetal slab 423.4 RHNC results for the projections occuring in the wall-solvent pair distribution function for a dipolar hard sphere fluid with u. = and p 0.7.In (a)-(d), the solid curves represent the present planar wall calculationsand the dotted curves are results obtained previously (Ref. 10) for a macro-sphere of diameter 30d8 573.5 The wall-solvent pair distribution function g5(O, z) for the pure dipolarfluid near an inert wall. The solvent parameters are as in Fig. 3.4 . . . 593.6 The asymptotic behaviour of (a) h° and (b) h2. The dipole momentand 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 60ix3.7 RHNC results for the projections occuring in the wall-solvent pair correlation function for a dipolar hard sphere fluid with= /‘ and p = 0.7.The solid curve represents results for an inert wall ( = 1). The long- andshort-dashed curves were obtained in the cavity approximation for = ccand 5, respectively. The long-dash-dot and short-dash-dot curves are theresults obtained with the superposition approximation for = cc and 5,respectively 623.8 The wall-solvent pair distribution function gws(&, z) for a pure dipolar fluidnear a dielectric surface with,= 5. The results in (a) and (b) are forthe cavity and full superposition approximations, respectively. The modelparameters for the solvent are given in Fig. 3.7 643.9 The wall-solvent pair distribution function g(O, z) for a pure dipolar fluidnear a conducting surface ( = cc). The results in (a) and (b) are forthe cavity and full superposition approximations, respectively. The modelparameters for the solvent are given in Fig. 3.7 653.10 Oscillations in the net potential drop across the jellium-vaciium interfacefor r 2.65a0 663.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 r8 = 2.0, 2.65 and 3.0a0,respectively. In (b), the dash-dot curve shows the self-consistent effectivepotential when no fluid is present and r3 = 2.65a0 The distribution ofelectrons 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 presenceof the fluid (r = 2.65a0) 68x3.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. Thedotted, solid and dashed curves are as in Fig. 3.11. The dash-dot curve isfor an inert surface (r3 oo) 693.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 asin Fig. 3.11 703.14 The wall-solvent pair distribution function g(O, z) near a jellium metalwall with r = 2.65a0 The solvent parameters are as in Fig. 3.7 713.15 The reduced electrostatic potential across the interface when r3 2.65a0(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 ofcontributions from the jellium (dotted curve) and the orientationally polarized fluid (dashed curve). The net potential across the jellium-vacuuminterface is represented by the dash-dot curve. The scale on the rightshows the electrostatic potential çb(z) in volts at 298K. (b) The electrostatic potential for z” > 0 on an expanded scale; curves are the same asin (a) 733.16 Comparison of the projections occuring in the wall-solvent pair correlationfunction for a dipolar hard sphere fluid near a jellium metal (r3 = 2.65a0)surface and a dielectric surface. The short-dash curve is the jellium metalresult. The long-dash and dash-dot curves are the results near a dielectricsurface using the cavity and superposition approximations, respectively.The solid curve is result near an inert surface 75xi3.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 puresolvent result and the dashed curve is for 0.945M CsF 773.18 Wall-ion pair distribution functions for 0.945M electrolyte solutions in contact 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 LiT; (---) Eq in EqEq;(_ —) Cs in CsT. In (b) the curves are: (—) Cl— in NaC1; (———) Eq inEqEq; (.-•—.) F— in CsF; (— —) 1 in CsI 793.19 Projections of the wall-solvent pair distribution function for NaC1 solutions in contact with a jellium slab. The solid curve is the pure solventresult. 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 813.20 Projections of the wall-solvent pair distribution function for RbF solutions in contact with a jellium slab. The solid curve is the pure solventresult. 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 823.21 The wall-solvent pair distribution function g(O, z) for the 0.945M NaC1solution in a dipolar solvent and near a jellium metal wall 833.22 The wall-solvent pair distribution function gws(O, z) for the 7.56M NaC1solution in a dipolar solvent and near a jellium metal wall 843.23 The wall-cation (—) and wall-anion (----) pair distribution functions forseveral 0.945M electrolyte solutions in contact with a jellium slab. Theparameters are as in Tables 3.2 and 3.3. In (h), the wall-cation distributionfunctions for RbF, CsI, CsF and ILi are shown on an expanded scale. . 873.24 An illustration of the distortions of the solvation shell due to the jelliuminduced interfacial polarization 88xii3.25 Wall-ion pair distribution functions for NaC1 solutions at 0.945M. Thedotted, solid and dashed curves are for r8 = 4.0a0, 2.65a0 and 2.0a0,respectively 893.26 The wall-ion pair distribution function for NaC1 solutions near a jelliumsurface. The curves are labeled as in Fig. 3.23 903.27 The wall-ion pair distribution function for RbF solutions near a jelliumsurface. The curves are labeled as in Fig. 3.23 913.28 The self-consistent mean-field electrostatic potential for 0.945M NaC1 incontact 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 thepure solvent. The solid line in (a) is the total electrostatic potential. Theparameters are as in Tables 3.2 and 3.3 934.1 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 dottedcurve is the HNCP result and the dashed curve is the result obtained usingthe 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) 1074.2 The wall-hard sphere distribution function gws(z; 1) for several wall-wallseparations 1. The parameters are as in Fig. 4.1 109xlii4.3 The (a) HNC and (b) RHNC contact density minus the bulk densityas a function of wall-wall separation 1. The crosses denotes values ofpg(d/2; 1)— Ps• The curves are the net HNC and RHNC pressurescalculated from Eq. (4.2.9) 1104.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 reduceddipoles moments of j = 0 (hard spheres), and = respectively. In (b), the HNC asymptotes calculated from Eq. (4.4.5) (smoothcurves) are compared with numerical HNCP results on an expanded scale 1114.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 areas in Fig. 4.4 1124.6 RHNC results for the projection g2 of the wall-particle pair distributionfunction for a dipolar hard sphere fluid confined between hard walls. Results are shown for several wall-wall separations 1. The parameters are asin Fig. 4.4 1134.7 Wall-solvent distribution functions for a Lennard-Jones fluid confined between Lennard-Jones 10—4 walls, Eq. (4.3.3) (T* = 1.21, p = 0.5925 andPsurf2 = 0.8). The wall-wall separations are (a) l = 2cr and (b) 1* = 3cr.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 correlationfunction of perturbation theory [124,125]. The circles are results frommolecular dynamics simulations [108] 114xiv4.8 Wall-solvent distribution functions for a Lennard-Jones fluid confined between hard walls for three state points and two wall-wall separations. Thesolid 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; andin (c) and (d) T* = 1.0 and p = 0.75 ± 0.01 1154.9 The grand canonical Monte Carlo average density (N)/1L2,betweenhard walls for a Lennard-Jones fluid in equilibrium with the bulk. Resultsfor the first state point (res/kBT = —3.29 and T* = 1.0) are shown in(a) and for the second state point (res/kBT = —3.84 and T* = 0.9)in (b). Crosses represent points that are thermodynamically stable andopen squares represent those that are metastable. Open circles are usedwhen the equilibrium and metastable fluids could not be distinguished.Except where illustrated, the standard deviations in 5 are smaller thanthe symbols. The inset shows an expanded view of the gas branch 1184.10 Grand canonical Monte Carlo density distributions for the equilibrium andmetastable 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 1214.11 The net pressure between hard walls. The symbols are the same as inFig. 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 inthe net pressure are smaller than the symbols 122xv4.12 The interaction free energy per unit area measured between hydrophobicsurfaces in water [128]. The solid curve is the conventional van der Waalsattraction, with a Hamaker constant of 2 x 102J, and the dashed curve isthe fitted mean-field result, Eq. (4.5.1) (A = 1.2 x 10’2Jm,lo = l2nm,and a = 4 x 104nm2) 125xviTable of SymbolsSymbol Descriptiona0 Bohr radius = O.529177A.b(12) Bridge function between particles of a and /3.b(12) Reference bridge function. In this thesis, either theHenderson-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 /3.c(r) Projection of the direct correlation function.C(r) x transform ofc.1(r). Also C,(r) when x = 0.e Elementary charge.dc, Diameter of species a.ffll Nonzero constant. Here, fmnl = (2m + 1)(2n + 1).g(l2) Rotationally and translationally invariant pair distributionfunction between particles of a and /3.g(1, 2) 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 distanceg(1; 2) Three-body distribution function defined after Eq. (3.3.5).g(l2) Projection of the pair distribution function.xviih 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 mnl mnHpj,;(r) x transform of h,j,;(r). Also H,;a(r) when x = 0.jn(r) Spherical Bessel function of order ri.kB Boltzmann constant.1 Separation between walls and length of simulation cell alongz axis.L Length of grand canonical Monte Carlo simulation box alongthe x and y axes.m Mass.Average number density of elementary charges in metal fromatomic cores.N Number of particles.N(l) iiww(l) per unit area.pw(0,z) Probability density that a particle of species /3 is at a distancez from a wall and at an angle 0 with respect to the surfacenormal vector.P Pressure. Specifically the net pressure P(l) = P(l)— Pb.Pt Internal pressure. Pressure at a surface due to fluid confinedbetween two surfaces.Pb Bulk fluid pressure.Pj(cos 0) Legendre polynomial of order 1.qa Charge on a particle of species a.xviiiElectric multipole moments of order mq.I of species o in spherical tensor notation and measured in the reference frame ofthe molecule.Electric multipole moments of order mqi of species o in spherical tensor notation and measured in the interparticle reference frame.r General vector. r1 is the vector denoting the position of particle labeled 1; r12 = r2 — r1.Unit vector i = r/r. When used as an argument for theWigner spherical harmonics, i = (q, 0, 0) denotes the Eulerangles.r3 Wigner-Seitz radius. r3 = (3/47rn)’/.Wigner generalized spherical harmonic.S(x) Polynomial. See page 30T 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).xixy y = 47rps/9kBTY(O, q) Spherical harmonic.z With reference to the partition functions and Monte Carlosimulations, z is the activity. Otherwise, z denotes a positionalong 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 =kBT in zu3 where z is the activity.Dipole moment of species a.E(, V, T) Grand canonical partition function.Number of species composing system.p Average bulk number density of species a.p(’)(1, 2,. ..,n) n-particle density.Pth Threshold density defined on page 104.Lennard-Jones collision parameter.xxZ2,112) Translational invariant basis function.One of the azimuthal Euler angles.cb(z) Scalar potential.x One of the azimuthal Euler angles.2,r12) 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 symbolsUsed to denote —m.h Fourier or Hankel transform of h.() Grand canonical ensemble average.xxiAcknowledgementIt is my pleasure to thank my academic supervisor, Prof. G. N. Patey, whose hasgiven me guidance, encouragement and independence in my studies over the last fiveyears.I would also like to thank the Natural Sciences and Engineering Research Council(NSERC) of Canada and the University of British Columbia Chemistry Department forfinancial support and Dr. Phil Attard, Dr. Masahiro Kinoshita, Dr. Glen Torrie, Dr.Dongqing Wei and Dr. Amalendu Chandra for numerous interesting discussion pertainingto my work.xxiiChapter 1IntroductionThe solid-liquid interface has received a great deal of interest in the past few decadesfrom both theoreticians and experimentalists. The principle reason for this interest isthat the properties of the individual phases at the interface can be significantly differentfrom those of the bulk. Furthermore, the interfacial properties depend on the nature ofboth phases and the properties of a liquid, for instance, can vary significantly at differentsurfaces. A familiar example at the macroscopic level are the wetting properties of waterat hydrophobic and hydrophilic surfaces. However, fundamental to an understanding ofany physical system is knowledge of the underlying microscopic “structure”.Statistical mechanics provides a means of calculating the structure and other observable properties of a given system based only upon a description of the microscopicinteractions and the state parameters. One of the principle tasks in any statistical mechanical theory is the description and development of microscopic models and two generalapproaches are often followed. One of these employs potentials that have been fitted todescribe the particle interactions in specific physical systems. Such potentials are oftenrestricted to a particular state of the system and may only model a small set of propertiesaccurately. Alternatively, simple model potentials which represent certain intrinsic properties of molecular interactions can be employed. Such an approach is directed towardsunderstanding the relationship between general features of intermolecular interactionsand observable properties. The simplest model potential which can reproduce an observable property often yields the most insight into the nature of a system.1Chapter 1. Introduction 2Particle distribution functions, in particular the singlet and pair distributions, are ofcentral importance to an understanding of the microscopic structure of fluids and solids.The singlet distribution function describes the probability density of a particle being ata given position and orientation. The pair distribution function describes the probabilitydensity of two particles simultaneously having given sets of positional and orientationalcoordinates. If we assume that the potential energy can be written as the sum of pairwisepotentials, all thermodynamical properties of a system can be calculated from the singletand pair distribution functions.Integral equation methods provide one approach for calculating pair distribution functions. Methods based on the Ornstein-Zernike equation are commonly used in liquid statestudies and will be employed here. The Ornstein-Zernike equation relates the pair distribution g(12) to the direct correlation function c(12) and may be considered as a definitionof c(12). In order to employ the Ornstein-Zernike equation, a second closure relation between g(12) and c(12) is necessary. An exact closure can be formally written but thisintroduces a third unknown function b(12). The Ornstein-Zernike equation together withan approximate closure form an integral equation theory. Approximate closures can beviewed as different approximations for the unknown function b(12).Computer simulations provide an alternative approach to the liquid state. In orderto understand the influence of a surface on a fluid, it is particularly useful to comparethe fluid properties at the interface with those of the bulk. Grand canonical MonteCarlo calculations are especially useful in this regard. Monte Carlo methods are basedon sampling system configurations. Starting with an initial configuration of particles,additional configurations are generated by translating, rotating, inserting and removingparticles according to their relative Boltzmann probabilities. The equilibrium distribution functions and other properties of the system are calculated as ensemble averages.The only inputs required by Monte Carlo and other computer simulation methods areChapter 1. Introduction 3the interaction potentials and state parameters. Hence, for a given model, computersimulation methods yield essentially exact results.In this thesis, the liquid side of the interface is of primary concern and we will studythe 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 tunnel-junction devices [1] and x-ray diffraction measurements [2]. A few limited results havebeen obtained for water at mercury and silver surfaces.There has been a considerable amount of theoretical work focussed on inhomogeneoussystems. These include dense fluids near simple surfaces [3] and electron density distributions at metallic surfaces [4]. However, these methods emphasize a detailed modeland theory of one phase and simplify the other phases enormously. For example, metalsurfaces have been studied primarily in vacuum or in dilute gases (emphasizing adsorption phenomena) [4]. When the liquid phase is of interest, detailed classical liquid statemodels and theories have been employed at inert and simple attractive walls.This thesis is divided into five chapters. Chapter 2 is devoted to the fundamentalsof liquid state theory and begins with a discussion of the properties of pair functions inhomogeneous and isotropic systems. The integral equation theories and grand canonicalMonte Carlo method are presented and the chapter concludes with a discussion of pairpotentials. Chapter 3 is devoted to molecular fluids near isolated planar surfaces andcontains an introduction to the relevant literature. A general reduction of the OrnsteinZernike equation for fluids near planar surfaces is presented and theories for molecularfluids near inert, dielectric and metallic surfaces are derived. Simple dipolar hard spherefluids and electrolyte solutions in contact with these surfaces are considered. The asymptotic behaviour of the pair distribution function for the pure dipolar fluid at long range isderived and compared with analytical continuum and exact asymptotic profiles. ChapterChapter 1. Introduction 44 is concerned with the structure of fluids confined between surfaces. General reductionsof the Ornstein-Zernike equation describing the pair correlation functions in these systems are presented. The grand canonical Monte Carlo method used is described in detailand a general method of simulating metastable states with Monte Carlo is included. Avariety of model fluids are studied between inert and attractive surfaces. The integralequation results are compared with simulations and a discussion concerning the applicability of integral equation methods to different interfacial models is given. A detailedMonte 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. Thepossible relevance of the Monte Carlo results to the experimentally measured attractionbetween hydrophobic surfaces is discussed. Finally, a summary of conclusions is given inChapter 5.Chapter 2Statistical MechanicsIn this chapter, the basic equilibrium statistical mechanical tools and theories usedin this thesis will be outlined.2.1 Two-body functions and their propertiesIn classical statistical mechanics, the grand canonical n-particle density is defined asNp(1, 2,. . . , n)=(N — n)! f e_U(12NTd n + 1) d(n +2). d(N) , (2.1.1)where the arguments 1, 2, ..., N denote the coordinates of particles labeled 1, 2,N, z is the activity, U(1, 2,. . . , N) is the N-particle potential, kBT is the Boltzmannconstant times the temperature and is the grand canonical partition function. Thepair distribution function will be of fundamental importance and is defined in terms ofthe pair density p(2)(1, 2) as(2)(1 2)g(1,2)p(’)(1)p(’)(2) (2.1.2)For a homogeneous and isotropic system of spherically symmetric particles, the pair distribution is simply a function of the interparticle separation r12 = r2 — r1 and describesthe probability density of finding a particle 2 at distance r12 from particle 1. Conversely,for a molecular system, g(1, 2) can be a function of as many as nine independent coordinates in a translationally invariant environment and twelve in an inhomogeneousenvironment. In order to work with these functions, it is necessary to first introduce themathematical tools needed to reduce them to a tractable form.5Chapter 2. Statistical Mechanics 6The translationally invariant pair function for a molecular system can be written asg(1,2) = g(1Z1,2r2),where 1 (qj,Oj,xj) represents the set of Euler angles denoting the orientation of molecule j and r12 = r2 — r1. In this form, the two-body functioncan be expanded in an orthogonal basis set which spans the anisotropic orientationalspace of the two particles. Specifically,g(1,2) = g(Z1jl2,r2)= > g,)1(ri2) T ,I(11,12,i12) , (2.1.3a)mnlwhere r12 = Iri2, i12 =r12/r and the basis functions areT),(1Z1,cz2,112) = ffl1R(1)R,(fZ2)R,0(i12) . (2.1.3b)Here, R(fZ) is a Wigner generalized spherical harmonic [5] and fmnl is an arbitrarynonzero constant chosen for mathematical convenience. When used as an argument of theWigner spherical harmonics, i denotes the Euler angles of r, (q, Or, 0). The eight indicesin Eq. (2.1.3) take on all integer values for which the Wigner generalized harmonics aredefined. That is,m,n,l = 0,1,2,... (2.l.4a)—m<<m and —m’<m (2.1.4b)—n v n and —n i/ n (2.1.4c)—l )‘ <1 . (2.1.4d)Provided the expansion is rapidly convergent and the sum can be truncated for relativelysmall values of m, n and 1, the one-dimensional coefficients are relatively easyto calculate.In this work, only translationally and rotationally invariant pair functions will be considered and, in general, these can be described by six (or fewer) independent coordinates.Chapter 2. Statistical Mechanics 7In 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 ofthe molecular pair (while keep the relative orientation and interparticle spacing fixed)leaves g(1, 2) unchanged. Mathematically, the concerted rotation is achieved by applying a rotational operator R(f) to each particle and to the interparticle vector 112 inEq. (2.1.3) [6]. This givesg(12) fmnlgnn,,(ri2) R,,(1Z)R1(mnl z)’xR11() RIA,l(f)R0(i,2 (2.1.5)where Eq. (A.13) has been used’. Integrating over all orientations and using Eq. (A.1O)yieldsg(12) = g’(r,2)mnl(f,,2,f,), (2.1.6a)mnlwhere the so-called rotational invariants are defined asg(r,2) = gsl(r,2), (2.1.6b)/I,V,)mnl( 2, r,2) = ( R(,)R(2R,0(,2) (2.1.6c)?I,v,), I) is the usual 3-j symbol [5] and, in order to emphasize that the pair functionsdepend only on the interparticle coordinates, the commas have been removed from thearguments of g(12). The contracted basis functions have five indices (one for each degreeof rotational freedom). These indices can take on the values given in Eq. (2.1.4); however,the 3-j symbols are nonzero only ifIm — iI 1 m + n. (2.1.7)properties of the Wigner functions are summarized in Appendix A.Chapter 2. Statistical Mechanics 8The angle-averaged distribution function is simply g°(r).The expansion (2.1.6) provides a general mechanism for describing two-body functions which only imposes the constraints that the functions are invariant to rotationsand translations. The pair functions of interest for molecular fluids must satisfy a fewadditional constraints [6]. Specifically, the two-body functions of interest here are alsoreal, invariant to the interchange of particles and invariant under the symmetry operations of molecules 1 and 2. As a result of these properties, many of the coefficientsbecome degenerate or vanish, and the number of unique coefficients is reduced. The firstconditions that the functions are real impliesrnnl*(— ( ‘m+n+1+#+v mnl(r’ 2 1g r1—\ ,where i = —ii, and is obtained by requiring that the product of two real rotationally andtranslationally invariant functions is real. The second conditions can be expressed asg(12) = g(21) . (2.1.9)Using this, the parity of the generalized spherical harmonics R0(i) = (_)lR0(—i) andthe symmetry of the 3-j symbols (see Appendix A) gives the identityg(r) = (_)m+ng1l(r) . (2.1.10)The molecular symmetry properties further reduce the number of unique coefficients andidentify those which must be zero (see Ref. [6] for a discussion). This reduction will bediscussed in Chapter 3 for dipolar particles.2.2 Integral equation theories2.2.1 The Ornstein-Zernike equationThe Ornstein-Zernike (OZ) equation is the basis for many of the integral equation theoriesof fluids and was first introduced in 1914 [7] to study density fluctuations in fluids nearChapter 2. Statistical Mechanics 9the critical point. The OZ equation relates the pair and direct correlation functions, hand c, respectively, through a convolution and is given in its most general form for aone-component molecular fluid ash(1,2) c(1,2)+--Jp(’)(3)h(1,3)c(3,2)d(3), (2.2.1)where p(’)(3) is the singlet density, h(1, 2) = g(1, 2) — 1 and fd(3) = fdr3d3 denotesthe integral over all positions and orientations of particle labeled 3. The OZ equation isoften viewed as the defining relation for the direct correlation function; however, c canalso 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 athand. Using this definition, and starting from the grand canonical ensemble of classicalstatistical mechanics, the OZ equation can be derived if we assume that all particlesinteract through a pairwise additive potential such that [3]NU(1,2,. ..,N) = u(i,j) . (2.2.3)i<jThe assumption that the fluid behaves according to classical mechanics implies thatA =2rh<<a, (2.2.4a)N rnkBTwhere A is the thermal de Broglie wavelength, h is Planck’s constant divided by 27r, mis the particle mass and a represents the average nearest-neighbour separation betweenparticles. For molecular fluids, the additional constraint thath2Orot = <<T (2.2.4b)21kBis required, where °rot is the characteristic rotation temperature and I is the molecularmoment of inertia. These conditions are satisfied for typical bulk fluids and, in particular,Chapter 2. Statistical Mechanics 10are 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 densityp(1). Eq. (2.2.1) is usually referred to as the inhomogeneous OZ equation since thecorrelation functions are solved in an external field.When the two-body correlation functions are invariant to translations and rotationsand the functional derivative in Eq. (2.2.2) is evaluated at q(l) = 0, the OZ equationhas the simpler form2h(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 commas have again been removed from the function variables to emphasize that they dependonly on the interparticle coordinates. Equation (2.2.5) is referred to as the homogeneousOZ equation and can not be solved for both h and c until a second closure relation isintroduced. The inhomogeneous OZ equation requires two additional relationships. Todate, many approximate closures have been proposed and some of these will be consideredin Sec. 2.2.3.The designation “homogeneous” and “inhomogeneous” are somewhat inappropriateand have apparently led to some confusion in the literature concerning the applicabilityof Eq. (2.2.5) [8,9]. Hence, it is useful to define these terms in the present context. Sincethe 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 ashomogeneous (or pairwise homogeneous) if the interactions and correlations between apair of particles depends only on the interparticle coordinates (i.e., translationallyand rotationally invariant). In this way, homogeneous also implies that the systemis pairwise isotropic and is valid whenever q(l) = 0; and2The homogeneous OZ Eq. (2.2.5) is also be obtained if particle 1 is assumed to be the source of q(l).Chapter 2. Statistical Mechanics 11inhomogeneous if the two-body correlations are translational and/or rotational variants.The external field in the inhomogeneous equation can often be defined in terms of the fieldof one or more particles. In this case, the distinction between a pairwise homogeneous anda pairwise inhomogeneous system is academic. In practice, the inhomogeneous equationnot only provides more information about the microscopic system structure (effectivelyas the three-body correlation functions), but is typically more accurate for a given closureapproximation. However, it is worthwhile to note that the inhomogeneous theories basedon Eq. (2.2.1) are significantly more complicated and time-consuming to solve. For thesystems that will be considered in this thesis, all species are modeled explicitly and thesimpler homogeneous Eq. (2.2.5) is used.Before proceeding, we generalize the homogeneous OZ equation to a multicomponentsystem of + 1 species to obtainh(12) = ca(12) + f ha(13) c(32) d(3), (a, = 0,1,. . . , ), (2.2.6)Chapter 2. Statistical Mechanics 12where p is the average bulk number density of species 7. In this work, the species labeledw, for wall, is considered to be infinitely dilute (Pw = 0) and will be used to represent asurface species. As a result, the coupled set of Eqs (2.2.6) is reduced to four smaller setsh(12) — c(12) Jh7(13) c(32) d(3), (2.2.7a)h(12) — c(12) = J h(13) c(32) d(3), (2.2.7b)haw(12) — c(12) = 2 f h(13) c(32) d(3), (2.2.7c)h(12) — c(12) = f2 f h(13) c(32) d(3) , (2.2.7d)where the labels (cr, /3 = 1, 2,.. .,now exclude the dilute species w. The first set ofEqs (2.2.7a) describe the correlation functions for the bulk fluid in the absence of a walland 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 remainingspecies 1, 2,. . . , c—and only one is required. The last equation describes the correlationsbetween two dilute w particles. The solution of Eq. (2.2.7a) is required as input toEqs (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 equationsAlthough we are primarily interested in the solution to Eqs (2.2.7b)—(2.2.7d), we requirethe bulk pair correlation functions before any of these equations can be solved. Hence, webriefly discuss the solution of Eq. (2.2.7a) for nonspherical particles. As we have seen, therotationally and translationally invariant two-body functions vary with six independentvariables and it is obviously not realistic to solve Eq. (2.2.7a) in its present form. Thereduction of the OZ equation for pairwise homogeneous and isotropic fluids to a set ofChapter 2. Statistical Mechanics 13coupled 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 reviewis given here.Since the OZ equation is in the form of a convolution, the spatial integral can beeliminated by taking the three-dimensional Fourier transform. This gives7a(12)—cyp(12)=J cv() y(32)d3, (2.2.8)where the Fourier transforms are defined asiai3(Zi, 2, k)= J e’’ ha(i, 12, r) dr , (2.2.9a)= feirc&(i,1z2,r)dr, (2.2.9b)and k is a reciprocal space vector. The k-space correlation functions can be expanded ina rotational invariant basis set yieldingJia(12)= E k) ‘(fZ1,2, 1(12) , (2 2 lOa)mnla(12) = k)mnl() (22 lOb)mnlj.LVThe coefficients and are related to the real space coefficients by the Hankeltransforms defined ask) = 47ri1 j2jj(kr)h(r)dr , (2.2.11)where ji(kr) is the spherical Bessel function of order 1. The Hankel transforms can beevaluated using Fourier transform methods [10,6,17,19]. Substituting these expansionsinto Eq. (2.2.8) and using properties of the rotational invariants (see Appendix A), thebulk OZ equation reduces to relations of the form=(2.2.12)-y 12flChapter 2. Statistical Mechanics 14where = h — c,‘712111 (i)m+n+ni 21 + 1 fmnilifninl2 f 12 11 l 1 /1211 l’\Limnn 2n1 + 1 fmnl mn nif 0 o oJ (2.2.13a)— 2l + 1 fmnhllfnlnl2— 2n1+1 f1mn1l’\ /n1rjl2’ /mnl”x ( (2.2.13b)= —x and { fl n } is a 6-j symbol [5]. If we use the definitionfmnl= /(2m + 1)(2n + 1), (2.2.14)referred to henceforth as Blum’s notation, and define the so-called x transforms as1,nnl)-mnlymn;x (k’ = c,(k) , (2.2.15a)iw;a/3 //mnl)mnlk) = (x o , (2.2.15b)then Eq. (2.2.12) simplifies considerably to giveBlumfmn;x(j ()x+iii Ffmn1;xfk+mn1;x(k] nin;x 1k (2.2.16)I. 1wl;a7 I wi;ay\ j 1i’;8 )ILv;a/3\ /-y 1’1Blumwhere = denotes the use of Eq. (2.2.14). These equations can also be convenientlyassembled into matrix form as shown in Refs [17], [18] and [10].2.2.3 The reference hypernetted-chain theoryThe OZ equation (2.2.5) relates the unknown pair and direct correlation functions. Inorder to make use of the relationship, a second closure condition relating h and c isrequired. Such a relation was simultaneously developed by several authors [20—24] around1960 and can be derived as a functional Taylor series expansion of Eq. (2.2.2) [3,25] togiveu(12)— b ‘12), (2.2.17)h(12) — c(12) = in [h(12) + + kBTChapter 2. Statistical Mechanics 15where u is the pair potential between particles of species c and 3. The function bafirepresents the bridge (or elementary) diagrams and consist of an infinite series of multi-centered integrals. With the exception of the first two- and three-centered integrals (toorder p2 and p3 in density) [26], the bridge functions are too complicated to evaluate.They are usually considered to be short ranged; decaying at long range as h2. However,for fluids near macroscopic surfaces it is expected that b may become very importantin certain cases (for example, as the equilibrium bulk phase approaches the coexistencecurve).In order to obtain a useful integral equation theory, an approximation of the exact relation (2.2.17) must be constructed. The approximation bai = 0 is called thehypernetted-chain (HNC) theory since the resulting expression for c consists of diagramscalled simple chains, netted chains and bundles. Another approximation is to replace b.owith those of a simpler reference system b, such as a hard sphere fluid. The resultingclosureha(12) — c(12) = ln [ha(12) + + u(12) — b(12), (2.2.18)is called the reference hypernetted-chain (RHNC) approximation. Since the bridge diagrams for hard spheres are also difficult to evaluate, several indirect methods of estimatingthe bulk hard sphere-hard sphere [26—30] and the hard wall-hard sphere [31,26,32] functions have been developed. Most of these methods are based on analytically fitting thepair distribution function. The bridge function is then calculated from Eq. (2.2.18) usingthe OZ definition for c(12). In the method of Attard and Patey [26], the hard spherebridge diagrams are calculated explicitly to second, b2(r), and third, b(r), order indensity. The bridge function is then approximated by the Padé formb(r) Ps WS\ / (2.2.19)1 — p5b(r)/bws (r)Chapter 2. Statistical Mechanics 16This equation is referred to as the HNCP approximation. Although the resulting bridgefunctions from the various methods can differ significantly, most of these methods givevery good estimates of the hard sphere-hard sphere and hard wall-hard sphere pair distribution functions. Unless otherwise stated, we use the Verlet-Weis [27] method generalizedto multicomponent fluids by Lee and Levesque [28] for the bulk bridge functions andeither the Henderson-Plischke fit [31] or the HNCP approximation for the wall-particlefunctions b.In order to make use of a closure for nonspherical particles, Eq. (2.2.17) must beexpanded in a rotationally invariant basis set. One method of achieving this is to rearrange the closure equation and take the partial derivative with respect to the interparticleseparation r12 [10]. After rearranging the result, this givesOca(12) = hafl(12)Owa,3(1 ) 18a,(12) +3(12), (2.2.20)Or12 kBT Or12 kBT Or12 Or12wherewa(12) = —kBTlng,3(12)= —kBTi(12) + ua(12) — kBTba(12) (2.2.21)is the potential of mean force. Integrating from an arbitrary radius r to oc and assumingthat ca(12) * ua(12)/kBT —* 0 as r12 —* oo (as is usually the case), one finds that1 Ow (12) (12)c(12)= kBT 112ha(12)‘rdr—+ ba(12) . (2.2.22)Expanding the two-body functions in rotational invariants using Eq. (2.1.6) and reducingyieldspmnl oc OWm2n2t2 (r’)c(r)=jj h’(r’) dr’m2flt.L2mnl ( ‘up,;af3r1 i jmn1— 1 m Y v;a/3T) ,Chapter 2. Statistical Mechanics 17wherefmn l fm2n212=(2m + 1)(2n + 1)(2l + 1)(_)m+n+1+1+1+2+M2ni1 n1 ‘1(m1 m2 m (ni n2 n (l 12 1X Tfl2 2 12 1 - 11 - 11 1 , (2.2.23b)\IL1 I2 /1) \i’1 V2 VI \000Jmnland131 J 3334 35 3618 39is a 9-j symbol [5]. Since the reference fluid consists of spherical particles in all currentapplications 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 coordinates have also been derived [33,34]. In this approach, the assumption that ca(12) _+ 0as i’12 —* oo is not necessary.2.3 Monte Carlo methodsThe grand canonical Monte Carlo method provides an alternative approach for liquidsin bulk or confined systems. The method basically consists of sampling configurationalphase space points occurring in the partition functionE(, V, T)=Z J -U(1,2 N)/kuTd(l) d(2)... d(N). (2.3.1)The distribution functions and observables,co 2N A(1 2(A)= N!f E(t,’çT) d(1)d(2)... d(N), (2.3.2)together with the associated statistical error, are then calculated as ensemble averages. Inprinciple, the Monte Carlo method is exact for a given model system and better estimatesof the observables can be obtained by sampling more phase space points.Chapter 2. Statistical Mechanics 18For nearly all systems of interest (N) c’-’ Avogadro’s number NA and a method ofreducing the number of independent particles N is necessary. The standard method ofdoing this while still considering a large system size is to introduce periodic boundaryconditions. In this approach, a given configuration of N particles in a finite volumeV is replicated periodically to form an infinite lattice. The N particles in the centralcell are free to translate over all space provided their replicate “images” follow identicaltrajectories.Even when (N) is relatively small, direct sampling of phase space is inefficient sincemost phase points make little or no contribution to the integral. In an alternative technique called importance sampling, phase points are chosen from a weighted distribution.The Metropolis method implements this technique by starting with an initial configuration and follows a Markov chain with a limiting distribution of states given byz V N e_U(1,2,N)PBTPiVT= N! E(t,V,T) ‘for spherical particles. The Metropolis method in the grand canonical ensemble consistsof a series of particle translations, insertions and deletions. In a translation from coil-figuration m to n, a particle is chosen at random and a new set of test coordinates areselected from within a volume dxdydz centered at the old position. The move is acceptedif the new configuration is lower in energy Un — Urn <0; otherwise, it is accepted with aprobabilityN N_Ufl/kBT/N!’= z e / .= e_(_Jm)T < 1 . (2.3.4)zNe_Um/IcBT/N!’’ —rm / .In either case, the resulting configuration contributes to the ensemble average. In adeletion step, a particle is again chosen at random. In an insertion step, a new setof 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 newChapter 2. Statistical Mechanics 19configuration of N — 1 or N + 1 molecules is accepted; otherwise, it is accepted with aprobabilityN NPn= C(UnUm)/kBT, (2.3.5)Pm zVfor a deletion andN— Z Ve__Um)T (2 3 6)— N+1 ‘for an addition.The principle task in a Monte Carlo simulation is the continual calculation of thepotential energy U in the Markov chain of configurations. When periodic boundary conditions are used, the potential is usually partitioned into short and long ranged contributions. The short-range contribution involves direct summation of the pair interactionsusing the minimum image convention where a particle i interacts with the nearest imageof particle j. In addition, the short-range contributions to U is typically truncated usinga spherical cut-off. The longer ranged contributions from particles beyond this truncationsphere can be included using various mean-field techniques [35].The grand canonical Monte Carlo calculation is usually carried out at constant reduced activity z3 = ehh/Tcr3/A ( has units of length and is an intrinsic propertyof the molecular model) instead of constant chemical potential . This is equivalent toholding it fixed but has the advantage that the deBroglie wavelength A does not have tobe defined. Using z, the excess chemical potential is given byitex = ititid= itres— kBTln3 , (2.3.7)where it’ = kBT ln za3 is often called the residual or configurational chemical potential,itid = kBT ln+ 3 in A is the chemical potential for an ideal gas and = (N)/V. AnotherChapter 2. Statistical Mechanics 20method, due to Adams [36,37,35], can also be used. Here the simulation is performed forconstant B, V and T whereB=lnzV= ex/kT+l_V (2.3.8)This method also has the advantage that the deBroglie wavelength does not have to bedefined. However, even when the chemical potential and activity are constant, B is afunction of V, the dimensions of the simulation cell.The Monte Carlo method has the advantage that it can provide an exact solutionfor a given model provided that a sufficient number of statistical sample configurationshave been averaged. However, Monte Carlo simulations can require several orders ofmagnitude more computer time than the integral equation theories.2.4 Models and pair potentialsIn order to use the integral equation theories, Monte Carlo method or any other statistical mechanical theory, it is necessary at some point to consider the interaction potential.Fundamentally, the potential is governed by the dynamic electronic and nuclear structure of the constituent molecules. At the pair level, the potential can be partitionedartificially into two basic components including a short-range repulsive term, US, arisingfrom electron and nuclear repulsion as the molecular charge distributions overlap, and along-range electrostatic term, el, to getup(l2) = u(12) + uj(l2) . (2.4.1)The short-range potential is the least-well characterized contribution to the pair potential and either fitted potentials or simple generic models must be used to represent it.Chapter 2. Statistical Mechanics 21Here, we are not interested in modeling specific systems and will employ some simplifiedapproximations for u including the hard-sphere potentialHS I oo, r<(d+d)/2,u(r) (2.4.2)(0, r>(da+d)/2,and the algebraic inverse power-law potentialzt(r) cx -— , (2.4.3)where n is typically chosen between 9 and 12. Although both forms are chosen in part formathematical convenience, they provide a reasonable approximation of us’ for “classicalfluids” and give good results for many properties [3,11—16,38].The electrostatic pair interaction between particles of species a and 9 in electrostaticunits (esu) is given byu(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 speciesa. In a point charge model, na(r) is simplyn(r) = (2.4.5)where j is the number of point charges q at r of species a. Equation (2.4.4) can beconveniently arranged in the form of a rotational invariant expansion [17],u(12) = u(rl2) (1M2,i2), (2.4.6a)mnlin,mnl;el (_)m6m÷ni (21 + 1)! Q;c Qi;1312— \j (2m)!(2n)! r1 2.4.appropriate for use in the integral equation theories. Here represents the electricmultipole moments of species a in spherical tensor notation and in the reference frameChapter 2. Statistical Mechanics 22of the molecule. For a rigid molecular model, the moments are constant and givenby= eJrma(r)R(I)dr, (2.4.7)where ia(r’) is the charge distribution measured in the molecular reference frame. Themultipole moments in the space-fixed, Q,, and molecule-fixed, reference framesare related by a rotational transformation where= R(f1)Q,, (2.4.8)andQ(a) = eJrmn (r)Rm (I)dr . (2.4.9)The average contribution to the pair potential from molecular polarizability can be included in a mean-field dipole-dipole contribution to the electrostatic potential. Whenthe molecules have no permanent multipole moments, the Lennard-Jones potential,u’(r) = 4e[(a)12 — (cT)6j, (2.4.10)where e is the depth of the potential well and a, the collision parameter, is the particle separation when uLJ(r) = 0, provides a reasonably good approximation of the pairinteraction.Chapter 3Fluids near planar surfaces3.1 IntroductionAs discussed in Chapter 1, the solid-liquid interface is of importance in many physical systems and has received a considerable amount of theoretical attention in recentyears [39—79]. The objective of much of this work is to develop a completely moleculardescription of the interface and, towards this goal, two fundamental approaches arisenaturally. One of these approaches uses quantum density-functional theory to study themetal-vacuum interface [80]. The second is based on classical statistical mechanical methods including integral equation theories [39—53] and computer simulations [54—61]. Thestatistical mechanical methods have been used to study liquids and electrolyte solutionsnear simple non-metallic surfaces. During the last decade, additional advances have beenmade towards combining the density-functional and statistical mechanical approaches,and developing more realistic descriptions of the metal-liquid interface. However, with theexception of work based on the formalism described in Sec. 3.2 [52,81], recent theoreticalmethods of the liquid phase structure near planar surfaces have been limited to dipolarsolvent models. Furthermore, these calculations have been confined to the mean spherical approximation (MSA) or linearized hypernetted-chain (LHNC) theory and have beensolved only near inert surfaces. The more robust RHNC theory has never been solved formolecular fluids near planar surfaces. Hence, before considering model potentials for thedielectric-fluid and metal-fluid interfaces, the reduction of the OZ theory will be given23Chapter 3. Fluids near planar surfaces 24for molecular fluids near planar surfaces.Among work dedicated to the study of interfacial models, there has been some contribution treating the surface as a dielectric continuum [60,61,9,64—66]. However, thesestudies have been limited to computer simulations of water-like fluids [60,61] and to fluidsof spherical ions [9]. A formal analysis of the long-range asymptotic profile of dipolarfluids near dielectric surfaces has also been given [62,63]. In addition, there has been aconsiderable amount of work devoted to the development of theories for the metal-liquidinterface [64—79]. In these studies, particular attention was given to the metal-electrolytesolution interface and the importance of explicitly including the metal side of the interface was clearly demonstrated. In these calculations, the metal was modeled as anelectron gas in a uniform background of positive charge. This so-called “jellium” modelwas then solved at varying levels of approximation using the general approach of theHohenberg-Kohn-Sham density-functional (DF) formalism [82,83]. However, we notethat the self-consistent method developed for the jellium-vacuum interface [80] was onlyapplied 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 liquid side at a molecular level, they did not properly couple the solution structure at theinterface to the electronic structure of the metal. That is, while the solvent and ion fieldswere allowed to influence the electron density distribution at the metal surface, the liquidstructure itself could not react to the field generated by the electron distribution in themetal.For polarizable interfaces, the principal problem hindering the development of fullyself-consistent theoretical methods arises from the fact that the forces involved are fundamentally many-bodied in nature. For example, the electrostatic interaction of a particlein solution with a polarizable surface depends in principle upon the coordinates of allChapter 3. Fluids near planar surfaces 25other particles in both the solid and liquid phases. Unfortunately, at present most statistical mechanical theories applicable to interfaces are only useful if formulated at thelevel of the translationally and rotationally invariant pair potential.The purpose of this chapter is to propose a nontrivial practical solution to this problem. We develop self-consistent mean-field theories for the dielectric-liquid and metal-liquid interfaces which reduce the many-body wall-particle interaction to an effective pairpotential 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 fora 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 bydifferent dielectric constants.The remainder of this chapter is divided into four sections. In Sec. 3.2, the reduction of the OZ equation is given for multicomponent molecular fluids near a planarwall. Sec. 3.3 is devoted to the wall-particle potentials and we begin by developing ageneral self-consistent mean-field theory for the dielectric wall-particle interaction. Aself-consistent mean-field theory is then presented for the metal-particle interaction. InSec 3.4, results are presented for pure dipolar fluids near inert, dielectric and metallicsurfaces. Results for model electrolyte solutions near inert and metallic surfaces arediscussed in Sec. 3.5.3.2 Reduction of the planar wall-fluid OZ equationThe Ornstein-Zernike equation for a multicomponent fluid was introduced in Sec. 2.2and the reduction of the bulk equations was briefly outlined. For molecular fluids incontact with a planar surface, the reduction of either Eqs (2.2.7b) or (2.2.7c) is required.Chapter 3. Fluids near planar surfaces 26If Eq. (2.2.7c) is written ash(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 clearthat the reduction of either equation is in principle equivalent. However, in practice froma numerical solution perspective, there are small differences in the reduction as discussedby Kinoshita and Harada [51] using the present formulation. In this section the reductionof 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 werefirst considered by Perram and White [84] and by Henderson et al. [85,86]. In theirpioneering work, the planar surface is obtained in the limit that a spherical particlebecomes infinitely large. Subsequently, Sullivan and Stell [87] provided a fast and efficientFourier transform method for solving the resulting wall-particle equation. Since thesurface is uniform, it has no orientational dependence and the rotationally invariantexpansion (2.1.6) reduces tog(01) = (3.2.la)= gw(fZi,roi)Chapter 3. Fluids near planar surfaces 27Figure 3.1: The coordinate system used in the reduction of the wall-particle Ornstein-Zernike equation.where we have used 0 to denotes the wall coordinates. The wall-particle correlation functions will be considered exclusively in a reference frame directed along the interparticlevector where r01 = z1 and z1 = 0 at the surface (see Fig. 3.1). Since corresponds tothe Euler angles ioi = (0,0,0), R,0(O) = &‘o and the basis set reduces considerably to‘n Onn, 0)= +R(q, 0, x) (3.2.lb)(_)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 spherical harmonics Y(x, 0) [5].2r21=(r2i,,4)1z1Chapter 3. Fluids near planar surfaces 28Starting with Eq. (2.2.7b) and expanding the correlation functions in rotational in-variants yields( )nfOnnOnn /77W/3(01)=(3.2.2)=( )fli fOni ni fm2 2 12 R(1Z)=1 i”i mnl /(2ni + 1)12 ‘2Im2 n2 12) j R, (12)R(1d1122 v2 2x I h° m2n1j . (3.2.3)Since the coefficients of the correlation functions depend only on interparticle spacing,the definition R,0(i’) = eiA’d,o(O) is used to reduce the spatial integral toJ h.(z2) m2n212 /C;.yi3iT21)R0(i21dr2— Jh0 / \ m2n1 I \ 112— ovi;w2)c7T2l)ao(O l ’2 , (3.2.4)where 921 is the polar angle of r21 (see Fig. 3.1). Subsequent application of the orthogonalproperties of the Wigner functions (A.9) yields relations of the formo;wl) =(_)fl+fla+i (2n + 1) In1 nOnn / Blum v’(2n1 + 1) o o o)x fh7(z2)C(r21) Pi(cos 0) dr2, (3.2.5)where Pi(cos 0) d0(0) is the Legendre polynomial of order 1 and we have used Blum’snotation [17] for the [See Eq. (2.2.14)].We now take the Fourier Transform of Eq. (3.2.5). However, since (or c) doesnot vanish in the limit z —* —oo, its transform diverges. This problem is easily solvedfor a hard wall interaction whereh,(Z,z) —1, z < d/2 , (3.2.6)Chapter 3. Fluids near planar surfaces 29and d is the diameter of species 3. For a soft wall interaction, the same conditionwill hold at a smaller separation z and is nearly always satisfied for z < 0; hence, forsimplicity and generality we useh(f,z) = —1, z <0 . (3.2.7)In terms of the expansion coefficients, Eq. (3.2.7) implies thath(z) = z <0 . (3.2.8)Taking advantage of this condition, Eq. (3.2.5) can be further simplified. Towards thisend, it is convenient to defineh(z) = h(z)H(±z) , (3.2.9)where 11(x) is the Heaviside step function(1, x>0,11(x) = (3.2.10)10, x<0.Equation (3.2.5) can now be expressed as the sumz) = 5z) + , (3.2.lla)wherep(_)fl+f1+M1) (nini)xJ h 4 (z2)c (r2i)Pi(cos O)dr2 . (3.2.llb)Here, the Fourier transform of (z) is now well defined since (z) —* 0 asz —* +oo and the remaining integral (z) can be evaluated in real space as shownbelow.Chapter 3. Fluids near planar surfaces 30If we started with Eq. (2.2.7c) in the reduction, Eq. (3.2.5) would be unchangedexcept that the function labels h and c inside the integral need to be switched, and adivision of the functions similar to those yielding Eq. (3.2.llb) is necessary. However, inthis 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. Nevertheless the direct correlation function does approach a limit and we can make use ofthe approximationc(fZ,z) Z < Zt < 0 , (3.2.12)where ITd/kBT = —(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 andc inside the integral switched. In this case, we emphasize that the equal sign is strictlyvalid only in the limit Zt — —°°.Let us now consider 5,(z). Expanding the Legendre polynomials in Eq. (3.2.llb)explicitly asPi(cosO)=acosiO=a(Z1Z2), (3.2.13)i=O r21where for convenience the coefficients a are given in Table 3.1, and using cylindricalOnn(—)coordmates, 7?Ov;w (z) reduces to00Onn(—) 1 2 Onn7OW/3 (zi) 21rp,j r Sn(zi/r)cou;.y(r)dr-y IzilIzi I—4ir6oH(—zi) p j r2c?(r) dr , (3.2.14)wherex—1, n=0,S(x)=a (xi+1—1) = (x2 — 1)/2, n = 1, (3.2.15)n=2,4,6,...Chapter 3. Fluids near planar surfaces 31Table 3.1: The first few coefficients of the Legendre polynomial expansion,Pi(x) =Pi(x)Po 1P1 1P2 -1/2 3/2P3 -3/2 5/2P4 3/8 -30/8 35/8P5 15/8 -70/8 63/8P6 5/16 105/16 -315/16 231/16P7 -35/16 315/16 -693/16 429/16and the particle labels have been dropped from the coordinates.Onn(+)The Fourier transform of iOv.w (z) is given by100________onn(+)(k) ikz1 Onn(+)710v;w/3 e 7lOii;w (z1)dziJ—00(2n+1) /n1nñ=(_)fl+fl1+niai (2ni+1) o oo)x I f ikz1 hofhfh(+) n1l0vi;wy (z2)cj,;y,(r21) Pi(cosO) dr2 dz1 . (3.2.16)Since in these equations k spans the reciprocal space of z, we define k = k to gete1 etl’21. Using the Rayleigh expansion [6,88]00eik.r= (21+1) ji(kr) Pi(cosO) (3.2.17)1=0and the identity f fdr1 f fdr21, together with the orthonormal properties ofthe Legendre polynomials we obtain=(_)fl+fl1+v1 (2n+1)7OL’;w13 )Y n1vil (2n1 + 1)n l jOnini+X ( 0 0 oJ (k) , (3.2.18)Chapter 3. Fluids near planar surfaces 32where00 Onn;(+)= I e’h (z) dz (3.2.19)Ov,w/Onn;(+)is the Fourier transform of h0; (z) and is the Hankel transform defined inSec. 2.2.2. Taking the inverse Fourier transform of Eq. (3.2.18) and using Eqs (3.2.11)and (3.2.14), the zi) coefficients for z1 > 0 can be expressed in the formP00Onn / ‘1ov;w*Z1) = 2irp7 / r2S(z1/r)cg3(r) drJz1-y(_)f11 (2n+1) (ninl\(2n1 + 1) 0 0(k) k)dk . (3.2.20)x J e1’Onini(+)Ovi ;wThe equations for z1 < 0 will also be useful in Chapter 4 and are given by Eq. (3.2.20)after replacing the first term by Eq. (3.2.14). When all indices are zero, it is easy toverify 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 withEq. (2.2.7c). The approximation (3.2.12) gives( Onn(+)Onn ( \ —CowZ)— { C0,;y (z), z > Zt (3.2.21)Z <ZtwhereOnn(+) OnnCOV;W (z) = Cov;w(Z)H(Z — Zt) . (3.2.22)Following the above analysis, we obtainP00z) —* 2irpJ r2S(z/r) h(r) dr-y(_)fl+fli+Vi I(2n+1) /n1nlnavil (2n1 + 1) 0 0P00X / e_i0fh71 kCoL,;.q3 ( ) h4y,y(k)dk , Zt , (3.2.23)for z> 0.Chapter 3. Fluids near planar surfaces 33Closure relations for the OZ equation were considered in Sec. 2.2.3 including a generalreduction in rotational invariants to yield Eq. (2.2.23) [10]. When c is a planar wall, onlycoefficients in the rotational invariant expansion of the closure relation Eq. (2.2.17) withrn = m1 = m2 = 0 contribute and Eq. (2.2.23a) reduces considerably to givec(z) Blum (—)“(2n+1) (71171271) (nln2n)V1 ‘2 VV1 V21 °° 8w°fl2fl2 ( ) Onn ( )x_jj h(z’) ‘‘ dz’ — tt w/3Z + bD(z) (3.2.24)where the property of the 9-j symbol [5](On1 li•i In1 11 00 2 12 = . 2 12 0 , 6fl16fl1116fl212 (3.2.25)(0 71 1 J (71 1 0 J /(2n + 1)(2n + 1)(2n2 + 1)has been used to reduce p tofOnin, fOn2nOnn — J J n2i Z’l+’2Poz— fOnn“ fl-i- ) ij6n11,6n212 (nin2n) (nin2n’ (3.2.26)+ 1)(2n2 + 1)(2n + 1) \“i V2 V \ 0 0 01(_y.(2n+1)(nhn2) (71171271). (3.2.27)3.3 Mean-field theories for polarizable surfacesIn this section, two mean-field approximations for the wall-particle potentials are considered. In one, the wall is modeled as a dielectric continuum and in the other as a jelliummetal.3.3.1 Dielectric surfacesIn the present model, the particles of the fluid (immersed in a vacuum) are near a semi-infinite dielectric continuum, characterized by its dielectric constant.The wall-particleChapter 3. Fluids near planar surfaces 34potential is given bysr I \ dielttwj3 UwZ) + U ‘where u(z) and u1 represent the short-range repulsive and electrostatic contributions,respectively. The functions and u1 depend on the orientational and positionalcoordinates of all fluid particles. In the dielectric model, when a charge q resides atdistance z from the dielectric (at z), an electric field results from both the charge andthe polarization of the surface. The field can be calculated using Maxwell’s equation ofelectrodynamics [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 —zinside the wall. The so-called image charge q’ which satisfies this condition is given by[89,90]= U q. (3.3.2)In the presence of the fluid, the polarization of the dielectric surface, and hencedepends upon the positions and orientations of all particles in the fluid. Here, we willdevelop a mean-field approximation which reduces the many-body potential to an effective wall-particle pair interaction suitable for use with the homogeneous integral equationtheories. However, it is worthwhile to note that an exact formulation of this many-bodyproblem could be obtained using the inhomogeneous equations when the dielectric surfaceis viewed as an external potential [9].In modeling the wall as a dielectric continuum, the mean potential generated by thesurface is calculated after considering the image-charge distribution of all particles in thefluid (see Fig. 3.2). The electrostatic interaction between a solvent particle and the wallconsists of a self-image (SI) term plus contributions from the images of all remainingparticles in the fluid (denoted 01 for other images). Therefore, one can writediel_ SIir \ i 01uw13 — tw/31,z1) 1- tL13 , .Chapter 3. Fluids near planar surfaces 35Figure 3.2: The coordinate system used in Sec. 3.3.1 for a dipolar fluid near a dielectricsurface showing two particle-image pairs. The images are shown to the left of the interfacewith a prime included in the labels.where z1 is the separation of the solvent particle from the surface and fi = (i, Oi, Xi)represents the Euler angles.The self-image term is a function of only the coordinates of a particle. If we denotethe image to be a second particle of species /3’, then u(O1) can be expressed in the formu(1i,zi) = ,—2zi) , (3.3.4)where u(fZi, f2,r12) is the electrostatic pair potential [Eq. (2.4.6)] and denotesthe orientation of the image 1’. Here, the self-image potential is reduced by a factor ofone half from the electrostatic pair potential to account in part for the work required topolarize the surface.In Eq. (3.3.3), the term u represents the instantaneous interaction of a molecule in2’Wall, Er12 — 2zFluidz2 2r121Chapter 3. Fluids near planar surfaces 36the fluid with the polarized surface for a given configuration of the remaining particles.We propose a mean-field approximation in which the instantaneous interaction betweena molecule and the surface is replaced by the result obtained after averaging over thepositions and orientations of the “other images”. This yields the effective 01 potentialu(i, z1) = J 1, r1 — 2z)g(l; 2) d(2) , (3.3.5)where g,(1; 2) is the probability density that particle 2 of species y has coordinates(12, r2), and hence 2’ has (fZ, r2 — 2z), given that particle 1 has coordinates (1Z, zi).The 01 potential is also reduced by a factor of one half to account for the work requiredto polarize the surface. Note that we distinguish the effective potential u(Zi, z1) fromthe instantaneous result u by specifying that it depends only on the coordinates of asingle particle.In order to proceed with the reduction of the wall-particle potential to a rotationalinvariant expansionu(O1) =we need to define the multipole moments of the images. To do this, consider the chargedistribution n(r) of a molecule of species ci in the space-fixed reference frame. Thedistribution of charge of the “image species” cV (denoted na’) is obtained after reflectingn,. through the surface at z = 0. Since the image of a vector r = (r, 0, so) is (r, r — t9, so),we have—= U ::) na(r,,so), (3.3.6)Chapter 3. Fluids near planar surfaces 37where 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_— w\ (a)(_)fl+V’+‘.(3.3.7)However, in the pair potential we require the multipole moments in the molecular-fixedreference frame. The moments in the space-fixed and molecule-fixed reference frames arerelated by a rotational transformation,= R7(’) Q2 , (3.3.8)Li1where 1’ denotes the orientation of the image. If a has orientation Z= (, 0, x) then=(q, ir — 0, r— x) and we can use Eq. (3.3.7) and the property [5](_)‘“ R() , (3.3.9)to obtainQp;a • (3.3.10)— \1 +The coefficients of the particle-image pair potential can now be written asmnl;el —) Om+n,l j (21 + 1)! Qi;13v;’(r) ( m’UILi;/3l= fmnl (2rn)!(2n)! r11= oh (1_—Ew)mnl;el\1+ujp;(r). (3.3.11)and the self-image potential is given by1 i’l—______mnl;el mu(01) =+(_)Vup;1313zi). (3.3.12)mnlpa,Chapter 3. Fluids near planar surfaces 38With 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 wall-solvent rotational invariant expansion appropriate for use with the wall-solvent RHNCtheory. Doing this givesuj((01) = u ‘(zi)8(O,1i,O) , (3.3.13a)— n+L.’ 3/2U’(zi)= 1 (1— :w) ( ) (2fl+1) A u11(2z) , (3.3.13b)m1nIil.1 ‘1whereA = (::‘;) [(-)‘ ( ) ( g)]. (3.3.13c)Turning now to u(01), the distribution function g(1; 2) used in Eq. (3.3.5) isdifficult to obtain exactly but can be approximated by the superpositiong(1,zi;fZ2r) . (3.3.14)If only the excluded volume contribution to the particle-particle distribution function isretained, one has the simpler approximationg(i,zi;1Z2,r) g-(2,z)H(ri — 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 particleparticle and wall-particle functions as in Eqs (2.1.6) and (3.3.13a), relationships of theformug°’(zi) = !i ( ::) (_)n+vFmnl Ju(\/r + 4z1z2)x [h(z2)+ 6o] H(r12 — d,)P,(cosO2)dr (3.3.16)Chapter 3. Fluids near planar surfaces 39are obtained. Here Pj(x) is the Legendre polynomial of order 1, cos O2 = —(z1 +z2)//r + 4z1z2 andF m+nfOnnfmnl/ 2m+1 (mnl\f°mm (2n+1) o oo)Blum m+m /2m+1 “mni”V 2n + 1 0 0 (3.3.17)Noting that one can write H(r12 — d) = 1 — H(d — ri2), a trivial rearrangement yields/1— 6w”.ç-Omm;OI__________nI —mnl;elUO;w/3 (z1)=irp7( + ) (_)n+um-yTi ,00 12 0z1 +z2)’[h (z2) + o] dz2 / x2Pj(x) dxxL-J J—1rz+d 1—lr Onn +60]dz2 f’xl_2p( )dx—Jz1-d(1 + 2) [07(2)JA ], (3.3.18)—mnl;el 1-H mnl;elwhere = r fl,;-y (r), x = cos°2’ and A = (z1 +z2)/d + 4z1z2. Using theorthogonality property of the Legendre polynomials,1I P(x)Fj(x)dx = (2j+ )6i3 (3.3.19)—1together with the expansion [88]1—2x12 =bP(x) , (3.3.20)i=0one can demonstrate that the first integral over Legendre polynomials is zero when 1 2(Note when 1 > 2, i 1 — 2 < 1 and the integral jx2Pi(x)dx = 0). The remainingintegrals are evaluated using the related expansionPi(x) = ax , (3.3.21)i=0where the coefficients a are listed in Table 3.1. Combining Eqs (3.3.18)—(3.3.21) givesOmm;OI______/ 1—jzi +d, Onn mn;cavUoj;13 (z1) = 7rYZpy (+ z1—cihoz,;w(z2)T,;fry (zi,z2)d , (3.3.22a)13-yChapter 3. Fluids near planar surfaces 40wherea (x+y)T;c(x, y) = Sm+ni()m mç (3.3 .22b)(1+ i — 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 functions and reducing with the use of Eqs (3.3.11), one obtains the general resultP00Omm;OI_______UO;W (zi) =(1— I Onn mmn;fuj0 (zi,z2)d ,(3.3.23a)whereF00 r m1nili;el ‘ 2mn;ful (zi,z2) r / r12dr Livi; (Vr12+4z1z2)m1n m2n212fL1Vj ?z’A’m2n212X 2M2/3Y()d,0(O12d,0(O2)] (3.3.23b)and£Onn= ()m+n+fll+LIl+ILfmlnhllfm2n2j (2m + 1)3/2 (m2m1 m) (n2 n1 n\fOmm, + 1 2 1 2 l/m1nl’1m2n\/rn2m1’ ,‘n2n1\, ) o) . (3.3.23c)Here, d() = R(0,O,0) and cos 012 = (z2 — zi)/r12. The numerical evaluation ofmmflfUlEqs (3.3.23) is time consuming. However, the functions (z1,z2) need to be calcuOmm;OIilated only once for a given model of the bulk fluid, and when is known, UO;W z1)can be obtained rapidly from Eq. (3.3.23a). We will discuss results for dipolar fluids inSec. 3.4.Chapter 3. Fluids near planar surfaces 413.3.2 Jellium surfacesIn this subsection, we consider a simple quantum mechanical model for a metallic surfaceand 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 phasesand hence represents a model for what is known as an “ideally polarizable electrode” inelectrochemistry.Model, potential and general theoretical approachWe represent the metal surface by the so-called “jellium model”, in which the atomiccores (including the nucleus and core electrons) are replaced by a uniform backgroundof positive charge [80,91]. This model is characterized by specifying the average numberdensity of the atomic cores n+, or equivalently, the Wigner-Seitz radiusr3 = (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 metalslab of finite width 2z (see Fig. 3.3) with the origin chosen at the center. This is themodel used by Gies and Gerhardts [91] to study the metal-vacuum interface and froma numerical perspective is a more convenient and accurate model than the semi-infinitewall of Lang and Kohn [80]. In particular, the boundary conditions are well defined inthe finite slab (either the wavefunction or its derivative is zero at the center of the slaband the function decays to zero at long range). In the Lang and Kohn method for semiinfinite 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 thebulk (plane wave) wavefunction. This matching condition makes it difficult to satisfycharge neutrality with sufficient accuracy.Chapter 3. Fluids near planar surfaces 42r0Figure 3.3: The coordinate system used in Sec. 3.3.2 for a dipolar fluid near a jelliummetal slab.It is well-known from early work on metallic surfaces [80,91] that the electrons ofthe metal spill several angstroms out from the surface and provide an electric field Ejelstrong enough to polarize the liquid near the interface. The liquid in turn polarizes themetal through a potential which depends in principle upon the position and orientationof all particles in the fluid. Hence, in general, E1 also depends on the instantaneousconfiguration of the liquid. Here, as in the dielectric continuum case, we will use amean-field approximation which reduces these inherently many-body terms to an effectiveone-dimensional interaction.2T21 =(r2i,6ci)1z1I Iz LChapter 3. Fluids near planar surfaces 43The interfacial structure (including the electron and molecular distributions) is obtained by the self-consistent solution of two coupled theories. The liquid phase is described by the RHNC approximation [92,51], and the metal side by the HohenbergKohn-Sham implementation of density-functional theory [82,83], as further developed byGies 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 thekth state, m is the mass of an electron, is the exchange and correlation potentialand Ue represents the instantaneous interaction energy of an electron in the field of thejellium and surrounding liquid for a given configuration of the particles in the fluid. Weuse the local density approximation to solve Eq. (3.3.25) with Wigner’s expression forthe correlation energy’ [93] such that[0.611+1474r8(r) + 23.4ao], (3.3.26)r5(r) [r5(r) + 7.8a0]where r3(r) = (4n(r)/3)/3, e is the magnitude of the charge of an electron,n(r) = 2 ‘T!k(r) 12 (3.3.27)is the electron density distribution and F is the Fermi energy. The factor of 2 inEq. (3.3.27) is to account for the electron spin degeneracy.In addition, we propose a mean-field approximation in which the instantaneous interaction of a particular electron with the external potential is replaced by the resultobtained 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 numerator7.8 should read 23.4.Chapter 3. Fluids near planar surfaces 44yields the mean-field electrostatic potentialUe(Z) = el(Z) + ect(Z) , (3.3.28)distinguished from the instantaneous potential Ue b specifying that it depends only onthe linear coordinate z. In Eq. (3.3.28), UJe1(Z) and u(z) denote the contributions to theu(z) due to the jellium background and fluid, respectively. The Schrödinger equation isthen simplified considerably and reduced to the one-dimensional form [91]h2 d2——-----i/,,(z) + [Ueff(Z)— €,] i/,(z) = 0 , (3.3.29)2mdzwhereUeff(Z) = Ue(Z) + u(z)= iel() + ut(z) + u(z) , (3.3.30)— —(k + kg), (3.3.31)= 2’re1 ‘‘IJk(r) , (3.3.32)and hk, hk are the x and y components of the electron momentum. The electrondensity can be reduced ton(z) = - (F-e) I (z) 12. (3.3.33)Using normalized wavefunctions, the Fermi energy can be determined from the chargeneutrality conditionn(z)dz = 2zn÷, (3.3.34)which together with Eq. (3.3.33) givesCF= (2nfl+Zw+ /, (3.3.35)Chapter 3. Fluids near planar surfaces 45where 11F is the number of eigenstates with energy e, < F•The potential contribution from the metal is obtained from the solution of Poisson’sequation (multiplied by the charge of an electron)= 4e2 [nW(z) — n(z)] , (3.3.36)where11, HI<z,W(z) = (3.3.37)10, HI>z.Integrating once and dividing by e, we get the electric fieldEiel(z)= 4e j [nW(z’) — n(z’)] dz’. (3.3.38)A second integration gives the interaction potential— uiel(O)= 4re2 j(z — z’) [nW(z’) — n(z’)] dz’= e[ZEe1(Z)— D(z)] , (3.3.39)whereD(z) = 4e j z’ [nW(z’) — n(z’)] dz’. (3.3.40)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 DFtheories, respectively.Consider first the mean scalar potential qt(z) = _uext(z)/e at z defined by= Jgwa(02) q(2,1)d(2) , (3.3.41)where q(2,1) is the scalar potential at a position z1 from the center of the slab due toa single molecule of species a with coordinates 2. Since Eq. (3.3.41) is in the form of aconvolution, Fourier transform methods can be used to calculate qt(z). However, sinceChapter 3. Fluids near planar surfaces 46a(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,-m(2, 1) =(_)m2m + 1 Q;amom, Or21) (3.3.42)mifmom m+1 i0 2,r21(this is equivalent to the pair potential expansion [Eq. (2.4.6)] when= toobtain(_)m+a g(z2)Pm(cos 021)eXt() =2m + 1 f m+1 dr2. (3.3.43)r21Using the Legendre polynomial expansion Eq. (3.2.13), the integral reduces toOmmfJ go;wz2, Pm(cosO2l) m r°° z21 0mmm+1 dr2 = 2ir a_______dz2r21 (m + i — 1) Jo iIm+i_1Y0;wa(Z2)= 27r ( a__0mm j/ m ){jzi 1(m + i — 1) 0+(_)m I 1g(z2)dz2] , (3.3.44)Jz1 Iz2iwhere we have used the property that (_)i = (_)m for nonzero coefficients ar (seeTable 3.1). When m 2, we havem a 1(m + = Xm 2Pm(X) dx = 0, m 2, (3.3.45)That is, a uniform planar distribution of quadrupole and higher order moments make nocontribution to qSt(z) and a test charge interacts only with wall-ion and wall-dipole pairdistributions functions. Hence, ext() is given in general asrext(1)= —47r j (zi — z’) L paqa2wa(z’)] dz’4(YPc;a j’ g.(z2)dz2, (3.3.46)3 c4twhere we have added the constantP00—27r I (zi — z’) [Pa92(z’)j dz’ — (YPa;a j g,(z2) dz2.JoChapter 3. Fluids near planar surfaces 47The interaction potential of a particle with the metal wall is derived in Appendix Band can be written asua(O1) = u a(zi)(O,i,O) , (3.3.47a)whereOnnel —n /2n + 11 0ieI(z)uov;wa(Z) Qv;a (3.3.47b)andiel() = _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 metallicsurfacesIn this section, we consider three model surfaces in contact with a fluid of dipolar hardspheres and in contact with solutions of charged hard spheres in a dipolar hard spheresolvent. The effects of salt concentration and ion diameter on the interfacial structureare examined. The long-range asymptotic behaviour of a pure dipolar fluid near aninert surface as implied by the RHNC theory is also derived and compared with classicalexpressions and exact results [62].3.4.1 Fluid modelThe solvent (s) is characterized by the dipole moment p and hard sphere diameter d5;the cation (+) and anion (—) are characterized by charges qj, and diameters 4. If wechoose 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 orChapter 3. Fluids near planar surfaces 48the correlation functions and, by symmetry, only coefficients in the rotational invariantexpansion with(3.4.1)are nonzero. For this choice of molecular reference frame = = ,u, is the onlymultipole moment needed to describe the solvent, and the solvent-solvent, solvention, and ion-ion, u±±(r), interactions are given by the pair potentials= uHS(r) + u(r) 2(12) , (3.4.2a)u±(12) u(r) + u±(r) iO1(12) , (3.4.2b)u±(12) = u(r) + ug?(r) , (3.4.2c)respectively, where the hard sphere potential u and multipole moment interactions aregiven by Eqs (2.4.2) and (2.4.6).The condition (3.4.1) also implies a significant reduction of the wall-solvent pair functions and the expansion (3.2.1) can now be written simply asgw(O1) = (_)?2gTh(zi)Pn(cos0i) , (3.4.3)where the subscripts (v = 0) have been dropped from the coefficients and O is theangle between ji and the surface normal vector . Since the wall-particle interactionsalso include a hard repulsive contribution, the wall-solvent, and wall-ion, pairpotentials can be written as(Do, IzH(zw+ds/2,‘u(01) = (3.4.4a)1. u(z), I z z + d/21°° Izl<z+d±/2,u±(01) = (3.4.4b)I z z + d/2where u(z) denotes the electrostatic contribution to the wall-particle potential definedfor the dielectric and metallic surface models and we have used z to denote the zChapter 3. Fluids near planar surfaces 49coordinate of the surface. For an inert wall, u(z) = 0. In the case of a jellium metalwall, u(01) is given by Eq. (3.3.47). For dielectric walls, we consider only the puredipolar hard sphere solvent andu(O1) u(01) + u(01) , (3.4.5)where the self-image term is given byu(O1) = ( w) (2) [2Po(cosOi) +P2(cosOi)]. (3.4.6)For dipolar fluids, one finds that the cavity approximation Eq. (3.3.22) implies thatu(01) = 0 . (3.4.7)This is only true for simple dipolar fluids. In general, there are contributions from theimages of other particles, even when the cavity approximation is employed. Using thesymmetry 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 yieldu(O1) HmuOrnm;OI(z) Pm(cos O) , (3.4.8)whereum;OI(z1)= 2Ps ( f°°g(z2)7;(zi, z2)dz2, (3.4.9a)and z2) has only the dipole-dipole contributionI’’(zi, z2) = r J r12 dr12 [u2;e1(\/r? + 4z1z2)m2n212 1z2 zilx g2fl2l2(ri)d,0(O12)d10(O2)j . (3.4.9b)Chapter 3. Fluids near planar surfaces 50In general, the properties of the bulk solution are completely determined by specifyingthe reduced parameters1LS/\/kBTd , (3.4.lOa)= pd , (3.4.lOb)q* ±q±/kBTd , (3.4.lOc)pd = p_d , (3.4.lOd)4 = d/d8 , (3.4.lOe)where Pa denotes the number density of species a. When we consider fluids near metallicsurfaces, it will also be necessary to specify the fluid temperature and the solvent diametersince the reduced elementary charge e* = e/V”kBTd is needed to define the electrostaticmetal-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* 8. For monovalent electrolytes with pa”, and q* fixed, the onlyremaining solution variables are the concentration and the ion diameters. The reducedion 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-raydensity measurements of ionic crystals [94] and were chosen to be consistent with earlierwork [11—13,47—49,95]. The ions labeled “Eq” have the same diameter as the solventparticles. Below we will also refer to a salt labeled “ILi” where the ion diameters are thesame as model Lii, but the charges are switched. The characteristic parameters for allsystems considered are summarized in Table 3.3. It should also be noted that the modelion 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 dipolarhard sphere solvent have previously been obtained using the RHNC theory [10,45] andChapter 3. Fluids near planar surfaces 51Table 3.2: Reduced ion diameters d±/d together with the examples and labels used inthis work. The solvent diameter d was taken to be 2.8A.d±/d cation anion0.68 Li0.84 Na F1.00 Eq Eq1.16 Rb C11.28 Cs Br1.44the reader is referred to these articles for details. A basis set study was also includedin Ref. [10] and it was found that truncating the basis at n, m 2 was sufficient tocalculate the pair distribution functions. In the present work, the basis set is truncatedwith n, rn < 4 or 6 to ensure accurate results. Results for n, m 4 are indistinguishablefrom those obtained with ii, m <6.All RHNC calculations were performed on a grid of 0.01d8 extending at least 212points from the surface. The bulk particle-particle reference bridge functions were obtained from the parameterization of Lee and Levesque for hard sphere mixtures [27,28];the one-component wall-solvent reference bridge functions from the parameterization ofHenderson and Plischke [31]. Unfortunately, wall-particle bridge functions are not available for hard sphere mixtures. Hence for electrolyte solutions near a wall, the necessarybridge functions were constructed by “shifting” the Henderson-Plischke results [31]. Thatis, we have used the Henderson-Plischke function (for a reduced density of 0.7) with thecontact position shifted to correspond to the appropriate ion diameter. We would expectthis to be a reasonable approximation for the present systems where the reduced iondensities and size differences are not very large.The density functional calculations for the jellium metal used a grid width of 0.02AChapter 3. Fluids near planar surfaces 52Table 3.3: Systems considered in the present work. In all cases the surface is in contactwith a solution of total reduced density p” = 0.7. The reduced solvent dipole momentand ionic charge are = and q* = 8, respectively. For the metal slab, we assumethat T = 298K and, unless otherwise stated, r3 2.65a0 where a0 = 0.529177A is theBohr radius.salt p p molaritynone 0.7 0 0LiT 0.675 0.0125 0.945NaC1 0.675 0.0125 0.9450.6 0.05 3.780.5 0.1 7.56EqEq 0.675 0.0125 0.9450.6 0.05 3.780.5 0.1 7.56RbF 0.688 0.006 0.4530.675 0.0125 0.9450.6 0.05 3.780.5 0.1 7.56CsF 0.675 0.0125 0.945CsI 0.675 0.0125 0.945ILi 0.675 0.0125 0.945and an infinite potential barrier was placed 8A out from the jellium surface. At thisdistance, the electron density distribution is essentially zero in the absence of the potentialbarrier and is unaltered by its presence.3.4.2 Asymptotic behaviour of h(01) near an inert wallIt is possible to determine analytically the long-range asymptotic behaviour of the wallparticle pair correlation function within the RHNC (or HNC) approximation when thebridge functions rapidly decay at long range. Near an inert wall, the RHNC equationsfor simple dipolar fluids require that only even-n projections occur in the pair expansionChapter 3. Fluids near planar surfaces 53Eq. (3.2.1) such that Eq. (3.2.20) leads to the simplified relationship= 2Ps jrS(zi/r)c(r)dr(2n+1) (ninl+PsN(2+1) 0 00xJhmnh(z)c(ri)Pi(cosO)dr, z1 >0. (3.4.11)Since c(01) —* —u(O1)/kBT as z —* oo, we would expect c(01) to be shorter rangedthan h5(01) and h(01) —* (01) in the same limit. Hence the required asymptoticlimit can be obtained from Eq. (3.4.11). Explicitly, for those projections which dominateasymptotically, it is not difficult to show that as z —* och(z) 2rp° jr2S(z/r)c(r)dr, (3.4.12)where= {i —p20(0)//2n +[i + pssTh0(0)//2Tj + 1] . (3.4.13)The OZ equation at k = 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 thelimiting behaviour of cz(r).Within the HNC/RHNC approximation, the solvent-solvent direct correlation function has the asymptotic limitc8(12) h(12) — u(12) , r —* oo. (3.4.14a)For the dipolar model the leading terms are determined by the projections h2(r) andu’2(r) which decay as r3, and Eq. (3.4.14a) yields1 2 112( )c(12) [h2(r)j2(12)]—kBT2(12), r —* oc . (3.4.14b)Chapter 3. Fluids near planar surfaces 54For mril 112 we obtain________) -‘-‘oo15 1( —“2I2pmnlc’(r) (47Tps)2 L 3Ys r6 r —+ oo, (3.4.15a)wherepmnl — I (ii2(f 2, 2 mnl*(ç1,2,I12)dcdfzf—00 (3.4. 15b).LJ00— f ‘(1Z1,2, ri2) U1 (ni, 2,12)dfzdfi’Also, we have used the exact limiting expression [96]h2(r)— 1) 1V3 47r,osEsY — r —* oo (3.4.15c)where e is the dielectric constant of the solvent and y = 4Kp/9kBT. From Eq. (3.4.15b)we find that all the B’ vanish except for B°°° — p022 — 1/5. Using Eqs (3.2.15), (3.4.12)00 ‘‘OO —and (3.4.15b) we find that as z —* co0001 (_—_1)212 (3.4.16a)_____I________WS/ 327rp5 L 3yc j z32201 F(_—_1)212h°22z I“/ 64rp L 3YEs ] . (3.4.16b)Equations (2.1.6) and (3.4.16) imply that within the HNC/RHNC approximation1 F( — 1)212 1 000 1 220P(cosO)), (3.4.17)h(Ofl —_____ __ _32Ps [ 3Ys ( +as z —* oo.It is interesting to compare the HNC/RHNC expression with the exact asymptotederived by Badiali [62] and with the result implied by the “image-dipole” potential ofclassical continuum theory. For comparison purposes it is convenient to express theresults as followshHNC(Ol) — i (S_1)2(S_121 4 9___ _F 000+ X220P(cosO)] , (3.4.18a)WS \ kBTz3 3y \. jXss01) I—1 (—1) 1lexact(It1v \ kBTz3 3y ) c(c + 1) 16 Z7(n)x P(cosO) , (3.4.18b)hdl(01) ( — 1) 1 4 2— kBTz3 + 1) +P(cos 0)] . (3.4.18c)WS \Chapter 3. Fluids near planar surfaces 55In Eq. (3.4.18b), n must be even and7(n) = n,0 + + a(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 thatthe cx(n, 0) contributions do not contribute in the HNC/RHNC case. Also, it is obviousfrom Eqs (3.4.18a) and (3.4.18b) that the HNC/RHNC “prefactor” differs from the exactresult. The classical continuum expression does not include x° contributions which isanticipated in the absence of a structured solvent.It is possible to show [62] that when averaged over orientations Eq. (3.4.18b) yields1 1 i 000(3 4 19)WS 16K c(c + 1) op z3 ‘which can be compared with the HNC/RHNC result given by Eq. (3.4.16a). Furthermore,using the expansion1 + 3y + 3y2(1 + b) + Q(p) , (3.4.20)where b is related to the second dielectric virial coefficient, together with the result1 3421XSS— 1+2B2(T)p+ ‘ . . )where 1T = is the isothermal compressibility and B2(T) is the second virial coefficient, it is possible to derive low density expansions for the h°(z) asymptotes. Explicitly,in the low density large z limit we obtain{aps + [4a2b — 2B(T)a — 2a]p + Q(p)} , (3.4.22a)8T {aps + [3a2b—2B2(T)a— a2]p + O(p)} , (3.4.22b)—8T {aPs + [a2b — a2]p + Q(p)} , (3.4.22c)where a = 4irz/9kBT. We note that all three results are identical only to linear orderin density.Chapter 3. Fluids near planar surfaces 563.4.3 Results for a pure dipolar fluidThe structure of dipolar fluids near inert planar and spherical surfaces has been the topicof several publications [54,45,50,51]. However, these systems have never been solved nearplanar surfaces using the full RHNC theory and it is useful to understand the structure ofthe fluid in this case before considering fluids near dielectric and metallic surfaces. Thelong-range asymptotic behaviour implied by the RHNC theory will also be considered.Pure solvent near inert surfacesThe four nonzero projections of the wall-solvent pair distribution function g(z) included in the numerical calculations are shown in Fig. 3.4 (note that z = (z— z)/dwhere z denotes the z coordinate of the wall). The projections for the dipolar fluid neara large inert sphere of diameter 30d [50] and with a hard sphere fluid near a planar inertwall are included. Only even-n projections occur in the expansion because symmetrydictates that gws(0, z) = gws(lr — 0, z) and because there is equal probability for a dipoleto be directed towards or away from the surface.The oscillation in the angle-averaged distribution function g° indicates ordering ofthe fluid perpendicular to the surface and arises from the inhomogeneous pressure exertedon a particle near the interface. In the case of dipolar hard sphere fluids, the contactvalue drops from the pure hard sphere value as a result of electrostatic interactions amongthe solvent particles and corresponds to a reduction in the bulk pressure Pb. This resultis anticipated from the contact value theorem which states for a hard wall interactionthat the contact density is proportional to the bulk pressurepg°(d/2) = Pb/kBT . (3.4.23)At contact, dipoles favour orientations parallel to the surface (i.e., cos 0 i—’ 0 and g2 <0),since this configuration lowers the interaction energy among particles in the dense surfaceChapter 3. Fluids near planar surfaces 57Figure 3.4: RHNC results for the projections occuring in the wall-solvent pair distributionfunction for a dipolar hard sphere fluid with j = and p = 0.7. In (a)-(d), the solidcurves represent the present planar wall calculations and the dotted curves are resultsobtained previously (Ref. 10) for a macrosphere of diameter 30d.5432100.40.200 1 2 3z*0—1—20—0.02—0.04—0.06—0.080 1 2 3z*0 1 2 3z* 0 1 2 3z*Chapter 3. Fluids near planar surfaces 58layer. Further out, as indicated by g2, the orientations oscillate between parallel andperpendicular 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 distribution functions are almost identical. Thus, for a simple dipolar fluid, the 30d macrosphereresults already approach the flat plate limit at short range. However at long range, thiscannot be true since the macrosphere-solvent (ms) projections h(r) and h2(r) (r isthe center-center separation) decay asymptotically as r6, whereas the correspondingwall-solvent projections decay as z3 as discussed in Sec. 3.4.2.The long-range behaviour of h°(z) and h2(z) is shown in Fig. 3.6. The RHNCasymptotes were obtained from Eq. (3.4.16) using the values = 49.24, x° 0.0738and x° = 1.349 determined from the solution of the RHNC equations for the bulkdipolar hard sphere fluid. The small discrepancy discernible in Fig. 3.6b between theanalytical and numerical RHNC asymptotes gives an indication of the accuracy of thenumerical solutions (i.e., i0). The classical continuum limits are apparent fromEq. (3.4.18c) and were evaluated using the RHNC dielectric constant. Also included inFig. 3.6a is the curve given by the formally exact Eq. (3.4.19) using RHNC values for thebulk fluid properties. It should be noted that this result is “exact” only in the sense thatthe formally exact expression is used. The true exact asymptote could only be found ifexact values of Es, 9Es/öP, and x° were known and this is not the case for the presentmodel. It can be seen from Fig. 3.6 that in the full RHNC result the oscillatory structurecharacteristic of molecular solvents persists until z 13d5. At larger separations theprojections obey the analytically determined asymptotic limits. For h°(z), we observethat both the RHNC and “exact” curves lie below the continuum result, but that theChapter 3. Fluids near planar surfaces 59Figure 3.5: The wall-solvent pair distribution function g(9, z) for the pure dipolar fluidnear an inert wall. The solvent parameters are as in Fig. 3.4g5(O,z)z/d•5_-__%.,_• cos9Chapter 3. Fluids near planar surfaces 60Figure 3.6: The asymptotic behaviour of (a) h° and (b) h2. The dipole moment anddensity 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 continuumasymptotes. The dash-dot curve included in (a) represents the exact asymptotic behaviorgiven by Eq. (3.4.19) and evaluated as discussed in the text.RHNC approximation overestimates the magnitude of the coefficient occurring in theasymptotic expression. From Fig. 3.6b, it is obvious that the RHNC result for h2(z)also lies well below the continuum curve. Unfortunately, unlike the case for hO;ct(z),there does not appear to be any way of obtaining h 2;et(z) in the absence of explicitresults for the three-body direct correlation function.Pure solvent near dielectric surfacesWe now consider the dipolar fluid near hard walls having dielectric constants of 5 and ccusing both the cavity approximation for wall-solvent pair potential,0.0—0.1—0.20.0—0.2—0.4—0.6—0.8—1.05 10 15 z’ 5 10 15 zu(01)=u(1Zi, zi), z > 0 (3.4.24)Chapter 3. Fluids near planar surfaces 61which amounts to retaining only the self-image contribution in the present model, andthe full superposition resultu(O1) = u(fi,z1)+ u(fZi,zi), 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 calculations are shown in Fig. 3.7 for both approximations and dielectric constants c. Theresults near an inert surface (c = 1) are also included. In general, results for the cavityand superposition approximations are significantly different whereas those obtained fora given approximation with = 5 and € = co are similar. Also, it is apparent fromFig. 3.7a that g°(z) does not have a strong dependence upon the value of.exceptvery near contact where the density is nearly double that of an inert surface. The largeincrease in g°(z) at contact is a direct consequence of the self-image potential whichis attractive for all molecular orientations. For theg2(z) projection, the superpositionand cavity approximation curves are significantly different even several diameters fromthe surface. Also note that the full superposition results at short range are similar tothose obtained for an inert wall. This is important because this projection makes a largecontribution to g(01). The curves obtained forg4(z) andg6(z) with the full superposition and cavity approximations are also plotted in Fig. 3.7 and essentially bracketthe inert wall results.The physical significance of the results obtained for different models and approximations is most easily understood by considering the full wall-solvent pair distributionfunction g(01) plotted in Figs 3.8 and 3.9 for = 5 and w = respectively. Inthe cavity approximation (Figs 3.8a and 3.9a), particles near the surface show relativelylittle orientational preference when = 5 and prefer orientations parallel to the surfaceChapter 3. Fluids near planar surfaces 6210286 0________—22—400 00.01.0—0.10.5—0.20.00Figure 3.7: RHNC results for the projections occuring in the wall-solvent pair correlationfunction for a dipolar hard sphere fluid with = and p = 0.7. The solid curverepresents results for an inert wall ( 1). The long- and short-dashed curves wereobtained in the cavity approximation for = cc and 5, respectively. The long-dash-dotand short-dash-dot curves are the results obtained with the superposition approximationfor = cc and 5, respectively.0008.0(a)I I_\I I I I04.001 2 3z* 1 2 3z*0 1 2 3z* 1 2 3z*Chapter 3. Fluids near planar surfaces 63normal when = oc. Further from contact, the particles shows an oscillatory structure similar to the case for an inert wall but with orientations parallel to more stronglyfavoured. 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. Itis clear (for orientations parallel to ) that the images of other particles act to opposethe influence of the self-image interaction and again orientations perpendicular to arepreferred. Further from the surface, the structure is also similar to that found for theinert 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 termalone is clearly insufficient to describe the interaction with a dielectric continuum surface.Pure solvent near jellium metal surfacesHere we consider the molecular and electronic structure of the metal-solvent interface andcompare results for several jellium densities. However, before proceeding, it is necessaryto consider the finite slab model used in the jellium calculations since many propertiesof the metal (such as the Fermi energy, the potential drop and the electron densitydistribution) vary as function of the width 2z. In Fig. 3.10, the potential drop LSacross the jellium-vacuum interface is plotted as a function of z, when r3 = 2.65a0 (AWigner-Seitz radius of 2.65a0, where a0 = 0.529177A is the Bohr radius, correspondsapproximately to the bulk electron density of Hg.) The potential oscillates with a periodof 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.1volts (approximately the order of magnitude of weak van der Waals interactions). Forlarger 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-infiniteChapter 3. Fluids near planar surfaces 64gws(O,z)z /dgws(6,z)cos 0cos 0Figure 3.8: The wall-solvent pair distribution function gws(O, z) for a pure dipolar fluidnear a dielectric surface with E = 5. The results in (a) and (b) are for the cavity andfull superposition approximations, respectively. The model parameters for the solventare given in Fig. 3.7.(a).- —--.(b)f10z/d•Chapter 3. Fluids near planar surfaces 65gws(6,z)1.010.0z Id8gws(O, z)cos 0cos 0Figure 3.9: The wall-solvent pair distribution function g(O, z) for a pure dipolar fluidnear a conducting surface ( = co). The results in (a) and (b) are for the cavity andfull superposition approximations, respectively. The model parameters for the solventare given in Fig. 3.7.(a).100(b)z.Chapter 3. Fluids near planar surfaces 663.32 3.323.28 3.283.24 3.243.2 3.28 29 30 31z(A)Figure 3.10: Oscillations in the net potential drop across the jellium-vacuum interfacefor r8 = 2.65a0j ellium-vacuum interface2.The Fermi energy and potential drop across the jellium-liquid interface show similarbehaviour. In particular the magnitude of the oscillations in the Fermi energy is i-.’ 0.016volts at narrow slab widths and 0.002 volts at the larger widths. In comparison to thesolvent effects considered below, the dependence on slab width is very small, particularlyat 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. Alsoin these calculations, the solvent particles have a diameter of 2.8A and the temperatureis taken to be 298K. The solvent diameter and the temperature are required in order tospecify the jellium-solvent interaction.The self-consistent electron density distributions for a jellium slab immersed in thedipolar fluid are shown in Fig. 3.lla for three values of r. The electron density shows thecharacteristic Friedel oscillations inside the metal and quickly decays to the bulk density2The potential drop across the semi-infinite jellium-vacuum interface was calculated using N.D. Lang’soriginal self-consistent program [801.I II I liii —: -Lçb (volts):119 10z(A)Chapter 3. Fluids near planar surfaces 67further in the slab. The magnitude and frequency of the oscillations depend stronglyon r5. Just outside the jellium surface, the density decays exponentially and drops to0.1% of the bulk value at 3A.It is interesting to compare the electron distribution at the jellium-solvent interfacewith the jellium-vacuum case. The inset in Fig. 3.lla shows the change in the electrondensity distribution resulting from the presence of the dipolar fluid. It can be seen thatthe solvent only slightly perturbs the distribution and that the maximum deviation isapproximately three orders of magnitude smaller than the bulk electron density. Alsofrom the inset, we see that the electron density spills further out of the metal due to thepresence of the solvent. The effective potentials obtained for the jellium-solvent systemsare plotted in Fig. 3.llb. The result for the jellium-vacuum interface with r = 2.65a0is also shown. The oscillations which occur just outside the jellium edge arise from thefluid structure. The effective potential is dominated by the exchange-correlation energy.The electrostatic contribution to Ueff will be considered in detail below (see discussionrelated to Fig. 3.28).We now consider the self-consistent fluid structure at the interface. The first fewnonzero projections of the pair correlation functions are shown in Fig. 3.12 and theresults are compared with those for an inert surface. (The inert surface can be thoughtof as a semi-infinite jellium metal with r = oo). A direct consequence of the interactionof the two phases is that the probabilities of a dipole pointing into and out of the surfaceare no longer equal. Hence, odd-ri projections now occur in the Legendre polynomialexpansion. Near the surface, the curves exhibit general3-dependent trends and rapidlyapproach the inert wall results at larger separations. From Fig. 3J2a, it is apparent thatnear the surface g°(z), and hence the local solvent density, increases with deceasingr8. We also note that the second peak in the function moves closer to the surface asr decreases, implying that the solvent packs more tightly for smaller values of r3. TheChapter 3. Fluids near planar surfaces 681 00.8—0.50.60.4 —10.2—1.50—15 —10 —5 0 5 z —15 —10 —5 0 5 zFigure 3.11: The self-consistent electron density (a) and effective potential (b) of ajelliumslab, 64A wide, immersed in the pure dipolar hard sphere fluid. The dotted, solid anddashed curves are for r5 = 2.0, 2.65 and 3.0a0, respectively. In (b), the dash-dot curveshows the self-consistent effective potential when no fluid is present and r3 = 2.65a0 Thedistribution 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(r8 2.65a0).projection g’(z) [Fig.3.12b] is a measure of the local orientational polarization of thesolvent. For inert surfaces, g1(z) = 0 since symmetry requires that the solvent has nonet polarization. Forg2(z) [Fig. 3.12c], the most striking feature is the change of signnear 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) is given simply by1 Onnp(O,z) = (_) O()Pn(cos) . (3.4.26)Orientational probability densities for four wall-solvent separations are compared inFig. 3.13. The probability density near an inert surface is included. For r3 = 2.65a0,the complete pair distribution function, in the interval 0.5d z — z 3.5db, is plottedChapter 3. Fluids near planar surfaces 691 2 3z* 0Figure 3.12: RHNC results for the projections of the wall-solvent pair correlation functionfor a dipolar hard sphere fluid in contact with a jellium metal. The dotted, solid anddashed curves are as in Fig. 3.11. The dash-dot curve is for an inert surface (r5 = oo).0—5—10—150011—10(b)__L1 2 3z*0.5 0.6108642010500.60.40.20—0.2—0.41 2 3z*1111111111 IIIg 10 (c)1111111 I 111110—1—21 2 3z*0 1 2 3z*Chapter 3. Fluids near planar surfaces2.01.51.00.50.00 cosO 12.01.51.00.50.0I III II II Ill III I Ip(, O.8d)700 cosO 12.01.51.00.50.0—1 0 cosO 12.01.51.00.50.0—1 0 cos9 1Figure 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.I I I(b)—1 —1I 111111 111111 III-p(6, 1.2d5) (c) ZI ii I 11111 I I iii ii II I I I I I I I I I I I I I I I I IEp(,1.5d8) (d)—l I I I I I IChapter 3. Fluids near planar surfaces 71Figure 3.14: The wall-solvent pair distribution function g(O, z) near a jellium metalwall with i’3 = 2.65a0 The solvent parameters are as in Fig. 3.7.gws(,z)10.00z/d e5 cos6Chapter 3. Fluids near planar surfaces 72in Fig. 3.14. The structure of the dipolar fluid near an inert surface was discussed earlier (see page 56 and Fig. 3.5). In contrast, near the jellium-liquid interface the solventorientations are highly asymmetric. At contact, the dipoles orient preferentially alongthe surface normal vector (i.e., the positive end of the dipole points away from the metalsurface) giving rise to significant surface induced polarization of the solvent. This is qualitatively consistent with calculations for dipolar monolayers [67,68]. As we move furtherfrom the surface the most probable dipole orientation oscillates. For example, at O.8dand 1.2d5 [Figs. 3.13b and 3.13c] the dipoles have rotated more towards the surface, butat 1.5d [Fig. 3.13d] they again prefer to align along the surface normal. We note that thesolvent structural features again decay rapidly with distance from the surface althoughthe correlations are a little longer ranged than in the inert case (see Fig. 3.14).Electrostatic potentials obtained for r3 = 2.65a0 are plotted in Fig. 3.15. The totalpotential across the jellium-liquid interface (solid curve) is the sum of the contributionsfrom the electron density (dotted curve) and from the solvent dipoles (dashed curve). Wesee that as expected the net potential drop across the interface is reduced by the presenceof the solvent. Structural features related to solvent granularity and orientational orderare also evident.The potential drop, /q = 4(z oo) — q!(z = 0), across the jellium-dipolar fluidinterface is given by=471-= S Sf g’(z’)dz’3 0+4ej 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 r3 values are given in Table 3.4. The potential drop across thejellium-vacuum interface is also included. We note that the results obtained for dipolarChapter 3. Fluids near planar surfaces100500321020100—10—200 2 473—5 0 5z* 6 zFigure 3.15: The reduced electrostatic potential across the interface when r3 = 2.65a0(a) The net potential for the jellium slab in a dipolar hard sphere fluid (parametersas in Fig. 3.11) is represented by the solid curve and is the sum of contributions fromthe jellium (dotted curve) and the orientationally polarized fluid (dashed curve). Thenet 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) Theelectrostatic 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 effectof varying the solvent diameter and/or density on the total potential drop is shown inTable 3.5. If the solvent diameter is reduced while holding p’ constant, the dipoles canapproach the jellium edge more closely and the potential drop decreases. Also, as wewould 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 r3 = 2.65a0 and ata continuum conducting surface using the cavity and superposition approximation areplotted again in Fig. 3.16 for comparison. The dielectric wall can only influence thesolvent’s preference to point either parallel to the surface or in and out from the surfacebut can not induce a net polarization. Hence, only even-n projections are shown inChapter 3. Fluids near planar surfaces 74Table 3.4: The self-consistent electrostatic potential drop across the jellium-vacuum andjellium-solvent interfaces for z = 32A. When present, the liquid consists of dipolar hardspheres with p = 0.7, = d = 2.8A and T = 298K.L4 (volts)j ellium-vacuum j ellium-solvent2.00 —6.787 —6.3012.65 —3.275 —2.9503.00 —2.318 —2.045Table 3.5: The self-consistent electrostatic potential drop across the jellium-solvent interface for r3 = 2.65a0 and z = 32A. In all cases the liquid consists of dipolar hardspheres with t = at T = 298K.p d8 (A) Liq!’ (volts)0.663 2.8 39.9 —2.9690.7 2.6 50.3 —2.8980.7 2.8 50.3 —2.9500.7 3.0 50.3 —3.0070.75 2.8 77.7 —2.940Fig. 3.16. In both approximations the contact density increases from the inert surfacevalue 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 allcases. In addition, we note that the jellium curves for even-n projections (short-dashcurves) resemble the inert wall results (solid curves) except very near contact where theyare most similar to the cavity approximation (long-dash curves).3.4.4 Results for electrolyte solutionsIn this section we will discuss electrolytes near inert and metallic surfaces.Chapter 3. Fluids near planar surfaces 7510286 04________—22—40 00.11.00.00.5—0.10.0—0.20 0Figure 3.16: Comparison of the projections occuring in the wall-solvent pair correlationfunction for a dipolar hard sphere fluid near a jellium metal (r8 = 2.65a0) surface and adielectric surface. The short-dash curve is the jellium metal result. The long-dash anddash-dot curves are the results near a dielectric surface using the cavity and superpositionapproximations, respectively. The solid curve is result near an inert surface.000 (a)8.0-\\4.0001 2 3z* 1 2 3z*1 2 1 2Chapter 3. Fluids near planar surfaces 76Electrolyte solutions in contact with an inert surfaceThe first three wall-solvent projections obtained for a O.945M CsF solution near an inertwall are shown in Fig. 3.17. The results for pure dipolar hard spheres are included. Ingeneral the solvent structure is only slightly modified by the electrolyte. As the anionand 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 polarized by the presence of an uncharged inert wall and the potential drop across theinterface is zero. However, if the ions are of unequal size the symmetry is broken and asmall surface-induced solvent polarization results. For such systems, dipoles will have apreference to point either into (g1 > 0) or out of (g’ < 0) the surface. For CsF at0.945M (Fig. 3.17) we see that g’ is positive near contact indicating that the dipolesprefer 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 polarization of the solvent. At higher concentrations, the decay to zero becomes more rapidas 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.945Mg1(d/2)is negative for Lii and NaC1 and positive for CsF. This is as one would deduce from thesign of the smaller ions; i.e., small positive ions located between the surface and thesolvent particles will direct the dipoles outward and small negative ions will do thereverse. For CsI, where both ions are larger than the solvent but Cs+ is smaller than1 (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’(d8/2) is not so easily predictedat high concentrations. For example, for NaC1 solutions g’(d/2) becomes positive atChapter 3. Fluids near planar surfaces 772.00.021.50.001.0—0.020.50.0—1.0Figure 3.17: Projections of the wall-solvent pair distribution function for inert walls. Theparameters are as in Tables 3.2 and 3.3. The solid curve is the pure solvent result andthe dashed curve is for O.945M CsF.0 1 23 4z* 0 1 2 30 1 2 3 4 zChapter 3. Fluids near planar surfaces 78Table 3.6: Contact values of the wall-solvent pair distribution function near an inertsurface.model molarity g° g°” g°22pure solvent 0 4.39 0 —2.58LiT 0.945 4.89 —0.073 —2.41NaC1 0.453 4.50 —0.017 —2.390.945 4.58 —0.024 —2.273.78 4.87 0.007 —1.927.56 5.23 0.114 —1.78EqEq 0.945 4.43 0 —2.18CsF 0.945 4.73 0.025 —2.37CsI 0.945 5.46 0.005 —3.06higher concentrations indicating that correlation effects become more important as theconcentration is increased.Apart from these relatively small polarization effects, the addition of salt does notstrongly influence the interfacial solvent structure. We see from Fig. 3.17 that the even-nprojections are only slightly influenced at 0.945M and in fact salt effects are not verysignificant even at 7.56M. The value of g°(d/2) (Table 3.6) is increased for all systemsconsidered, most significantly for the 7.56M NaC1 and 0.945M CsI (where both ion diameters are larger than d) solutions. Also,g2(d/2) decreases in all cases except for0.945M CsI. The more relevant quantity for comparisong(d/2)/g°(d/2), which isproportional to (P2 (cos 0)) per particle, decreases in all cases indicating that in the presence of salt the dipoles have a reduced tendency to align parallel with the surface. Thisis obviously consistent with the polarization effects discussed above.The wall-cation and wall-anion distribution functions for a few salt solutions at0.945M are shown in Fig. 3.18. At 0.945M the anion and cation distributions are nearlyindependent of the coion diameter; for example, the distributions of Cs in CsF (notshown) and CsI are similar. Also note, that in the simple dipolar solvent the distributionChapter 3. Fluids near planar surfaces 795 54 43 32 21 10 00 1 2 3z* 0 1 2 3z*Figure 3.18: Wall-ion pair distribution functions for O.945M electrolyte solutions in contact with an inert surface. The parameters are as in Tables 3.2 and 3.3. In (a) the curvesare: (—) Na in NaC1; (----) Li in Lii; (----) Eq in EqEq; (— —) Cs in CsI. In (b) thecurves 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 significantfeature evident in Fig. 3.18 is a sharp increase in the contact values as the ion diameteris increased. This comes about simply because the smaller ions (which can approachthe solvent particles more closely) are better solvated than the larger ones in the dipolarsolvent.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 dielectricconstants of solution are also included. If the ion pairs are of equal diameter but oppositecharge (e.g., RbF and NaC1), IL4 has the same magnitude but opposite sign. We notethat although iI and ILci range between 0 - 146mV, the ion and dipole contributionsare always opposite in sign and cancel such that qS is always less than 5OmV. Aswe would expect, ILq generally increases with the difference in ion size. Further, Lqpasses through a minimum at ‘S-’ 0.945M. It is also evident that ILq is not simply relatedChapter 3. Fluids near planar surfaces 80Table 3.7: The potential drop in volts across inert surface-electrolyte solution interfaces.The dielectric constant is also included. In all cases p = 0.7, = /ä and q* = 8.model molarity c Lqpure solvent 0 50.3 0 0 0LiT 0.945 30.23 —0.060 0.048 —0.012NaCl 0.453 37.22 —0.034 0.036 0.0020.945 30.64 —0.032 0.029 —0.0033.78 17.23 —0.068 0.042 —0.0267.56 11.47 —0.146 0.099 —0.047EqEq all 0 0 0RbF 0.453 37.22 0.034 —0.036 —0.0020.945 30.64 0.032 —0.029 0.0033.78 17.23 0.068 —0.042 0.0267.56 11.47 0.146 —0.099 0.047CsF 0.945 30.94 0.042 —0.036 0.006CsI 0.945 35.36 —0.019 0.014 —0.005ILi 0.945 30.23 —0.060 0.048 —0.012to the dielectric constant of solution. For example, although f decreases with increasingsalt concentration, Lq does not.Electrolyte solutions near a metal surfaceWe begin by considering the first three projections of the metal wall-solvent pair distribution 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 included. The full pair distribution function for 0.945M and 7.56M NaCl are shown inFigs 3.21 and 3.22, respectively. We note that in all cases the general structure is qualitatively similar to that of the pure solvent. However, the salt does induce significantperturbations which depend on both the cation and anion diameters. The g’ and g2projections are of most interest here.Chapter 3. Fluids near planar surfaces 812.0 0.51.5 0.0—0.51.0—1.00.50.0—1.0Figure 3.19: Projections of the wall-solvent pair distribution function for NaC1 solutionsin contact with a jellium slab. The solid curve is the pure solvent result. The dotted anddashed curves are for O.945M and 7.56M solutions, respectively. The parameters are asin Tables 3.2 and 3.3.0 1 23 4z* 0 1 2 30 1 2 3Chapter 3. Fluids near planar surfacesFigure 3.20: Projections of the wall-solvent pair distribution function for RbF solutionsin contact with a jellium slab. The solid curve is the pure solvent result. The dotted anddashed curves are for 0.945M and 7.56M solutions, respectively. The parameters are asin Tables 3.2 and 3.3.820.50.0—0.5—1.02.01.51.00.50.0—1.00 1 2 3 4 z* 0 1 2 30 1 2 3 4 zChapter 3. Fluids near planar surfaces 83Figure 3.21: The wall-solvent pair distribution function g5(O, z) for the O.945M NaC1solution in a dipolar solvent and near a jellium metal wall.gws(O,z)10.0z /d5-j.cosChapter 3. Fluids near planar surfaces 84Figure 3.22: The wall-solvent pair distribution function g(O, z) for the 7.56M NaC1solution in a dipolar solvent and near a jellium metal wall.o.o gws(O, z)11z Id5-.-cos 6Chapter 3. Fluids near planar surfaces 85The g’ projections show the strong dipolar polarization induced by the metal surface. The large negative values at contact indicate that the dipoles strongly prefer toorient perpendicular to the surface. From Figs 3.19 and 3.20 we see that this effect ismagnified as the salt concentration is increased. Further from the surface g’ oscillatesand approaches zero. For the pure solvent these oscillations are centered about zeroand quickly decay. At 0.945M, g’ is longer ranged and oscillates about negative valuesfor both NaC1 and RbF solutions. At 7.56M the solvent is more strongly ordered nearcontact, but g’ is highly screened and decays very rapidly. The g2 projections exhibit qualitatively similar trends and become more positive at small separations as theconcentration is increased. These wall-solvent results will be discussed further after weconsider the wall-ion distributions.The short-range potential arising from the jellium is attractive to anions and repulsiveto cations. Thus if the solvent were not explicitly present (e.g., a primitive model) wewould 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 mostobvious features, common to all solutions considered, are the large cation peak near themetal surface and a major anion peak located a full solvent diameter from contact (nearz = 1 + (d_/d3)/2). This structure comes about because of the strong ordering of thesolvent particles near the metal surface and can be understood by considering the iondipole interactions. In bulk solution (see Fig. 3.24), the dipoles in the first solvationshell prefer to be directed towards an anion and away from a cation. As an ion is movedtowards the surface, an increasingly strong force tends to orient the dipoles normal tothe surface. For anions the orienting fields of the surface and of the ion act more or lessconstructively on dipoles located between them for separations d8 + d_ /2. (Note thatthe surface field is very short ranged and also partially screened on the far side of theanion.) Closer to the surface the fields act destructively and destabilize the solvationChapter 3. Fluids near planar surfaces 86shell. These observations explain all but the very short-range behaviour of the aniondistributions shown in Fig. 3.23. At small separations the direct metal-anion potentialis strongly attractive and consequently the anion density increases rather sharply atcontact. 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 ‘‘ O.6d3) the polarization of the dipoles creates a new region of stability (see Fig.3.24). However, at still smaller separations the direct metal-cation interaction begins todominate 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 followthe same trend as in Fig. 3.18 for an inert surface (i.e., the contact density increasessharply for larger ions as a result of the weakened ion-solvent interaction). We also notethat cation “dewetting” due to the repulsive interaction with the metal occurs only forsmaller ions (d 1.2d), an indication of the range of the jellium potential. The anionpeak at “-‘ d3 + d_/2 is relatively insensitive to the ion diameter, but increases in heightas the cation density near the surface increases.Cation and anion distribution functions obtained for different values of r3 are shownin Fig. 3.25. We observe that the qualitative features are similar for the different valuesof r3, but that the structure sharpens considerably as r3 is decreased from 4.0a0 to 2.0a0This reflects the increasing degree of orientational order in the solvent contact layer asthe metal-dipole interactions increase in strength.The results discussed above show that at O.945M the ion-solvent interactions largelydetermine the ion-wall distribution functions. The metal orders the solvent which in turndirects the ions. The direct ion-metal and ion-ion interactions tend to be dominated bythe ion-solvent forces. The only important exception to this occurs for small ions nearcontact where the direct ion-metal interaction sharply repels cations and attracts anions.Chapter 3. Fluids near planar surfacesFigure 3.23: The wall-cation (—) and wall-anion (---.) pair distribution functions forseveral O.945M electrolyte solutions in contact with a jellium slab. The parameters areas in Tables 3.2 and 3.3. In (h), the wall-cation distribution functions for RbF, CsI, CsFand ILi are shown on an expanded scale.8711111 liii 11111 111111 I I(a)LiI.I I I I I I j I I I I I I I I Ii II II •iI I i i i i I i i i I i i i I i i i420202020III\,IIIIIIii;IiiiIII I I Iii I I I I i i I i i I i20202020:L’>’IIIIIIIIIIIIIIIIIIIIIIII0 1 2 3I 1111111 I 11111(h)—CsF —RbFI I I I —I—1 I I I1 2 z1000 1 2 3 4z* 0Chapter 3. Fluids near planar surfaces 88G00eØBulkFigure 3.24: An illustration of the distortions of the solvation shell due to the jelhum-induced interfacial polarization.Results for higher concentrations (up to 7.56M) of NaC1 and RbF are shown in Figs 3.26and 3.27. Although screening and other concentration dependent effects can be seen inthe 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 withthe solvent structure illustrated in Figs. 3.19 and 3.20. The excess anion concentrationnear z — 1.5 induces an additional polarization of the solvent and hence g1 is negativein this region and approaches zero from below. Further, for 0.945M RbF (d_ = 0.84d3)g’ is more negative than in the case of 0.945M NaC1 (d_ = 1.16db). For 7.56M RbFeChapter 3. Fluids near planar surfaces 895H 54- 432z 21- 10_I 00 1 2 3z* 0Figure 3.25: Wall-ion pair distribution functions for NaC1 solutions at O.945M. Thedotted, solid and dashed curves are for r5 = 4.0a0, 2.65a0 and 2.0a0, respectively.the region near z’ = 2 is also of interest. Here g’ has a positive peak indicating thatdipoles prefer to orient towards the surface. This peak is presumably associated with thestrongly bound first solvation sphere of the small F— ion.The potential drops across the metal-solution interface are given in Table 3.8. Thecontributions from the dipoles, ions, and jellium, j, are included. Thenet contribution from the electrolyte solution is also included and denoted by LSd. Weobserve that the presence of salt increases zg,j by about a factor of three above thepure solvent value. However, this increase is largely canceled by zj such that the totalsolution contribution to remains near that obtained for the pure solvent. Also, we notethat 4dj is roughly an order of magnitude larger than corresponding values for an inertwall (see Table 3.7). The partial cancellation of /qf’d and qj serves to emphasize thatall effects must be taken into account in a self-consistent manner if accurate theoreticalresults are to be obtained. The ion and solvent contributions are strongly coupled andone cannot focus separately upon one or the other as is frequently done in double layertheories. Several other “trends” are evident from Table 3.8. For example, we see thatg (a)Il1 2 3z*Chapter 3. Fluids near planar surfaces 904 43 32 21 1o 0-j-0 1 2 3 4z*3210Figure 3.26: The wall-ion pair distribution function for NaC1 solutions near a jelliumsurface. The curves are labeled as in Fig. 3.23.Lf is not strongly concentration dependent. Also, the change in Lq with respect tothe jellium-vacuum interface, tends to decrease with increasing concentration,although for RbF a minimum occurs at r’.i 3.78M. Furthermore, metal induced chargeasymmetry 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 potentialdrop 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 NaClsolution are plotted in Fig. 3.28. Similar results were obtained for the other systemsconsidered. As in the pure solvent case (see Sec. 3.4), the jellium contribution decaysvery rapidly on the solution side of the double layer and is more or less negligible forz 0.6. However, it is the rapid decay of this potential which gives rise to the strongelectric field orienting the dipole in the contact layer. The separate contributions of theI I I I I I I I I I I I I I I I I I I(c) 7.56Mr’ii;iiiiiiIi,ii:0 1 2 3 4z*6420Figure 3.27: The wall-ion pair distribution function for RbF solutions near a jelliumsurface. 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 differgreatly from the pure solvent curve also shown in the plot. This demonstrates again thestrong coupling of the ionic and dipolar potentials emphasized above.ostaChapter 3. Fluids near planar surfaces 918 86 64 42 20 00 1 2 3 4z*0 1 2 3 4z*Chapter 3. Fluids near planar surfaces 92Table 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 r5 = 2.65a0model molarity LS4 /q’d1 LLqfvacuum 1.0 0 0 0 —3.275 —3.275 0pure solvent 0 50.3 0.412 0 0.412 —3.362 —2.950 0.325LiT 0.945 30.23 0.955 —0.550 0.405 —3.375 —2.969 0.306NaC1 0.945 30.64 1.047 —0.635 0.412 —3.375 —2.964 0.3113.78 17.23 1.191 —0.836 0.355 —3.384 —3.028 0.2477.56 11.47 1.013 —0.716 0.297 —3.379 —3.082 0.193EqEq 0.453 37.57 1.003 —0.574 0.429 —3.373 —2.944 0.3310.945 31.08 1.187 —0.774 0.413 —3.382 —2.969 0.3063.78 18.17 1.310 —0.925 0.385 —3.387 —3.002 0.273RbF 0.453 37.22 1.199 —0.768 0.431 —3.374 —2.943 0.3320.945 30.64 1.376 —0.956 0.420 —3.382 —2.962 0.3133.78 17.23 1.438 —1.008 0.430 —3.412 —2.982 0.2937.56 11.47 1.299 —0.852 0.447 —3.421 —2.974 0.301CsF 0.945 30.94 1.524 —1.099 0.425 —3.387 —2.962 0.313CsI 0.945 35.36 2.011 —1.592 0.419 —3.396 —2.976 0.299ILi 0.945 30.23 1.520 —1.084 0.436 —3.393 —2.957 0.318Figure 3.28: The self-consistent mean-field electrostatic potential for O.945M NaCl incontact with a metal slab (z = 32A). The long-dash-short-dash, dash-dot and long-dashcurves are the jellium, dipole and ion contributions, respectively. The short-dash curve isthe sum of the dipole and ion contributions and the dotted curve is the dipole contributionobtained 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 3. Fluids near planar surfaces10050093.IIiIIIIIIIIIIIIIIIIII.!(b)—--321040200—20—40—2 —1 0 1 2 z 0 1 2 3 4z*Chapter 4Fluids between planar surfaces4.1 IntroductionIn chapter 3, the effect of surface properties on the interfacial fluid structure was examined. This represents one approach towards understanding fluids in inhomogeneousenvironments. 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, phasebehaviour and pressure of the fluid. In these confined systems, integral equation theorieshave been used primarily to study primitive model electrolytes (spherical ions in a dielectric medium) [109—112]. Lennard-Jones fluids have also been considered near an isolatedinert surface [119] using homogeneous integral equation theories and between soft planarwalls using the inhomogeneous theories [113,114]. In addition, there has been an extensive 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] havealso been investigated.In this chapter we consider fluids confined between planar walls and study the separation-induced fluid structure and fluid-mediated force between surfaces. General reductions of the wall-wall and wall-solvent OZ equations are obtained for molecular fluidsconfined between planar surfaces. In addition, a general method is given for simulating metastable fluids beyond coexistence. Results are discussed in two parts. In the94Chapter 4. Fluids between planar surfaces 95first, a survey of various integral equation methods for different model fluids betweensurfaces is given. These results are compared with computer simulations to verify theaccuracy of the HNC and RHNC theories for different systems. In the second part, grandcanonical Monte Carlo is used to simulate subcritical Lennard-Jones fluids between inertwalls as a function of the wall-wall separation and to study capillary evaporation. Thefluid, which is liquid in the bulk, evaporates at small separations. These results shedlight on separation-induced phase changes. The possible relevance of these results to theexperimentally measured attraction between hydrophobic surfaces is discussed.The remainder of this chapter is organized as follows. The reduction of the planarwall-wall OZ equation for molecular fluids is presented in Sec. 4.2.1. In Sec. 4.2.2, amethod of calculating the wall-particle pair distribution function for a fluid confinedbetween planar walls (using the homogeneous OZ equation) is presented. In Sec. 4.3, weelaborate on the grand canonical Monte Carlo method and apply it to a fluid of sphericallysymmetric particles between planar walls. Here, general methods of including long-rangeinteractions for Lennard-Jones particles and of simulating metastable fluids are included.Results for hard sphere, Lennard-Jones and dipolar hard sphere fluids using the HNCand RHNC theories are presented in Sec. 4.4 and compared with simulation results. Ageneral discussion of the applicability of the integral equation theories near surfaces andbetween surfaces is presented. Results for equilibrium and metastable simulations ofLennard-Jones fluids near coexistence are given in Sec. 4.5. These results are comparedwith a simple theory and experimental data.Chapter 4. Fluids between planar surfaces 964.2 Theory4.2.1 The interaction free energy between two walls separated by a molecularfluidIn this subsection, we consider correlations between planar parallel walls immersed in amolecular 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 by1 and 2 respectively. For smooth surfaces, the wall-wall pair functions vary only withthe separation 1 and we use a coordinate system in which z1 = 0 (see Fig. 3.1) at the leftsurface and 1 = 0 at wall-wall contact. The reduction of the OZ relation follows from theanalysis of the wall-solvent OZ equation presented in Sec. 3.2.Using the identity c(32) =c7(23) [see Eq. (2.1.9)], the OZ equation is first writtenin the equivalent formh(l) — c(l)=J h(O, z3)c7(O, (z3 — l)) d(3). (4.2.1)Expanding the wall-solvent pair correlation function as in Eq. (3.2.1) and the wall-solventdirect correlation functions asc7(23) = c(l — Z3)—)= (_)nf ( ) c(l - Z3) R(, (4.2.2)Chapter 4. Fluids between planar surfaces 97where Eq. (A.8b) has been used, we get’fOmm\ /Onn\0mrn0nno ) o o8ir2 mj..7W0mm / ‘ Onn /p00 i00 [00x / dx / dy hOi;w.yIZ3) Co,,;w..y(l — z3) dz3J—oo J—oo J—oox(—)JRm(c *(ç3)dfz3. (4.2.3)o, 3i-oiApplying the orthogonality property of the Wigner spherical harmonics, Eq. (A.9), thisreduces immediately toOn \12,°mn jiui(1) = ( o)i,(2n+1)a0 ico p00iOnn I \ Onn ‘7x / dx / dy 10;w73)COz,;wyt — z3) dz3 . (4.2.4)J—oo J—oo 1—00Here we notice that = h — c,, diverges with the surface area. Hence, usingEq. (2.2.21), we introduce the potential of mean force per unit area asw(l)F(l)= f dxfdy= U(l)— kBTB(l) — kBTN(l) (4.2.5)where is the direct wall-wall pair potential per unit area, is the wall-wall bridgefunction per unit area and- l)N(l)— fdxfdy(_)n_v= (2n+1) L ooo]2=Ffonn(0flN00xJ hOni ‘ Onn /op;wy(Z3) COv;wl — z3) dz3 . (4.2.6)-00For hard walls, the bridge function B(l) 0 when 1 1 (i.e., when the walls areseparated 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 incorrectsign on the lower index of homm in Eq. (4.15) of that reference but has no effect on the results presentedO ;wyfor dipolar fluids.Chapter 4. Fluids between planar surfaces 98hg7(z) and c(z) for positive and negative values of z. For we haveh7(z) = nO&,O, Z < 0 . (4.2.7)Onn Onn Onn(+) Onn(—) Onn(+)For we use c = h — where 770v;wy (z) + iOu;w7 (z). 770v;wy (z) isobtained from the inverse Fourier transform of Eq. (3.2.18) for both positive and negativeOnn(—)values of z and 1lOuw7 (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 surroundingthem. When the walls are at infinite separation, no net force acts upon them; at smallseparations, a net force results from the pressure of bulk fluid pushing the walls intocontact and from fluid confined between the walls pushing them apart. A positive valueof F corresponds to a repulsion between surfaces. The excess potential of mean forceper unit area is given byF(l) = F(l) — U(l) , (4.2.8)and represents the fluid-mediated contribution to F(l). Notice that, since h,(z) andc,(z) —* 0 as z —* oc, we have F(l) and F(l) —* 0 as 1 — oc, and neither the bulkfree energy nor the surface energies contribute to F(l) or F(l). The net pressurebetween the walls, mediated by the fluid, is given asdF (P(l) = — . (4.2.9)Specifically, the net pressure is the normal component of the internal pressure at one walldue to fluid P(l) minus the bulk pressure Pb,P(l) = P1(l) — Pb . (4.2.10)Chapter 4. Fluids between planar surfaces 994.2.2 The wall-solvent-wall correlation functionsThe wall-wall OZ equation only describes correlations between surfaces. The three-bodywall-solvent-wall correlation functions are also of interest and can be defined from theinhomogeneous OZ equationh(l,3) — c(l,3) = ps(4) h(l,4) 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 definesa space-fixed coordinate frame. In this approach, the inhomogeneous direct correlationfunction c6(4, 3) and the singlet density distribution ps(4) near an isolated wall are required as input, and a new solution is obtained for all wall-wall separations 1. Thehomogeneous wall-particle OZ equation can be modified to obtain similar results withrelatively 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, z3)is replaced with the wall-particle-wall potential,uw,3(3; 1) = u(O, z3) + u(3,0, 1 — z3) , (4.2.12)and the resulting equations are solved for a given separation of the walls 1. The relevantOZ equation is now written ashp(3; 1) — c(3; 1) = 2 Jh7(4; 1) 43) d(4), (4.2.13)where the two walls are treated as a single confining particle and the coordinate systemdefined in Sec. 4.2.1 has been used2. The only modification in the reduction of Eq. (2.2.7b)2Comparing Eqs (4.2.11) and (4.2.13), we see that h(3; 1) h(l, 3), c(3; 1) c(1, 3) and theinhomogeneous singlet density p(3) psgw(13) in an exact theory. Hence, it is interesting to point outthat substitution of the results from the homogeneous calculations into Eq. (4.2.11) provides a relativelysimple means of calculating the inhomogeneous correlation functions.Chapter 4. Fluids between planar surfaces 100(see Sec. 3.2) arises because h,3 = —1 inside both surfaces,hg(z; 1) = z < 0 or z > 1, (4.2.14)and in place of the definitions for and we useh(3; 1) = h,(3; 1) W(z3 — 1/2) , (4.2.15)h(3; 1) = W(z3 — 1/2) — 1, (4.2.16)where(1, z<l/2,W(z) =. (4.2.17)l0, z>l/2.For z > 0, relations of the formz; 1) = 2irpyf r2S[(z + z)/r] c(r) drZyq-f-Z+2rp7jr2Sn[(zw— z)/r] c(r) dr[ (2n + 1) (n1 n 12rnl‘ (2n1 + 1) 0 0 0x j°° e(k; 1) k) dk , (4.2.18)are obtained. Before proceeding, we must consider closure relations for Eq. (4.2.18). Forspherical particles, the best method is direct application of Eq. (2.2.17) for some approximation of the bridge functions. For molecular fluids between two walls, the approach ofFries and Patey [see Eq. (3.2.24)] is no longer convenient since u(3; 1) —* oc as z3 —÷ +ooand the assumption that c(3; 1) —* 0 as z1 —+ oc is no longer valid. In this case, wefollow an analysis similar to Caillol [33,34] and take the derivative of Eq. (2.2.17) withrespect to the Euler angles of the fluid particle In general, this givesJ2h(12) = J2w(12), (4.2.19)Chapter 4. Fluids between planar surfaces 101orJ2c(12)= h(12)kBTJ2w(12) —-J2u(l2) +J2b(12), (4.2.20)where J2 represents a linear combination of partial derivative operators with respect tof2. Choosing J2 as the step operator (see Eq. (3.45) of Ref. [5])8 1 \ 81== ie cot +O+ ij (4.2.21)8X2 sinand applying J to the wall-particle basis functions yieldsJ±R,1(2)= /n(n + 1) — v(v + 1) R(i2 . (4.2.22)Expanding the closure Eq. (4.2.19) in rotational invariants for the wall-particle pairfunctions and applying the step operator yieldsJ±h(12) = f0ThTh (°)Vn(n+1)_v(v+1)h(z)R±i(00= —/n(+1)—vv±1)fl2xf°’’ (Oni fOfl2fl2 (On2 n201Oninixg0,1;(z) w(z)R1(2)R21(fZ2) . (4.2.23)Using either Eq. (A.9) or Eq. (A.11) then gives relationships of the formh0v;w(z) = —(_)fl+Thl+fl2+V(2n + 1)n2(n + 1) — v2(v + 1)Blum___________________________________kBT /n(n + 1) — v(v ± 1)fl2In1 n2 n\ In1 n2 n‘ Onnxo 1±1v+1 1oz>O. (4.2.24)A similar relation can be obtained starting with Eq. (4.2.20). When n = 0, the projection1000Oow3 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 forcez) = e_()hhJT — 1 (4.2.25)Chapter 4. Fluids between planar surfaces 102for the choice (0, 0, 0). This givesh(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 OZequation.4.3 Grand canonical Monte Carlo method for Lennard-Jones particlesThe grand canonical Monte Carlo method was introduced in Sec. 2.4. In this section, itis considered in detail for a fluid of Lennard-Jones particles confined between two planarwalls. The computer simulations are performed in a central cell of dimensions L x L x 1 inwhich periodic boundary conditions are imposed along the x and y axes and two planarwalls are located at z = 0 and z = 1. The solvent-solvent interactions are given by theLennard-Jones potentialu(r) = 4€[(u)12 — (7)6], (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 additivepotential [see Eq. (4.2.12)]. The wall-solvent interaction is modeled predominantly bythe hard wall potentialI co, z</2,u(z) = (4.3.2)(0, z>a/2.However, we also consider the Lennard-Jones 10—4 potentialI \1O ,2___ ___u(z) =2Psur g + u/2)—(z + u/2) (4.3.3)where Psin is the surface density. Equation (4.3.3) represents the interaction energy of aLennard-Jones particle in the fluid (at a separation z from the surface) with a uniformChapter 4. Fluids between planar surfaces 103distribution of Lennard-Jones particles in a plane centered at z = —u/2 [122]. To get theparticular 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 fluidsuggests a cylindrical truncation of u5(r) instead of the usual spherical truncation usedfor homogeneous bulk fluids. In this way, the potentialI ss(v”5 + z2) , S < Stu(r) = u(s,z) = (4.3.4)10, S>St,is defined where st L/2 is the cylindrical cut-off radius. The cylindrical truncationscheme has the advantage that the mechanical properties of the fluid can be correctedfor long-ranged structural effects in a straightforward manner.The average configurational energy of the central cell,U=u(s, z) + u(z; 1)) + U, (4.3.5)i j.<iwhere zj zj — z, sj =—s and the angular brackets denote a grand canonicalaverage, is calculated using the minimum image convention. The long-range correctionU is given byU= 2Jv 1J2p)(ri,r;l) {uss(Vs?2+z?)_u(si2zi)] , (4.3.6)where the integral over r1 is restricted to the volume V = L21 of the central cell (f r2in unrestricted) and p)(r1,r2; 1) is the inhomogeneous solvent-solvent pair distributionfunction between two planar surfaces separated by a distance 1. In Eq. (4.3.6), thetruncation cylinder is centered on particle 1 such that s = 0 and only values of s12 > .scontributes.For the inhomogeneous fluid, the pair distribution function is required in order toevaluate the correction U. If the truncation radius s is chosen sufficiently large, spatialChapter 4. Fluids between planar surfaces 104correlations can be ignored such thatp(r1,r2; 1) = p(’)(z1;1) p(’)(z2;1) g)(r,r2; 1)p(’)(zi; 1) p()(z2;1), s12 > s , (4.3.7)where p(’)(z; 1) is the singlet distribution function. Further, if p(1)(z; 1) is assumed to beuniform then the assumption used by Scheon and coworkers [107] is obtained. A slightlybetter approximation isp’>(z; 1) po(z; 1) (4.3.8)where, for example, po(z; 1) is the grand canonical Monte Carlo singlet distribution function when U = 0. This is the approximation used in the present calculations. Thecorrection to the energy can then be evaluated either in integral form,U j po(zi; 1)1(zi) dr1, (4.3.9a)or as an ensemble average,(4.3.9b)where1(zi)= f po(z2;1) [u(s + z?2) —u8(s2,z12)] dr2. (4.3.10)Explicitly, using the asymptotic limit of the Lennard-Jones potential, u(r) ‘-.‘ —4(u/r)6,one obtainsU _2(j po(z2; ( + _ Z)2) dz2) . (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 millionconfigurations on a grid of /z = 0.Olu. This tail correction is included when attemptingparticle insertions or deletions.Chapter 4. Fluids between planar surfaces 105The contact value theorem [3] is used to calculate the internal pressure and relatesthe 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, P1(l) was obtained from the average contact value at the two walls,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 =kBT in zo3, as discussed in Sec. 2.3. The equilibrium state of the system is obtained afterattempting 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 calculations. It is particularly useful for studying inhomogeneous systems and for studyingphase behaviour. In this work, we are interested in the pressure between walls as a function of separation at constant chemical potential, and further, in this pressure near phasetransitions. Conventional grand canonical Monte Carlo is suited for finding the equilibrium coexistence point of two phases, but it can be adapted to study the metastablefluid beyond the coexistence point, and to find the spinodal point itself. Consider a fluidclose to the liquid vapour coexistence curve where a phase transition can be inducedby small variations in V or T. In this region, there exists two minima in the freeenergy as a function of density at constant res, V, and T. The global minimum corresponds to the equilibrium state and a second minimum corresponds to the metastablestate. In general, for metastability to exist, the two minima are separated by a freeenergy barrier. The disappearance of this barrier corresponds to the spinodal point. Themetastable 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 106Since the two branches (the equilibrium state, and the metastable state) correspond todifferent densities, one can confine the fluid to the metastable branch by imposing athreshold volume-averaged density, 1öth, and by rejecting any attempted particle insertions or deletions that would cross this threshold. By monitoring fluctuations in thenumber of particles during the simulation, one can ascertain whether the average densitycorresponds to the metastable state, at the local minimum, or whether it correspondsto the threshold density. The spinodal point in which cavitation occurs is approximatedby the minimum separation at which a local metastable minimum is observed; the spinodal point in which condensation occurs is approximated by the maximum separation inwhich a local metastable minimum is observed.4.4 Results for pure fluids between planar wallsThis section is concerned with the structure of simple fluids and dipolar hard spherefluids confined between walls and the solvent-mediated wall-wall pressure. The pressurebetween hard walls immersed in a hard sphere fluid [calculated from Eq. (4.2.9)] isshown in Fig. 4.1 as a function of the reduced wall-wall separation 1* = l/d > 1. Here,results obtained from the HNC theory, RHNC theory and grand canonical Monte Carlosimulations [123] are compared for a reduced hard sphere density of p = 0.681. TheRHNC calculations were carried out using the Verlet-Weis [27] and Henderson-Plischke[31] parameterization for the hard particle-particle and hard wall-particle bridge functions, respectively. The direct hard wall-wall bridge functions are zero when 1* > 1 andare not required. The integral equation results generally show little dependence on theapproximation used. At narrow separations, the pressure is oscillatory with a period ofabout one hard sphere diameter. The oscillations arise from particle correlations betweenthe solvent layers at the surfaces.Chapter 4. Fluids between planar surfaces 107100.55 0.00 —0.5—1.01 2 3 41*5Figure 4.1: Net pressure between hard walls in a hard sphere fluid at a density of= 0.681. (a) The solid curve is the HNC result, the almost coincident dotted curveis 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 scalewith 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 Carlosimulation results of Karlström [123] on an expanded scale (note that the bulk pressure of3.695kBTd has been subtracted from the internal pressure data of Karlström in orderthat the net pressure approaches zero as z —* oc). Although noise in the simulationdata makes it difficult to assess the effects of including bridge functions, the RHNCresults generally show slightly better agreement with simulations. This system of hardspheres between hard walls has also been solved using the inhomogeneous OZ equationtogether with the Percus-Yevick (PY) approximation [113]. The homogeneous RHNCresults are about as accurate as the inhomogeneous results. The largest disagreementbetween 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 gws(z2;1) calculated from Eq. (4.2.13) areIII II 1111111 I I IIP(l)d/kBT (a)I I I I I I I I I2 3 4 5 6 71*Chapter 4. Fluids between planar surfaces 108shown in Fig. 4.2 for several wall-wall separations. The results for 1* = 10 are compared with grand canonical Monte Carlo results [123]. Again the agreement with simulation results is good except very close to the wall where the HNC theory overestimatesgw8(ds/2; 1). At 1* = 3.1 and 1* = 10, the distribution function can also be comparedwith the inhomogeneous PY results [113] (not shown). In both cases, the HNC resultsare about as accurate as the inhomogeneous PY theory except near contact. In an exacttheory, the contact theorem [see Eq. (3.4.23)] tells us that pg(d/2; 1) is equal to the internal 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 andcompare the results with the net pressure results from Eq. (4.2.9). These results showthat while the contact density in the HNC theory is too high, variations in the contactdensity due to confinement of the fluid are proportional to the net pressure and the contact density for hard spheres near a hard wall can also be used to obtain the pressurebetween the surfaces.The net pressure between hard walls mediated by a dipolar hard sphere fluid was alsocalculated using Eq. (4.2.9) for the RHNC theory. In this case, the RHNC calculationwas carried out using the Verlet-Weis parameterization [27] and HNCP approximation[26] for the solvent-solvent and wall-solvent bridge functions, respectively. Figure 4.4shows that the results are qualitatively similar to the hard sphere fluid (Fig. 4.1) andthat 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 ofN(12) between planar walls immersed in a dipolar fluid can be determined. For adipolar fluid near an inert wall, h5(12) decays asymptotically as [see Eq. (3.4.17)]h5(12) —* 32p[53_2]2(x + X 5P2(cosO)) z —+ (4.4.1)where = ‘r is the isothermal compressibility. In the HNC/RHNC approximation,Figure 4.2: The wall-hard sphere distribution function gws(z; 1) for several wall-wall separations 1. The parameters are as in Fig. 4.1.the wall-solvent direct correlation function has the limitc(21) h(21) —u8(21)/kBT, (4.4.2)and, hence, decays as 1/z6 far from an inert wall. In addition, since c(21) —* asz —* cc, we make use of the approximationc(21) Z 1 (4.4.3)kBT F(_l)2 2 I —dz—L 3YEs j .1110964200Chapter 4. Fluids between planar surfaces642010864201234 z*5 0 1 2 3 4z*5108641086420:‘ I : :‘ I1.5Ljd: : :0 lz2 0 lz Ozl2010864200 1 2z*3to obtainF(l)Chapter 4. Fluids between planar surfaces 1100.5 0.50 0—0.5 —0.52 3 4 5 6z’Figure 4.3: The (a) HNC and (b) RHNC contact density minus the bulk density as afunction of wall-wall separation 1. The crosses denotes values of pg(d/2; 1)— Ps• Thecurves are the net HNC and RHNC pressures calculated from Eq. (4.2.9).32l2 [ i)2]2 , 1 - oo, (4.4.4)where y = 47rpS/9kBT. The asymptotic limit of the net pressure, plotted in Fig. 4.4b,is given byP(l) = dF(l) (4.4.5)The first two projections of the wall-solvent distribution functions gws(z2; 1) are shownin Figs 4.5 and 4.6 for several wall-wall separations. Simulation results are available fordipolar hard spheres near isolated surfaces [54] and show good agreement with RHNCtheories.The distribution functions for a fluid of Lennard-Jones particles confined betweentwo Lennard-Jones 10—4 walls at separations of 1* = l/o = 2 and l = 3 are shown inFig. 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 reducedIII 1111111111 1(11111-1 (b)Z_i I -- I— Iz>j7-I I i i i I i i i I I I I i2 3 4 5 6z*Chapter 4. Fluids between planar surfaces 1110.50.40.20.00.0—0.2—0.5—0.4—0.6—1.08 9 10 11 12 lFigure 4.4: HNCP pressure between hard walls in dipolar hard sphere fluids at a densityof p’ = 0.7. In (a), the solid, dashed and dotted curves are for reduced dipoles moments of0 (hard spheres), = and u = respectively. In (b), the HNC asymptotescalculated from Eq. (4.4.5) (smooth curves) are compared with numerical HNCP resultson an expanded scale.temperature of T* = kBT/e = 1.21. The molecular dynamics simulations3were carriedout using a spherically truncated and shifted Lennard-Jones potential(u’(r) — nLJ(3.5a), r <3.5au(r)=(4.4.6)0, r>3.5u.In the HNC theory for Lennard-Jones particles, the bulk direct correlation functionc(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 correlationfunction was calculated using both the full Lennard-Jones, Eq. (2.4.10), and the truncated Lennard-Jones, Eq. (4.4.6), potentials, respectively. Using Eq. (4.4.6), the HNCresults 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,3The state parameters in the molecular dynamics simulation [108] were chosen such that the confinedfluid is in equilibrium with the bulk having an average density 0.5925.111111111111111111111111- P(l)d/kBT x io (b) --- 1 - — .—. —It —-t St. -LiI\J- —‘IIIIIIIIIIIIIIIIIIIIIIIII2 345 61*Chapter 4. Fluids between planar surfaces 1126 64 42 2o 006 64 4 42 2 2o 0 00 1 2z*3 0 lz*2 0 lz*2Figure 4.5: RHNC results for the angle-averaged wall-particle pair distribution functiong°(z; 1) for a dipolar hard sphere fluid confined between hard walls. Results are shownfor 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 between hard walls calculated from the HNC theory, RHNC theory with the HNCP hardsphere bridge function and grand canonical Monte Carlo simulations. The grand canonical Monte Carlo simulation method is described in Sec. 4.3 including a correction forlong-range Lennard-Jones interactions. Results for a wall-wall separation of 1* = 5 areconsidered for three dense bulk state points at reduced temperatures of T* = 1.0, 1.35and 2.0. For the low temperature state point, the distribution function for a wall-wall0 4z*561 2 3 4z*5Chapter 4. Fluids between planar surfaces 1130 1 2 3 4z*5—I I I I I I I I I—÷1111111÷o i2Figure 4.6: RHNC results for the projection g2 of the wall-particle pair distributionfunction for a dipolar hard sphere fluid confined between hard walls. Results are shownfor 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 inTable 4.1.In general, the HNC and RHNC results show relatively little dependence on thestate parameters, whereas the grand canonical Monte Carlo results clearly show that theLennard-Jones fluid is in the liquid phase for the first two state points with 1* = 5 (Figs4.8a and 4.8b) and in the gas phase for the third state with 1* 5 (Fig. 4.8c). In addition,results for the third state point (Figs 4.8c and 4.8d) indicates that the phase transitionis a function, not only of the state parameters, but also of the wall-wall separation. The0—1—2—31*= 50—1—2—30—1—2—3:1 I I I I I I I I I I I I I I I I I I I 1 II-I0 1 2 3 4z*5II lilt Illf—I I i I i I—l-f-I—l-f-I0—1—2—3I I I I I IIl*=2.5 —i-II I I 11111 I0—1—2—30 1 2z*3 0 lz2lilt III IChapter 4. Fluids between planar surfaces 1146 64 42 2o 00.0 0.5 1.0 1.5z*2.0 0Figure 4.7: Wall-solvent distribution functions for a Lennard-Jones fluid confined betweenLennard-Jones 10—4 walls, Eq. (4.3.3) (T* 1.21, p = 0.5925 and Ps2 = 0.8). Thewall-wall separations are (a) 1* = 2 and (b) 1* = 3u. The solid and dotted curves arethe 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 bulkdirect correlation function of perturbation theory [124,125]. The circles are results frommolecular dynamics simulations [108].grand canonical Monte Carlo results are characteristic of a drying transition as a fluidwith attractive interactions near an isolated “solvophobic” surface approaches the bulkcoexistence curves [102,119]. The phase diagram for the Lennard-Jones fluid [126] showsthat the density and temperature considered in Fig. 4.8c are for a state point close tothe coexistence curve.The failure of the HNC theory to account for dewetting at an isolated solvophobicsurface is well-known [119,8] and is a result of the approximation used in the closure. Itis not surprising that the RHNC theory also fails in this case since the hard sphere bridgefunctions are used as an approximation of the state dependent hard wall—Lennard-Jonesbridge function. At high temperatures, the strength of the Lennard-Jones attraction is1 2 z*3Chapter 4. Fluids between planar surfaces 115432103210o 1 2 3 4z*5o 1 2 3 4z*5321032100 1 2 3 4z*50 1 2 3 4z*5Figure 4.8: Wall-solvent distribution functions for a Lennard-Jones fluid confined betweenhard walls for three state points and two wall-wall separations. The solid curves areresults from grand canonical Monte Carlo simulations. The dashed and dotted curvesare the HNC and HNCP results. In (a) T* = 2.0 and p = 0.710 + 0.003; (b) T* = 1.35and p = 0.70 ± 0.01; and in (c) and (d) T* = 1.0 and p = 0.75 + 0.01._i11Fj11IIJIIIJIIIijiiii_g(z;5u) (a)1I I I I I I II J I II J Illij liiig(z;5u) (b)I I I I Iliii 1111111111111111111- g(z;5a) (c)I— I I —— I I —— I I ——I I —— I I —— I —— f\ /\ f\/-‘ / \/ ‘‘ I —/ \ I --IIIIIIIIIIIIIIIIII11_III I I I I I J I I I I I I I I I- gws(z; lOu) (d)I I I I I I I I I I IChapter 4. Fluids between planar surfaces 116Table 4.1: Details of the grand canonical Monte Carlo simulations for a Lennard-Jonesfluid confined between hard walls. The last column labeled “operations”, denotes thetotal number of particle moves, insertions and deletions included in the simulation. (Thegrand canonical Monte Carlo simulations for the bulk were performed in a central cell9u x 9u x 9a.)Operationsi/u L/u St/J T* preS/kT x1065 15 3.5 1.0 0.751 + 0.007 -3.287 1010 7 3.5 1.0 0.751 + 0.007 -3.287 105 7 3.5 1.35 0.700+0.006 -1.941 105 7 3.5 2.0 0.710+0.003 1.525 10attenuated and the bridge functions, defined as integrals over Mayer-f functions,f(12) = e_ 12)/kT — 1, (4.4.7)approach the hard wall results and better results are obtained. For narrow wall-wallseparations, even at relatively high temperatures, the HNC and RHNC approximationare less reliable. For the Lennard-Jones fluid near an attractive 10—4 wall, it is worthwhileto note that the state point considered in Fig. 4.7 (T* = 1.21 and p = —.5925) is alsovery close to the coexistence curve yet the integral equation theories show relatively goodagreement with simulations.4.5 Cavitation of a Lennard-Jones fluid between hard wallsAs pointed out in the last section, the HNC and RHNC theories, using hard wall-hardsphere bridge functions, are inadequate to describe Lennard-Jones fluids between hardwalls for state points near the coexistence curve. In this section, grand canonical MonteCarlo is used to investigate the phase behaviour of this system for two bulk state pointsnear the coexistence curve. In addition to studying the equilibrium states of the fluid,Chapter 4. Fluids between planar surfaces 117we also consider the metastable liquid and gaseous branches. Tables 4.2 and 4.3 containdetails of the simulations. The threshold density th, defined in Sec. 4.3, is included for themetastable states. Table 4.2 is for the confined Lennard-Jones fluid with T* = =1.0 and ,t1e8/kBT = —3.29, which corresponds to a bulk having excess chemical potential,eX/kBT = —3.00; density, Ps3 = 0.75 ± 0.01; and pressure, Pbu3/kBT = 0.44 + 0.05.Table 4.3 is for the second state point, T* = 0.9 and u/kBT = —3.84, which correspondsto ex/kBT = —3.60, p8a3 = 0.79 ± 0.01, PbJ3/kBT = 0.4 ± 0.1. Again, the grandcanonical Monte Carlo simulations for the bulk liquids were performed in a central cubiccell using periodic boundary conditions. The Lennard-Jones pair potential was truncatedspherically at half the cell length and the potential energy was corrected assuming auniform distribution of particles beyond the truncation sphere.Figure 4.9a shows the average density between the walls as a function of wall separation for the first state point. The equilibrium phase points, denoted by crosses, show thata high density fluid is stable at large separations, and that a distinct low density fluid isstable at small separations. The two fluids coexist at a separation of ‘-. 6.Ou. It was possible to follow the metastable liquid branch, denoted by open squares, into separations assmall as 4.Ou. Hence this separation approximates the spinodal point at which cavitationoccurs. Although decreasing, the average density of the metastable liquid remains clearlygreater than the density of the gas. Moreover, this density is larger than the thresholddensity (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 toseparations of 15o- and 20u. Although there is a significant L dependence in the largerseparation result, it is clear that the spinodal point for condensation is at very largeseparationsFigure 4.9b shows the average density for the second, higher density, lower temperature, state point. Here, the gas phase is also stable at narrow separations, 1 5 6.5a.Chapter 4. Fluids between planar surfaces00.80.60.40.20-x XXX5L:J10 15118C20 i*Figure 4.9: The grand canonical Monte Carlo average density = (N)/1L2, betweenhard walls for a Lennard-Jones fluid in equilibrium with the bulk. Results for the firststate point (&/kBT = —3.29 and T* = 1.0) are shown in (a) and for the secondstate point (‘/kBT = —3.84 and T* = 0.9) in (b). Crosses represent points that arethermodynamically stable and open squares represent those that are metastable. Opencircles are used when the equilibrium and metastable fluids could not be distinguished.Except where illustrated, the standard deviations in are smaller than the symbols. Theinset shows an expanded view of the gas branch.I I I I I0.80.60.40.2XXI I I I I I I I I I I I I(a)XXI I I I I I I I I I I5 10 15 20 *I I I I I I I I I I I I I I I I I: (b)—0 . 11111 II I I I0 IC —E 0.04- -- 0.02— -— 111111111110 20 40 60XX 0I I I I I I I I I I IChapter 4. Fluids between planar surfaces 119Table 4.2: Details of the grand canonical Monte Carlo simulations for a Lennard-Jonesfluid 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 cell9u x 9o x 9, has (pex)/kBT = —3.00, PsU3 = 0.75 ± 0.01 and Fbo3/kBT 0.44 ± 0.05.)Operationsi/u L/u St/U iiu3 PIflt(l)cx3/kBT x1063.0 7 3.5 0.0493 + 0.0002 0.048 + 0.002 105.0 10 5.0 0.0621 + 0.0005 0.051 ± 0.002 405.0 10 5.0 0.47 ±0.01 0.29 ±0.01 0.40 405.2 10 5.0 0.50 +0.01 0.31 ±0.03 0.40 205.5 10 5.0 0.065 +0.001 0.052+0.002 305.5 10 5.0 0.551 ±0.007 0.36 ±0.02 0.45 205.8 10 5.0 0.572 ±0.006 0.38 +0.02 0.48 206.0 10 5.0 0.068 ± 0.001 0.052 + 0.002 306.0 10 5.0 0.588 ±0.006 0.38 ±0.02 306.5 10 5.0 0.607 +0.006 0.40 ±0.03 307.0 10 5.0 0.624 + 0.007 0.405 + 0.02 308.0 10 5.0 0.643 ±0.007 0.41 ±0.03 3010.0 10 5.0 0.668 ±0.006 0.42 ±0.03 6010.0 15 5.0 0.084 +0.002 0.052±0.002 0.10 3015.0 15 5.0 0.095 ±0.001 0.053+0.03 0.10 3020.0 15 5.0 0.0988±0.0002 0.041±0.03 0.10 30At large separations and in the bulk, the liquid phase is stable and the gas phase ismetastable. For separations 6.5u 5 1 51u, the grand canonical Monte Carlo methodcould not distinguish the equilibrium and metastable states after at least 20 million configurations, and thus the coexistence separation was not determined. Compared to thestate point shown in Fig 4.9a, the spinodal separation for cavitation is slightly smaller at4.9u and although the coexistence separation is not known, it is at least slightly largerat 1 6.5cr.The density profiles for the equilibrium and metastable liquid branches are shown inFig. 4.10 using the midplane between the two walls as the origin. The sharp decrease inChapter 4. Fluids between planar surfaces 120Table 4.3: Details of the grand canonical Monte Carlo simulations for the second statepoint (T* = 0.9 and res/kBT = —3.84). (The grand canonical Monte Carlo simulationfor the bulk, performed in a central cell 7o’ x 7o x 7o-, has ex = 3.60, Ps°’3 = 0.79±0.01and Pb0’3/kBT = 0.4 ± 0.1.)Operationsl/ L/ St/0’ Pt(l)cr3/kBT Pth0’3 x1064.9 10 5.0 0.496 ±0.024 0.25 ±0.02 0.40 305.0 10 5.0 0.0287 ± 0.001 0.023 ± 0.002 305.0 10 5.0 0.562 ±0.01 0.31 ±0.02 0.40 205.2 10 5.0 0.591 ±0.01 0.34 ±0.02 0.40 205.5 10 5.0 0.619 +0.01 0.36 ±0.02 0.42 205.8 10 5.0 0.634 +0.01 0.37 ±0.02 0.42 206.0 10 5.0 0.0294 ± 0.0001 0.023 ± 0.002 306.0 10 5.0 0.644 ±0.007 0.38 ±0.02 0.45 206.5 10 5.0 0.662 ±0.006 0.39 ±0.02 207.0 10 5.0 0.674 ±0.005 0.39 ±0.02 208.0 10 5.0 0.695 ± 0.007 0.41 + 0.02 2010.0 10 5.0 0.718 +0.006 0.41 ±0.02 2021.0 16 8.0 0.0320 ± 0.0001 0.023 ± 0.002 3051.0 10 5.0 0.0323±0.0002 0.023+0.002 30Chapter 4. Fluids between planar surfaces—4 —2 0 2 4z* —4 —2 0 2 4z*Figure 4.10: Grand canonical Monte Carlo density distributions for the equilibrium andmetastable 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 thewalls is at z = 0.the density adjacent to the walls is characteristic of a liquid close to drying. The densityprofiles show the oscillations that are characteristic of a liquid, and their amplitude isslightly 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 netpressure must approach zero at large separations. In the simulations, the net pressureat large separations is equal to the difference in two nearly equal quantities (the totalpressure 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 simulationcell, and also from the different potential truncations used (spherical in the bulk andcylindrical 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 onthe liquid branch rather that the bulk pressure obtained by direct simulation. This onlyshifts the curve by a few percent.In general, Fig. 4.11 shows that the force between the walls is attractive, and on10.80.60.40.20I I I I I I I I I I I I I Ipg5(z;l) (a)l*=567810.80.60.40.201211111 I 111111111 IIpg(z;l) (b)— LO—78liii I I I liiiChapter 4. Fluids between planar surfaces 1220 0—0.2 —0.2—0.4 —0.42 4 6 8 101*Figure 4.11: The net pressure between hard walls. The symbols are the same as inFig. 4.9. The solid curve is the conventional van der Waals attraction, Eq. (4.5.3), andthe dashed curve is the fitted mean-field result, Eq. (4.5.2). The fitting parameters forthe first state point (a) are A0/kBT = 0.15, cw2 = 0.02, and l = 3.Ou. For the secondstate point (b), AJ/kBT = 0.13, a2 = 0.02, and 10 = 3.0g. For the gas branch, thestandard deviations in the net pressure are smaller than the symbols.any one branch the magnitude decays monotonically with increasing separation. It ispossible that the pressure on the liquid branch is oscillatory due to molecular packingbetween 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 isa discontinuous jump in the pressure as the liquid evaporates to the gas, in qualitativeagreement with earlier results [104,116]. At small separations, the force is dominated bya large attraction. Along the gas branch, the total pressure in the slit is very small, andhence, the net pressure is essentially the negative of the bulk pressure. This represents alarge 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 toAttard [127]. In this theory, the fluid at the midplane between the surfaces is assigned aneffective chemical potential l(l) equal to the chemical potential of the bulk liquid minus2 4 6 8 101*Chapter 4. Fluids between planar surfaces 123the excess mean-field contribution from fluid displaced by the presence of the walls. (Thefluid 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 chemicalpotential. At a smaller separation, it will equal the bulk spinodal value which defines thespinodal separation 10. Close to the spinodal separation, the potential of mean force perunit area, ‘(l), and pressure, P(l), between walls decay asw(l) “ i1e__b0)1, 1 .‘ 10, (4.5.1)P(l) 1‘ lo, (4.5.2)where A and a are positive parameters of the theory. Strictly speaking, the pressurebetween 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 simpletheory gives a rather good representation of the data for the liquid branch. In bothcases the values of lo are somewhat smaller than those estimated directly from the grandcanonical Monte Carlo simulations. This is probably due to the simplicity of the mean-field theory which completely neglects the structure of the fluid adjacent to the walls.Asymptotically, the net pressure between hard walls mediated by a Lennard-Jonesfluid, is given by [121]P(l) = 2irpu6 (4.5.3)This is the pressure that arises from the van der Waals forces and is included in Fig. 4.11as solid curves. It may be seen that the van der Waals contribution greatly underestimatesthe attractions for the metastable liquid branch. This is unambiguous at separations nearthe spinodal point and is expected to hold at larger separations.Figure 4.12 shows experimental measurements of the potential of mean force perunit area between two hydrophobic surfaces in water [128]. The measurements wereChapter 4. Fluids between planar surfaces 124obtained at standard temperature and pressure. In Fig. 4.12, the solid curve representsthe 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 102J. This value is an upper limit for a hydrocarbon-water-hydrocarbon or a mica-water-mica system. We see that the measured attraction issignificantly larger than the van der Waals force and increases sharply as the surfaces arebrought into contact. The dashed line represents the mean field result, with parametersA = 1.2 x 10’2J/m,‘0 = l2nm, and a = 4 x 104nm2. The measured force is verywell fitted by this choice of parameters. At a narrow separation (1 10 ‘-. l2nm), theattraction between surfaces exceeds the counter force of experimental apparatus and thesurfaces jump into contact. When the surfaces are subsequently separated, cavitation canpersist in some cases for separations on the order of micrometers [99], suggesting thatwater is metastable between these surfaces at relatively large separations. Currently, theorigin of this long-ranged attraction between hydrophobic surfaces is unknown. However,comparison of the experiments and simulations with the mean-field theory, suggest thatthe measured attraction might result from a separation-induced phase transition.The simulation results show that cavitation will occur for water at subcritical statepoints near the coexistence curve and between hydrophobic surfaces at sufficiently narrowseparation. In the simulation results of a Lennard-Jones fluid between hard walls at aparticular bulk state point showed a cavitation separation of about 4 diameters, whereaswater at standard temperature and pressure between hydrophobic walls apparently cayitates at separations of about 40 water diameters.Chapter 4. Fluids between planar surfaces 125I 111111111 liiw(l) (mJ/m20---ctJ0.05: d--—0.1—0.15 --I I I I I I I I I I I I I I I I I20 40 60 80 100 120 1 (nm)Figure 4.12: The interaction free energy per unit area measured between hydrophobicsurfaces in water [128]. The solid curve is the conventional van der Waals attraction,with a Hamaker constant of 2 x 102J, and the dashed curve is the fitted mean-fieldresult, Eq. (4.5.1) (A = 1.2 x 10’Jm,to = l2nm, and = 4 x 104nm2).Chapter 5Summary and ConclusionThis thesis has been concerned with nonspherical particles near planar surfaces andwith both spherical and nonspherical particles between planar parallel surfaces. Generaltheoretical methods for obtaining the wall-particle pair distribution function for multicomponent fluids near isolated inert, dielectric and metallic surfaces were presented.These were used to study simple dipolar hard sphere fluids and simple model electrolytesolutions near such surfaces. General methods for finding the fluid-mediated net pressure between planar walls and for calculating the wall-particle correlation functions fornonspherical particles between surfaces have also been derived. Fluids of hard spheres,dipolar hard spheres and Lennard-Jones particles were considered between inert andLennard-Jones surfaces. The theories were used to investigate the fluid-mediated netpressure between surfaces and the effects of confinement on the fluid structure. Comparison of HNC and RHNC results with grand canonical Monte Carlo calculations andmolecular dynamics simulations allowed us to establish the validity of these integralequation theories for certain systems. Finally, a grand canonical method for studyingmetastable 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 measuredattraction between hydrophobic surfaces.General reductions of the OZ equation and the HNC/RHNC closure relation weregiven for multicomponent fluids of nonspherical particles near an isolated wall. Theequations are presented in a formalism similar to that previously derived for bulk fluids126Chapter 5. Summary and Conclusion 127[6,10—13,17,18], and provide a relatively efficient method of calculating the structureof molecular liquids near planar walls. In addition, fully self-consistent theories arederived for multipolar particles near dielectric and metallic surfaces. In these theories, themany-body interfacial potentials are reduced to an effective mean-field pair interactionappropriate for use in integral equation theories.For a pure dipolar solvent near an inert surface, analytical expressions for the asymptotic behaviour of the wall-solvent pair correlation function at large wall-solvent separations were obtained. Within the RNC approximation, the resulting equations werecompared with the classical continuum (image-dipole) expression and with an exact limiting result obtained by Badiali [62]. It is shown that all three limiting expressions agreeonly to linear order in density. The HNC and continuum theories give different secondorder terms and both are inexact. Numerical solutions of the RHNC equations wereobtained for a dense system and these indicate that solvent structural effects cease to beimportant and the z3 asymptotic limiting behaviour is reached at about 11-13 solventdiameters. Comparison with the exact asymptotic result shows that the RHNC approximation overestimates the magnitude of the coefficient multiplying z3 in the limitingexpression.At a semi-infinite dielectric slab, the interaction energy of a fluid particle with thepolarizable material is approximated by a uniform mean-field pair potential. Formalexpressions have been obtained for general multipolar models. For the metal-liquid interface, a metal slab is modeled at the jellium level and the electron density distributionacross the interface is obtained from density functional theory. The liquid and metal interact electrostatically and the resulting electron density distribution and metal-particledistribution functions obtained are completely self-consistent. This self-consistency is aprinciple feature of the present theory for the liquid-metal interface and is a significantChapter 5. Summary and Conclusion 128improvement upon earlier treatments based upon the mean spherical closure approximation 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 andmetallic surfaces. Symmetry requirements do not allow the simple dipolar solvent to bepolarized by an inert surface. This remains true if the dissolved anions and cations areof the same size. If the ions are of unequal size a net solvent polarization arises fromthe asymmetric density distributions of anions and cations across the interface, but thiseffect is not very large for simple ions. For inert walls, the ion distributions are governedprincipally 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 whichare large have a high contact density. This holds independent of the sign of the ion chargeand 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 thepure dipolar case, it was shown that images of other solvent particles make a large contribution to the interaction between a particular particle and the wall. If these contributionsare ignored and one includes only the simple self-image potential then rather dramaticstructural changes occur as the dielectric constant of the wall is increased. However, theimage contribution from all remaining fluid particles acts to largely cancel the self-imageterm and, as a result, the qualitative solvent structure near the surface is not stronglyinfluenced 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 inertwalls.At a metal surface the solvent and ion distributions differ fundamentally from theinert and dielectric wall cases. The metal induces a very strong short-range electric fieldand a net polarization of the solvent distribution results. This was most dramatic in theChapter 5. Summary and Conclusion 129layer 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 fordipolar monolayers [67,68] and contrasts with the inert wall case where the dipoles preferto lie perpendicular to the surface normal [54,50].The polarization of the solvent at the surface plays an important role in determining the ion distribution. This gives rise to a “layered” structure with anions located atd3 + d_ /2 (measured from the surface) and cations nearer the surface. These qualitative structural features were found for all solutions considered and were not particularlysensitive to ion size or salt concentration. It is worth noting that the direct ion-jelliuminteraction is attractive for anions and repulsive for cations and hence would imply an “inverse” order with more anions than cations near the surface. Thus the observed ion-walldistributions are to a large extent determined by the ion-solvent forces and are greatlyinfluenced by the jellium-induced solvent order. These calculations emphasize that thesolvent response to the metal potential is of crucial importance and must be included indouble layer theories.The electron density distribution of the metal is also altered by the presence of thefluid and extends a small distance further from the jellium edge in the presence of theliquid. The electrostatic potential drop across the interface can be written as the sumof jellium, ion and dipole terms. In all cases, it was found that the jellium contributionincreases due to the presence of the fluid while, as expected, the net electrostatic potentialdrop across the interface is reduced. Compared with the vacuum values the potentialreductions for the pure solvent range from 6—10% (i.e., 0.42—0.24 volts) as r3/a0 is variedfrom 2.0—3.0. It is interesting to observe that these values do not differ greatly from thosegiven by the monolayer model [67,68]. In the presence of ions, the jellium contributionto the potential drop is not greatly altered from the pure solvent value. On the otherhand, the ion and dipolar contributions are relatively large and are sensitive to bothChapter 5. Summary and Conclusion 130salt concentration and ion size. However, these terms cancel to give a net solutioncontribution comparable in magnitude with that of the pure solvent. This demonstratesthat the ion and dipole contributions are very strongly coupled and should never betreated as independent quantities in theories of the double layer.Fluids between planar surfaces were also considered using the HNC and RHNC theories and grand canonical Monte Carlo methods. A general reduction of the wall-wallOZ equation was presented for nonspherical particles between walls. This allows thewall-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 reducedfor fluids of nonspherical particles. Previously this wall-solvent OZ equation had onlybeen considered for spherical particles [109]. These equations are used to investigate theeffects of confinement on the net pressure between surfaces, fluid structure and phasebehaviour.For simple hard sphere fluids between hard walls, the net pressure between surfacesas a function of the wall-wall separation was obtained from the HNC and RHNC theories. 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 grandcanonical Monte Carlo simulations [123] and with results obtained from the inhomogeneous PY theory [113]. The small discrepancies between the RHNC and grand canonicalMonte Carlo results were comparable to those between the inhomogeneous PY and grandcanonical Monte Carlo results. HNC results for the hard sphere distribution function between hard walls were also obtained for a few wall-wall separations in the range of 1.2—10solvent diameters. The distribution function at a separation of 10 solvent diameters wasin 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 resultsfrom the inhomogeneous PY theory. While the HNC contact densitiesp5g(d/2; 1) areChapter 5. Summary and Conclusion 131incorrect, the contact theorem was found to still give good results for the net pressurebetween hard walls.Dipolar hard sphere fluids between hard walls were considered using the RHNC approximation with HNCP bridge functions. The net pressure was calculated for two valuesof the reduced solvent dipole moment. It was found that the solvent-mediated pressureis similar to the hard sphere-mediated pressure and the net pressure is smaller for largerdipole moments. The wall-solvent structure was also examined for several wall-wall separations.The distribution functions for fluids of Lennard-Jones particles confined between wallswere also calculated from the HNC and RHNC theories. The bulk direct correlationfunction needed in the wall-solvent integral equation theories was obtained from perturbation theory [124,125]. These results were compared with computer simulations. TheLennard-Jones parameters for the wall-particle pair potential were equal to the valuesused for the particle-particle pair interaction. For a Lennard-Jones fluid confined betweenLennard-Jones 10—4 plates, the HNC distribution functions showed very good agreementwith simulations. It is worthwhile noting that these results were obtained for a bulkLennard-Jones state point near the coexistence curve.When the Lennard-Jones fluid was confined between hard walls, grand canonicalMonte Carlo results indicate that the wall-solvent distribution function is highly sensitiveto the bulk state parameters and wall-wall separation. In particular, it was found that thefluid, which is liquid in the bulk and at large wall-wall separations, evaporates at narrowseparations. This separation-induced phase transition arises from the loss of attractivematerial (displaced by the inert walls) and corresponds to a decrease in the interactionenergy per particle relative to the bulk. The HNC and RHNC results showed relativelygood 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 failChapter 5. Summary and Conclusion 132to describe solvent dewetting at a hard solvophobic wall.Grand canonical Monte Carlo was used to study confined Lennard-Jones systems indetail. For one subcritical state point, (T* = 1.0, p 0.75), the coexistence separationfrom grand canonical Monte Carlo was at l = 6. The metastable branches on either sideof coexistence were also explored. Spinodal cavitation was found at a separation of 1* = 5and spinodal condensation occurred near = 15. The net pressure between walls was ofparticular interest. As the wall-wall separation was decreased, the pressure was found tobe attractive and increase monotonically on the liquid branch as the spinodal cavitationseparation was approached. As previously reported [115], a discontinuous increase in thepressure 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 betweensolvophobic surfaces, separation-induced cavitation is certain to occur in the limit thatthe state parameters approach the bulk coexistence curve. The separations at whichcoexistence and spinodal points occur depend upon the specific system and the distanceof the bulk state point from coexistence. Although experimental results for aqueoussystems at standard temperature and pressure indicate a much longer ranged attraction(measurable to lOOnm), both experiments and simulations show qualitatively similarbehaviour. Our simulation results suggest that a separation-induced phase transitionmight be the origin of the experimentally measured attraction between hydrophobicsurfaces immersed in water.Bibliography[1] J. D. Porter and A. S. Zinn. Ordering of liquid water at metal surfaces in tunneljunction 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, London, 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: Thermodynamic 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 certainintegral 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 withimage 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 approximationfor fluids of nonspherical particles, a general method with applications to dipolarhard spheres. J. Chem. Phys., 82:429—440, 1985.[11] P. G. Kusalik and 0. N. Patey. The solution of the reference hypernetted-chainapproximation for water-like models. Mol. Phys., 65:1105—1119, 1988.[12] P. G. Kusalik and 0. N. Patey. On the molecular theory of aqueous electrolytesolutions. I. The solution of the RHNC approximation for models at finite concentration. J. Chem. Phys., 88:7715—7738, 1988.133Bibliography 134[13] P. G. Kusalik and 0. N. Patey. On the molecular theory of aqueous electrolyte solutions. II. Structural and thermodynamic properties of different models at infinitedilution. J. Chem. Phys., 89:5843—5851, 1988.[14] P. 0. Kusalik and G. N. Patey. On the molecular theory of aqueous electrolytesolutions. IV. Effects of solvent polarizability. J. Chem. Phys., 92:1345—1358, 1990.[15] P. G. Kusalik. The structural, thermodynamic and dielectric properties of electrolyte 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, Universityof British Columbia, Vancouver, Canada, 1990.[17] L. Blum. Invariant expansion. II. The Ornstein-Zernike equation for nonsphericalmolecules 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 modelfor 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. NumericalRecipes: 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 calculation of the pair correlation function. Physica, 25:792—808, 1959.[21] E. Meeron. Nodal expansions. III. Exact integral equations for the particle correlation 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 classical 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 highdensity and overlap regions. Mol. Phys., 67:431—438, 1989.[31] D. Henderson and M. Plischke. Pair correlation function in a fluid with densityinhomogeneities: 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 OrnsteinZernike equation for hard spheres near a hard wall. a rapid method of numericalsolution 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 molecular fluids. Chem. Phys. Lett., 121:347—350, 1985.[34] J. M. Caillol. A method for solving the hypernetted-chain approximation for molecular fluids. Chem. Phys. Lett., 156:357—362, 1989.[35] M.P. Allen and D. J. Tildesley. Computer Simulations of Liquids. Oxford University, 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 solutions. III. A comparison between Born-Oppenheimer and McMillan-Mayer levelsof description. J. Chem. Phys., 89:7478—7484, 1988.[39] D. J. Isbister and B. C. Freasier. Adsorption of dipolar hard spheres onto a smoothhard 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 fordipoles 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 chargedsurfaces: Ion-dipole mixtures. J. Chem. Phys., 73:2949—2957, 1980.[42] L. Blum and D. Henderson. Mixtures of hard ions and dipoles against a chargedwall: The Ornstein-Zernike equation, some exact results, and the mean sphericalapproximation. 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 thesolid-electrolyte solution interface. I. Structure of a hard sphere ion-dipole mixturenear 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: Adiscussion 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 anelectrical double layer: Reference hypernetted-chain (RHNC) results for solventstructure 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 anelectrical double layer: Reference hypernetted-chain results for ion behaviour atinfinite dilution. J. Chem. Phys., 89:3285—3294, 1988.[49] 0. M. Torrie, P. G. Kusalik, and G. N. Patey. Molecular solvent model for anelectrical double layer: Reference hypernetted-chain results for potassium chloridesolutions. J. Chem. Phys., 90:4513—4527, 1989.[50] 0. M. Torrie, A. Perera, and G. N. Patey. Reference hypernetted-chain theory fordipolar hard spheres at charged surfaces. Mol. Phys., 67:1337—1353, 1989.[51] M. Kinoshita and M. Harada. Numerical solution of the RHNC theory for fluidsof non-spherical particles near a uniform planar wall: A study with dipolar hardspheres as an example system. Mol. Phys., 79:145—167, 1993.[52] M. Kinoshita and M. Harada. Numerical solution of the RHNC theory for water-like 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 inhomogeneous dipolar fluid. Mol. Phys., 68:903—915, 1989.Bibliography 137[54] D. Levesque and J.-J. Weis. Orientational structure of dipolar hard spheres near ahard 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 dipolarfluid 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 atan 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 uncharged, unpolarized surface. J. Chem. Phys., 86:4162—4170, 1987.[58] S. B. Zhu and G. W. Robinson. Structure and dynamics of liquid water betweenplates. J. Chem. Phys., 94:1403, 1991.[59] Zixiang Tang, L. E. Scriven, and H. T. Davis. Structure of a dipolar hard spherefluid 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 doublelayer 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 waterbetween 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 behaviorof 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. Roleof 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 propertiesof 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 tothe 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 polarizableinterface 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 forthe liquid-ionic solution interface. J. Electroanal. Chem., 158:253—267, 1983.[70] Wolfgang Schmickler and Douglas Henderson. The interphase between jellium anda 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 crystalface of an FCC metal/electrolyte interface. J. Chem. Phys., 82:2825—2829, 1985.[72] Wolfgang Schmickler and Douglas Henderson. The interphase between jellium anda 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 thedouble layer. J. Electroanal. Chem., 150:347—353, 1983.[74] J. W. Halley, B. Johnson, D. L. Price, and M. Schwalm. Quantum theory of thedouble 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 includingsolvent 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 couplingat 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 ofan 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 anideally 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 surfaceenergy. Phys. Rev. B, 1:4555—4568, 1970.Bibliography 139[81] M. Kinoshita and D. R. Berard. Numerical methods for analyzing the bulk andsurface-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 correlation effects. Phys. Rev., 140A:1133—1138, 1965.[84] J. W. Perram and L. R. White. Structure of the liquid/vapour and liquid/solidinterfaces. Discuss. Faraday Soc., 59:29—37, 1975.[85] D. Henderson, Farid F. Abraham, and John A. Barker. The Ornstein-Zernikeequation 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-Zernikeequation 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. Structurenear 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 profilesat 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 fluidsof 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: Statistical 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 between 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 betweensolvophobic 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 incontact with a (100) Lennard-Jones wall and its relationship to idealized fluid/wallsystems: 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 ofadsorption 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 surfacetransition. 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. structureof 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, liquidfilled 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 hypernettedchain 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 results for the force between plates. J. Chem. Phys., 93:1386—1398, 1990.[112] E. GonzJez-Tovar, M. Lozada-Cassou, and W. Olivares. Electrokinetic flow inultrafine slits: Ionic size effects. J. Chem. Phys., 94:2219—2231, 1991.[113] Roland Kjellander and Sten Sarman. On the statistical mechanics of inhomogeneous fluids in narrow slits. An application to a hard sphere fluid between hardwalls. Chem. Phys. Lett., 149:102—108, 1988.[114] Sten Sarman. The influence of the fluid-wall interaction potential on the structureof 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 fluidinterfaces 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 forfluids 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: Layeringtransitions from a continuum theory. J. Chem. Phys., 89:4412—4423, 1988.[118] Ziming Tan and Keith E. Gubbins. Selective adsorption of simple mixtures in slitpores: 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 freeenergy between planar walls in dense fluids: An Ornstein-Zernike approach withresults 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 simulations 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 theequilibrium structure of simple liquids. J. Chern. Phys., 54:5237—5247, 1971.[125] Hans C. Andersen, John D. Weeks, and Davis Chandler. Relationship between thehard-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 fluidbetween hard walls, and the possible relevance to the attraction measured betweenhydrophobic surfaces. J. Chem. Phys., 98:7236—7244, 1993.[128] H. K. Christenson, J. Fang, B. W. Ninham, and J. L. Parker. Effects of divalentelectrolyte 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 AProperties of the Wigner symbols, rotation matrices and rotationalinvariants.In deriving much of the theory presented in this thesis, the orthogonality and symmetry properties of the Wigner 3-j symbols, rotation matrices and rotational invariants areused extensively. For convenience, many of the properties of these functions are listedhere. All of these relationships can be found in Ref. [5J or derived from relationshipstherein.The matrix R associated with rotating a rank-rn spherical tensor has the elementsR(f) (known as the Wigner generalized spherical harmonics or the symmetric-topwavefunctions 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 anasymmetric particle and can be written asR4() = e’d (0)& , (A.l)where the functions d(0) are real. The index m takes on all positive integer and halfinteger values; and t’ are given by= —rn, —rn + 1, ..., rn — 1, m, (A.2a)= —m, —rn + 1, ..., m — 1, m. (A.2b)Since the rotation matrices are unitary (RRt = I, where Rt is the complex transposeof R and I is the identity matrix), we have= [RtA,(x,o,)j* = [R(,o,x)j*, (A.3)143Appendix A. Properties of the Wigner symbols, ... 144where R is the complex conjugate of R, and the elements R must satisfy the sumrules’R(q, 0, x)R1>,,(qf, 9, x) d(0) d,A,I(0) = (A.4)A’ A’= d,A(0)d,,A(0) = A’ . (A.5)A AThe functions d(0) have several symmetry properties which will be of use. These ared(0) = d,(0) (A.6a)(—)‘d(0) (A.6b)= d1(—0) = (—)‘d(—0) (A.6c)d(r — 0) = ()m+LLdm(0) (A.6d)where j = —p. Using Eq. (A.6b), we see thatR() = (—)‘R(f) . (A.7)In addition,dL(O) = (A.8a)d(+r) = Hm±IL’6, (A.8b)A few useful integrals, sums and products of the Wigner spherical harmonics follow:87r2fR() dfz = 6mn6,’a,”5,w , (A.9)2m+12Imn i fmnlf = 87r , (A.1O)1m n 1 (mn 1) Rl*() (A.n)R(f)R1() = (2l -- 1),IA’A‘mn 1” frnnl\RA(fZ) = (2l + 1),I-I’ Lx R() R,(f) , (A.12)‘The sum rules given in Eqs (3.77) and (3.78) of Zare [5] are incorrect.Appendix A. Properties of the Wigner symbols, ... 145where ) is a Wigner 3-j symbol.Applying two rotations in succession and in the order R(f1) then R() is equivalentto a single rotation R(f),= R(fZ)= R11()R,,L(Zl) . (A.13)IL”Using the unitary property, the inverse expression isR(fZ)[ft(1i)] == [R(1)]. (A.14)Ii”These relations are useful when analyzing the pair functions for molecular fluids (seeAppendix B).The 3-j symbols have a few useful symmetry properties. Exchanging two column isequivalent to multiplication by (_)m+n+l. For example,(m fl— m+n+l (m 1 m+n+l (n m l = m+n+l (1 n m (A 15V)Two such exchanges gives(mn l — (n 1 m — (1 A 16V)VIL)ILV)In addition, (:) = m+n+lThe orthogonality properties of the 3-j symbols are(inn l” (mn l’’\—6ll’5)\’ A 18jt—(21+1) ‘ ( )(2l +1) (: ) ( , = (A.19)Appendix A. Properties of the Wigner symbols, ... 146If, in Eq. (A.18), the sum over,\ is included, we getfrnnl\ frnni’\=6. (A.20)A general expression for evaluating the Wigner symbols, as well as FORTRAN code, canbe found in Appendix A of Ref. [5].The orientational space of one asymmetric particle with respect to another is spannedby the rotational invariantsmnl(, 2, p21)mni Rm ()R,(2R,0(21) . (A.21)/W i/The complex conjugate of mnl ismnl*(, 2,r21) = mni (m fll Rrn* ()R(2R0(21) . (A.22)IL’1)’’ /2 V X) L1/tA few useful integrals over the rotational invariants areJ ‘(1, 2,r21) 4minili*(f 2,r21)df1 dfZ2/L1111 \512ir6 (fmnl)2(2n + 1) (2m + 1) (21 + 1)8mm161hh111 (A.23)J Jrnfll(j2) mmfuh1(21\ m2z2l2(12) dfZ1 dfZ2/L1L’l \ ) /.42L’2= 5127r6 fmnlfmlnhllfm2n2l2(m m1 m2 1 2 / i 11 i m m1 m22 rj j1i 1L2) ( 2) 2) { 11 12 } .(A.24)Since the rotational invariants are independent of the orientation of the pair, the integration over the 121 is usually omitted from the orthogonality relation (A.23) and,consequently, the normalization constant is reduced by a factor of 87r2.Appendix BThe electrostatic interaction potential between a molecule and a uniformsurface charge distribution.In this appendix, the metal wall-particle pair potential is reduced to a rotationalinvariant expansion. Consider the static charge distribution n(r) of a molecule of speciesnear a planar polarizable surface having a charge distribution n(z). The electrostaticinteraction energy between the surface and molecule is given byup(O1) = e2JJfldrdt, (B.1)where 0 and 1 denote the coordinates of the wall and molecule 3, respectively, andZt = t . We assume that the charge distribution of the molecule is localized near thecenter of mass and further that the surface and molecular charge distributions do notoverlap. Choosing an origin such that the molecule is at z1 and r = z1 + r’, we can usethe Taylor series expansion1 1 00 m— ri = t’ — r’ = (t)m+l Pm(cos ), (B.2)where t’ = t — z, t’ = 1t9, r’ = Ir’I and is the angle between the vectors t’ and r’.Using the spherical harmonic addition theorem,Fi(cos7) = R0(’)R,i’) , (B.3)“I147Appendix B. The electrostatic pair potential... 148Eq. (B.2) and the identity drdt = dr’dt’, Eq. (B.1) reduces tou(O1) = e2 J(r)mnp(r + ziê)R0(’) dr’ f n(zi + Zt’) R(’) dt’mp’(t)m+l= 2e f dz’ ft’ dt’ n(zi + Zt’)mi’(j!)m+l Pm(CO5O) (B.4)where Eq. (2.4.9) has been used to defines the molecular multipole moments andcos 8 = z’/t’ = (zi —z1)/t’ denotes the polar angle of t’. Writing the multipole momentsin the molecular frame [using Eq. (2.4.7)], we getu,3(O1) = n(zt)t_zim+1Pm050t , (B.5)miwhere 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 thecharge distribution n(z),= f dt. (B.6)z’When m = 1, we getn(z)dt f n(zt) cosOdt = —E(zi) (B.7)OZ It—ziI = t—z12where E(zi) is the electric field normal to the surface. In general,I n(zt) 1 am n(zt)= (B.8)Thus, we have1 Om(z)R(1Z1 (B.9)u(01)= t9zror, using the rotational invariant expansion Eq. (3.2.1),0mm— (_)mI2m + 1 —m 8m(zi) (B.10)uj;w(01)— f°mm m! Q1L;/3 OzAppendix CNumerical considerationsIn this appendix, the computer algorithms used to obtain the wall-particle pair anddirect correlation functions are outlined. The computer programs used were written fora general multicomponent fluid composed of particles of C2 or higher symmetry. Thealgorithm used to solve the fluid structure near an isolated surface and between twosurfaces 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 long-range pair functions in the integral equation theories has been described elsewhere [15].The basic algorithm used to solve the wall-particle correlation functions is summarizedbelow.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 andPlischke [31] or Attard and Patey [26].Onn(—)3. Calculate 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(—1, z<d/2,h(z) = (C.1)10, z>d/2,is sufficient for the systems considered.Onn(+) Onn(+)6. Calculate the Fourier transforms of the projections hv;w (z), h0; (7. Calculate (k) using Eq. (3.2.18).149Appendix C. Numerical considerations 150.-Onn(+) Omm(+)8. Inverse Fourier transform 7lOv•w (k) to get iOu.w (z).9. Calculate z) and the projections of potential of mean force, w(z) [seeEq. (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) fromWth(z) + (1 — W) [q,3(z) + c(z)] , (C.2)where h(z), and c(z) are the most recent estimates and W isa mixing parameter in the interval 0 W < 1.12. Repeat steps 6—12 until converged. Convergence was obtained when the pair correlation function was constant to five significant figures.Step 11 is characteristic of the Picard iteration strategy [52]. A faster convergence strategy has recently been developed by Kinoshita and Harada [52]. When the wall-particlepotential depends on the surface and fluid structure in a self-consistent manner, an additional iteration cycle is required. At a dielectric surface, the algorithm is summarizedbelow.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), usingeither 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 wallparticle 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 usingthe iteration strategies of Kinoshita et al. [52]. When the jellium metal is immersed influid, the basic algorithm is summarized below.Appendix C. Numerical considerations 1511. Calculate the wall-particle pair and direct correlation functions at an inert wall asdefined above.2. Calculate the electron density distribution, electrostatic potential and effective potential at a jellium-vacuum interface using Eqs (3.3.33), (3.3.39) and (3.3.30) withext()= 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 functionsusing 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 potential for the jellium metal using ut(z).7. Repeat steps 3—7 until the electron density distribution and fluid structure areself-consistent.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Structure and properties of fluids near inert and metallic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Structure and properties of fluids near inert and metallic surfaces Bérard, Daniel Robert 1994
pdf
Page Metadata
Item Metadata
Title | Structure and properties of fluids near inert and metallic surfaces |
Creator |
Bérard, Daniel Robert |
Date Issued | 1994 |
Description | General reductions of the Ornstein-Zernike equation are given for multicomponent molecular fluids near an isolated planar surface and between planar surfaces. This allows integral equation approximations such as the hypernetted-chain (HNC) and reference hypernetted-chain (RHNC) theories to be solved numerically. Theoretical methods for treating multipolar particles near inert, dielectric and metallic surfaces are considered. The metal is represented using a jellium model together with density-functional (DF) theory and the two phases interact electrostatically. Mean-field theories which reduce 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 surfaces, for the electron density distributions. Explicit results are given for dipolar hard sphere fluids and for electrolyte solutions. The wall-induced fluid “structure”, electrostatic 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 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. |
Extent | 2515363 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-06-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0061608 |
URI | http://hdl.handle.net/2429/8750 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 1995-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
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
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0061608/manifest