Adsorption of hard spheres: structure and effective density according to the potential distribution theorem

We propose a new type of effective densities via the potential distribution theorem. These densities are for the sake of enabling the mapping of the free energy of a uniform fluid onto that of a nonuniform fluid. The potential distribution theorem gives the work required to insert a test particle into the bath molecules under the action of the external (wall) potential. This insertion work W_ins can be obtained from Monte Carlo (MC) simulation (e.g. from Widom's test particle technique) or from an analytical theory. The pseudo-densities are constructed thusly so that when their values are substituted into a uniform-fluid equation of state (e.g. the Carnahan-Starling equation for the hard-sphere chemical potentials), the MC nonuniform insertion work is reproduced. We characterize the pseudo-density behavior for the hard spheres/hard wall system at moderate to high densities (from \rho^*= 0.5745 to 0.9135). We adopt the MC data of Groot et al. for this purpose. The pseudo-densities show oscillatory behavior out of phase (opposite) to that of the singlet densities. We also construct a new closure-based density functional theory (the star-function based density functional theory) that can give accurate description of the MC density profiles and insertion works. A viable theory is established for several cases in hard sphere adsorption.


Introduction
Inhomogeneous fluids pose, due to their vast varieties, challenges to the liquid-state theories, and are of capital importance in many industrial processes. The nonuniformity in the fluid can be produced by a one-body potential: such as the surface potential from the interfaces separating two or three coexisting phases (solid/liquid/vapor phase transitions), electrical fields acting on ionic species (the electric double layer), or magnetic fields on ferromagnetic liquids. The conventional classical density functional theory (DFT) employed to deal with such systems in the last half century has been based on the Hohenberg-Kohn [1] formulation of a grand potential, Ω, that is expressed in terms of an intrinsic Helmholtz free energy functional F [ρ (1) w ] (IHFEF): (1) w ] + d rρ (1) where ρ (1) w is the singlet (nonuniform) density, w(r) the external (1-body) potential, and µ 0 the chemical potential of the bulk fluid. At equilibrium, the grand potential is minimized with respect to the singlet ρ (1) w , and it yields the Euler-Lagrange equation δβΩ δρ (1) δβF [ρ (1) w ] δρ where β = 1/(kT ): k is the Boltzmann constant, T is the absolute temperature. The IHFEF can be separated into two parts: the ideal part F id (which is known: F id [ρ] = d rρ (1) w ( r)[ln ρ and the excess part F ex . Thus F [ρ (1) w ] = F id + F ex . Recently, we have proposed a star-function based density functional theory [2,3] (s-DFT) which expresses F ex succinctly as − βF ex [ρ] = −βF ex 0 + d rC (1) 0 ( r)δρ w ( r) + 1 2 d rd r C (2) 0 ( r, r )δρ w ( r)δρ w ( r ) + S * w (1. 3) where S * w is the star function [4] defined as the "primitive" of the bridge function B w S * w = d r δρ w ( r) γ w ( r) 0 ( r, r ) is the pair direct correlation function (2 − DCF ). F ex 0 is the bulk fluid free energy. Subscripts 0 and w indicate uniform and nonuniform fluid properties, respectively. δρ w = (ρ (1) w − ρ 0 ) is the difference between the nonuniform density ρ (1) w and the uniform density ρ 0 and γ w is the indirect correlation. We note that if the bridge function B w can be obtained exactly, the star function S * w will consequently be exact from equation (1.4), thus the excess free energy F ex in equation (1.3) will also be exact. The Euler-Lagrange (EL) equation (1.2) can be expressed as where C (1) w is the nonuniform singlet direct correlation function (1 − DCF ). Alternatively, it is possible to show [2] via functional expansions that equation (1.5) can be written as (1.6) Thus, we have the equality [3] Equations (1.3) and (1.6) constitute the basis of the s-DFT theory.
In this article, we shall explore two essential aspects of nonuniform fluids: (1) the effective density that enables the mapping between the uniform-fluid free energies and the nonuniform-fluid free energies; (2) the closure relation between the bridge function and other correlation functions. The first task is mediated through the potential distribution theorems [5,6]; the second is based on a successful uniform liquid theory [7] developed in recent years. Historically, there have been interchanges and cross-fertilization of uniform and nonuniform liquid theories [8,9]. In this paper we test the new formulations on hard spheres adsorbed on hard planar walls.
In weighted-density density functional approaches (WDA) [10,11] the "mapping" of the free energy is done through some weighted density: ρ. The weighted density is obtained from the true density ρ (1) w via a convolution integral with a weighting function, ω, as the kernel: i.e. ρ( r) = d r ρ (1) w ( r )ω( r, r ). We shall, however, develop an effective density without going through this weighting procedure. This is made possible by the potential distribution theorem (PDT), since PDT naturally lends itself to yield the work of insertion (the free energy of an inhomogeneous system).
The PDT is a fundamental theory of statistical mechanics [5,6,[11][12][13]. It is concerned with the work required to insert test particles in an equilibrium ensemble representing a uniform fluid or a nonuniform fluid. In an earlier paper [14] we have generalized the PDT from the chemical potential, to the cavity functions, and to higher-order correlation functions. The PDT has previously been extended to nonuniform systems [5,6]. However, here we distinguish two types of test particles. The first type is for a test particle subject to the wall force as well as to the forces of the bath molecules. This is the commonly studied type. The second type is for a test particle free of the wall force. Its interaction with the wall is intentionally removed (similar to the cavity function 33601-2 in homogeneous fluids). The PDT of the second kind is directly related to the Euler-Lagrange equation [3]. The PDT establishes a connection among (i) the insertion work W ins , (ii) the bulk fluid chemical potential µ 0 , and (iii) the singlet distribution function ρ (1) w . To illustrate the PDT, we carry out calculations on the system of hard spheres on a hard wall (HS/HW), for which MC data are available (e.g. from Groot et al. [15]).
The mapping between the uniform and the nonuniform systems can be constructed as follows: given an equation of state for the uniform fluid (which is utilized to calculate the uniform chemical potential βµ 0 ), we solve for the hypothetical densities ρ pseudo (z) that will reproduce the insertion works W ins (z) obtained from the MC, i.e. βµ 0 (ρ pseudo (z)) = βW ins (z). This density profile is an artificial construct and is called the pseudo-density ρ pseudo (z). Its sole role is to reproduce the nonuniform free energy via the bulk equation of state. For instance, in the HS/HW system, the bulk fluid equation is taken to be the Carnahan-Starling (CS) [16] equation which is known to be highly accurate. It yields the uniform hard-sphere chemical potential βµ 0 . If we know the insertion work βW ins (z) (from MC) in HS/HW, a series of values of ρ pseudo (z) can be generated so that the equality βµ HS 0 = βW ins (z) holds. The ρ pseudo (z) obtained thusly will enable the mapping of the free energies ( βµ HS 0 ) of the uniform fluid to those (βW ins ) of the nonuniform fluid. In this respect, the pseudo-densities perform precisely the same function as the weighted densities ρ(z) in the WDA. The difference is that our approach is based on the PDT, an exact theory, not on approximate theories (such as the Percus-Yevick equation [17]).
Should the bridge function B w be known, one can obtain the fluid structure at the wall via the Euler-Lagrange equation (1.6) [2,3]. B w as a correlation function is well defined in terms of infinite series of cluster diagrams [18] or functional expansions [19]. However, direct evaluation from its definition is numerically intractable. But it can be inverted from machine data (as in reverse engineering [20]) if we consider equation (1.6) as the defining equation [19]. The inverted B w can be given as a set of numerical data [3], say, in the form of tabulated entries; or it can be expressed as a function B w = f (r) for some analytical function f (e.g. a polynomial). Another approach, and a more appealing one is to formulate a liquid theory (à la mode of the Percus-Yevick equation [17]) that relates B w to other correlation functions B w = f (γ) , γ being a well-defined quantity in liquid state theory. This in uniform liquid theory is called the closure relation. This relation may or may not exist, because B w is in fact a functional [3,21] B w = Φ[ρ (1) w , w] of the singlet density ρ (1) w and/or the external potential w (Φ[.] being the functional sign). Thus, it may not be simply related to any one correlation function. For uniform liquids, the determination of the function B = f (γ), γ being some correlation, is already a major task (examples such as the Rogers-Young [22] (RY) closure or the Martynov-Sarkisov [23] (MS) closure). We have proposed a theory upon resummation of the functional expansions of B w in [3]. This resummation seems to work well for a number of pair potentials. For nonuniform fluids, the PY and hypernetted-chain [24] (HNC) closures have been tried and shown not to work well in the past [25]. And this prompted the proposition that closurebased approaches might not be suitable for nonuniform systems. In this work, we shall propose a zero-separation type (ZSEP) closure derived from the study of uniform liquid theory [7]. This closure will be tested here on selected cases of the HS/HW systems. Once the singlet density is obtained from equation (1.6) given B w , one can use the Euler-Lagrange equation (1.5) to obtain the 1-DCF C (1) w . The 1-DCF is simply related to the insertion work as Thus, we can access the insertion work through (1.5), using the MC data [15] on ρ (1) w as input. The thusly obtained βW ins is considered as the MC-derived insertion work. In addition, we shall construct a theoretical method for obtaining the insertion work.
Section 2 gives the potential distribution theorems for nonuniform fluids with general interaction forces. Section 3 shows the calculations of the effective densities based on the MC-derived insertion works for the hard spheres/hard wall system. The mapping of the nonuniform free energies is mediated through an accurate equation of state by Carnahan and Starling. We characterize the behavior of ρ pseudo and compare with the commonly used weighted densities found in literature. In section 4 we test the new closure theory (ZSEP) on the same systems. We compare the theoretical structures and free energies with the MC data. In section 5 we draw the conclusions.

Potential distribution theorems for nonuniform systems
In this section, we shall introduce the potential distribution theorems for two types of test particles. One type of test particle is the well-known one, i.e. it is subject to interactions with the bath molecules as well as to forces of the wall. The second type is the one that interacts only with bath molecules and not with the wall. The latter result will be the PDT that we shall use in this article. We start with the averages in a canonical ensemble. (This can also be done in the grand canonical ensemble [5,6]). The N-body system consists of N fluid molecules and they border on one side at a solid surface. The Hamiltonian H N is written as where r N ≡ ( r 1 , r 2 , . . . , r N ) is a shorthand for the N -vector of the positions of N particles, and p N is the N -momenta vector, K N is the kinetic energy, while is the total potential energy between fluid particles (assumed to be pairwise additive), and is the sum of one-body energies arising from the external (wall ) potential w. This is the potential energy that drives the inhomogeneities in the system. Next, we consider an (N-1)-body system: the sum of pair energies becomes and the sum of one-body energies A test particle is introduced to the (N-1)-body system as particle with label "N" located at the distance r N . This test particle can interact with the other N − 1 bath molecules through the pair potential u (2) (r i , r N ), i < N , as well as with the wall through the one-body potential w( r N ). This particle is designated the type-1 test particle. The excess potential energy Ψ w ( r N ) due to the presence of type-1 test particle is On the other hand, the type-2 test particle, by choice, does not interact with the wall (w( r N ) = 0). The excess potential energy Ψ −w ( r N ) is (subscript "-w" stands for "without the wall interaction") (2.7)

33601-4
The potential distribution theorem for the first-type test particle (interacting with w) is obtained from the (N-1)-body canonical ensemble average of Ψ w ( r N ) Furthermore, from the definition of the singlet density ρ (2.10) Equation (2.8) has been obtained earlier [5,6]. Note that Q N is the N-body configurational integral; Λ is the de Broglie wavelength; and The potential distribution theorem for the type-2 test particle (without w) is obtained instead from the (N-1)-body ensemble average of the wall-less excess potential The insertion work W ins required to insert a type-2 test particle with the excess potential 14) The quantity W ins can be directly simulated via the MC method using the Widom particleinsertion technique [5,26], or from the EL equation (1.5). W ins can also be considered as the intrinsic work of insertion.

Effective density and work of insertion
For the hard sphere system, we adopt the MC data by Groot et al. [15] as the basis of calculation. The number densities chosen for the hard spheres are ρ * = 0.5745, 0.715, 0.758, 0.813, and 0.9135. The hard spheres are adsorbed on a planar hard wall stretching over the x-y plane. Thus, the fluid-fluid pair potential u (2) is u (2) (r) = ∞, r σ, u (2) (r) = 0, r > σ.
w (z) from MC data of Groot et al. [15]. The brown line is the insertion work βWins produced by the FMT theory. Z ≡ distance from the wall. Unit of length is taken to be the hard-sphere diameter σ. The horizontal line is at ρ * = 0.5745 (a guide for the eye). The pseudo-density ρ pseudo (z) oscillates weakly around the bulk value and in opposite periodicity as compared to ρ (1) w (z). The insertion work βWins also oscillates synchronously with ρ pseudo (z) but oppositely to ρ The wall potential w on the fluid particles (in the z-direction perpendicular to the wall) is Thus, the closest approach of the center of a hard sphere to the wall is at z = σ/2 (we shall use the hard sphere diameter σ as unit of length below). For ρ * = 0.5745, the pure hard-sphere fluid excess chemical potential βµ 0 is 4.897, as calculated from the CS equation. The singlet density profile ρ (1) w (z) has been obtained by Groot et al. [15] via MC simulation and is displayed in figure 1 (see the diamond symbols ♦). The insertion work βW ins (z) is identified as −C The insertion work is plotted as triangles ( ). We remark that both the singlet density ρ (1) w (z) and the insertion work βW ins (z) are "real" quantities, i.e. they are measurable quantities for the system of interest. We next formulate the pseudo-densities as mentioned earlier in section 1. We ask at what fictitious densities ρ pseudo (z) for a uniform HS fluid (which follows the CS equation), will we obtain a chemical potential value βµ 0 that is equal to the MC nonuniform insertion work βW ins (z) where η = π 6 ρ pseudo (z)σ 3 . In figure 2, we show this mapping via a diagram. Note that Z HS is the non-ideal compressibility factor from the CS equation. The pseudo-density ρ pseudo (z) thus calculated is also plotted (as asterisks * ) in figure 1. We observe that the pseudo-density stays fairly flat and close to the bulk value of 0.5745, oscillating only weakly up and down with respect to the horizontal line. The oscillations are out of phase with respect (opposite in period) to the singlet density ρ  w (z) has a contact value ≈ 2.296 (z = σ/2). It oscillates vigorously with pronounced peaks and valleys. We note that the insertion work βW ins (z) also oscillates with trends similar to the pseudo-density: when ρ pseudo (z) is up, βW ins (z) is also up, and when ρ pseudo (z) is down, βW ins (z) is also down.
A question naturally arises: how does this pseudo-density ρ pseudo (z) compare with the commonly used weighted densities, say, n α (z) in the fundamental measure theory (FMT) of Rosenfeld [27]. We   selected two weighted densities n 0 (z) (green line) and n 3 (z)/V (blue line) from FMT for comparison in figure 3 (there are in fact six weighted densities n α (z) in FMT (α = 0, 1, 2, 3, V 1, V 2) based on the geometries of hard spheres). The basis set of weights are composed of three functions: ω 2 is based on the Dirac Delta function, thus generating n 0 , ω 3 is based on the Heaviside function, generating n 3 , and ω V 2 is a vector function (see Rosenfeld [27]). The other linearly-dependent weighting functions are related to this basis set by scaling through spherical geometries. The first three weighted densities are related by the relations: n 0 (z) = n 1 (z)/R = n 2 (z)/S (R and S are the radius and surface area, respectively, of a hard sphere), while n 3 (z) is normalized by the spherical volume V : n 3 (z)/V .
In figure 3 we notice that the magnitudes of both n 0 (z) and n 3 (z)/V are close to the pseudodensity ρ pseudo (z) and far from the singlet density ρ (1) w (z). There are visible differences between ρ pseudo (z) and n 0 (z) as well as n 3 (z)/V . This is not surprising, since these coarse-grained densities are derived from different theories (the FMT vs. the PDT). ρ pseudo (z) is close to n 0 (z) only near the wall and at long ranges. The cusp in n 0 (z) (due to the delta weighting) is missing from ρ pseudo (z).
Values of n 3 (z)/V are quite distinct from both n 0 (z) and ρ pseudo (z). The observations also apply to other two weighted densities owing to scaling: n 1 (z)/R and n 2 (z)/S. On the whole, the weighted densities and the pseudo-density are commensurate with each other (of similar magnitudes) even though they came from entirely different origins.
To make further comparisons, we also calculate the insertion work arising from the FMT. The excess free energy F ex in FMT according to literature [27]) is expressed as: βF ex = d rΦ(n α ( r)) where Φ is a free energy density. This quantity is a function of the weighted densities n α and itself arises from a uniform fluid model (such as the compressibility equation from the Percus-Yevick equation [17] (PY c ) or variations thereof [28]). The FMT insertion work can be evaluated as the functional derivative [28] δF ex /δρ (1) w (z), thus in Fourier space This FMT insertion work (W ins (FMT)) at ρ * = 0.5745 is plotted in figure 1 (the brown line). It quite closely follows the simulated W ins (MC) (triangles). This is to be expected, since both the FMT and the PDT-based theories are designed to match the nonuniform free energy. The FMT free energy lies somewhat above the PDT values, because the PY c solution for the chemical potential is higher than the MC value for hard spheres. w (z) from MC data of Groot et al. [15]. The brown line is the insertion work βWins produced by the FMT theory. The horizontal line is at ρ * = 0.715. The pseudo-density ρ pseudo (z) oscillates weakly around the bulk value and in opposite periodicity as compared to ρ (1) w (z). The insertion work βWins also oscillates synchronously with ρ pseudo (z) but oppositely to ρ (1) w (z).
Similar calculations are made for HS/HW at a higher density ρ * = 0.715. Figure 4 shows the singlet density ρ In figure 5 ρ pseudo (z) is compared with n 0 (z) (green line) and n 3 (z)/V (blue line). Again, we see that at a few points, n 0 (z) and ρ pseudo (z) match closely, while n 3 (z)/V has larger differences from ρ pseudo (z) at this higher density. The comparisons show that the PDT approach is clearly different from the fundamental measure theory. FMT was constructed based on the PY c theory (or modifications thereof) in order to map the compressibility free energy (see also later improved versions [29]) to the nonuniform free energy. In our approach, we use the PDT to map the CS free energy to the nonuniform free energy, and we use MC data to obtain the insertion work. The objectives are the same; but the theoretical bases are different.  To further examine high-density hard spheres, we next examine the case ρ * = 0.813. The results are plotted in figure 6. Similar oscillatory behavior is observed for the pseudo-density and the insertion work (reduced to 20% in order to fit in the graph). The phases of oscillations of the singlet ρ (1) w (z) and ρ pseudo (z) (and βW ins as well) differ by a half period π. Figure 6. The effective (pseudo) density ρ pseudo (symbol * ; from PDT) and the insertion work βWins (symbol ∆; from MC)) at bulk hard sphere density ρ * = 0.813 for the HS/HW system. The insertion work is scaled down to 20% in order to fit in the graph. Symbol (♦) ≡ the singlet density profile ρ  w (z). The insertion work βWins also oscillates synchronously with ρ pseudo (z) but oppositely to ρ

Construction of a closure theory
It has been recognized that the DFT for the structures of nonuniform fluids should first be based on valid uniform-fluid properties as inputs. In particular, the structure of nonuniform fluids depends critically on the accuracy of the pair direct correlation functions 2-DCF of the uniform fluid that must be supplied beforehand. We have developed in the last decades an accurate theory [7] for the uniform hard sphere fluid: the zero-separation-theorem-based closure (ZSEP). We can thus easily generate accurate 2-DCF C (2) 0HS for use in the present study of nonuniform hard-sphere fluids.

33601-9
We shall look at the five densities where the MC data are available for the HS/HW systems: namely ρ * = 0.5745, 0.715, 0.758, 0.813, and 0.9135. The ZSEP closure for the uniform bridge function B 0HS is of the form [7] (caret indicates a function) where α, φ, and ζ are adjustable parameters, and γ * = γ + ρ f 2 . f is the Mayer factor of the repulsive WCA (Weeks-Chandler-Andersen) potential. For details, see references in [7]. The adjustable parameters α, φ, and ζ are determined using structural consistencies (the zero-separation theorem and the contact-value theorem) and thermodynamic consistencies (e.g., the pressure consistency, and Gibbs-Duhem relation). The uniform 2-DCF's C (2) 0HS obtained are depicted in figure 7. They have all satisfied closely the consistency relations imposed. For example, in the case of the high density ρ * = 0.90, the C 0HS is compared with the MC-generated data [30] in figure 8. Diamond 0HS (r) of uniform hard spheres generated by the ZSEP closure [7] for the bulk density ρ * = 0.90. The lines are from ZSEP theory, and the symbols (♦) from MC data of Groot [30]. The inset shows the magnified view of the non-zero tail and discontinuity in C   symbols denote the MC data, and the solid line is the 2-DCF from the ZSEP(4.1). We see very close agreement of the two curves. The zero-value C  while MC gives 1.05. This is in obvious contrast to the PY-based theories that would produce C (2) 0HS (r) = 0 for r > 1 (an erroneous result). The same observation applies to all five present cases ( figure 7).
For the nonuniform HS/HW case, we developed a closure similar to ZSEP [3]: where α , φ , and ζ are adjustable parameters for the nonuniform B w . sgn(.) is the sign function. γ w is the indirect correlation function defined as We found from practice that "renormalization" of the indirect correlation is not needed for this nonuniform fluid. We also note that the expression (4.2) is modeled after the ZSEP equation (4.1) for uniform fluids [7]. In earlier studies [2,3] we have shown that this closure was effective for adsorption of high-density Lennard-Jones fluids. The similarity between the uniform closure and the nonuniform closure used represents an ongoing effort in the transfer of successful theories for uniform liquids to nonuniform liquids.  Equations (1.6), (4.2), and (4.3) are coupled. The convolution was solved by standard numerical methods [3] in bipolar coordinates. Trapezoidal rule was used in integration. Picard's iterations with mixing (relaxation) on input and output iterates were employed. The grid of numerical integration was ∆r = 0.01σ, and total grid number N = 2048. Thus, the maximum z-distance reached is 20.48σ. Cauchy's absolute convergence was enforced for the γ w -function with convergence criterion δ = 0.0001. There were three free parameters to be determined in equation (4.2): α , φ , and ζ . We needed three conditions. The first condition was the hard-wall sum rule: the contact density ρ (1) w (σ/2) = P/(kT ). In practice, we set φ to a default value of 0.4. α and ζ were adjusted in tandem until the wall rule was satisfied (see table 1). In the future, other sum rules [31] will be tested and incorporated. Figure 9 shows the result of application of ZSEP. The density varies from the moderate ρ * = 0.5745 to the highest ρ * = 0.9135. The symbols are the MC data [15], the lines are from the ZSEP theory. We observe overall superb agreement for densities ρ * from 0.5745, 0.715, 0.758 to 0.813. At ρ * = 0.813, our results compare well with previous accurate theoretical attempts (e.g. the White-Bear version of FMT [29]). The ingredients for success are from (i) the high accuracy of the input uniform 2-DCF's; and (ii) the effectiveness of the ZSEP closure (4.2). For the highest density ρ * = 0.9135 there is some deterioration of the theoretical curve at the second peak. The heights of the peaks (of ZSEP and MC) are similar, however the predicted location is displaced to a longer distance. A closer reading of the original MC source [15] reveals that at this high density (packing fraction η = 0.4783), there are clusters of molecules forming hexagonal structures for the 33601-11   [15]. Lines ≡ ZSEP results. We discern close agreement between the ZSEP curves and the MC curves except at the highest density (0.9135) where the second peak of the theoretical curve is displaced to a larger z-value. (See discussions in the text). first layer of particles (figure 6 of [15]), and reduced normal-direction diffusion coefficient for the surface molecules: all indications of two-dimensional quasi-solid-fluid transition. Note that the fluid limit of homogenous hard spheres is at about η ≈ 0.4948. For inhomogeneous hard spheres this may occur earlier. This puts stringent demand on a closure theory which is not intended to predict solid formation.
Next, we calculate the insertion work via equation (2.14) using ZSEP-density ρ

Conclusions
The two major questions that we set out to explore in this study are (i) how does the effective density based on the potential distribution theorem look like and how does it behave? (ii) Can we obtain the structure of hard spheres adsorbed on the hard wall through a closure theory for nonuniform fluids? We formulate two potential distribution theorems: one for the insertion work of the usual test particle that interacts with the wall, and the other for the test particle that ignores the wall. It is the latter insertion work that corresponds to the singlet direct correlation and as a consequence the second potential distribution theorem becomes equivalent to the Euler-Lagrange equation.
The pseudo-densities ρ pseudo , due to their construction, are much smoother quantities than the singlet densities. They oscillate around the bulk value with periods out-of phase with respect to the nonuniform ρ (1) w . The insertion work is derived from the MC data on HS/HW (Groot et al., 1987 [15]) via the Euler-Lagrange equation (with the identification βW ins = −C (1) w ). Its behavior for several densities (ρ * = 0.5745, 0.715, and 0.813) is determined. βW ins oscillates in phase as the pseudo-density, but out of phase with respect to the singlet density ρ (1) w . The density functional theories need accurate uniform pair direct correlation functions C 0HS as input. Earlier we have formulated an accurate closure theory for the 2-DCF (the ZSEP theory) of uniform hard spheres. It is computationally cheap now to generate these functions: it took less than a minute of CPU time on a PC computer to produce an accurate DCF. Here we propose a new closure theory for the nonuniform hard spheres equation (4.2). Its form is inspired by the uniform fluid ZSEP equation. This bridge function is used in the star-function based density functional theory equation (1.6) to generate the singlet density profiles. These ZSEP densities are compared with the MC data (figure 9) at five densities ρ * = 0.5745, 0.715, 0, 758, 0.813, and 0.9135. Except for the highest density, the ZSEP gives accurate nonuniform density profiles. This shows that the closure-based density functional theory can perform reasonably well for hard spheres when a suitable closure (bridge function) has been adopted. We shall have further developments on this in the future.
These theoretical predictions for the singlet densities are used to calculate the insertion works (figure 10). They also compare favorably with the MC data.
In summary, we have proposed and evaluated a new type of effective densities via the potential distribution theorem. We characterize their behavior for the hard spheres/hard wall system at moderate to high densities (up to ρ * = 0.9135). The free energies (the insertion works W ins ) of the nonuniform system are also calculated. They show oscillations opposite in phase to those of the singlet densities. We also construct a new closure-based density functional theory that can give accurate reproduction of the computer simulated densities and insertion works.
The newly developed star-based density functional theory can be extended to fluids interacting with soft/attractive potentials. This task is under active investigation.