Within the extended Capillary Wave Theory (ECWT), to extract the bending modulus of a liquid surface from the total structure factor of the interfacial region requires to separate the capillary waves (CW) signal from a non-CW background. Some years ago, Höfling and Dietrich (HD), working in the strict grazing incidence limit qz = 0, proposed a background that combines the liquid and vapor bulk structure factors in the amounts set by Gibbs’s plane. We contrast that proposal with Molecular Dynamics (MD) simulations of the Lennard-Jones model analyzed with the Intrinsic Sampling Method (ISM). The study is extended to qz ≠ 0, to test the stronger consistency requirements of the ECWT and the experimental conditions; it shows a good MD-ECWT matching although we need some fine tuning over HD proposal. Then, the agreement with the ISM result for the surface bending modulus is good and that provides an interpretation, in terms of the molecular layering at the liquid edge, for the fluctuating surface represented by the CW signal in the surface structure factor, both for MD simulations and surface diffraction experiments.
I. INTRODUCTION
X-ray and neutron diffraction experiments may test the structure and fluctuations of liquid surfaces.1,2 With the beams on a specular plane, the total (T) structure factor ST(q‖, qz) over the interfacial region depends on , the wave vector modulus on the interface plane, and on qz normal to it. The surface reflectivity (the signal peak for q‖ = 0) is proportional to , with the Fourier transformed derivative of the density profile ρ(z) across the interface, from vapor (v) to liquid (l) bulk
That provides some experimental access to ρ(z) through the surface factor |Φ(qz)|2, as that presented in Fig. 1. For q‖ > 0, the diffuse scattering by density fluctuations, represented by ST(q‖, qz), is also proportional to |Φ(qz)|2 although with entangled dependence on qz and q‖.
The MD surface factor (full line) at T* = 0.763, and its (low qz) best Gaussian fit (dashed line), expected for a step-like intrinsic density profile (IDP). Inset, black line: mean density profile ρ(z); red line: IDP showing the molecular layering with the optimal ISM parameter n0 = 0.75, i.e., the integral of the main (z ≈ 0) IDP peak. All in LJ unit σ = 1.
The MD surface factor (full line) at T* = 0.763, and its (low qz) best Gaussian fit (dashed line), expected for a step-like intrinsic density profile (IDP). Inset, black line: mean density profile ρ(z); red line: IDP showing the molecular layering with the optimal ISM parameter n0 = 0.75, i.e., the integral of the main (z ≈ 0) IDP peak. All in LJ unit σ = 1.
The thermal capillary waves (CW) make ρ(z) and Φ(qz) to depend on the sampled area set by the section of the incident beam or by the Molecular Dynamics (MD) box size, A = L2 on the plane. Only in the strict (experimentally inaccessible) grazing incidence limit (qz → 0), the reflectivity (for q‖ = 0) and ST(q‖, qz) (for q‖ ≥ 2π/L) may be calculated in terms of the density gap between the coexisting bulk phases, Φ(0) = Δρ ≡ ρl − ρv, that is independent of A.
The Capillary Wave Theory (CWT) analyzes the CW effects, with the intrinsic density profile ρin(z) and surface factor |Φin(qz)|2, associated with the sampling of the interface over a small area (assumed of molecular size, but often not explicitly set) and referred to the instantaneous mean position of the surface over that area. The theory describes a fluctuating intrinsic surface (IS)
with CW amplitudes . The early CWT version3–7 predicted the mean square amplitudes from the surface tension γo and temperature (kT = β−1), up to an upper wavevector (q‖ ≤ qu) at the (assumed sharp) end of the mesoscopic CW fluctuations. The extended (E) version8–10 ECWT uses a wavevector dependent surface tension with a surface bending modulus K,
as in Fig. 2. Formally, any K > 0 avoids the ultraviolet catastrophe of the theory (if the CW spectrum with γo were extended to unlimited large q‖) so that K and qu could be regarded as alternative empirical ways to set an effective molecular top to the CW spectrum.11,12 However, the theory predicts13 a CW contribution to ST(q‖, qz) proportional to the mean square amplitude . Therefore, , or at least the bending modulus K, would be physical characteristics of liquids surfaces, accessible to experiments and MD simulations. Nevertheless, experimental,11,12,14–16 theoretical,8,9,17–24 and computer simulation attempts25–28 to characterize (3) led to rather confusing results.
The T* = 0.763 LJ extended surface tension γ(q‖) for several choices of the ISM parameter n0. Horizontal black dashed line: thermodynamic surface tension γo cut at qu = 0.7, to get the Gaussian fit in Fig. 1 from a step IDP. Inset: red line with squares, the bending modulus K for each choice of n0; the black square with error bars gives the optimal ISM choice n0 = 0.75 ± 0.05, with K = 0.28 ± 0.05; the purple box KHD, with HD background that would be consistent with ISM for n0 ≈ 1. All in LJ units ϵ = σ = 1.
The T* = 0.763 LJ extended surface tension γ(q‖) for several choices of the ISM parameter n0. Horizontal black dashed line: thermodynamic surface tension γo cut at qu = 0.7, to get the Gaussian fit in Fig. 1 from a step IDP. Inset: red line with squares, the bending modulus K for each choice of n0; the black square with error bars gives the optimal ISM choice n0 = 0.75 ± 0.05, with K = 0.28 ± 0.05; the purple box KHD, with HD background that would be consistent with ISM for n0 ≈ 1. All in LJ units ϵ = σ = 1.
At low q‖ and , the CW contribution to ST(q‖, qz) diverges11,12,29 as . For low qz (i.e., at or near grazing incidence), this CW contribution is large at the wavevectors (q‖ ∼ 10−3 nm−1) accessible with the size (L ∼ 1 μm) of the incident x-ray beams. The MD box size gives access only to q‖ ≥ 2π/L, with CW contributions to ST(q, qz) that, although not as large as in experiments, allow to test the CWT prediction for γ(0) = γo.
In principle, the mesoscopic description of the liquid surface as a smooth surface (2) could be formally extended up to wavelengths of molecular size (σ), i.e., q‖σ ≲ 2π. However, even with the best tuned definition for in terms of the molecular positions,30,31 MD simulations show that the ECWT assumptions are not quantitatively accurate beyond q‖σ ≈ 2. That typical range of theoretical validity is large enough to get γ(q‖) significantly different of γo and it should be enough to get the surface bending modulus in (3). Thus, K(T) would be a physical characteristics of the liquid surface, which quantifies the smooth matching between the mesoscopic and the molecular descriptions, for wavelengths of a few molecular diameters. The main difficulty is that for q‖σ ≳ 1 the CW contribution to ST(q‖, qz) is similar (or even lower) than the scattering from the bulk-like liquid region sampled at the interfacial region. The intensity of that liquid bulk signal is set by the penetration of the incident beam in experiments, while in MD we may control directly the depth of liquid and vapor bulk sampled to get ST.
Höfling and Dietrich (HD)25 presented an extensive and thoughtful analysis of this problem, with large MD simulations for the (truncated) Lennard-Jones (LJ) model. ST(q‖, qz), calculated per unit area, has a background (bg) Sbg(q‖, qz) with the contribution from non-CW sources. To make the CW contribution, H(q‖, qz) ≡ ST(q‖, qz) − Sbg(q‖, qz), independent of the depth sampled in the bulk phases, the background has to include contributions of the usual (per particle) structure factors in the bulk phases, Sl(q) and Sv(q), functions of the 3D wavevector modulus and combined as
with the nominal numbers of liquid Nl and vapor Nv molecules in the sampled volume VT.
In the strict qz = 0 limit, the ECWT gives Werthein’s relationship4,32,33 between γ(q‖) and H(q‖, 0),
Although inaccessible to experiments, this limit is very useful to simplify the analysis of MD simulations. From it, the value of γ(0) = γo is got independently of the assumed background,10,19,34 but any change Sbg → Sbg + ΔSbg would shift K by .
Within their qz = 0 analysis, HD25 assumed a background , with Gibbs’s (G) plane to set , and with an interfacial (but not-CW) contribution δSHD. This was done assuming a virtual (CW-less) interface formed by the molecules that, along the MD, happen to be at one side of a plane that cuts through the bulk phase, as an open boundary.35 Near the triple point temperature, T ≈ Tt, HD got a very flat γHD(q‖), with marginally negative KHD(Tt). At higher T, they got 0 < βKHD(T) ≲ 0.1, but all the reported functions γHD(q‖) bent down and went below γo for q‖σ ≳ 2. Thus, HD results would require to reintroduce an (undetermined) upper cutoff q‖ ≤ qu, as the early CWT versions, failing to the ECWT expectations for a raising γ(q‖) to provide a smooth molecular top to the CW spectrum. Nevertheless, HD estimation for the surface bending modulus K, as the behavior of γ(q‖) around q‖σ ≈ 1, could still be considered the result of a quite sensible assumption for Sbg(q‖, 0).
In MD, we may get γ(q‖) by a very different method. Instead using only the (experimentally accessible) pair correlations described by ST(q‖, qz), we may take advantage of the full N-particle structure over the interfacial region, since we know the instantaneous positions of all the particles. The Intrinsic Sampling Method (ISM)30,31,36 and similar algorithms37–40 define the instantaneous IS shape (2), to describe the edge of the percolating liquid cluster at any MD snapshot. The sampling of the IS fluctuations in MD allows us to get an intrinsic density profile (IDP) ρin(z), which gives a different view of the interfacial structure, much sharper for the molecular details than the usual mean density profile ρ(z). Moreover, the same analysis for a large set of MD configurations gives the mean square CW amplitudes, and therefore γISM(q‖) so that the ECWT assumption relating ρ(z) to ρin(z) and γ(q‖) may be quantitatively tested.
The aim of this paper is to compare the two routes to extract γ(q‖) within the ECWT. The ISM is obviously restricted to MD, but it gives a useful intuition for how the molecular details in the definition of (2) are related to the choice of the non-CW background Sbg(q‖, qz) in the total structure factor ST(q‖, qz) of surface diffraction experiments. We extend the analysis to qz ≠ 0 because the Φ(0) = Δρ limit misses the ECWT self-consistency between the CW blurring in ρ(z) and the mean square CW amplitude described by γ(q‖). Here and in the supplementary material, we discuss the claims of enhanced CW [γ(q‖) < γo] that had been made from the analysis of experimental results14,15 and interpreted as an effect of long-range dispersion forces.8,9
II. MD SIMULATIONS AND THEIR LINK TO CWT
NVT-MD simulations were done with LAMPPS package,41 for the LJ model truncated at rc = 4.4σ and a Nosé–Hoover thermostat. We give all our results in LJ units ϵ = σ = 1 and explore temperatures from T* ≡ kT/ϵ = 0.678 (just above the triple point) to 0.848, well below the critical one (, with our rc).42,43
N = 176 440 particles in a box with transverse (XY) size L = 41.28σ and Lz = 6L, form thick liquid slabs with the coexisting densities (ρv, ρl) and the surface tension (γo) given in Table I. Notice that our rc is larger than used by Höfling and Dietrich,25 which produces some differences in the liquid–vapor coexistence data. The statistical averages ⟨⋯⟩, over the two interfaces, are represented here with the liquid at the z > 0 side. 5000 configurations, separated by 5000 MD steps and taken after 2 × 106 MD steps for thermal equilibrium, were used to get the density profile
and, for qx,y = 2πnx,y/L (with integer nx,y, as set by the periodic boundaries) and any qz, the structure factor
The index i in (6) and (7) runs over the NT particles in the chosen region (za ≤ zi ≤ zb) of the MD box, with volume VT = L2(zb − za), and the limits inside the vapor [ρ(za) = ρv] and liquid [ρl(zb) = ρl] bulks. The index j in (7) is allowed to run far enough beyond the limits imposed to the index i as in open boundaries.35 The bulk structure factors (per particle) Sv(q) and Sl(q), used in (4) as functions of the 3D wavevector q, are obtained with qz = 0 (i.e., q = q‖) but taking both za and zb within the same bulk phase, and replacing A by NT in the denominator of (7).
MD results for the liquid–vapor surfaces studied here, in units of the LJ parameters ϵ and σ. Thermodynamic data: reduced temperature T* = kT/ϵ, vapor (ρv), and liquid (ρl) densities and surface tension γo. ISM: optimal 2D density of the outermost liquid n0 and surface bending modulus KISM. Results with the non-CW background of Höfling and Dietrich (HD): surface bending modulus KHD and the approximate equivalent of the virtual surface term in terms of the shift from Gibbs’s plane value . This work: best fitting value of the shift from Gibbs’s plane and surface bending modulus K.
. | ISM . | HD . | This Work . | |||||||
---|---|---|---|---|---|---|---|---|---|---|
T* . | ρvσ3 . | ρlσ3 . | γ0σ2/ϵ . | n0σ2 . | K/ϵ . | δnlσ2 . | K/ϵ . | Δnlσ2 . | K/ϵ (qz = 0) . | K/ϵ (qz ≠ 0) . |
0.678 | 0.0019 | 0.8428 | 1.02 | 0.80 ± 0.05 | 0.20 ± 0.05 | 0.04 | 0.0054 ± 0.01 | 0.20 ± 0.02 | 0.26 ± 0.05 | 0.27 ± 0.08 |
0.763 | 0.0053 | 0.8051 | 0.84 | 0.75 ± 0.05 | 0.28 ± 0.05 | 0.04 | 0.056 ± 0.01 | 0.27 ± 0.03 | 0.33 ± 0.05 | 0.30 ± 0.09 |
0.848 | 0.0117 | 0.7650 | 0.66 | 0.70 ± 0.05 | 0.33 ± 0.07 | 0.05 | 0.05 ± 0.01 | 0.43 ± 0.04 | 0.35 ± 0.07 | 0.10 ± 0.3 |
. | ISM . | HD . | This Work . | |||||||
---|---|---|---|---|---|---|---|---|---|---|
T* . | ρvσ3 . | ρlσ3 . | γ0σ2/ϵ . | n0σ2 . | K/ϵ . | δnlσ2 . | K/ϵ . | Δnlσ2 . | K/ϵ (qz = 0) . | K/ϵ (qz ≠ 0) . |
0.678 | 0.0019 | 0.8428 | 1.02 | 0.80 ± 0.05 | 0.20 ± 0.05 | 0.04 | 0.0054 ± 0.01 | 0.20 ± 0.02 | 0.26 ± 0.05 | 0.27 ± 0.08 |
0.763 | 0.0053 | 0.8051 | 0.84 | 0.75 ± 0.05 | 0.28 ± 0.05 | 0.04 | 0.056 ± 0.01 | 0.27 ± 0.03 | 0.33 ± 0.05 | 0.30 ± 0.09 |
0.848 | 0.0117 | 0.7650 | 0.66 | 0.70 ± 0.05 | 0.33 ± 0.07 | 0.05 | 0.05 ± 0.01 | 0.43 ± 0.04 | 0.35 ± 0.07 | 0.10 ± 0.3 |
To get HD background (extended to qz ≠ 0), we take (7) with both limits for za ≤ zi ≤ zb within the liquid (or vapor), but the restriction za ≤ zj is applied also to the second particle in the pair. That corresponds to use the plane z = za to create the virtual (CW-less) interface assumed by Höefling and Dietrich. Equivalently, the surface excess background is minus the result of (7) but taking only the pairs with zi ≤ za and zj > za for any za well within the bulk phase.
The main dependence of ST(q‖, qz) on the shape of volume VT goes into (4), to leave H ≡ ST − Sbg independent of that choice. The sharp boundaries at z = zl,v used in MD to define VT could be changed into an exponential decay toward the liquid bulk25 to represent the adsorption of the x-ray beams in diffraction experiments. That would change ST and Sbg equally to leave the same H unless the exponential decay is set to be too sharp.25 Therefore, we waive this question here and we present MD results without any damping over the sampled slab.
The CWT assumes that an intrinsic density profile (IDP), ρin(z), is locally shifted by the IS corrugations to produce the mean density profile as the average
with a Gaussian probability distribution , set by the mean position of the interface and the mean square CW amplitude
The density profile (6) and the surface factor (1) get their size dependence from the lower limit (q‖ ≥ 2π/L) in the wavevectors sum. That dependence with L is eliminated in the IDP and in the intrinsic surface factor
Notice that these generic predictions come both from the earlier or the extended versions of the theory, since they share the central hypothesis that the IDP is locally shifted (without any damping or deformation) to follow the corrugations of the IS. The early (CWT) and modern (ECWT) versions differ essentially in our expectations for the range of validity of the hypothesis. The original CWT was developed to describe the low q‖ effects that set the dependence of ρ(z) with the area and the divergence of ST(q‖, qz). The assumption that the CW amplitude is uncorrelated with the intrinsic density operator28 is very robust for low q‖. The theory may be applied without paying attention to the precise definition for or to the molecular details in ρin(z). Moreover, the Gaussian fluctuations of are well described by the macroscopic γo. The upper cutoff q‖ ≤ qu was imposed to formally avoid the ultraviolet catastrophe of the theory, but without any expectation to quantify the CW fluctuations at very short wavelengths. Thus, a typical early CWT analysis of our MD data would use the simplest step function to represent the IDP, i.e., constant Φin(qz) = Δρ carrying no information on the liquid interface at molecular scale. That would imply from (9) a Gaussian Φ(qz), which in Fig. 1 seems a fair approximation to the MD result (1). With the macroscopic γo in Eq. (9), the MD Gaussian width (at T* = 0.763) sets quσ ≈ 0.7; which may be used to accurately predict how ρ(z) and Φ(qz) change when we increase the size of the MD box. The CW divergence in ST(q‖, 0) would be fairly well represented by Wertheim’s (5) with γo. However, nothing peculiar is observed in the MD results for ST(q‖, qz) when q‖ goes across that (CWT) upper bound. That value of qu cannot be taken as a quantitative description for the physical top of the CW at molecular scale, and it is just an empirical fit to keep the theory free of unphysical divergences.
In contrast, the expectation of the ECWT is to push the validity of the CWT hypothesis to the highest possible q‖ to get a quantitatively accurate match between the mesoscopic and the molecular views. To achieve that, it matters a lot how do we define the continuous surface from the discrete set of molecular positions. The shape of the intrinsic profile has to be obtained consistently with that definition, and the Gaussian fluctuations of for large q‖ require to use a function γ(q‖) that goes beyond its γ(0) = γo limit. Then, the self-consistency requirement [(9) and (10)] becomes trivial only for qz = 0, with Φin(0) = Φ(0) = Δρ. Its validity for qz ≠ 0 provides a direct test for the assumptions of the mesoscopic theory, which fail if we push q‖ beyond a value that depends on the IS definition and at most about q‖ ≲ 2.
The qz ≠ 0 theoretical prediction for the CW contribution to ST(q‖, qz) goes beyond Wertheim’s (5), that Bedeaux and Weeks13,29 (BW) showed to be just the first term in a qz power series,
with the height–height correlation, , for any two points on the surface plane, and in (9). This SBW(q‖, qz) is the theoretical (ECWT) prediction for H(q‖, qz) ≡ ST(q‖, qz) − Sbg(q‖, qz), but we keep the different notation to remark that H(q‖, qz) reflects our choice for the non-CW correlation background, while SBW(q‖, qz) is a strong theoretical prediction in which the joint dependence on q‖ and qz comes from their separated dependence in γ(q‖) and Φ(qz). Our goal of H(q‖, qz) ≈ SBW(q‖, qz), could only be obtained with consistent results for γ(q‖) and Sbg(q‖, qz), that get the best for the ECWT assumptions. In contrast, if we restrict ourselves to the qz = 0 case, Wertheim’s prediction (5) would give directly a function γ(q‖) for any choice of Sbg(q‖, 0), with no touchstone to test the theoretical consistency.
The rules to define the IS shape from the molecular positions are obviously open to choice.34 Given a rule, the statistical average over MD configurations gives the IDP, , and the function γ(q‖) that quantifies the surface fluctuations.30 The ISM30,31,36 is a method based on a percolation analysis and an iterative procedure to locate the particles selected as surface pivots, representative for the outermost layer in the liquid. These pivots determine the instantaneous shape of as a Voronoi tessellation31 that Fourier transformed gives the CW amplitudes . The selection of the surface pivots is crucial to push the ECWT to molecular scales. The main parameter in the ISM is the 2D density of those pivots, n0, i.e., the integral of the first peak in the IDP (inset in Fig. 1).
With early insight, Stillinger44 foresaw that the molecular packing in a dense liquid should produce a strong layering effect from the outermost liquid layer of the IDP toward the inner bulk, as it produces spherical shells in the pair distribution function of the bulk liquid. The ISM focuses on this molecular layering structure to find the best ECWT-MD link, i.e., the optimal value for n0 to describe the outermost liquid layer. Other physical properties (kinetics, dynamics,31 and even heat transport45) have confirmed that the molecular layering is indeed a main characteristic of liquid surfaces (at least for T/Tc ≲ 0.8). The ISM may find its clearest view through the optimal value of n0 that depend on the temperature and the molecular interactions.31
Figure 2 presents γ(q‖) at T* = 0.763, over the range 0.5 ≤ n0 ≤ 1 for the ISM parameter. The fits to (3) have common γo but different surface bending modulus (see the inset) that cover the range 0.6 ≳ K ≳ 0. HD background gives, through (5), a low KHD that could be reached in ISM with n0 ≈ 1. That n0 is well above the ISM optimal choice, n0 = 0.75 ± 0.05, that gives the best representation of the intrinsic layering structure (Fig. 1 inset). It is only after such selection of the optimal n0 (and within the error bars associated with that choice) that we get the ISM prediction for the surface bending modulus KISM(T*) given in Table I, always well above KHD(T*). Moreover, it is only with the optimal n0 that the ECWT-MD matching is extended to wavelength of just a few molecular diameters (i.e., up to q‖ ≈ 2), with quantitative accuracy for the ECWT predictions [(8)–(10)] that link ρ(z), and γISM(q‖).
The MD snapshot in Fig. 3 (bottom) gives a pictorial idea on how the choice of n0 is reflected in the IS. The coarse shape of , (i.e., the surface undulations with longer wavelengths) is quite robust, similar for the two choices of n0 presented in the MD snapshot. For these low q‖ values, the statistical analysis over thousands of MD configurations extrapolates to γISM(0) ≈ γo, in agreement with the macroscopic surface tension obtained by the virial method. However, at the shorter wavelengths that determine the surface bending modulus K in (3), the choice of n0 becomes relevant. From a single MD snapshot, we could hardly determine which value of n0 gives the best IS representation of the interface, but their statistical analysis, along the full MD trajectory, makes clear that the view of the surface layering is much better with n0 = 0.75 (leading to KISM = 0.28 ± 0.05) than with n0 = 1, as needed to recover KHD ≈ 0.
Bottom: Typical MD snapshot for a narrow slice (Δy = 2σ) across the liquid surface for the T* = 0.763 LJ surface. The dots give a few cuts of the ISM intrinsic surface within the slice: black dots with the optimal n0 = 0.75 that gives K = 0.28; red dots with n0 = 1 and K ≈ 0 (see Fig. 2). Top: snapshot for a virtual (CW-less) surface, as used by Höfling and Dietrich25 to define the background .
Bottom: Typical MD snapshot for a narrow slice (Δy = 2σ) across the liquid surface for the T* = 0.763 LJ surface. The dots give a few cuts of the ISM intrinsic surface within the slice: black dots with the optimal n0 = 0.75 that gives K = 0.28; red dots with n0 = 1 and K ≈ 0 (see Fig. 2). Top: snapshot for a virtual (CW-less) surface, as used by Höfling and Dietrich25 to define the background .
What we propose here is to take the MD-ISM results for γISM(q‖) and Φ(qz) into SBW(q‖, qz), in (11). Then, we run in reverse the problem studied by Höfling and Dietrich, i.e., we get the background that is consistent with the molecular definition for . This exercise had been done28 for the pair correlation function G(z, z′, q‖) and its corresponding BW series GBW(z, z′, q‖), with real space variables z and z′ that are Fourier transformed as z − z′ to get ST(q‖, qz).46 The result for G − GBW was compared with a direct estimation of the non-CW background from an intrinsic pair correlation obtained directly from the MD-ISM analysis,28 as a check for the internal consistency of the ISM-BW link. However, the function G(z, z′, q‖) is not experimentally accessible. Our present study compares with HD proposal , and it lead to a molecular interpretation of the function γ(q‖) that is extracted from surface diffraction experiments.
III. STRUCTURE FACTOR IN OBLIQUE PLANES
An experimental way to extract the non-CW background in ST(q‖, qz) is to keep the same relative geometry for the incident and diffracted beams but to rotate their plane an angle ϕ out of the reflection plane. If , then most of the CW contribution would be cancelled, , while keeps unchanged the contribution of the bulk phases to ST(q‖, qz) that depends on the 3D wavevector modulus q = q′.
At low qz, the CW signal with may still be important and may be a poor approximation to H(q‖, qz). However, if Sbg(q‖, qz) depends mainly on the 3D q, ΔST should be similar to the equivalent difference between the BW predictions .
Figure 4 compares ΔST(q‖, qz; ϕ) and ΔSBW(q‖, qz; ϕ) as functions of qz and with ϕ = π/4. We present the lowest transverse wavevector in our MD box, q‖ = 0.152 [panel (a)], and a larger q‖ = 0.456 [panel (b)]. The lowest qz > 0 accessible with the periodic boundaries for each q‖ is that which makes , reducing to a half the CW factor in (11). As qz increases, grows and its CW contribution vanishes. Therefore, we expect a smooth trend from ΔST ≈ H/2 for the lowest qz, to ΔST ≈ H for large qz. The excellent agreement between the MD result for ΔST and the mesoscopic prediction ΔSBW supports the accuracy of BW theory and the hypothesis that Sbg(q‖, qz) may be approximated by a function of the 3D wavevector q.
Difference ΔST(q‖, qz; ϕ) (blue squares) between the MD specular structure factor and that in a plane rotated by ϕ = π/4 for the T* = 0.763 LJ interface; the red lines are the equivalent BW predictions ΔSBW(q‖, qz; ϕ). Panel (a): q‖ = 0.152; panel (b): q‖ = 0.456. We use LJ units σ = 1.
Difference ΔST(q‖, qz; ϕ) (blue squares) between the MD specular structure factor and that in a plane rotated by ϕ = π/4 for the T* = 0.763 LJ interface; the red lines are the equivalent BW predictions ΔSBW(q‖, qz; ϕ). Panel (a): q‖ = 0.152; panel (b): q‖ = 0.456. We use LJ units σ = 1.
IV. ISM-MD CONSISTENT BACKGROUND
Figure 5 presents the MD result (7) for ST(q‖, qz) at T* = 0.763, and the CW contribution SBW(q‖, qz), built as (11) with γISM(q‖) and Φ(qz). The data are presented for several q‖ ≥ 0.152 and as functions of the 3D q. For q ≲ 0.5, ST(q‖, qz) reflects the larger CW contribution for smaller q‖, but for larger q, the bulk-like contributions dominate and for q ≳ 3 all the (q‖, qz) points merge in a single line. Notice that as q‖ grows the decay of the CW contributions is dominated by |Φ(qz)|2 rather than by q‖, consistently with (10). The relative weight of CW vs bulk contributions depends on the shape of the volume used in Eq. (7), that we take here to give NT/A ≈ 13σ−2 (about 15σ depth in both the liquid and vapor sides). Assuming that the CW contribution is represented by H(q‖, qz) = SBW(q‖, qz), we get the difference , as the ISM-consistent background, which is close to be a function of the 3D q over its full range, with all the values of (q‖, qz) merged in a narrow band.
MD results (7) for ST(q‖, qz) (full circles, with dashed lines). The full lines give the CW contribution SBW(q‖, qz) from (11) using γISM(q‖) and the MD density profile ρ(z). The empty squares are the non-CW background . The results, at T* = 0.763 and in LJ units (σ = 1), are presented as functions of the 3D wavevector q, for different values of q‖ = 2π/L = 0.152 (black), 0.304 (red), 0.456 (blue), and 0.608 (orange). The inset shows a detailed study of the background contribution. Green line: average of the over q‖; red line: HD proposal .
MD results (7) for ST(q‖, qz) (full circles, with dashed lines). The full lines give the CW contribution SBW(q‖, qz) from (11) using γISM(q‖) and the MD density profile ρ(z). The empty squares are the non-CW background . The results, at T* = 0.763 and in LJ units (σ = 1), are presented as functions of the 3D wavevector q, for different values of q‖ = 2π/L = 0.152 (black), 0.304 (red), 0.456 (blue), and 0.608 (orange). The inset shows a detailed study of the background contribution. Green line: average of the over q‖; red line: HD proposal .
The inset in Fig. 5 shows a clear gap between our ISM background and HD proposal. Notice that the original HD surface term , that reflects a virtual (CW-less) interface, has to be extended to qz ≠ 0 but still it may be obtained from the bulk structure factors, or directly from the MD sampling (7). Although strictly is not just a function of the 3D q, in practice (within the error bars of the MD results), its contribution may still be assimilated to a small change of Nl with respect to Gibbs’s value in Slv(q), keeping fix Nl + Nv = NT, with the values of given in Table I. A larger change Δnl of that effective number of liquid-like particles may be used to close the gap between and . The value of the best fitting Δnl depends on T* (see Table I) and it is only a fraction of the (ISM) liquid monolayer, i.e., equivalent to shift Gibbs’s plane by a fraction of a molecular diameter.
As an alternative to the above discussion, we could keep HD background proposal (i.e., , Δnl = 0), and change the ISM parameter n0 for the LJ surface, to achieve . However, with the insight of the N-particle structure in MD, we know that the mesoscopic view of the liquid surface given by would be spoiled if we take n0 = 1, instead of its optimal value n0 ≈ 0.75. Figure 3 shows that the sensible HD proposal, to represent the interface as a virtual cut across the liquid bulk located exactly at Gibbs’s plane, should not be taken too strictly. The local correlation structure at the interface (bottom panel) may likely differ from that at the virtual interface (top panel), with an effect equivalent to the required change Δnl in the liquid-like contribution to the Slv(q) background.
Figure 6 gives a more detailed comparison of ST, and for a lower temperature (T* = 0.678). We use now the on-plane q‖ as variable for qz = 0 (as HD) and for qz = 1.0622. At low q‖, when ST(q‖, qz) and SBW(q‖, qz) are large and nearly equal, we have to assume a larger error bar for their difference , particularly for qz = 0. However, the results clearly suggest changing HD background proposal allowing for a Δnl ≠ 0. Sticking to Gibbs’s plane to set Nl not only changes the non-CW background at low q‖ (therefore pushing KHD well below KISM) but it fails to recover Sbg(q‖, qz) = ST(q‖, qz) at large q‖ when the CW contribution has to become negligible. The difference ST(q‖, qz) − Sbg(q‖, qz) ≈ 0.1 produces a sharp drop of γ(q‖) for q‖ ≳ 2, as reported in HD results.25 The accuracy of our proposed fit is not perfect, particularly for qz ≠ 0 and low q‖. However, the useful range to estimate the surface bending modulus is 0.3 ≲ q‖ ≲ 1.2, and in that range the best fitting Slv(q) (with the Δnl in Table I) gives an evident improvement over the original HD proposal. Keeping HD virtual surface term in the fit makes too little difference to be worth, since it is absorbed by a small change in Δnl.
Non-CW background (green squares and lines), at T* = 0.678, for qz = 0 (top) and qz = 1.0622 (bottom). The results, in LJ σ = 1 units, are presented as functions of the on-plane wavevector q‖. Red lines: HD background (with , ). Full orange lines: fitted Slv(q) with . Dashed orange lines: similar fits (with Δnl − δnl) if we include the HD (virtual) surface contribution . Blue circles and lines: the total ST(q, qz) to show the CW relevance for these q‖.
Non-CW background (green squares and lines), at T* = 0.678, for qz = 0 (top) and qz = 1.0622 (bottom). The results, in LJ σ = 1 units, are presented as functions of the on-plane wavevector q‖. Red lines: HD background (with , ). Full orange lines: fitted Slv(q) with . Dashed orange lines: similar fits (with Δnl − δnl) if we include the HD (virtual) surface contribution . Blue circles and lines: the total ST(q, qz) to show the CW relevance for these q‖.
The supplementary material provides evidence that the description of Sbg(q‖, qz) ≈ Slv(q), applied here to the LJ liquid surface, is also very accurate for a different (SA28) molecular interaction model.
V. BEST FITTING BULK-LIKE BACKGROUND
In Sec. IV, we took advantage of the ISM that uses the positions of all the particles (as given by MD) and produces a γISM(q‖) optimized to give the best mesoscopic description of the liquid surface. Then, we got the non-CW background that, within the ECWT, is consistent with the MD results for pair correlation structure ST(q‖, qz). The good news is that is not far from the proposal of Höfling and Dietrich25 , which combines the structure factors of the bulk liquid and vapor phases; although it requires some fine-tuning in the effective contribution of each phase in (4), rather than taking for granted Gibbs’s plane values for Nl and Nv.
In surface diffraction experiments, the estimation of the non-CW background has to be done directly from the total ST(q‖, qz). In the analysis of MD simulations, the use to fit , as we have done in Sec. IV, requires the prior application of the ISM to get γISM(q‖). It would make little sense, as a practical method, to use γISM(q‖) to extract from ST a γ(q‖) that should be very similar to the input. In this section, we show how γ(q‖) may be obtained solely from the pair correlations in ST(q‖, qz) and that it comes to be similar to γISM(q‖). To achieve that, we take advantage of the rapid decay of the CW contribution SBW(q‖, qz) with increasing qz, so we fit Slv(q) ≈ ST(q‖, qz) over the 3D q range in which the CW contribution is negligible. Then, the best fitting value of Δnl is used to approximate the non-CW background Sbg(q‖, qz) ≈ Slv(q) over the lower values of q‖ that give access to γ(q‖).
Figure 7 presents the fit ST(q‖, qz) ≈ Slv(q) at T* = 0.678, over the range 2.5 ≲ q ≲ 4. At that q-range, there is still a visible CW contribution for the lowest q‖. However, the data for q‖ ≳ 0.6 collapse into a single line that we may approximate by Slv(q) with Δnl = 0.20 ± 0.02. The range of uncertainty reflects the error bars in the MD results and the precise range of q used for the fit. Although we fit at larger q than in Sec. IV, the estimated backgrounds are very similar, which supports Slv(q) to represent Sbg(q‖, qz). Notice that to use still higher q for the fit gives less accuracy for Δnl because Slv(q) becomes very steep.
MD ST(q‖, qz) at T* = 0.678, as function of the 3D wavevector q. Full circles: q‖ = 0.152 (black), q‖ = 0.304 (red), blue: q‖ = 0.456 (blue) q‖ = 0.608 (orange), q‖ = 0.912 (green), and q‖ = 1.216 (cyan). Red line: Höfling–Dietrich background,25 , with Gibbs’s surface to set Nl. Orange line: best fitting to Slv(q), done for the highest q‖ over the range 2.5 ≤ q ≤ 4., with Δnl in Table I. The bottom panel gives an enlarged view of the relevant region for the fit, with the estimation for the error bars in Δnl (dashed orange lines).
MD ST(q‖, qz) at T* = 0.678, as function of the 3D wavevector q. Full circles: q‖ = 0.152 (black), q‖ = 0.304 (red), blue: q‖ = 0.456 (blue) q‖ = 0.608 (orange), q‖ = 0.912 (green), and q‖ = 1.216 (cyan). Red line: Höfling–Dietrich background,25 , with Gibbs’s surface to set Nl. Orange line: best fitting to Slv(q), done for the highest q‖ over the range 2.5 ≤ q ≤ 4., with Δnl in Table I. The bottom panel gives an enlarged view of the relevant region for the fit, with the estimation for the error bars in Δnl (dashed orange lines).
Then, in our MD, we may use qz = 0 Wertheim’s equation (5) to get γ(q‖) from H(q‖, 0) = ST(q‖, 0) − Slv(q). In Fig. 8 we present the results of this method at two T*, to compare with γISM(q‖) and γHD(q‖). Table I gives K ± ΔK as extracted by this method. The error bars are an important part of the physical information given by γ(q‖). It is obvious that we cannot expect a unique and exact IS to describe the molecular positions in a MD snapshot. Apparently, mild changes in the IS shape (as in Fig. 3) produce important changes in the bending modulus K. The ISM optimal value for n0 has its own error bar, since it comes from a quite robust but soft combination of different physical aspects (surface structure, kinetics, dynamics, …), sampled over many configurations along the MD trajectory.31 Therefore, we present in Fig. 8 the results with K = KISM ± ΔKISM (as given in Table I), that bracket those obtained from ST(q‖, qz). Notice that the uncertainty in γ(q‖) grows with q‖, but still it leaves a clear agreement between the estimates for the surface bending modulus given by the ISM and by the present method of analysis for ST(q‖|, qz). The results with the original HD proposal give much lower KHD ≈ 0. Our cutoff of the LJ potential is not the same as used by HD;25 we refer here to the results of HD method with our MD data.
γ(q) from H(q, 0) through Eq. (5), at T* = 0.678 (top) and T* = 0.848 (bottom). Full blue circles: background fitted to Slv(q) = ST(q‖, qz), with error bars from those for the best fitting Δnl in Table I. Blue dashed lines: HD background. Full red lines, ISM results for the optimal n0, with red dashed lines for the error bars in KISM, as given in Table I.
γ(q) from H(q, 0) through Eq. (5), at T* = 0.678 (top) and T* = 0.848 (bottom). Full blue circles: background fitted to Slv(q) = ST(q‖, qz), with error bars from those for the best fitting Δnl in Table I. Blue dashed lines: HD background. Full red lines, ISM results for the optimal n0, with red dashed lines for the error bars in KISM, as given in Table I.
The results in Table I show a clear trend in the dependence with the temperature. While γo decays with increasing T*, the surface bending modulus K increases. Thus, the mean square CW amplitude for low q‖ grows with the temperature, but the CW fluctuations become restricted to a narrower range of q‖. This tendency was already observed (although with their lower K) by Höfling and Dietrich,25 who compared the length with the correlation length in the liquid bulk. At the low T* explored here, we get that ℓ is just a fraction of the molecular diameter. Near Tc, when the bulk correlations are dominated by the compressibility, rather than by the molecular coordination shells, the ISM becomes less discerning. Therefore, the critical behavior of K(T) or ℓ(T) is out of our present analysis.
VI. ON THE ANALYSIS OF EXPERIMENTAL SURFACE DIFFRACTION DATA
For the analysis of experimental surface diffraction data,14,15 the separation of the non-CW background from the CW signal, H = ST − Sbg, has to be done essentially as for MD data in Sec. V, with a fit ST(q‖, qz) ≈ Slv(q) over q-range where the CW effects are negligible. We cannot use Gibbs’s plane to get and (taking into account the exponential decay from the X-rays adsorption25) since the experimental reflectivity data do not give the location of ρ(z) relative to the signal in ST. Moreover, the strict qz = 0 limit is never reached in experiments so that the ECWT analysis has to be done with the full BW series (11) that may be summed as
To extract γ(q‖) from equating H(q‖, qz) = SBW(q‖, qz) is not as simple as in the qz = 0 case (5), since to get (12) for each value of q‖ we need the full function γ(q‖). As a practical method, we may truncate (3) at the order, leaving K as the only free parameter. Given γo and |Φ(qz)|2, we look for the value of K that best fit H(q‖, qz) ≈ SBW(q‖, qz) over the available range for qz, with . That would be our best estimation for K, extracted from Φ(qz) and ST(q‖, qz) at qz > 0, as it would be accessible to experimental data. We follow this method, to get the results in Fig. 9, which should be compared with those in Fig. 8 for the same MD data analyzed in the qz = 0 limit. The estimates for K ± ΔK (Table I) are compatible, but the error bars become larger particularly at the highest T*, from the more convoluted analysis for qz ≠ 0. The use of Gibbs’s plane to determine Sbg (i.e., HD method extended to qz ≠ 0) becomes even worse than for qz = 0.
Equivalent at qz ≠ 0 for the (qz = 0) analysis presented in Fig. 8 at the same T*. The CW signal H = ST − Sbg is represented as , following Eq. (14), to get γ(q‖) except for a constant factor (including the size of the MD box). Blue circles with error bars: results with our best fitting background. Red lines: the same representation (and box size) for SBW using , with K ± ΔK as a free parameter, to get the mean value and to bracket the (blue) error bars (see Table I). Dashed blue lines: results with HD background (extended to qz ≠ 0).
Equivalent at qz ≠ 0 for the (qz = 0) analysis presented in Fig. 8 at the same T*. The CW signal H = ST − Sbg is represented as , following Eq. (14), to get γ(q‖) except for a constant factor (including the size of the MD box). Blue circles with error bars: results with our best fitting background. Red lines: the same representation (and box size) for SBW using , with K ± ΔK as a free parameter, to get the mean value and to bracket the (blue) error bars (see Table I). Dashed blue lines: results with HD background (extended to qz ≠ 0).
Under experimental conditions, the continuous limit q‖L ≫ 1 may be used to simplify (12) as
where behaves as log(x/xo) for large distances, with a molecular-like length parameter xo that reflects the deviation of γ(q‖) from γo. The upper limit for the integral in (13) should be treated as a smoothed far boundary for the oscillatory decay of Bessel function Jo(q‖xo). Thus, with the intrinsic surface factor (10), the dependence on the system size and the boundary conditions disappears in (13). Following Sinha et al.,29,47 if L−1 ≪ q‖ ≪ σ−1, a still simpler approximation to (13) may be obtained as
with the gamma function and the exponent , which carries all the dependence on γ(q‖). For low qz, η(q‖) ≪ 1 and . Therefore, a fairly accurate and simple way to get γ(q‖) out of H ≡ ST(q‖, qz) − Slv(q) is to keep just the γ(q‖) in the denominator of (14) and use otherwise ηo instead of η(q‖) [i.e., γo rather than γ(q‖)].47 This is the usual procedure for the analysis of the experimental diffraction data,14,15 although it is not accurate for the analysis of MD simulations, since condition L−1 ≪ q‖ ≪ σ−1 becomes too narrow and the continuous limit (13) is not valid for the required values of q‖.
Enhanced CW, i.e., functions γ(q‖) with a minimum well below γo, for q‖ ∼ 1 nm−1, were reported14,15 from x-ray surface diffraction data. The key difference between the analysis used in those reports and our analysis here is that the non-CW background of the liquid bulk was assumed to be a constant and interpreted as proportional to the compressibility, i.e., to the limit Sl(0) rather than keeping the dependence Sl(q) over the range used to extract H(q‖, qz) and γ(q‖).
In Fig. 10, we mimic that procedure to interpret our MD results for ST(q‖, qz) with a constant background . Although its value is assumed to represent the q = 0 limit of the non-CW signal (times the sampled depth of liquid), it has to be fitted to the values of ST(q‖, qz) with q ≳ 2, where the CW contribution may be considered negligible. Moreover, we expect , since otherwise the assumed CW signal would become negative. Therefore, in Fig. 10 (top), we fix a value of q‖ and take as the minimum of ST(q‖, qz), which happens around q ∼ 1.5 − 2.5, in between the q = 0 compressibility peak and the main (q ≈ 2π) peak of Sl(q). Now, we may use that (assumed constant) background to obtain H and γ(q‖), as we have done with Sbg(q) ≈ Slv(q), or with . We take the simplest qz = 0 case for this analysis, with the results shown in Fig. 10 (bottom). The function γ(q‖) goes through an apparent enhanced CW range γ(q‖) < γo because for low q we are missing the compressibility peak in which . That peak (of bulk like correlations) is therefore taken as part of the CW contribution H0. Hence, the (spuriously) large fluctuation of the interface is interpreted as a too low surface tension for those q‖ values. As q‖ increases, approaching the point where , γ(q‖) has to raise back [since H0(q‖, 0) goes to zero], producing a shape qualitatively similar to that reported as enhanced capillary waves in some experimental studies.14,15
Top panel: Mimic with our MD results for ST(q‖, qz) of the experimental procedure to assume a constant non-CW background , which (formally) represents the bulk (q = 0) compressibility, but in practice is fitted at q large enough to discard any CW contribution in ST, as explained in the text. Bottom panel: the apparent enhanced CW signature, i.e., γ(q) < γo, for the two choices (orange and cyan lines) of the background (same colors in top panel). The same ST(q‖, qz) interpreted in this work with Sbg(q) = Slv(q) (black circles and line) gives a positive bending modulus, i.e., γ(q‖) ≥ γo, similar to the ISM result (red line).
Top panel: Mimic with our MD results for ST(q‖, qz) of the experimental procedure to assume a constant non-CW background , which (formally) represents the bulk (q = 0) compressibility, but in practice is fitted at q large enough to discard any CW contribution in ST, as explained in the text. Bottom panel: the apparent enhanced CW signature, i.e., γ(q) < γo, for the two choices (orange and cyan lines) of the background (same colors in top panel). The same ST(q‖, qz) interpreted in this work with Sbg(q) = Slv(q) (black circles and line) gives a positive bending modulus, i.e., γ(q‖) ≥ γo, similar to the ISM result (red line).
VII. DISCUSSION AND CONCLUSIONS
At the never still surface of a liquid, the capillary wave (CW) spectrum goes continuously from molecular to macroscopic sizes, without any gap to set the statistical physics separation between micro and macro-states. The lack of a macroscopic limit (until Earth gravity sets it at millimeter range) produces the size dependence of the density profile and the q‖ → 0 divergence of the surface correlation, as observed in MD simulations and surface diffraction experiments.5 The ECWT, the mesoscopic theory to describe these effects, uses an extended wavevector dependent surface tension γ(q‖) to match the molecular and mesoscopic descriptions. Bedeaux–Weeks (BW) series (11) is the ECWT prediction for the CW contribution to the structure factor H(q‖, qz) ≡ ST(q‖, qz) − Sbg(q‖, qz), but to extract γ(q‖) from ST(q‖, qz) has been difficult and controversial.
As pointed by Höfling and Dietrich (HD),25 the problem is the non-CW background Sbg(q‖, qz), that gives a regular (non-diverging at q‖ = 0) but important contribution to ST(q‖, qz). The thermodynamic limit γ(0) = γo is independent of Sbg, but the physical description of the surface fluctuations with wavelengths of a few molecular diameters, i.e., the molecular top to the CW spectrum within the ECWT, is described by the surface bending modulus K in (3), which depends on Sbg.
HD proposed a background, , defined on physically sound terms, and they used it to evaluate a γHD(q‖) (and KHD) from MD simulations for the (truncated) LJ liquid surface. They took the qz = 0 limit that simplifies the analysis, and Gibbs’s plane to set the Nl (liquid-like) and Nv (vapor-like) contributions of the particles with bulk-like correlations sampled within ST(q‖, qz), together with a virtual (CW free) interface, as a cut through the bulk phase. We build here on that HD proposal with three points of support:
We take advantage of the ISM, a method based on the full (N-particles) structure of the interface to get γISM(q‖) consistently with the other crucial ingredient of the ECWT, the intrinsic density profile ρin(z). The method can only be applied to MD simulations, but its comparison with the analysis of the pair-correlations in ST gives a molecular interpretation for the function γ(q‖) extracted from x-ray surface diffraction experiments.
Our study covers qz ≠ 0, which requires a more complex analysis but unveils the strong consistency requirements of the ECWT, linking the mean (size dependent) density profile ρ(z) with ρin(z) and γ(q‖).
We relax HD use of Gibbs’s plane to set Nl and Nv in the definition of the background.
The latter point comes out as a consequence of the first two, and it heals the dropping of γHD(q‖) that goes well below γo for q‖ ≳ 2/σ, a behavior that would produce unphysical predictions if that γHD(q‖) were used within the ECWT. The fall of γHD(q‖) at large q‖ comes from the HD underestimation of the background, which fails to recover the full structure factor ST when the CW effects should vanish. That may be achieved with a minor (but crucial!) modification of Nl and Nv from the values given by Gibbs’s plane.44 Notice that HD setting is a natural and reasonable choice, but not a physical requirement as strong as getting ST = Sbg for large q‖. The latter requirement may be achieved with a small change (just a fraction of a monolayer) in Nl.
A serendipity of that large q‖ requirement is that the behavior of γ(q‖) changes also over the CW range at lower q‖. That takes the bending modulus to a value much larger than KHD and quite close (within the error bars) to the ISM value. This is a most welcome agreement, not only because it gives support to the results of both methods, based on very different analysis of MD simulations, but also because it provides a molecular image for the fluctuating mesoscopic surface that describes the edge of the liquid through the function γ(q‖), as extracted from the surface structure factor.
The ISM and its estimation of the bending modulus KISM are based in the definition of an intrinsic surface that optimizes the intrinsic layering structure at the liquid surface, as it had been foreseen by Stillinger long time ago.44 Getting similar values for K from the pair-correlation structure implies that the surface bending modulus K may be regarded as soft but robust (physically measurable, although with inherent error bars) property of liquid surfaces. It is soft because we cannot expect a fully determined value of K, since it gives a mesoscopic description of the liquid surface that would always be open to our choice for the mathematical shape that we assign to a MD snapshot, or to the non-CW background that we use to extract H(q‖, qz) and γ(q‖) from ST(q‖, qz); and it is robust because we may get similar values through very different ways, asking for the best fitting bulk-like background or asking for the best view of the intrinsic layering structure. Therefore, the raising γ(q‖) (positive K) obtained from the CW signal H(q‖, qz) = ST(q‖, qz) − Slv(q) with the best fitting bulk-like background, should be associated with the fluctuating surface that best represents the local displacements of the (highly correlated) molecular layers at the edge of the liquid.
This interpretation gives a tight physical link between the mesoscopic and the molecular views of the interface because it provides a smooth high q‖ end to the CW spectrum in H(q‖, qz). The surface layering is produced by the short-ranged excluded volume interactions in the dense liquid and it represents precisely what the background Sbg has to take out of ST. Therefore, we may understand why taking out that bulk-like contribution is equivalent to define the intrinsic surface that better follows that local layering structure at the edge of the liquid, as the ISM does. Notice that in our analysis for ST, as in general for the applications of the ISM, we cannot take the liquid too close to the critical temperature. For T/Tc ≳ 0.85, the 3D correlation length associated with q ≈ 0 fluctuations becomes comparable or larger than the correlation length for layering (at q ≈ 2π) and the use of that molecular structure to identify a physical intrinsic surface becomes troublesome. The behavior of γ(q‖) approaching the critical liquid–vapor region could only be afforded with the parallel exploration of molecular scale definitions for the IS and of approximations for the non-CW background when the fluctuations in the liquid bulk are dominated by the q = 0 compressibility peak.
Finally, we have shown that the accuracy of the background Sbg(q‖, qz) ≈ Slv(q) (fitted with the Δnl parameter) goes much beyond the assumption of a constant background, Sbg(q‖, qz) ≈ Sl(0)Nl/A, which was used in the interpretation of x-ray diffraction experiments14,15 with the apparent result of strongly enhanced CW (i.e., γ(q‖) well below γo) and attributed there to the nonanalytic behavior of the 1/r6 tails in van der Waals interactions. A fit to ST(q‖, qz) over the region of its minimum (q‖ ≈ 2 − 3) goes below the maximum of at q = 0. Therefore, those bulk-like fluctuations are spuriously interpreted as enhanced CW since they would correspond to γ(q‖) < γo. The ISM analysis for LJ surfaces, as a function of the cutoff distance for the potential tail,48 observed the dispersion force effects in γISM(q‖), as generically predicted by Napiórkowski and Dietrich,8 but its quantitative effect is small49 and the non-analytic behavior induced in γ/q‖) appears to be well below what could be extracted from ST(q‖, qz).
Our computer codes for the ISM30,31 and for the present analysis of ST(q‖, qz) have been and are available by request to the corresponding author. Access to the files with the MD data used for our present analysis could also be provided on reasonable request.
SUPPLEMENTARY MATERIAL
The supplementary material provides evidence that the description of Sbg(q‖, qz) ≈ Slv(q), applied here to the LJ liquid surface, is also very accurate for a different (SA28) molecular interaction model. In light of this result, we discuss enhanced CW [γ(q‖) < γo].
ACKNOWLEDGMENTS
We acknowledge the support of the Spanish Secretariat for Research, Development, and Innovation under Grant Nos. PID2020-117080RB and PGC2018-096955-B-C44 and from the Maria de Maeztu Programme for Units of Excellence in R&D (Grant No. CEX2018-000805-M).
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Jose Hernández-Muñoz: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Pedro Tarazona: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Enrique Chacón: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
REFERENCES
Theoretical analyses have often been done in terms of the so-called local structure factor, S(q‖, z), which integrates G(q‖, z, z′) over z′ and keeps the real space dependence on z. This function is not experimentally accessible, and it misses the z − z′ dependence that is crucial to describe the molecular layering effects in ST(q‖, qz).