On Surface Properties of the One-component Plasma

. We consider a plasma of point ions in the presence of a non-uniform neutralising background. This background. the source of an external field, may have some of its parameters (density, form of surface profile, etc) modified, as long as the total charge is maintained. By considering such modifications in the context of the density-functional formalism for the ions, we prove sum rules giving the first and second moments of the ion density p(z) in terms of other properties (bulk pressure and temperature derivative of surface tension). The Poisson-Boltzmann functional is considered in detail. We show that the first and second momenr conditions on p(z) are verified. We calculate p(z) exactly for this system, and also perform variational calculations: comparison shows the importance of respecting the asymptotic behaviour of p(z). Variational calculations have been performed, using the density-functional formalism in the square-gradient approximation. for systems with plasma parameter r from 1 to 10. For r > 3, important oscillations appear in the profile, as shown by recent Monte Carlo calculations. The profiles calculated variationally also show increasing oscillations. but are not in good agreement with the Monte Carlo results. The surface energies are poor even for r = 1 showing the inadequacy of the square-gradient expansion for this system.


Introduction
There has been a growing interest in the past few years in the surface properties of Coulombic systems such as molten salts, electrolyte solutions and liquid metals. The one-component plasma (OCP) is a simple and useful prototype of such Coulombic systems, whose bulk properties are now well established over a wide domain of temperature and concentration. The density profile and other surface properties, such as surface energy and surface tension, have recently received great attention (Ballone et al 1981a, Evans and Hasegawa 1981, Jancovici 1982, Badiali and Rosinberg 1982 and Monte Carlo computations are now available (Badiali er af 1983). In the case where the OCP is in contact with an impenetrable wall, comparison has been made with predictions of analytical methods such as the mean spherical approximation and the hypernetted chain equation (Badiali er a1 1983). In this paper we calculate, for small values of the plasma parameter, the density profile, the surface tension, and the surface energy of the OCP in a non-uniform background. This is done using the free-energy density functional formalism (FDFF). In 9 2 we give several general thermodynamic relations for systems of this type which can be derived quite simply in the FDFF. They follow from the fact that @ 1983 The Institute of Physics certain parameters of the background can be modified arbitrarily. Two of these are the classical analogues of the sum rules derived by Vannimenus and Budd (1974) for the jellium model of a metal surface: one is for the surface potential and the other for the rate of change of surface tension with bulk density.
These theorems are valid for any functional of the free-energy and this is illustrated in 0 3 by the behaviour of the exact solution of the Poisson-Boltzmann equation for the semi-infinite OCP. In 0 4, we consider variational calculations in the FDFF; instead of solving the equation for the profile directly, one may seek an approximate solution in the form of a parametrised trial function. We discuss the implications of our theorems for such an approach. The variational calculations are done by minimisation of the free-energy density functional in the square-gradient approximation. We calculate the surface density profile, the surface tension and the surface energy. Comparison with Monte Carlo results and concluding remarks are given in § 5 .

General theorems
We consiser a classical fluid of point ions with charge Ze and bulk density p b at temperature T , embedded in a non-uniform neutralising background of density n ( z ) such that n(z) + n (bulk density, n = zpb) for z + -cc and n ( z ) + 0 for z ---$ + x . This system is characterised by the plasma coupling parameter r = ( Z e ) * / a k T where k is the Boltzmann constant and a = (4np1,/3)-';~ is the ion sphere radius. Following the pioneering work of Hohenberg and Kohn (1964) and Mermin (1965) one can show that the free energy of the system can be obtained from the functional where a separation has been made (Evans and Sluckin 1980) between Coulombic and non-Coulombic contributions to the free energy. q ( r ) is the electrostatic potential related to p(r) and n(r) by the Poisson equation. The functional G [ p ] is a unique functional of p(r) and pis the position-independent chemical potential of the particles.
The minimum of R[ p] is the grand potential of the system and the equilibrium density satisfies the Euler-Lagrange equation (Evans and Hasegawa 1981) We now consider general theorems which are consequences of this approach. Since this model is quite similar to the uniform background model in the inhomogeneous electron-gas theory, some of these will be analogous to the sum rules derived in this case by Vannimenus and Budd (1974). However, one here considers changes in the free energy of ions instead of changes in the internal energy of electrons.
Consider an infinitesimal modification &(r) of the background profile such that the total number N of particles is unchanged. The modification 6n(r) will cause a change in ionic density dp(r) and possibly a shift in the chemical potential 6p. To first order, the change in the grand potential is where Nis the total number of ions, N = J p ( r ) dr.
Using the equilibrium condition (2) and the fact that J 6p(r) d r = 0 we get SR = -e d r q(r)Sn(r) -Nap. (4) In addition, at the minimum of R[ p] one has where Fo is the free energy for the homogeneous N-particle system, A is the surface area, and ythe surface tension. For any change holdingA and Nconstant, we have One can now compare this expression for SQ with the expression of the form ( 4 ) to deduce interesting relations.
Consider first the case of the semi-infinite OCP: the neutralising background of our system is a slab of density nb extending from z = -L to z = 0 with the faces normal to the z axis each having area A (i.e. in the limit Heaviside function). The area is supposed to be large enough so all properties are uniform in the x and y directions. Following Budd and Vannimenus we consider the change SQ associated with stretching the slab so that it extends from z = -L to z = SL while holding A , T and N constant. It is straightforward to show that the electrostatic term in (4) becomes where the new background density is nb -Snb such that to first order LSnb = nbSL, and q ( -L / 2 ) is the potential in the centre of the slab. On the other hand, equation (6) becomes (with A = 2A since there are two surfaces) where 6pb is Snb/Z and (ay/apt,)~ is the derivative of surface tension with bulk ion density. Since

P b ( a F O / a P b ) T ' NP/pb
where P i s the bulk thermodynamic pressure, we get where ALpb = N has been used.
have Then, inserting (7) in (4), and comparing the expression for SR with that of (8), we

L e [ q ( O ) -q ( -L / 2 ) ] -2 e r -L'2
Each member consists of a 'volume' term (proportional to L ) and a surface term, which must be separately equal, so that where it has been assumed that q(-L/2) tends toward its asymptotic value q ( -= ) faster than 1/15. Equation (9a) has been already given by Ballone et a1 (1981b) and provides an exact relation between the potential difference [ q ( O ) -q(x ) ] , which is a surface property, and the bulk quantity P . This relation can also easily be derived by a direct integration of the first equation of the BGY hierachy. As noted by Ballone et a1 (1981a) the total potential drop for this system, [ q( =)q( -=)I, is divergent because of the z -2 asymptotic behaviour of p(z) in the absence of a finite background density for z > 0 (the density p. of the vapour phase is strictly 0 in our system; in a real system, where f i is small but not vanishing, the total potential drop of course remains finite).
On the other hand, equation (9b) is a new result and can also be written, after insertion of the Poisson equation and partial integration, as This equation, at least theoretically, provides a way of calculating the surface tension directly from the density profile p ( z ) . It relates the change of ywith bulk density to the second moment of the ionic distribution. By comparison, equation (9a) is a condition on the first moment since integration of the Poisson equation gives f 0 Note that the relations (9) are valid when n ( z ) is a step function. However, from (4) and (6), similar results can be derived when the background profile has a more general shape.
In addition to (lo), there exists a relation between (dY/apb)T and y itself. One can show easily by scaling arguments that the surface tension can be written where & = m , with a representing the parameters a l . . . am specifying the background profile. For a step profile one has simply so that, using the definition of r, one finds by direct differentiation (13a) Consider now the case of an arbitrary background profile n ( z ) depending on a parameter a: (13a) no longer holds, but, differentiating (12a) with a constant, we find Let us consider the effect on the surface tension of a change of the parameter a, keeping A and Nconstant. According to (6) and (4) an ( The effect of a non-Coulombic interaction between the ions and the background can also be considered. Suppose that we have in the RHS of (1) the extra term where A represents a set of parameters {Al . . . An}. If A is modified, the condition of minimisation of R and (6) gives (n is unchanged) For example, if A is the radius of the Ashcroft pseudopotential, U = Ze2r-' e(), -r), aFo/aA is easily calculated and one finds Finally, the effect of a non-Coulombic external potential can be investigated; it adds to R a term drp(r)Vext(r; i b l . . . An).

I
If the parameters A, can be varied, e.g. to change the 'softness' of a wall, we have The relations (9) and (13)-(15) are exact regardless of the form of the functional, so do not provide a test of its validity. Note that P and ( d Y / d p b ) T are supposed to be calculated with the same functional as the moments of the profile in (9) and (10). To test the functional, one requires an alternative route to the quantities considered. such as a calculation of ( a y / d p b )~ from the Hamiltonian of the system that the functional is supposed to describe. It may also be noted that the theorems require for their proof an arbitrarily modifiable external field. For a metal, where both electrons and ions can rearrange simultaneously, such theorems do not exist. Indeed, in this case, one has to consider a functional which depends both on p(z) and n ( z ) (Evans and Hasegawa 1981). R must be a minimum with respect to these two profiles, hence 6R = 0. Equation (3) no longer holds and the theorems cannot be derived.
In the remainder of this work we are principally concerned with the relations (9) and (13). We shall use them as a test of a correct minimisation of the functional Q .

The local approximation and the Poisson-Boltzmann equation
The local approximation consists of taking G [ p] in equation (1) to be given by where g ( p ( r ) ) depends on p ( r ) and represents the free-energy density of a homogeneous OCP of density p(r). The surface tension may be defined as the surface excess free energy After integration by parts and use of the Euler-Lagrange equation (2) we get Thus y can be calculated without explicit reference to the nature of g ( p ) ; p(z) is, however, determined from (2), in which g plays a role. When n ( z ) is a step, the first integral in (18) vanishes.
The simplest local approximation, which should be valid for low densities, is obtained when one supposes that g(p(r)) in (16) is given by the free energy of a non-interacting gas. Then one gets , Ballone et a1 1981a immediately from (2) where the origin of the potential is at z = -C C . Equation (19), together with the Poisson equation, yields the Poisson-Boltzmann equation which, if n ( z ) is a step, is where p' ( z ) = dp( z)/dz, and K = (4xpb Z2e2/kT)"* is the reciprocal Debye-Huckel length.
In terms of the reduced length U = KZ and p ( U ) = p( u ) / p b , (20) simply reads For U > 0 this equation can be solved analytically (Ballone er a1 1981a) to give The value of p ( 0 ) can readily be obtained after some simple manipulations. From (21) we have Using (22) With this value it is clear that equation (9a)  In fact, (9a) could have been used to find p( 0). Let us now consider the verification of equation (96). Within the Poisson-Boltzmann approximation, the surface tension for a step background profile is simply which can be written, with the help of equation (19), as By integrating equation (23) we get where we have used equation (25) and the overall electroneutrality condition. Hence Therefore we get, using the definition of K, which clearly satisfies equation (9b).
For z < 0, the profile can be obtained only numerically. A method for generating it is given in Appendix 1. By inserting our solution into (29) we get y (compare Ballone et a1 (1981a)) in the form of (126): The excess surface energy U, is given by the thermodynamic relation us= y -T(ay/aT)v= y -T(ay/aT),= -2y+ 3 p b ( a y / a~b )~. (32) The use of (dy/aT), instead of (dy/dT)Vis discussed after (36). Equation (13a) has been used in deriving the second part of equation (32). It follows from the form (12b) which can be also obtained by the direct definition of the surface energy.
More generally equation (9a) can be used in order to get the exact value of p ( 0 ) in the local approximation, for any free energy density g( p ) . According to equation (2) we have and hence, using (9a), (dg/dp),=o = P -Z e q ( O ) = pp/pb.
If g is known for the homogeneous system as a function of p or r, both sides of the equation are calculable. Then one inserts the bulk density in the expression g ( p ) / p and finds the density for which dgldpis equal to it. This is the value of p ( z ) for z = 0.
For example, when 1 S r S 160, Slattery et a1 (1980) have given an expression which fits very well the Monte Carlo results for the free-energy density of the bulk OCP g(p)/pkT = -0.89752 r + 3.78176 r1'4 -0.71816 Y1' + 2.19951 In r -3.30108. In the weak-coupling limit 0.1 S r S 1 we have derived (36a) g(p)/pkT = -0.64986 -0.33676 r -0.19797 r2 + 0.04929 r3 + 3.012285 In r (36b) by numerical integration for the internal energy (for 0.1 S r s 0.6 we have used the HNC results and for 0.6 6 I' s 1 the Monte Carlo results). With the use of equations (36) in equation (35) we find, for a given value of r b = T(z+ -m ) , T(z = 0) and hence p(0). Results are shown in table 1. The decrease of p(0) with r, corresponding to a displacement of the ion profile into the jellium and implying oscillations for z < 0, is to be noted. As noted after equation (32), the thermodynamic formula for U, in terms of yrequires the temperature derivative of y with V constant, whereas we have taken the derivative at constant p b , To justify this, we note that we consider a system of N positive ions in a volume V ( N -=, V+ =, but NIVfinite) and in the presence of a neutralising background whose rigidity is maintained by some external constraint. When the temperature Tvaries the background is not affected and hence the bulk value nb remains unchanged. Since electrostatic considerations require that po = Znb, pb is also constant. It may be noted that in the recent Monte Carlo computations (Badiali er a1 1983) the thermodynamic limit of the finite physical system is supposed to be taken as follows: Consider N ions in a spherical volume V = 4nR3/3, at the centre of which is a sphere of radius Ro containing neutralising charge at density nb, so 4xRinb/3 = Z N . Then, as V and N become infinite, the ratio Ro/R is to be held constant. With the constancy of N/V this implies the constancy of ZN/Ri and hence n b , On the other hand, if the profile of the neutralising background is not a step, ambiguities arise.
In the local approximation it is easy to verify that y -T(dy/dT),, correctly gives U, since the surface tension is given by (17). Now g ( p ) , the free-energy density of the homogeneous system of density p, is related to the energy density of such a system u ( p ) by 4 P ) = -"/dT).
(37) Therefore y -T(aY/aT), = 1 [ u ( p ) -u(pb) e(-z>l dZ + I j i Z Pn , q ( z ) dZ which is clearly U,. Although there is a T dependence through p, terms in ( d p / d T ) (6R/6p) vanish because the minimisation process makes SsZ/Sp vanish.

The square-gradient approximation
The simplest improvement to the local approximation is to add a square-gradient term to the expression for G [ p ] : In the limit where the ionic density varies slowly and exhibits only small departure from the bulk value p, it can be shown (Evans 1979) that the quantityg2(p(r)) is determined by the second moment of the direct correlation function c ( r ) of the homogeneous OCP. More precisely where a ( p ) is the coefficient of q2 in the expansion of the non-Coulombic part of c ( q ) , the Fourier transform of c ( r ) . Following the scheme used by Evans and Sluckin (1981) to calculate c ( q ) it can be shown that where K o ( r ) is a dimensionless parameter determined by fitting the isothermal compressibility of the OCP. Using the fits of (36) we have (40) &-) = -0.0246 + 0.89802 r + 1.31982 r2 -0.59154 r3 &r) = 2.3934 -2.0484 r-34 -0.32916 r-j4 + (i.6008/1-) 0.1 s r s 1 1 s r s 160.
Although the fits (36b) and (40) are strictly valid for > 0.1 we shall use them even for r < 0.1, since we can verify from (38) that the region of very small r (when p ( z ) e p b ) does not give any significant contribution to the free energy when rb 2 1.
The truncated gradient approximation is now fully specified and the problem is to solve the non-linear differential Euler- Lagrange equation ( 2 ) . This is a difficult task and one usually prefers to assume some parametrised form for the profile and choose the parameters to minimise the surface excess free energy per unit area (R -Rb)/A.
We first considered a three-parameter class of trial functions  (22) and (25)) p ( z ) = ~( K Z + (2e)''2)-2 for z > 0 and the first part of (41) for z s 0, demanding that p be continuous, have continuous slope, and satisfy electroneutrality; the remaining free parameter is chosen variationally. The profile obtained is much better (figure 1) and y/apbkT = -2.157 (3F)-'''. This shows the importance, at least for small r, of taking into account the z-' asymptotic behaviour. We then consider the second class of functions: where A , B and zo are calculated as functions of (Y, p and y in order to satisfy continuity of p and p' at z = zo and the electroneutrality condition, so (42) is, like (41), a threeparameter trial function. The optimised value of p may be compared with the exact asymptotic value r (since K ' = 3 r if we take the ion radius as a length unit). We now use these trial functions with the more general density functional including gradients (cf equations ( l ) , (36) and (38)). We show in Appendix 2 that in this case too the parametrisation (42) correctly describes the asymptotic behaviour of p(z) for z > 0.
The surface tension may be written Y = -2 ( z p -n ) y ; ( z ) c l z + j [ g ( p ) -g ( P b ) e ( -Z ) l d z i jgZ(p)(vP)'dz where yes is the electrostatic contribution, yl comes from the non-gradient term and yl, from the square-gradient correction. Minimised density profiles are given in figures 2 and 3, and numerical results in tables 2 and 3. We have verified that the trial functions z l o (41)givehighervaluesof y(forinstance, y/apbkT = -0.806forr = 1). Thispoint shows that the class of function (42) remains better far from the Poisson-Boltzmann limiting case (this is not a general statement, as we shall see in the case r = 10). As expected, the square-gradient contribution yl, increases with r . For I-= 1 and r = 2 , it gives a negligible or very small contribution for the surface tension but, because g ( p ) for the OCP is different from the Poisson-Boltzmann expression k T [ p ln(A3p) -11, surface tensions are quite different from the Poisson-Boltzmann values. Similarly, U, is quite different from -+y(see below for calculation of U,). For r = 3, oscillations become significant in the profile and yI, is no longer negligible. For = 5 the gradient term becomes 80% of the non-gradient one and twice the value of y, as important oscillations in the profile appear. These effects can also be seen from the value of the profile at z = 0, in comparison with the local approximation results given in table 1.
On the other hand, when r increases the behaviour of p ( z ) for z > 0 differs more and more from the Poisson-Boltzmann behaviour, as we can see from the variational parameter p which takes the values 1.097, 2.475, 4.225 and 8.209 for r = 1, 2, 3, 5 instead of the asymptotic value p = r (at the same time zo remains in the vicinity of One can compare our surface tensions with those calculated by Ballone er a1 (1981a) who have solved the differential equation (2) in the local approximation, using for g ( p ) a mean spherical approximation result. As shown in table 2, the agreement is satisfactory as long as there is no significant contribution coming from the gradient correction. This comparison tends to prove that our minimisation process is correct (at least for r < 3) which means that the class of trial functions (42) permits a good approximation of the exact profile.
We can also use other tests. The first one is the balance between the various contributions yes, y L , x, , to the surface tension in the density functional formalism.
Starting from (43) and making explicit the relations between the different contributions to y implied by the differential equation (2), it is possible to show (Evans and Sluckin 1980) that Comparing with (43) we get The results of table 2 satisfy (45) quite well. We must note, however, that (45) is a necessary but not sufficient condition for a good minimisation since a single-parameter profile (for instance a single exponential form) can be easily shown to satisfy (45), although this one-parameter profile is not at all the exact solution. As indicated in 0 2, the two sum rules (9) can also be used as a test of minimisation since they are valid for any approximation of Q [ p ] as long as the true minimum is obtained. The sum rules are checked in table 3. To evaluate the left-hand side of (9b) we require the isothermal derivative of y with respect to the bulk density p b . The derivative ( d~/ d p b )~ (or (dy/dT),) has been calculated by changing pb (or T ) by a small amount, while holding constant the parameters of the profile (since SQ/Sp = 0 for the equilibrium density profile, variations of the optimal parameters accompanying infinitesimal changes in pb or T can be neglected in first order). The accuracy of the numerical differentiations can be checked by the relation (13a) which gives a thermodynamical link between y , ( d y / d T ) , and ( d~/ d p b )~, validwhenever R isminimised within a givenclass of functions.  From table 3 we see that, while sum rule (sa) is rather well satisfied, the agreement is much worse for the sum rule (9b). We have noted that these sum rules are related respectively to the first and second moment of the profile for z < 0, so (9b) is more sensitive to the asymptotic behaviour for z < 0. Because of this sensitivity, it is unfortunately not possible to calculate (dy/dfi)Tfrom the Monte Carlo profile p ( z ) obtained by Badiali et a1 (1983). Note that in the Poisson-Boltzmann case, the precision on the sum rule (9b) was only 8% while the minimised value of ywas better than 0.1%.
Unlike equation (45), sum rules (9) can be shown not to be satisfied by the minimisation of a single exponential type profile. show that within statistical errors this dependence is negligible. On the other hand, the discrepancies between Monte Carlo and density functional results are not due to inadequacies in the trial functions since the Monte Carlo profile can be satisfactorily fitted within the class (42). The problem indeed seems to be that the truncated gradient expansion itself is incorrect even for this small value of r. As noted by Alastuey and Levesque (1983) for the two-dimensional OCP, the non-local terms neglected in the square-gradient expansion (38) may have the same order in I as the first terms considered.

Comparison with Monte Carlo results and concluding remarks
We have also done the minimisation work for r = 10 keeping the same functional and the two classes of functions (41) and (42). The class (41) appears to give the lowest minimum for the free-energy (the asymptotic tail for z > 0 is now making small contributions) but disagreement is important with Monte Carlo results as we can see from figure 4 and table 4. It is possible to increase the oscillations by multiplying the gradient term by an adjustable parameter: If the value of k is chosen in order to get the Monte Carlo results for the surface energy, we find k = 0.535. We see in figure 4 that the agreement with the Monte Carlo profile becomes satisfactory now. But if the background profile becomes of exponential form (characterised by a width A) we observe again important deviations in the surface energies (table 4). This fact may show the inadequacy of the ad hoc correction, although there may be also a problem in the definition of the thermodynamic limit for the Monte Carlo computations when the background profile is not a step (Badiali et a1 1983).
Anyway the adjustable parameter k is r dependent.
This failure of the truncated gradient expansion, even when no oscillations are present in the profile, is of great importance since this approximation has been used in z l o  other Coulombic systems such as liquid metals in order to obtain the surface tension. We see that the discrepancies between theory and experiment, which have been attributed to an oversimplified treatment of the electron-ion interaction (Evans andHasegawa 1981, Goodisman and, may come also from the truncated gradient expansion itself. A l . We give the values of the first aj calculated from various values of k, as well as the slope at 0, calculated as The exact value of this quantity, from the solution for U > 0, is -4(2e)-3'2 = -0.3155.
The convergence is quite satisfactory. We can then use the solution to compute the integral needed for the dimensionless part of the surface tension y, equations (27) and (28). Thus becomes -1.08206 using the last line of table A l , and y = -2.1641 apbkT (3r)-".

Appendix 2. Asymptotic behaviour in the gradient approximation
When the square-gradient approximation (38) is used for G [ p ] , the Euler-Lagrange equation (2)  and we know (see equation (22)) that the asymptotic behaviour of the solution is p ( z ) z -~, We argue that we have the same algebraically decay for the solution of (A2.4). Indeed when p ( z ) -Y 2 kT In p( z ) / p ,q ( z ) -In z and it is easy to check that the two gradient terms have a faster decay, in Y 4 I 3 .
We can also look at the asymptotic behaviour of p ( z ) in the bulk phase, i.e. when z+ -=. Differentiating (A2.1) twice in order to introduce the Poisson equation and linearising with respect to h ( z ) = p ( z )p b , it is clear that we get a linear differential equation with constant coefficients. The general solution is then (A2.5) p ( z )constant em cos( yz + q ) .