Density functional calculations for liquid metal surfaces

Various forms for the ion density profiles of Cs and Na have been tested in the density functional formalism. The large oscillations expected for Coulombic systems are largely suppressed by the pseudopotential. While oscillatory profiles give some improvement in surface tension, the values obtained are still about twice the experimental ones. It is concluded that the treatment of the pseudopotential or the gradient expansion of the functional is responsible for the unsatisfactory results obtained.


Introduction
The local density functional formalism (Evans 1979) , extended to Coulombic systems by Evans and co-workers (Telo da Gama and Evans 1979, Telo da Gama et a1 1980, Evans and Sluckin 1981, can provide a general and powerful description of the thermodynamics of such systems. For an inhomogeneous system, it can yield surface properties and density profiles for the different species present. The difficulty lies in finding a good approximation to the density functional. For the vapour-liquid interface of a simple liquid metal, Evans and Hasegawa (1981) proposed and tested a local density functional which ( a ) includes the long-range Coulombic interactions exactly, ( b ) uses a truncated density gradient expansion for the interacting system, (c) treats equivalently electrons and ion cores (allowed if one invokes the pseudopotential formalism for the non-Coulombic part of the electron-ion interaction and treats the pseudopotential in first order). The functional is chosen to yield the correct results for the homogeneous system (no gradients). Since the functional gives the free energy of the inhomogeneous system in terms of the density profiles, minimisation of it yields the surface tension as well as the profiles. Instead of solving the resulting differential equation directly, one can insert parametrised forms for the profiles and minimise the functional with respect to the parameters.
Using one-parameter exponential forms for both ionic and electronic profiles, Evans and Hasegawa (1981) found surface tensions about twice as high as the experimental values for the alkali metals. They argued that this was due to the first-order treatment of the pseudopotential. However, since the first-order theory has successfully been applied to the surfaces of solid metals (Lang 1973), it seems useful to examine the effect of the choice of the ion profile in the present theory, before turning to a more exact one. Since the variational principle is the minimisation of the surface free energy, calculated 0022-3719/83/061143 + 10 $02.25 0 1983 The Institute of Physics 1143 surface tensions are necessarily too high, and can only be improved by the use of more flexible trial functions for the profiles. The work reported here involves the use of oscillatory profiles for the ions; it is known that substantial oscillations exist in the density profiles for related Coulombic systems (Badiali er a1 1983). Experimental results (Lu and Rice 1980) may also indicate oscillatory profiles for the liquid metal surface. For the surface of asolid metal, it is also clear (Lang 1973) that there is a substantial displacement of the electronic profile relative to the ionic profile, with important physical consequences. Such a displacement is not allowed by the exponential profiles so far used for the variational calculations. It is conceivable that this lack is responsible for the high surface tensions obtained.
Trial functions for the ionic profiles, allowing for both displacement and oscillation, have been tested for this density functional. These calculations are presented in § 2 . The improvement in surface tension is not substantial. The profiles obtained imply that oscillations are not very important. In P 3 we ascertain, by changing the pseudopotential radius, that at least part of the reduction of oscillations compared with other systems is due to the pseudopotential, which removes part of the Coulombic attraction between electrons and ions. Another reason for our finding smaller oscillations is the width of the determined electronic profile, as the calculations of this section show. Our conclusions about the validity of the functional considered, and about the accuracy of the trial functions used to minimise it, are summarised in § 4.

Improved ionic profiles
The density functional proposed for the liquid metal may be written with the various contributions defined below, using atomic units. The electrostatic contribution is Yes = 4 j j d z d r ' [ p i ( z )pe(z)j[pi(Zt)pe(Z'11 / r -1 -1 and depends on both the ion density profile pi and the electron density profile pe (both functions of the coordinate z normal to the interface). The electron density profiles must satisfy -(3/4n)(32ne) li3 -0.0474 -0.0155 In (3r~he) "7 ( 2 ) the second integral being the contribution of bulk densityn,. It includes the conventional terms for electronic kinetic, exchange, correlation, and inhomogeneity energies. The non-Coulombic part of the ion-electron pseudopotential contributes where we use a pseudopotential of Ashcroft form: -r-' + Vp with I The values used for Rc are given in table 1. The contributions ye[ and yps are commonly used in the density functional theories for metal electrons (Lang 1979). The remaining terms in yinvolve the ions alone. In terms of the plasma parameter r = (R,kT)-', where 4npiRs3/3 = 1, Evans and Hasegawa (1981) proposed two terms. The first is based on the free energy of the homogeneous system: where r b is r for p, = ne, li = -0.89752, 6 = 3.78176, c = 0.71816, d = 2.19951, and 2 = 3.30108. The second is proportional to the square of the density gradient: where the value of xo is calculated from the Monte Carlo results for the compressibility. Evans and Hasegawa (1981) used a value of xo = 1.53, except that several different values were tried for Na. We have modified the calculation by calculatingxo as a function of local density rather than taking a constant value corresponding to bulk fluid. For the large densities considered here, corresponding to r b % 1, this makes only a minor change.
In table 1, we give the parameters obtained by minimising y for the trial function used by Evans:

I
The values of &and E found by Evans and Hasegawa (1981) are given for comparison. The various energy contributions are tabulated in table 2, as well as the potential drop across the interface: For the profiles of (6), A V = 4 n~, ( d -~ -E -~) .

RC
Contributions to ywith non-oscillatory profile (dyn cm-') A more complicated form for the ionic profile, which permits oscillations, is For v = 0, (7) reduces to (6). There are five parameters (A, B , CY, v and 20) related by the three conditions: continuity of pion at zo, continuity of dpion/dz at zo, and electroneutrality : Using the first two to eliminate A and B , we may rewrite the third equation as follows: Thus, for any choice of CY and v one determines zo, and then A and B from the other conditions. For certain values of cyand v, equation (9) has more than one solution, giving rise to a variety of profiles. Some of these are shown in figure 1 for CY = v = 1; they correspond to zo = -0.13409, -0.78538, and -0.83602, the last being indistinguishable from the second on this plot. A profile with zo large and negative is suggestive of the situation for a solid surface, in which the ions are arranged in planes and the electron density is fairly constant at the bulk value until the last ionic plane, outside of which it decreases gradually. To the extent that this structure is maintained in the liquid, one would have an oscillatory ion density profile and an electron density profile centred several a. to the right of the last ion peak (Allen and Rice 1978). One could expect that in order to get a good surface tension from a variational calculation, a profile of this kind would be required. Nevertheless, our experiments with function (7) always showed that the lowest v was obtained when the root zo closest to 0 was chosen.
Accordingly, we considered the related trial function The continuity at zo is assured by A = 1 -B , the continuity of the slope at zo implies A a = BP and hence A = p/( a + p), and the electroneutrality condition may be written a p 2 -2 -2 ' " = ( p + a ) ( 2 + 2 ) .
Unlike the corresponding condition (9) for the function ( 7 ) , this admits only a single solution for 20, so that the present function is easier to work with systematically. Although it does not allow the large displacement of ionic from electronic profiles allowed by (7) it does include oscillations. On the other hand, the possibility of p unequal to &gives it additional flexibility. Note that the new function also includes pori of (6) as a special case when Y equals zero. There are now three free parameters in the ionic density profile, and one in the electronic density profile. The minimisation of ywith respect to the four parameters was performed using a gradient search method. The results for sodium and caesium are given in table 2. It can be seen that, although there is noticeable improvement (lowering) of the surface tension, it is by no means substantial enough to give satisfactory agreement with experiment. The profiles obtained for the ions are somewhat oscillatory, and appear quite different from those of the simpler exponential form, since the values of v obtained by minimising the surface tension (see table 3) differ from zero. Nevertheless, the total potential drop across the interface AV is not much changed by the introduction of the additional variational parameters (see tables 2 and 3). We take this as evidence that we have come close to finding the profile that minimises the surface tension, i.e. that the introduction of additional variational parameters would not produce much further change. (These parameters might allow oscillations in the electronic profile, which could decrease yes by allowing the electrons to follow the ions, at the expense of an increase in the electronic gradient contribution, but yes is a small contribution to y.) If our calculated values for yare close to the true minima of the density functional, we must conclude that the large difference between calculated and experimental surface tensions is due to problems with the functional itself, or the physical model on which the functional is based. Similarly, the narrowness of the profiles (compared with typical interionic spacings) and the large magnitude of AV would be due to flaws in the model. The aspect of the model which seems the most questionable is the use of a pseudopotential in first order for the ion-electron interaction. This suggests testing the functional for a model system without a pseudopotential, and work along these lines is in progress.
An apparent problem with the profiles found has to do with the lack of large oscillations. In Coulombic systems with no pseudopotential, one finds (Badiali et a1 1983) highly oscillatory profiles, except for very small bulk densities. We consider in the next section the reasons for the absence of strong oscillations in our calculations.

Importance of oscillations
It should be noted that the present trial function is not capable of giving very large oscillations. For fion to exceed 2ne, the value of A must be greater than unity. Since A = @/(a + p) with a and / 3 positive, this is impossible. One can show however that smaller oscillations are easily obtained and their absence certainly implies that the minimisation of the functional employed does not require substantial oscillations.
One explanation of the lack of oscillations is the presence of the pseudopotential. The Ashcroft (empty core) pseudopotential corresponds to the removal of the Coulombic attraction for small distances ( r < R,), and large oscillations do not occur in the profiles for non-Coulombic systems. To test the hypothesis that the pseudopotential removes the oscillations, we performed variational calculations for decreasing values of R,. These calculations also serve to show the effect of the choice of R,, which after all is a semi-empirical parameter of the model, on the calculated surface tension. The results for several choices of R, are given in table 4. The four parameters a, p, y and E were determined variationally for each value of R,. In table 4, we also give the maximum value of the ion density functions, which indicates the importance of oscillations around unity. It is evident that a decrease in R, leads tc an increase in the strength of the oscillations, It is also evident that an unreasonably large increase in R, would be needed to produce a surface tension in accord with experiment.
The contribution yps is large and negative (see table 2). It becomes less negative for the oscillatory profile. However, the resulting increase in y is more than compensated by the decrease in yll, producing a smaller surface tension. It appears that ylI is the term in the functional responsible for the existence of the oscillations. It may be noted that yll (gradients) and yes (electrostatic) change only slightly. A second reason for the lack of large oscillations seems to be the electron density profile. This is shown by the results of calculations (see table 4) in which we fix &instead of changing it variationally. Larger values of E produce larger oscillations in the ion density profile although the effect of increasing E on the oscillations is smaller than the effect of increasing R,. Of course, using &larger than the variational value increases the surface tensions, due to the increasing electronic contributions. An oscillatory electronic profile could also favour oscillations in the ionic profile, while (if electronic and ionic oscillations were in phase) decreasing the magnitude of AV. Note, though, that the Monte Carlo calculations give oscillations in the ion profile for non-oscillatory electron profiles (Badiali etal 1983).

Conclusions
Minimisation of the surface free energy, as given by the Evans-Hasegawa functional, yields surface tensions too high by a factor of 2. Experiments with more flexible trial functions for the ion profiles lead to relatively little improvement, in agreement with these authors' conjecture. Ions and electrons being treated similarly, one can expect that more flexible trial functions for the electrons would likewise not greatly improve calculated surface properties, as appeared from previous calculations (Badiali etal 1981) on surface energies. One could expect some improvement in the electrostatic terms if the electrons were allowed additional flexibility, but these terms are not very important. The problem seems thus to be with the truncated gradient expansion for kinetic and correlation energies or the first-order treatment of the pseudopotential for the ionelectron interactions. For the former, it has been suggested (Evans and Hasegawa 1981) that the form of the present functional is not appropriate when the profile is highly oscillatory. Other functionals, representing a partial summation of terms in the gradient expansion, are possible, but considerably more difficult to use. Assuming the same first-order pseudopotential and other interactions as in the present work, we have performed calculations (Badiali et a1 1981, Amokrane et a1 1981 for the alkali metals in which an attempt was made to model the interionic correlations, leaving out the ion-electron correlations. Electronic profiles were determined assuming steps for the ionic profiles (this theory does not provide a route to the determination of ionic profiles). Surface energies were in qualitative accordance with experiment for the alkalis, but surface tensions were much too low (sometimes negative, although absolute errors were comparable with those of the present work), suggesting that ion-electron correlations in the interface (i.e. a better treatment of the pseudopotential) needed to be considered. The results of table 4 indicate that, had those calculations used less abrupt ion profiles, higher (and positive) surface tensions would have been obtained. The main difference between the present and former approaches is in the treatment of interionic correlations. The density functional calculations and the step-function calculations both contain the errors associated with a first-order pseudopotential theory, so differences between their results certainly reflect this. However, calculations which neglected interionic correlations (Badiali eta1 1981) gave about the same surface tensions as calculations (Amokrane er a1 1981) incorporating them in an approximate way. Then the vastly different surface tensions of the density functional model suggest that interionic correlations are not well treated. In the density functional calculations these are found in yll and yl,. It may be noted that incorporates the response function or direct correlation function of the homogeneous systems of various densities, including unphysical ones. Further, the homogeneous systems are necessarily electrically neutral at every point, so a treatment based on them would seem suspect for regions of space where ions and electron densities are very different. This may also be why highly oscillatory profiles were not allowed by the present functional.
We conclude that the high surface tensions result, at least in part, from the way in which the density functional treats interionic correlations. Probably the first-order treatment of the pseudopotential is also a significant source of error. Going to higher order would require a fundamental change in the formalism, which treats electrons and ions in an equivalent manner. A better description of interionic correlations could be incorporated into the present formalism through a different choice of density functional ( y,l and y,,). Density functional calculations for a system with no pseudopotential should help in understanding the problems with the choice used here.