Hard spheres at a planar hard wall: Simulations and density functional theory

Hard spheres are a central and important model reference system for both homogeneous and inhomogeneous fluid systems. In this paper we present new high-precision molecular-dynamics computer simulations for a hard sphere fluid at a planar hard wall. For this system we present benchmark data for the density profile $\rho(z)$ at various bulk densities, the wall surface free energy $\gamma$, the excess adsorption $\Gamma$, and the excess volume $v_{ex}$, which is closely related to $\Gamma$. We compare all benchmark quantities with predictions from state-of-the-art classical density functional theory calculations within the framework of fundamental measure theory. While we find overall good agreement between computer simulations and theory, significant deviations appear at sufficiently high bulk densities.


Introduction
Hard spheres are a central model in statistical physics and act as a reference system for simple and complex fluids. The bulk phase behavior of the hard-sphere system is simple because the system is athermal -there is no energy scale to be compared to the thermal energy of k B T , where k B is the Boltzmann constant and T is the absolute temperature. Therefore, the phase behavior is controlled solely by the particle density ρ = N /V -or equivalently by the packing fraction η = ρ * π/6, where the reduced density is defined as ρ * = ρσ 3 , with σ being the hard-sphere diameter. For packing fractions η < 0.492, the equilibrium thermodynamic phase is fluid, while above this value of η, hard spheres are capable of forming a crystal.
In general, the structure of simple liquids is dominated by packing effects generated by the steep and short-ranged repulsive part of the interatomic interaction, which can be very well approximated by a hard-sphere interaction. The usefulness of the hard-sphere system as a reference is not limited to bulk systems, but also extends to inhomogeneous fluids. For example, the system studied here, i.e., a hard-sphere fluid confined by planar hard walls, is one of the simplest model systems for a fluid-solid interface. The wall is modelled by the external potential V ext (r) = V ext (z) = ∞, z < σ/2, 0, otherwise, (1.1) where z is the distance normal to the wall. Close to the wall, the fluid develops an inhomogeneous structure, which can be described through the ensemble averaged density profile ρ(r). At sufficiently low bulk densities, when spontaneous symmetry breaking due to freezing can be ruled out, the equilibrium density profile possesses the same spatial symmetry as the external potential, so that we can assume that the density profiles depend only on z, that is, ρ(r) = ρ(z). From the density profile as a function of packing fraction, one can determine a number of thermodynamic interfacial properties, including the excess adsorption Γ, the excess volume v ex = ρΓ and the wall surface free energy γ. Recently, high precision results for the excess adsorption Γ and the wall surface free energy γ at various bulk densities have been obtained using molecular-dynamics simulation [1,2]. Due to their high precision, these data are well suited as benchmarks to enable the testing of theoretical predictions. In this paper we present the results for the density profile ρ(z), from these high precision simulations, together with comparisons with state-of-the-art classical density functional theories (DFT) for hard-sphere fluidsnamely DFT formulations based on fundamental measure theory (FMT) [3,4].

Simulation
Density profiles of the hard-sphere fluid at a planar hard wall were determined in the moleculardynamics (MD) simulation using the algorithm of Rapaport [5]. The walls were placed normal to the zaxis at a distance of about 65σ apart. The x-y cross section was approximately square with a side length of about 50σ. Periodic boundary conditions were employed in the x and y directions. To measure the density profiles, ρ(z), the system was divided along the z-axis into bins of width 0.02σ. The size of the simulation box was the same in all simulations while the number of spheres varied from about 8 000 for the lowest bulk reduced density of about 0.052 to about 150 000 for the highest reduced density of about 0.938. Systems at 17 different densities were simulated. At each density, we performed 50 independent runs starting from well equilibrated initial states. The data for the density profiles were averaged over the runs and over the two walls. The 95% confidence intervals were estimated from the scatter in the data for independent runs. We also simulated systems of smaller size at a range of densities and determined that the size of the systematic error due to the finite system size was much smaller than the statistical errors. In this manuscript, we include only a few examples of the density profiles. The complete tabulated density profiles can be found in the Supplementary Material [6]. MD results for the excess volume, v ex , were obtained from the density profiles using the procedure described in the Supplement to reference [7]. Near the freezing density, the hard-sphere fluid at a hard wall is known to exhibit a pre-freezing transition [8,9]; however, the transition is known to have a significant nucleation barrier due to line tension [10]. We have examined the structure of our hard-sphere fluid at and near the wall surface and find no evidence of crystalline ordering that would indicate that prefreezing has occured, even at the highest packing fractions; therefore, we are confident that we are examining a fully fluid-wall interface.

Density functional theory
In density functional theory (DFT), there exists a functional Ω[ρ] of the density distribution ρ(r) of the form [11] is the functional of the intrinsic Helmholtz free energy, V ext (r) is the external and µ the chemical potential. It can be shown that the functional Ω[ρ] is minimized by the equilibrium density profile ρ 0 (r) and that it reduces to Ω -the grand potential of the system in equilibrium, i.e., Ω = Ω[ρ 0 ] [11]. These properties can be employed in order to obtain the inhomogeneous structure of a fluid in an external potential within the same framework as thermodynamic quantities. From the variational principle of DFT we obtain the equilibrium density profile: From the density profile ρ 0 (r), in general, or ρ 0 (z) in the present study, one can directly compute the excess adsorption Γ via where ρ is the bulk density, and the wall surface free energy γ via Here, A is the area of the wall, which is assumed to be infinite and L = V /A is the extension of the system in z-direction. In order to make Γ and γ well-defined quantities, it is necessary to define the volume V in equations (2.3) and (2.4), i.e., one must define the dividing interface at which the system and the wall are separated [12]. Here, we use the location of the actual hard wall as a dividing interface. As a consequence, the excess adsorption Γ for a hard-sphere fluid turns out to be negative and the surface free energy γ positive. Because the wall surface free energy γ and the excess adsorption Γ are related by Gibbs's adsorption theorem This quantity, which is related to the interfacial adsorption by is useful because it provides a convenient route to the determination of γ directly from the density profile as was illustrated in reference [2].
In order to use DFT, one must specify the functional F [ρ] of the intrinsic Helmholtz free energy in equation (2.1). The intrinsic Helmholtz free energy can be split into with an exactly know ideal gas contribution and an excess (over the ideal gas) contribution F ex [ρ], which contains all the information about interparticle interaction. In equation (2.10) β = 1/(k B T ), and λ is the thermal wavelength.
For the system of interest, a hard sphere fluid, fundamental measure theory (FMT) [3,4] provides an accurate approach for the excess free energy functional. Within FMT, the excess free energy is written as follows: (2.11) where Φ, the excess free energy density, is a function of weighted densities n α (r) [3,4]. The details of FMT, including the definition of the weighted densities can be found in a recent review [4]. Here, we employ three different versions of FMT: (i) the original Rosenfeld functional [3], (ii) the White-Bear version of FMT [13,14], and (iii) the White-Bear version of FMT Mark II [15]. The main difference between these three versions of FMT are the equations of state underlying the functionals: the Percus-Yevick (PY) compressibility equation for the Rosenfeld functional, the Manssori-Carnahan-Starling-Leland (MCSL) [16] pressure for the White Bear version and a recently proposed, somewhat more consistent generalization [17], of the Carnahan-Starling (CS) pressure for the White-Bear Mark II functional.

Density profiles
We begin by presenting density profiles of a hard-sphere fluid in contact with a planar hard wall at selected values of the bulk density. The density profiles from the MD simulations act as benchmark data for a comparison with the results from different versions of FMT.
It is interesting to note that, for a fluid in contact with a planar hard wall, the density closest to the wall, the so-called contact density ρ c , is fixed by the contact theorem ρ c = ρ(z = σ + /2) = βp, (3.1) where p is the bulk pressure. Equation (3.1) is satisfied by the results of FMT [4].
It is well known that the CS pressure, which underlies the one-component White-Bear and White Bear Mark II versions of FMT, is more accurate, as measured by comparison to computer simulations, than the PY compressibility pressure, which underlies the original Rosenfeld functional. At low bulk density, the difference between the CS and PY equation of state is small, but at sufficiently high bulk fluid densities, the PY pressure significantly overestimates the pressure of a hard-sphere fluid relative to the computer simulation results. Hence, it is to be expected that very close to the wall, where ρ(z) is strongly influenced by the contact theorem [equation (3.1)], the density profiles obtained from the White-Bear versions of FMT will be closer to the simulation results than those calculated using the Rosenfeld functional.
In figures 1-4, we show density profiles at four different values of the bulk density. In figure 1, the bulk density ρ * ≈ 0.305, corresponding to a bulk packing fraction of η ≈ 0.159, is relatively low. All versions of FMT predict the density profiles that agree very well over the whole range of z with that obtained from computer simulations. At the wall, as expected, the White Bear versions of FMT are more accurate than the Rosenfeld functional, but the difference is very small, as can be estimated from the difference between the PY and the CS equations of state (about 0.3%). The packing effects close to the wall are small and decay fast.

23001-4
Hard spheres at a planar hard wall ( figure 4), respectively, the agreement between FMT and computer simulations is still good very close to the wall, but clearly less satisfying for the oscillatory structure -we observe that the FMT significantly underestimates the height of the density peaks about one particle diameter away from the wall. Especially at ρ * ≈ 0.938, a density close to bulk freezing, one can see that the height of the second density peak is also underestimated by FMT.
Since all versions of FMT seem to have a problem with the second peak in the density profile at high fluid density, it seems likely that it is not due to the precise form of the excess free energy density Φ employed, but rather the structure and range of the weight functions that can only approximate the complicated integrals of the virial expansion of the free energy [20][21][22].

Surface free energy γ
From the equilibrium density profile it is straightforward within DFT to compute the wall surface free energy from equation (2.4). Because various versions of FMT result in slightly different density profiles, as shown in figures 1-4, the wall surface free energy obtained from different versions of FMT are expected to differ by a small amount.
The MD simulation results were obtained, as in reference [2], using equation (2.8) and integrating the excess volume, v ex , with respect to pressure. In that work, the numerical error in the integration was minimized by subtracting from the integrand the expression for the excess volume from Scaled Particle Theory (SPT) and then adding back, after integration, the exact expression for the SPT. For the data used in reference [2] this was sufficient to ensure that the numerical integration error was smaller than the reported statistical error. However, this is not sufficient for the high precision simulations reported here.
To calculate γ from the current MD results, instead of subtracting the SPT expression from the integrand, we subtract the corresponding expression from a recent high-accuracy parameterization of γ, presented in reference [23]. This results in a numerical error that is smaller than the reported statistical error. A table of numerical values of γ is presented in the Supplementary Information [6].
In  At higher values of the packing fraction, η 0.3, one can observe that the FMT results slightly overestimate the wall surface free energy compared to the MD results [23]. In the inset, we highlight the region of high values of η, from which one can clearly see that the wall surface free energy obtained by the original Rosenfeld functional (red line) overestimates the values from the simulations the most, while the prediction from the White-Bear version of FMT (green line) is closest to the MD data.
Recently, a new parametrization of the simulation data for γ was given by us [23], that is accurate in the whole range of fluid densities. Especially at higher densities, an additional term with a high power in η is required in the parametrization in order to account for all the simulation data. Such a term is not reflected in the excess free energy density Φ of FMT.

Adsorption Γ and excess volume v ex
As discussed in section 2, the excess adsorption Γ can be calculated via two different routes: (i) by integrating the density profile, equation (2.3), or (ii) using the adsorption theorem, equation (2.5). We have confirmed that our DFT and MD results are consistent, in that they yield the same values of Γ via both routes. The MD data for Γ can be found in the Supplemental Information [6].
In figure 6, we show the excess adsorption as a function of the packing fraction η. At low values of η, the DFT results are indistinguishable from each other and agree very well with the MD data. In contrast to the wall surface free energy γ, shown in figure 5, FMT results start to deviate from the MD data already at a moderate value of η 0.25. For η 0.3, the results from different versions of FMT begin to deviate from each other and for η 0.4, there is a significant difference between the simulation results and those from FMT. While the excess adsorption increases towards less negative values for η 0.4 in the simulation results, the FMT results display a local minimum and then decrease further for η → 0.5. It should be noted that the presence of the local minimum is a consequence of the current choice of volume definition. If one uses another commonly used volume definition in which the position of the wall is coincident with the center of a hard sphere in contact with the wall then the local minimum disappears; however, the underestimation by the DFT of Γ at high packing fractions remains.
The excess volume, v ex as a function of η is shown in figure 7 -see the Supplemental Information volume has a direct relationship to the surface free energy through equation (2.8) and it is through the integration of this equation that the MD simulation values for γ shown in figure 5 are calculated [2]. Because v ex is directly calculated from the density profile using equation (2.6), it provides a direct link between the density profile and interfacial thermodynamics, as measured by the surface free energy.
In fact, the overestimation of v ex by the FMT at high packing fractions is dominated by a significant underestimation, at high η by the FMT, of the second and third peaks of ρ(z) ( figure 3 and 4). Because the pressure is a monotonously increasing function of η, the overestimation of v ex at high packing fractions by the FMT leads directly, through equation (2.8) to the observed overestimation of γ by the FMT.

Discussion
The overall agreement between the benchmark simulation data for the density profiles, the wall surface free energy and the excess adsorption of hard sphere fluid at a planar hard wall and the results obtained using accurate FMT DFT calculation is very good. We have found, however, that systematic deviations develop at sufficiently large packing fractions η. First, small deviations develop at around η 0.3. They become more serious for η 0.4 and show up in all quantities that we have studied here, i.e., the density profiles, the surface free energy and the excess adsorption. This was also observed in a study of confined hard-sphere fluids [24].
We find that for the dividing surface employed in our study, the actual hard wall position, the DFT results of the surface free energy γ is less sensitive to the deviations in the density profiles than the excess adsorption Γ. While the deviation in γ between the DFT and the simulation results is small, even at high packing fractions, the DFT results for Γ show far larger relative deviations compared to the benchmark simulations. At the moment it is not clear whether these deviations can be reduced by changing the form of the excess free energy density Φ -the functional form of the empirical parametrization of γ and Γ [23] could give a hint -or if they are due to the range of weight functions employed by FMT.