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 *q*_{z} = 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 *q*_{z} ≠ 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 *S*_{T}(*q*_{‖}, *q*_{z}) over the interfacial region depends on $q\Vert =[qx2+qy2]1/2$, the wave vector modulus on the interface plane, and on *q*_{z} normal to it. The surface reflectivity (the signal peak for *q*_{‖} = 0) is proportional to $|\Phi (qz)|2/qz2$, 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* |Φ(*q*_{z})|^{2}, as that presented in Fig. 1. For *q*_{‖} > 0, the diffuse scattering by density fluctuations, represented by *S*_{T}(*q*_{‖}, *q*_{z}), is also proportional to |Φ(*q*_{z})|^{2} although with entangled dependence on *q*_{z} and *q*_{‖}.

The thermal capillary waves (CW) make *ρ*(*z*) and Φ(*q*_{z}) to depend on the sampled area set by the section of the incident beam or by the Molecular Dynamics (MD) box size, *A* = *L*^{2} on the $x\u20d7=(x,y)$ plane. Only in the strict (experimentally inaccessible) grazing incidence limit (*q*_{z} → 0), the reflectivity (for *q*_{‖} = 0) and *S*_{T}(*q*_{‖}, *q*_{z}) (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}(*q*_{z})|^{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 $\xi \u0302(q\Vert \u20d7)$. The early CWT version^{3–7} predicted the mean square amplitudes $\u27e8|\xi \u0302(q\Vert )|2\u27e9$ from the surface tension *γ*_{o} and temperature (*kT* = *β*^{−1}), up to an upper wavevector (*q*_{‖} ≤ *q*_{u}) at the (assumed sharp) end of the *mesoscopic* CW fluctuations. The *extended* (E) version^{8–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 *q*_{u} could be regarded as alternative empirical ways to set an effective molecular top to the CW spectrum.^{11,12} However, the theory predicts^{13} a CW contribution to *S*_{T}(*q*_{‖}, *q*_{z}) proportional to the mean square amplitude $\u27e8|\xi \u0302(q\u20d7\Vert )|2\u27e9$. Therefore, $\gamma (q\Vert )\u2261kT/[q\Vert 2\u27e8|\xi \u0302(q\u20d7\Vert )|2\u27e9A]$, 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 attempts^{25–28} to characterize (3) led to rather confusing results.

At low *q*_{‖} and $\eta o\u2261qz2/(2\pi \beta \gamma o)<2$, the CW contribution to *S*_{T}(*q*_{‖}, *q*_{z}) diverges^{11,12,29} as $|\Phi (q\u2009z)|2/(\beta \gamma oq\Vert 2\u2212\eta o)$. For low *q*_{z} (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 *S*_{T}(*q*, *q*_{z}) 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 $\xi (x\u20d7)$ 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 *S*_{T}(*q*_{‖}, *q*_{z}) 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 *S*_{T}.

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. *S*_{T}(*q*_{‖}, *q*_{z}), calculated per unit area, has a background (bg) *S*_{bg}(*q*_{‖}, *q*_{z}) with the contribution from non-CW sources. To make the CW contribution, *H*(*q*_{‖}, *q*_{z}) ≡ *S*_{T}(*q*_{‖}, *q*_{z}) − *S*_{bg}(*q*_{‖}, *q*_{z}), 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, *S*_{l}(*q*) and *S*_{v}(*q*), functions of the 3D wavevector modulus $q=[q\Vert 2+qz2]1/2$ and combined as

with the nominal numbers of liquid *N*_{l} and vapor *N*_{v} molecules in the sampled volume *V*_{T}.

In the strict *q*_{z} = 0 limit, the ECWT gives Werthein’s relationship^{4,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 *S*_{bg} → *S*_{bg} + Δ*S*_{bg} would shift *K* by $\Delta K\u2248\beta \gamma o2\Delta Sbg/\Delta \rho 2$.

Within their *q*_{z} = 0 analysis, HD^{25} assumed a background $Sbg\u2009HD=Slv\u2009G+\delta S\u2009HD$, with Gibbs’s (G) plane to set $Nl\u2009G=N\u2009T\u2212Nv\u2009G=(N\u2009T\u2212V\u2009T\rho v)\rho l/\Delta \rho $, and with an interfacial (but not-CW) contribution *δS*^{HD}. 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* ≈ *T*_{t}, HD got a very flat *γ*^{HD}(*q*_{‖}), with marginally negative *K*^{HD}(*T*_{t}). At higher *T*, they got 0 < *βK*^{HD}(*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*_{‖} ≤ *q*_{u}, 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 *S*_{bg}(*q*_{‖}, 0).

In MD, we may get *γ*(*q*_{‖}) by a very different method. Instead using only the (experimentally accessible) pair correlations described by *S*_{T}(*q*_{‖}, *q*_{z}), 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 algorithms^{37–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 *S*_{bg}(*q*_{‖}, *q*_{z}) in the total structure factor *S*_{T}(*q*_{‖}, *q*_{z}) of surface diffraction experiments. We extend the analysis to *q*_{z} ≠ 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 results^{14,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 *r*_{c} = 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 ($Tc*>1.2$, with our *r*_{c}).^{42,43}

*N* = 176 440 particles in a box with transverse (*XY*) size *L* = 41.28*σ* and *L*_{z} = 6*L*, form thick liquid slabs $(\u223c3L)$ with the coexisting densities (*ρ*_{v}, *ρ*_{l}) and the surface tension (*γ*_{o}) given in Table I. Notice that our *r*_{c} 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 $(dt=510\u22123\sigma m/\u03f5)$ and taken after 2 × 10^{6} MD steps for thermal equilibrium, were used to get the density profile

and, for *q*_{x,y} = 2*πn*_{x,y}/*L* (with integer *n*_{x,y}, as set by the periodic boundaries) and any *q*_{z}, the structure factor

The index *i* in (6) and (7) runs over the *N*_{T} particles in the chosen region (*z*_{a} ≤ *z*_{i} ≤ *z*_{b}) of the MD box, with volume *V*_{T} = *L*^{2}(*z*_{b} − *z*_{a}), and the limits inside the vapor [*ρ*(*z*_{a}) = *ρ*_{v}] and liquid [*ρ*_{l}(*z*_{b}) = *ρ*_{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) *S*_{v}(*q*) and *S*_{l}(*q*), used in (4) as functions of the 3D wavevector *q*, are obtained with *q*_{z} = 0 (i.e., *q* = *q*_{‖}) but taking both *z*_{a} and *z*_{b} within the same bulk phase, and replacing *A* by *N*_{T} in the denominator of (7).

. | ISM . | HD . | This Work . | |||||||
---|---|---|---|---|---|---|---|---|---|---|

T*
. | ρ_{v}σ^{3}
. | ρ_{l}σ^{3}
. | γ_{0}σ^{2}/ϵ
. | n_{0}σ^{2}
. | K/ϵ
. | δn_{l}σ^{2}
. | K/ϵ
. | Δn_{l}σ^{2}
. | K/ϵ (q_{z} = 0)
. | K/ϵ (q_{z} ≠ 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}/ϵ
. | n_{0}σ^{2}
. | K/ϵ
. | δn_{l}σ^{2}
. | K/ϵ
. | Δn_{l}σ^{2}
. | K/ϵ (q_{z} = 0)
. | K/ϵ (q_{z} ≠ 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 $Sbg\u2009HD(q\Vert ,qz)$ (extended to *q*_{z} ≠ 0), we take (7) with both limits for *z*_{a} ≤ *z*_{i} ≤ *z*_{b} within the liquid (or vapor), but the restriction *z*_{a} ≤ *z*_{j} is applied also to the second particle in the pair. That corresponds to use the plane *z* = *z*_{a} to create the *virtual* (CW-less) interface assumed by Höefling and Dietrich. Equivalently, the surface excess background $\delta Sbg\u2009HD(q\Vert ,qz)$ is minus the result of (7) but taking only the pairs with *z*_{i} ≤ *z*_{a} and *z*_{j} > *z*_{a} for any *z*_{a} well within the bulk phase.

The main dependence of *S*_{T}(*q*_{‖}, *q*_{z}) on the shape of volume *V*_{T} goes into (4), to leave *H* ≡ *S*_{T} − *S*_{bg} independent of that choice. The sharp boundaries at *z* = *z*_{l,v} used in MD to define *V*_{T} could be changed into an exponential decay toward the liquid bulk^{25} to represent the adsorption of the x-ray beams in diffraction experiments. That would change *S*_{T} and *S*_{bg} 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 $P(\xi (x\u20d7))$, set by the mean position of the interface $\u27e8\xi (x\u20d7)\u27e9\u2261\xi \u0302o$ 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 *S*_{T}(*q*_{‖}, *q*_{z}). The assumption that the CW amplitude $\xi \u0302(q\Vert )$ is uncorrelated with the intrinsic density operator^{28} is very robust for low *q*_{‖}. The theory may be applied without paying attention to the precise definition for $\xi (x\u20d7)$ or to the molecular details in *ρ*_{in}(*z*). Moreover, the Gaussian fluctuations of $\xi \u0302(q\Vert )$ are well described by the macroscopic *γ*_{o}. The upper cutoff *q*_{‖} ≤ *q*_{u} 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}(*q*_{z}) = Δ*ρ* carrying no information on the liquid interface at molecular scale. That would imply from (9) a Gaussian Φ(*q*_{z}), 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 *q*_{u}*σ* ≈ 0.7; which may be used to accurately predict how *ρ*(*z*) and Φ(*q*_{z}) change when we increase the size of the MD box. The CW divergence in *S*_{T}(*q*_{‖}, 0) would be fairly well represented by Wertheim’s (5) with *γ*_{o}. However, nothing peculiar is observed in the MD results for *S*_{T}(*q*_{‖}, *q*_{z}) when *q*_{‖} goes across that (CWT) upper bound. That value of *q*_{u} 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 $z=\xi (x\u20d7)$ 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 $\xi \u0302(q\Vert )$ 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 *q*_{z} = 0, with Φ_{in}(0) = Φ(0) = Δ*ρ*. Its validity for *q*_{z} ≠ 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 *q*_{z} ≠ 0 theoretical prediction for the CW contribution to *S*_{T}(*q*_{‖}, *q*_{z}) goes beyond Wertheim’s (5), that Bedeaux and Weeks^{13,29} (BW) showed to be just the first term in a *q*_{z} power series,

with the height–height correlation, $S(|x\u20d7\u2212x\u20d7\u2032|)=\xi (x\u20d7)\xi (x\u20d7\u2032)\u2212\xi (x\u20d7)\xi (x\u20d7\u2032)$, for any two points on the surface plane, and $S(0)=\Delta \u2009CW2$ in (9). This *S*_{BW}(*q*_{‖}, *q*_{z}) is the theoretical (ECWT) prediction for *H*(*q*_{‖}, *q*_{z}) ≡ *S*_{T}(*q*_{‖}, *q*_{z}) − *S*_{bg}(*q*_{‖}, *q*_{z}), but we keep the different notation to remark that *H*(*q*_{‖}, *q*_{z}) reflects our choice for the non-CW correlation background, while *S*_{BW}(*q*_{‖}, *q*_{z}) is a strong theoretical prediction in which the joint dependence on *q*_{‖} and *q*_{z} comes from their separated dependence in *γ*(*q*_{‖}) and Φ(*q*_{z}). Our goal of *H*(*q*_{‖}, *q*_{z}) ≈ *S*_{BW}(*q*_{‖}, *q*_{z}), could only be obtained with consistent results for *γ*(*q*_{‖}) and *S*_{bg}(*q*_{‖}, *q*_{z}), that get the best for the ECWT assumptions. In contrast, if we restrict ourselves to the *q*_{z} = 0 case, Wertheim’s prediction (5) would give directly a function *γ*(*q*_{‖}) for any choice of *S*_{bg}(*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, $\rho \u2009in(z)=\u27e8\u2211i\delta (z\u2212zi+\xi (x\u20d7i))\u27e9/A$, and the function *γ*(*q*_{‖}) that quantifies the surface fluctuations.^{30} The ISM^{30,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 $z=\xi \u2009ISM(x\u20d7)$ as a Voronoi tessellation^{31} that Fourier transformed gives the CW amplitudes $\xi \u0302\u2009ISM(q\u20d7\Vert )$. 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, *n*_{0}, i.e., the integral of the first peak in the IDP (inset in Fig. 1).

With early insight, Stillinger^{44} 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 *n*_{0} to describe the *outermost liquid layer*. Other physical properties (kinetics, dynamics,^{31} and even heat transport^{45}) have confirmed that the molecular layering is indeed a main characteristic of liquid surfaces (at least for *T*/*T*_{c} ≲ 0.8). The ISM may find its clearest view through the optimal value of *n*_{0} that depend on the temperature and the molecular interactions.^{31}

Figure 2 presents *γ*(*q*_{‖}) at *T** = 0.763, over the range 0.5 ≤ *n*_{0} ≤ 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 *K*^{HD} that could be reached in ISM with *n*_{0} ≈ 1. That *n*_{0} is well above the ISM optimal choice, *n*_{0} = 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* *n*_{0} (and within the error bars associated with that choice) that we get the ISM prediction for the surface bending modulus *K*^{ISM}(*T**) given in Table I, always well above *K*^{HD}(*T**). Moreover, it is only with the optimal *n*_{0} 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*), $\rho \u2009in\u2009ISM(z)$ and *γ*^{ISM}(*q*_{‖}).

The MD snapshot in Fig. 3 (bottom) gives a pictorial idea on how the choice of *n*_{0} is reflected in the IS. The coarse shape of $z=\xi \u2009ISM(x\u20d7)$, (i.e., the surface undulations with longer wavelengths) is quite robust, similar for the two choices of *n*_{0} 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 *n*_{0} becomes relevant. From a single MD snapshot, we could hardly determine which value of *n*_{0} 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 *n*_{0} = 0.75 (leading to *K*^{ISM} = 0.28 ± 0.05) than with *n*_{0} = 1, as needed to recover *K*^{HD} ≈ 0.

What we propose here is to take the MD-ISM results for *γ*^{ISM}(*q*_{‖}) and Φ(*q*_{z}) into *S*_{BW}(*q*_{‖}, *q*_{z}), in (11). Then, we run in reverse the problem studied by Höfling and Dietrich, i.e., we get the background $Sbg\u2009ISM(q\Vert ,qz)\u2261S\u2009T(q\Vert ,qz)\u2212S\u2009BW(q\Vert ,qz)$ that is consistent with the molecular definition for $z=\xi \u2009ISM(x\u20d7)$. This exercise had been done^{28} for the pair correlation function *G*(*z*, *z*′, *q*_{‖}) and its corresponding BW series *G*_{BW}(*z*, *z*′, *q*_{‖}), with real space variables *z* and *z*′ that are Fourier transformed as *z* − *z*′ to get *S*_{T}(*q*_{‖}, *q*_{z}).^{46} The result for *G* − *G*_{BW} 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 $Sbg\u2009ISM(q\Vert ,qz)$ with HD proposal $Sbg\u2009HD(q\Vert ,qz)$, 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 *S*_{T}(*q*_{‖}, *q*_{z}) 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 $q\Vert \u2032=[q\Vert 2+qz2\u2061sin2\u2061\varphi ]1/2\u226bq\Vert $, then most of the CW contribution would be cancelled, $S\u2009BW(q\Vert \u2032,qz\u2032)\u226aS\u2009BW(q\Vert ,qz)$, while $qz\u2032=qz\u2061cos\u2061\varphi $ keeps unchanged the contribution of the bulk phases to *S*_{T}(*q*_{‖}, *q*_{z}) that depends on the 3D wavevector modulus *q* = *q*′.

At low *q*_{z}, the CW signal with $q\Vert \u2032$ may still be important and $\Delta S\u2009T(q\Vert ,qz;\varphi )\u2261S\u2009T(q\Vert ,qz)\u2212S\u2009T(q\Vert \u2032,qz\u2032)$ may be a poor approximation to *H*(*q*_{‖}, *q*_{z}). However, if *S*_{bg}(*q*_{‖}, *q*_{z}) depends mainly on the 3D *q*, Δ*S*_{T} should be similar to the equivalent difference between the BW predictions $\Delta S\u2009BW(q\Vert ,qz;\varphi )\u2261S\u2009BW(q\Vert ,qz)\u2212S\u2009BW(q\Vert \u2032,qz\u2032)$.

Figure 4 compares Δ*S*_{T}(*q*_{‖}, *q*_{z}; *ϕ*) and Δ*S*_{BW}(*q*_{‖}, *q*_{z}; *ϕ*) as functions of *q*_{z} 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 *q*_{z} > 0 accessible with the periodic boundaries for each *q*_{‖} is that which makes $q\Vert \u2032=2q\Vert $, reducing to a half the $1/q\Vert 2$ CW factor in (11). As *q*_{z} increases, $q\Vert \u2032$ grows and its CW contribution vanishes. Therefore, we expect a smooth trend from Δ*S*_{T} ≈ *H*/2 for the lowest *q*_{z}, to Δ*S*_{T} ≈ *H* for large *q*_{z}. The excellent agreement between the MD result for Δ*S*_{T} and the mesoscopic prediction Δ*S*_{BW} supports the accuracy of BW theory and the hypothesis that *S*_{bg}(*q*_{‖}, *q*_{z}) may be approximated by a function of the 3D wavevector *q*.

## IV. ISM-MD CONSISTENT BACKGROUND

Figure 5 presents the MD result (7) for *S*_{T}(*q*_{‖}, *q*_{z}) at *T** = 0.763, and the CW contribution *S*_{BW}(*q*_{‖}, *q*_{z}), built as (11) with *γ*^{ISM}(*q*_{‖}) and Φ(*q*_{z}). The data are presented for several *q*_{‖} ≥ 0.152 and as functions of the 3D *q*. For *q* ≲ 0.5, *S*_{T}(*q*_{‖}, *q*_{z}) reflects the larger CW contribution for smaller *q*_{‖}, but for larger *q*, the bulk-like contributions dominate and for *q* ≳ 3 all the (*q*_{‖}, *q*_{z}) points merge in a single line. Notice that as *q*_{‖} grows the decay of the CW contributions is dominated by |Φ(*q*_{z})|^{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 *N*_{T}/*A* ≈ 13*σ*^{−2} (about 15*σ* depth in both the liquid and vapor sides). Assuming that the CW contribution is represented by *H*(*q*_{‖}, *q*_{z}) = *S*_{BW}(*q*_{‖}, *q*_{z}), we get the difference $Sbg\u2009ISM(q\Vert ,qz)\u2261S\u2009T(q\Vert ,qz)\u2212S\u2009BW(q\Vert ,qz)$, 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*_{‖}, *q*_{z}) merged in a narrow band.

The inset in Fig. 5 shows a clear gap $Sbg\u2009ISM>Sbg\u2009HD$ between our ISM background and HD proposal. Notice that the original HD surface term $\delta Sbg\u2009HD(q\Vert ,0)$, that reflects a virtual (CW-less) interface, has to be extended to *q*_{z} ≠ 0 but still it may be obtained from the bulk structure factors, or directly from the MD sampling (7). Although strictly $\delta Sbg\u2009HD(q\Vert ,qz)$ 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 *N*_{l} with respect to Gibbs’s value $Nl\u2009G$ in *S*_{lv}(*q*), keeping fix *N*_{l} + *N*_{v} = *N*_{T}, with the values of $\delta nl\u2261(Nl\u2009G\u2212Nl)/A$ given in Table I. A larger change Δ*n*_{l} of that effective number of liquid-like particles may be used to close the gap between $Sbg\u2009ISM$ and $Sbg\u2009HD$. The value of the best fitting Δ*n*_{l} 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=Nl\u2009G$, Δ*n*_{l} = 0), and change the ISM parameter *n*_{0} for the LJ surface, to achieve $Sbg\u2009ISM\u2248Sbg\u2009HD$. However, with the insight of the *N*-particle structure in MD, we know that the mesoscopic view of the liquid surface given by $z=\xi \u2009ISM(x\u20d7)$ would be spoiled if we take *n*_{0} = 1, instead of its optimal value *n*_{0} ≈ 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 Δ*n*_{l} in the liquid-like contribution to the *S*_{lv}(*q*) background.

Figure 6 gives a more detailed comparison of *S*_{T}, $Sbg\u2009ISM$ and $Sbg\u2009HD$ for a lower temperature (*T** = 0.678). We use now the on-plane *q*_{‖} as variable for *q*_{z} = 0 (as HD) and for *q*_{z} = 1.0622. At low *q*_{‖}, when *S*_{T}(*q*_{‖}, *q*_{z}) and *S*_{BW}(*q*_{‖}, *q*_{z}) are large and nearly equal, we have to assume a larger error bar for their difference $Sbg\u2009ISM$, particularly for *q*_{z} = 0. However, the results clearly suggest changing HD background proposal allowing for a Δ*n*_{l} ≠ 0. Sticking to Gibbs’s plane to set *N*_{l} not only changes the non-CW background at low *q*_{‖} (therefore pushing *K*^{HD} well below *K*^{ISM}) but it fails to recover *S*_{bg}(*q*_{‖}, *q*_{z}) = *S*_{T}(*q*_{‖}, *q*_{z}) at large *q*_{‖} when the CW contribution has to become negligible. The difference *S*_{T}(*q*_{‖}, *q*_{z}) − *S*_{bg}(*q*_{‖}, *q*_{z}) ≈ 0.1 produces a sharp drop of *γ*(*q*_{‖}) for *q*_{‖} ≳ 2, as reported in HD results.^{25} The accuracy of our proposed fit $Sbg\u2009ISM(q\Vert ,qz)\u2248Slv(q)$ is not perfect, particularly for *q*_{z} ≠ 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 *S*_{lv}(*q*) (with the Δ*n*_{l} in Table I) gives an evident improvement over the original HD proposal. Keeping HD virtual surface term $\delta Sbg\u2009HD(q\Vert ,qz)$ in the fit makes too little difference to be worth, since it is absorbed by a small change in Δ*n*_{l}.

The supplementary material provides evidence that the description of *S*_{bg}(*q*_{‖}, *q*_{z}) ≈ *S*_{lv}(*q*), applied here to the LJ liquid surface, is also very accurate for a different (SA^{28}) 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 *S*_{T}(*q*_{‖}, *q*_{z}). The good news is that $Sbg\u2009ISM(q\Vert ,qz)$ is not far from the proposal of Höfling and Dietrich^{25} $Sbg\u2009HD(q\Vert ,qz)$, 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 *N*_{l} and *N*_{v}.

In surface diffraction experiments, the estimation of the non-CW background has to be done directly from the total *S*_{T}(*q*_{‖}, *q*_{z}). In the analysis of MD simulations, the use $\Delta nl\u2261(Nl\u2009G\u2212Nl)/A$ to fit $Slv(q)\u2248Sbg\u2009ISM(q\Vert ,qz)$, 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 *S*_{T} 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 *S*_{T}(*q*_{‖}, *q*_{z}) and that it comes to be similar to *γ*^{ISM}(*q*_{‖}). To achieve that, we take advantage of the rapid decay of the CW contribution *S*_{BW}(*q*_{‖}, *q*_{z}) with increasing *q*_{z}, so we fit *S*_{lv}(*q*) ≈ *S*_{T}(*q*_{‖}, *q*_{z}) over the 3D *q* range in which the CW contribution is negligible. Then, the best fitting value of Δ*n*_{l} is used to approximate the non-CW background *S*_{bg}(*q*_{‖}, *q*_{z}) ≈ *S*_{lv}(*q*) over the lower values of *q*_{‖} that give access to *γ*(*q*_{‖}).

Figure 7 presents the fit *S*_{T}(*q*_{‖}, *q*_{z}) ≈ *S*_{lv}(*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 *S*_{lv}(*q*) with Δ*n*_{l} = 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 *S*_{lv}(*q*) to represent *S*_{bg}(*q*_{‖}, *q*_{z}). Notice that to use still higher *q* for the fit gives less accuracy for Δ*n*_{l} because *S*_{lv}(*q*) becomes very steep.

Then, in our MD, we may use *q*_{z} = 0 Wertheim’s equation (5) to get *γ*(*q*_{‖}) from *H*(*q*_{‖}, 0) = *S*_{T}(*q*_{‖}, 0) − *S*_{lv}(*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 *n*_{0} 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* = *K*^{ISM} ± Δ*K*^{ISM} (as given in Table I), that bracket those obtained from *S*_{T}(*q*_{‖}, *q*_{z}). 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 *S*_{T}(*q*_{‖|}, *q*_{z}). The results with the original HD proposal give much lower *K*^{HD} ≈ 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.

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 $\u2113\u2261K/\gamma o$ 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 *T*_{c}, 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* = *S*_{T} − *S*_{bg}, has to be done essentially as for MD data in Sec. V, with a fit *S*_{T}(*q*_{‖}, *q*_{z}) ≈ *S*_{lv}(*q*) over *q*-range where the CW effects are negligible. We cannot use Gibbs’s plane to get $Nl\u2009G$ and $Nv\u2009G$ (taking into account the exponential decay from the X-rays adsorption^{25}) since the experimental reflectivity data do not give the location of *ρ*(*z*) relative to the signal in *S*_{T}. Moreover, the strict *q*_{z} = 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*_{‖}, *q*_{z}) = *S*_{BW}(*q*_{‖}, *q*_{z}) is not as simple as in the *q*_{z} = 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 $q\Vert 2$ order, leaving *K* as the only free parameter. Given *γ*_{o} and |Φ(*q*_{z})|^{2}, we look for the value of *K* that best fit *H*(*q*_{‖}, *q*_{z}) ≈ *S*_{BW}(*q*_{‖}, *q*_{z}) over the available range for *q*_{z}, with $\gamma (q\Vert )=\gamma o+Kq\Vert 2$. That would be our best estimation for *K*, extracted from Φ(*q*_{z}) and *S*_{T}(*q*_{‖}, *q*_{z}) at *q*_{z} > 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 *q*_{z} = 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 *q*_{z} ≠ 0. The use of Gibbs’s plane to determine *S*_{bg} (i.e., HD method extended to *q*_{z} ≠ 0) becomes even worse than for *q*_{z} = 0.

Under experimental conditions, the continuous limit *q*_{‖}*L* ≫ 1 may be used to simplify (12) as

where $\Sigma (x)\u22612\pi \beta \gamma o(S(0)\u2212S(x))$ behaves as log(*x*/*x*_{o}) for large distances, with a molecular-like length parameter *x*_{o} 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 *J*_{o}(*q*_{‖}*x*_{o}). 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 $\eta (q\Vert )=qz2/(2\pi \beta \gamma (q\Vert ))$, which carries all the dependence on *γ*(*q*_{‖}). For low *q*_{z}, *η*(*q*_{‖}) ≪ 1 and $qz2\Gamma (\eta (q\Vert )/2)\u22484\pi \beta \gamma (q\Vert )+O(qz2)$. Therefore, a fairly accurate and simple way to get *γ*(*q*_{‖}) out of *H* ≡ *S*_{T}(*q*_{‖}, *q*_{z}) − *S*_{lv}(*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 reported^{14,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 *S*_{l}(0) rather than keeping the dependence *S*_{l}(*q*) over the range used to extract *H*(*q*_{‖}, *q*_{z}) and *γ*(*q*_{‖}).

In Fig. 10, we mimic that procedure to interpret our MD results for *S*_{T}(*q*_{‖}, *q*_{z}) with a constant background $Sbg0$. 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 *S*_{T}(*q*_{‖}, *q*_{z}) with *q* ≳ 2, where the CW contribution may be considered negligible. Moreover, we expect $S\u2009T(q\Vert ,qz)\u2265Sbg0$, since otherwise the assumed CW signal $H0\u2261S\u2009T\u2212Sbg0$ would become negative. Therefore, in Fig. 10 (top), we fix a value of *q*_{‖} and take $Sbg0$ as the minimum of *S*_{T}(*q*_{‖}, *q*_{z}), which happens around *q* ∼ 1.5 − 2.5, in between the *q* = 0 compressibility peak and the main (*q* ≈ 2*π*) peak of *S*_{l}(*q*). Now, we may use that (assumed constant) background to obtain *H* and *γ*(*q*_{‖}), as we have done with *S*_{bg}(*q*) ≈ *S*_{lv}(*q*), or with $Sbg\u2009HD(q\Vert ,qz)$. We take the simplest *q*_{z} = 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 $Sbg(q)>Sbg0$. That peak (of bulk like correlations) is therefore taken as part of the CW contribution *H*^{0}. 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 $S\u2009T(q\Vert ,0)=Sbg0$, *γ*(*q*_{‖}) has to raise back [since *H*^{0}(*q*_{‖}, 0) goes to zero], producing a shape qualitatively similar to that reported as *enhanced* capillary waves in some experimental studies.^{14,15}

## 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*_{‖}, *q*_{z}) ≡ *S*_{T}(*q*_{‖}, *q*_{z}) − *S*_{bg}(*q*_{‖}, *q*_{z}), but to extract *γ*(*q*_{‖}) from *S*_{T}(*q*_{‖}, *q*_{z}) has been difficult and controversial.

As pointed by Höfling and Dietrich (HD),^{25} the problem is the non-CW *background* *S*_{bg}(*q*_{‖}, *q*_{z}), that gives a regular (non-diverging at *q*_{‖} = 0) but important contribution to *S*_{T}(*q*_{‖}, *q*_{z}). The thermodynamic limit *γ*(0) = *γ*_{o} is independent of *S*_{bg}, 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 *S*_{bg}.

HD proposed a background, $Sbg\u2009HD$, defined on physically sound terms, and they used it to evaluate a *γ*^{HD}(*q*_{‖}) (and *K*^{HD}) from MD simulations for the (truncated) LJ liquid surface. They took the *q*_{z} = 0 limit that simplifies the analysis, and Gibbs’s plane to set the *N*_{l} (liquid-like) and *N*_{v} (vapor-like) contributions of the particles with bulk-like correlations sampled within *S*_{T}(*q*_{‖}, *q*_{z}), 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*S*_{T}gives a molecular interpretation for the function*γ*(*q*_{‖}) extracted from x-ray surface diffraction experiments.Our study covers

*q*_{z}≠ 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

*N*_{l}and*N*_{v}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 *S*_{T} when the CW effects should vanish. That may be achieved with a minor (but crucial!) modification of *N*_{l} and *N*_{v} from the values given by Gibbs’s plane.^{44} Notice that HD setting $Nl=Nl\u2009G$ is a natural and reasonable choice, but not a physical requirement as strong as getting *S*_{T} = *S*_{bg} for large *q*_{‖}. The latter requirement may be achieved with a small change (just a fraction of a monolayer) in *N*_{l}.

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 *K*^{HD} 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 *K*^{ISM} are based in the definition of an intrinsic surface $z=\xi \u2009ISM(x\u20d7)$ 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 $z=\xi (x\u20d7)$ that we assign to a MD snapshot, or to the non-CW background that we use to extract *H*(*q*_{‖}, *q*_{z}) and *γ*(*q*_{‖}) from *S*_{T}(*q*_{‖}, *q*_{z}); 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*_{‖}, *q*_{z}) = *S*_{T}(*q*_{‖}, *q*_{z}) − *S*_{lv}(*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*_{‖}, *q*_{z}). The surface layering is produced by the short-ranged excluded volume interactions in the dense liquid and it represents precisely what the background *S*_{bg} has to take out of *S*_{T}. 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 *S*_{T}, as in general for the applications of the ISM, we cannot take the liquid too close to the critical temperature. For *T*/*T*_{c} ≳ 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 *S*_{bg}(*q*_{‖}, *q*_{z}) ≈ *S*_{lv}(*q*) (fitted with the Δ*n*_{l} parameter) goes much beyond the assumption of a constant background, *S*_{bg}(*q*_{‖}, *q*_{z}) ≈ *S*_{l}(0)*N*_{l}/*A*, which was used in the interpretation of x-ray diffraction experiments^{14,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/*r*^{6} tails in van der Waals interactions. A fit to *S*_{T}(*q*_{‖}, *q*_{z}) over the region of its minimum (*q*_{‖} ≈ 2 − 3) goes below the maximum of $Sbg\u2009ISM(q\Vert ,qz)$ 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 $\u223c1/r6$ 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 small^{49} and the non-analytic behavior induced in *γ*/*q*_{‖}) appears to be well below what could be extracted from *S*_{T}(*q*_{‖}, *q*_{z}).

Our computer codes for the ISM^{30,31} and for the present analysis of *S*_{T}(*q*_{‖}, *q*_{z}) 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 *S*_{bg}(*q*_{‖}, *q*_{z}) ≈ *S*_{lv}(*q*), applied here to the LJ liquid surface, is also very accurate for a different (SA^{28}) 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 *S*_{T}(*q*_{‖}, *q*_{z}).