Structure of the LiquidiVapour and Liquid/Solid Interfaces BY JOHN W. PERRAM AND LEE R. WHITE* Department of Applied Mathematics, Institute of Advanced Studies, Australian National University, Canberra, ACT. 2600, Australia Received 8th January, 1975 A theoretical method of determining the density profile p(')(z) at liquid/vapour and liquid/solid interfaces is developed, based on a procedure of Helfand, Frisch and Lebowitz. The density profile for a given fluid is shown to be related to the value of the pair distribution function &)(r) for a two component fluid in the limit that the molecular radius of the second component becomes infinitely large (i.e. macroscopic) and its density vanishes. The density profile for hard spheres against a hard interface is obtained in the P.Y.approximation. The extension to the liquid/vapour interface is discussed and the free-interfacial density profile for a fluid with an attractive inverse sixth power pair potential is displayed. A comparison is made with other modern theories. 1. INTRODUCTION The structure of fluids in inhomogeneous regions has been intensely studied. lml The principal methods of investigation being quasi-thermodynamic arguments, 1-6 integral equations 7-1 and computer simulations. l2-I7 Fundamental to all three approaches is the singlet density p ( l ) (z) which varies with some co-ordinate z in the interfacial region. All three approaches have been described by Croxton,18 but a brief description of each is given here. 1.01 :8 -6 -4 -2 0 2 4 6 8 Z FIG. 1 .-Typical monotonic density profile obtained by the constant chemical potential approach.T/Tc = 0.50 (A), 0.80 (B), 0.98 (C). The quasi-thermodynamic theories are familiar to workers in colloid science and electrochemistry, although the crude formulae used there have been considerably refined. Fundamentally, a local chemical potential p(z) = p[p(z), z] is defined. It is a function of the local density p(z) and is usually written p = p' + kTln a[p(z)] + u(z) 2930 STRUCTURES OF INTERFACES where po is at some reference state, a[p(z)] is the local activity, usually calculated from one of the hard sphere equations of state,lg* 2o and u(z) is the interaction of an atom at z with the inhomogeneity (and hence is a functional of ~ ( 2 ) ' ) . By requiring p(z) to be constant, a non-linear integral equation for p(z) has been ~ b t a i n e d , ~ ' ~ which has been solved to yield p(z).For all theories of this type, p(z) is a monotonic function of z,18 as shown in fig. 1. The starting point for most of the integral equation theories is the first order Born-Green equation, relating p'l' ( r ) to the (inhomogeneous) pair correlation function pC2) (rl, r,), viz. kTVlp'l' (u)+ JVl@(lrl - - U ' ~ I ) ~ ' ~ ' (rl, r2)drz = 0. (1 -2) Various closures 3* * * l1 permitting solution of this equation have been suggested. Toxvaerd writes Pt2) (r1, r2) = aP(2) 0.12lp(l) (0 + (1 - 4P'2' (r12lp(l) @,I) (1 -3) where p ( 2 ) (r121p) is the radial distribution function for the unstable liquid at uniform density p . With a further approximation on pc2) (rlp), eqn (1.2) is solved to yield a monotonic profile for p'l) (r) (see fig.1). The approximations used are difficult to assess. Nazarian has used two alternative closure approximations (1.4) and p'2' (rl, r2) = p&m (rl2)e(z1 +z2)+p$%oUR (r12)@-(zl +z2)] I 2 FIG. 2.-Oscillatory density profiles obtained by solving an approximate BBGY equation for p(') (2). I, Nazarian approximation (1.4) ; 11, Nazarian approximation (1.5) ; 111, Croxton and Ferrier " bootstrap " function solution.J . W. PERRAM AND L. R . WHITE 31 and solved (1.2) for a liquid argon potential (90 K) to obtain strong density oscilla- tions as shown in fig. 2 for both closures. Croxton and Ferrier * have used an unusual closure relation where gCz) (rl, v,) is replaced by its bulk liquid value and a function @(z2) is introduced multiplying the potential <D(r12) in (1.2).The function is chosen so that it effectively creates a free surface by cutting off the interaction between par- ticles as a function of the z component of the second particle. Croxton and Ferrier apply this procedure to a liquid argon (at its triple point) potential to obtain the mild density oscillation as shown in fig. 2 which should disappear as the critical point is approached. Although both methods give density oscillations, the method of Nazarian shows much larger effects. A controversy arises. Does p{l) (z) exhibit an oscillatory nature near the interface and if so how is such behaviour modified by the approach to the critical temperature ? As has frequently been the case,21* 22 when approximate statistical theories conflict, the best way of testing them is by some type of computer simulation.This has proved invaluable in the study of homogeneous fluids,18 but removal of one of the periodic boundary conditions means that this is not as straightforward as it might seem. The first study l3 was a molecular dynamics simulation of a two dimensional interface of particles confined to the surface of a cylinder. A random field was applied to particles at the base and the transition profile generated. Again, density oscillations were reported. Another method is to simulate a thin film,15 so that two interfaces are generated. Liu l7 has recently “ created ” a single interface by placing an attractive potential at one end of the computer “ box ”.In this way, condensation is forced on the system and again an oscillatory profile is obtained. However, Rowlinson l 6 reports much less pronounced oscillations. We present here a theory in which we aim to unify the various approaches to the liquid surface. The method is based on an idea of Helfand, Frisch and Lebowit~.,~ Consider a two component system where the interaction between type 1 particles is via a pair potential u(r12) and between type 1 and 2 particles by a potential @(z) where z = Ir, -r2( - B’ where 0’ is the hard sphere radius of type 2 particles. In such a system the thermodynamics is completely specified by a knowledge of the distribution functions g$$), &), g\$) (z) all three of which are functions of the number densities p and p’ of the two particle types.If we now take the limit u’--+co while taking p‘--+0 then we have effectively created an interface at z = 0 of particles of type 1 with an interaction potential @(z) with the interface. One has then for the density dis- tribution of type 1 particles away from the interface This method produces density oscillations as well. where gC2) (rlJ is the bulk distribution function at density p of particles of type 1. If we took as the interfacial potential @(z) = co if z < 0 then the limiting procedure of eqn (1.6) would yield the density profile for the liquid against a hard impenetrable interface with which it interacts via the potential $(z). Provided 4(z) has a short range repulsive component (i.e. the liquid particles have = 4(z) if z > 0 (1.8)32 STRUCTURE OF INTERFACES some sort of repulsive core, it is obvious that such an interface, with p"' (2) given by (1.6), will exhibit density oscillations.These oscillations will be smoothed by softening the repulsive core but will be little affected 2 5 by changes in the attractive part of 4(z)* The preceding argument does not enable us to say much about the liquid/vapour interface, however, and it is at this interface where the controversy exists. Suppose that type 2 particles are vacuum bubbles of diameter o', assumed large. Let us consider a liquid particle near the bubble. The force of interaction of the particle with the interfaces is always, for real fluids, in the direction of the bulk liquid. This may be interpreted as a repulsive force between the bubble and the particle.The potential energy of interaction @(z) is plotted schematically in fig. 3, together with its Boltzmann factor. @(z) must be constant as Izl-+ co the difference @( - co) - @(a) is the molecular heat of vaporization and is a function of density and temperature. 0 2 3 FIG. 3.-A schematic representation of the potential of interaction @(z) of an atom in a liquid with the rest of the liquid as it moves through the liquid/vapour interface (A). The associated Boltzmann factor exp[- p<D(z)] is also displayed (B). C represents the heat of vaporization. A soft potential @(z) of the form shown in fig. 3 would lead to a smoothing of the density profile. As the critical temperature is approached from below the heat of vaporization will tend to zero and the will be a temperature (5 T,) at which the density profile becomes monotonic.With the above approach, we can see how the density profile is a function of the softness of the interfacial interaction potential. The choice of a hard potential @(z) or a closure approximation which amounts to such a choice will give rise to strong density oscillations. Correspondingly, a closure which amounts to the choice of a soft potential @(z) will give rise to monotonic or, at least, weakly oscillatory density profiles. The approach outlined above can be applied to the pair distribution function pC2) (zl, z2, r12). If &)I (zi, z2, r12) is the triplet correlation function in the mixture between a type 2 particle and two type 1 particles, thenJ . W.PERRAM AND L . R . WHITE 33 In the mixture, the triplet correlation function g$$)l can be replaced by the super- position approximation 2 5 Therefore, in the limit, where we have invoked the results (1.6) and (1.7). Thus the result (1.11) for the pair correlation function at an interface which was suggested by Green 26 and used by Berry and Reznek 27 can be seen as an approximation of the same order of mag- nitude as the superposition approximation in two component mixtures of very unequal size. Further since the approximation (1.1 1) is now couched in the familiar terms of a superposition approximation, one can systematically correct the approxi- mation in a statistically rigorous way.28* 29 The degree of approximation can also be tested in an average way by invoking the relation between neighbouring orders of the hierarchy of distribution functions. Here the appropriate relation is 9131’1 (21, 22, Y l 2 ) = 9122) (2119122) (22)gW (r12).PC2’ (21, 22, r12) = P‘l’ (Zl>P“’ ( z M 2 ’ ($12) (1.10) (1.11) - 1 = 1 P‘l’ (22) W2’ (21, 22, G2)- w 2 (1.12) which should be obeyed for all zl. Thus f ( z J = -J P(1)(z2)c9(2)(% 2 2 , r 1 2 ) - 11 dr2 V may be plotted and its deviation from unity used as a measure of the degree of approximation involved in g(2) (zl, z2, r12). In the next section, we exhibit density profiles for a number of interfaces in the Percus-Yevick approximation. 2. PERCUS-YEVICK EQUATION FOR THE DENSITY PROFILE We consider an M+ 1 component system of particles, with numbers N j of type j , j = 1, . . .M+ 1, enclosed in a volume V.Each particle has a core of diameter X,, so that the distances of closest approach Rl, are assumed to be R,, = +(& + Rj). It is our intention to allow &+I to become much larger than other R,, at the same time allowing pM+1 = N M + l/Vto tend to zero. This creates an interface in the form of a macroscopic spherical cavity. The particles interact through potential functions Uafi (r) of the form Uafi(r) = +a r < Ra, = u$’(r) r > Rap We investigate various forms for the Uafi (r). The total correlation functions h.8 ( r ) and direct correlation functions cap ( r ) are connected by the Ornstein-Zernike equations 2o M+1 y= 1 ha,fi(r) = ca,p(r)+ c P y J ca,y(~s~)hy,fi(~r-s~> ds (2.1) which, when supplemented by the Percus-Yevick approximation ha,fi(r) = - 1 +~2fi(r)/ll -exp(uap/WI (2.2) allow (in principle) the calculation of all the correlation functions.For the case when u$) ( r ) = 0, these equations have the happy feature of an analytic solution.30 Let us focus attention on h ~ + l , j . If R M + ~ S R,, then the distribution function 59-B34 STRUCTURE OF INTERFACES g M + 1 , j = 1 + hM+ l , j , when multiplied by p j = N,/Y is a measure of the variation in density as we move away from our central (M+ 1)-type particle. The question as to when our ( M f 1) particle becomes macroscopic is a matter of semantics, but we shall study this by considering a range of values of up to a factor of 40 times other particle diameters. A very convenient numerical method has been given 24 for the numerical evaluation of the functions hij for hard sphere mixtures, and there is no doubt that h M + l , j ( r ) will be oscillatory.It is however, reasonable to ask whether these oscillations are a result of the hardness of the cores. Fortunately, methods have been given 31 for the easy modification of the h a p ( r ) to account for the softness of the core. An effective hard sphere diameter RZj is defined by R,*B = (1 - exp[ - uaB/kT]) dr. (2.3) So" Now for all potentials, although the functions g a b (Y) may be discontinuous, the function Ydls(d, defined by Y a p ( r ) = gus (Y) ex~[ua~lkTl Y,*B(d = -c,*,<d r < K s (2.4) is continuous. For effective hard spheres = s,*a(r> r > R$. Then, gives good estimates of the correlation functions as a function of distance.Fra. 4.-The mixed radial distribution function glz(r) (i.e., density profile) for a system of hard spheres (po3 = 0.85) against a hard sphere (p' = 0) with u' = Ma.J . W. PERRAM AND L. R. WHITE 35 From the above, we can see how surface density oscillations can arise, even when the j, N+ 1 exclusion is quite soft (of the order of 1-2 f i j values). If Rj < RM+l then RT, M+l % Aj, M+f very closely, and so the functions y , y* are essentially the same. Thus outside the range of the soft potential the form of the functions has (r), hzp (r) are the same. We first examine the limiting process of eqn (1.6) for a system of hard spheres, bulk density po3 = 0.85, against an impenetrable interface. We obtain, numerically, the solution of eqn (2.1) in the Percus-Yevick approximation for a mixture of hard spheres at the above density and hard spheres of radius o’(o/o’ 4 1) at zero number density.In fig. 4, the value of hI2 (2) is plotted for ~ / o ’ = 1/40. The results for a/o’ = 1/25, 1/20, 1/15 are scarcely distinguishable from the exhibited curve. We can examine the limiting process in more detail however. 4.0 I 1 I L 1 0 0. I 0 . 2 (ula’) FIG. 5.-The peak heights of the first three maxima in the density profile for a hard sphere system as a function of u/u’. A, first maximum ; B, second maximum ; C, third maximum. In fig. 5, the height of the first, second and third peaks are plotted as a function of (ole'). For the mixed hard sphere system it is easily shown that where n 3 q = -pa . 6 With pa3 = 0.85 we see that lim h1,(a/2) = 5.137.p’+O U’-+ Q) A linear extrapolation of the first peak height for the points corresponding to c/o’ = 1/40, 1/25, 1/20, 1/15 passes through c/o’ = 0 at h12 (1) = 5.14. Similar extra- polations for the second and third peak heights also seem indicated due to the linearity in the region o/o’ < 0.1. Thus very accurate numerical estimates of the limit (1.6) may be obtained by linear extrapolation from two or three h12 (2) curves calculated with (c/o’) < 0.1.36 STRUCTURE OF INTERFACES 3. THE LIQUID/VAPOUR INTERFACE We combine here the idea of the self consistent energy profile 4(z) mentioned in the Introduction with the equations (2.3)-(2.5) which give a prescription for the modi- fication of interfacial density. Suppose the small particles interact via the potential ull(r) = +a r -= Rll = -&(R11/r)6 r > R11 and suppose that the density profile pl(r) is given by Pl@) = Pls221’(r) then the energy of interaction of a particle 1 at r with every other particle is n Then +(P) may be defined as 4(r) = V(r)- V(m).(3.2) (3.4) If we replace u; by 4 in eqn (2.3) to calculate Rzl and use (2.5) to determine g$i), we obtain vW) + +(r) = p 1 J y~(rz>exp[ - BM,)I Id 1 1(1r - r2 ~)g\?(ir - r2 I) (3.5) vz which is the self-consistent equation for the energy profile +(T). The key question is however whether or not iteration of the equations produces a convergent answer, and what is the nature of p l ( r ) derived from it. ( z l 4 FIG. 6.-The density profile p(l)(z)/p~ at the liquid/vapour interfacelfor a system of particles (po3 = 0.85) with a hard sphere core plus an attractive interaction -E(o/r)6, for ElkT = 1.5.J . W.PERRAM AND L. R. WHITE 37 As a rough guess, we would suspect that at low temperatures, as measured by e/kT, stable solutions could be obtained, because the self-consistent potential will be quite hard. Such is indeed the case, for fig. 6 shows p l ( r ) in the vicinity of the interface for EIkT = 1.5 and density plR:l = 0.8. (The results of Lee, Barker and Pound, although for a different potential, correspond roughly to e/kT = 1.4 and pRfl = 0.81). We note that the density exhibits pronounced oscillations. We have also studied i$kT = 1.25, at the same density, and find similar, but less pronounced oscillations. T. L. Hill, J.Chem. Phys., 1952, 20, 141. I. W. Planer and 0. Platz, J. Chem. Phys., 1968,48, 5361. S . Toxvaerd, J. Chern. Phys., 1972,57,4092. S. Toxvaerd, J. Chem. Phys., 1971, 55, 3116. €3. W. Ninham and J. Mahanty, J. C. S. Faraday I& 1975,71. G. H. Findenegg, personal communication. J. G. Kirkwood and F. P. Buff, J. Chem. Phys., 1949,17, 338. * C. A. Croxton and R. P. Ferrier, J. Phys. C., 1971,4,1909; Phil. Mag., 1971,24,489. !a C. A. Croxton and R. P. Ferrier, J. Phys. C., 1971,4, 1921 ; Phil. Mag., 1971, 24,493. lo J. W. Perram and L. R. White, Nature, to be published, l1 G. M. Nazarian, J. Chem. Phys., 1972,56, 1409. l2 S. Toxvaerd, Mol. Phys., 1973, 26, 91. l3 C. A. Croxton and R. P. Ferrier, J. Phys. C., 1971, 4, 2447; Phys. Letters, 1971, 35A, 330. l2 S. Toxvaerd, Mol. Phys., 1973, 26, 91, l4 J. D. Bernal, Proc. Roy. SOC. A, 1964, 280,299. l5 J. K. Lee, J. A. Barker and G. M. Pound, J. Chem. Phys., 1974,60, 1976. l6 J. S. Rowlinson, personal communication. l7 K. S . Liu, J. Chem. Phys., 1974,60,4226. l9 L. Tonks, Phys. Rev., 1936,50,955. 2o H. Reiss, H. L. Frisch and J. L. Lebowitz, J. Chem. Php., 1959, 31, 369. 21 R. 0. Watts, Mol. Phys., 1974, 28, 1069. 22 A. Ralman and F. H. Stilfinger, J. Chem. Phys., 1971, 55, 3336 ; 1972,57, 1281. 23 E. Helfand, H. L. Frisch and J. L. Lebowitz, J. Chem. Phys., 1961,34, 1037. 24 F. Kohler, J. W. Perram and L. R. White, Chem. Phys. Letters, 1974, in press. 25 S. A. Rice and P. Gray, The Statistical Mechanics of Simple Liquids (Interscience, New York, 26 H. S . Green, Handbook of Physics, 1960,10,79. 27 M . V. Berry and S. R. Reznek, J. Phys. A, 1971,4, 77. 28 J. S. Rowlinson, Mol. Phys., 1964, 6, 591. 29 M. J. D. Powell, Mol. Phys., 1964, 7 , 591. 30 R. J. Baxter, J. Chem. Phys., 1970,52, 4559. 31 H. C. Anderson, D. Chandler and J. D. Weeks, J. Chem. Phys., 1972, 56, 3812. C. A. Croxton, Liquid State Physics-A Statistical Mechanical Introduction (Cambridge U.P., 1974). 1965).