Layering and capillary waves in the structure factor of liquid surfaces

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‖=[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 $|Φ(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)|², 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)|² 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² on the $x⃗=(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)|², 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 $ξ̂(q‖⃗)$⁠. The early CWT version^3–7 predicted the mean square amplitudes $⟨|ξ̂(q‖)|2⟩$ from the surface tension γ_o and temperature (kT = β⁻¹), 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¹³ a CW contribution to S_T(q_‖, q_z) proportional to the mean square amplitude $⟨|ξ̂(q⃗‖)|2⟩$⁠. Therefore, $γ(q‖)≡kT/[q‖2⟨|ξ̂(q⃗‖)|2⟩A]$⁠, 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 $ηo≡qz2/(2πβγo)<2$⁠, the CW contribution to S_T(q_‖, q_z) diverges^11,12,29 as $|Φ(q z)|2/(βγoq‖2−ηo)$⁠. For low q_z (i.e., at or near grazing incidence), this CW contribution is large at the wavevectors (q_‖ ∼ 10⁻³ nm⁻¹) 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 $ξ(x⃗)$ 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)²⁵ 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‖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 $ΔK≈βγo2ΔSbg/Δρ2$⁠.

Within their q_z = 0 analysis, HD²⁵ assumed a background $Sbg HD=Slv G+δS HD$⁠, with Gibbs’s (G) plane to set $Nl G=N T−Nv G=(N T−V Tρv)ρl/Δρ$⁠, 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.³⁵ 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

NVT-MD simulations were done with LAMPPS package,⁴¹ 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 = 6L, form thick liquid slabs $(∼3L)$ 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,²⁵ 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−3σm/ϵ)$ and taken after 2 × 10⁶ 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²(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.³⁵ 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).

To get HD background $Sbg HD(q‖,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 $δSbg HD(q‖,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²⁵ 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.²⁵ 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(ξ(x⃗))$⁠, set by the mean position of the interface $⟨ξ(x⃗)⟩≡ξ̂o$ 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 $ξ̂(q‖)$ is uncorrelated with the intrinsic density operator²⁸ is very robust for low q_‖. The theory may be applied without paying attention to the precise definition for $ξ(x⃗)$ or to the molecular details in ρ_in(z). Moreover, the Gaussian fluctuations of $ξ̂(q‖)$ 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=ξ(x⃗)$ 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 $ξ̂(q‖)$ 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⃗−x⃗′|)=ξ(x⃗)ξ(x⃗′)−ξ(x⃗)ξ(x⃗′)$⁠, for any two points on the surface plane, and $S(0)=Δ CW2$ 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.³⁴ Given a rule, the statistical average over MD configurations gives the IDP, $ρ in(z)=⟨∑iδ(z−zi+ξ(x⃗i))⟩/A$⁠, and the function γ(q_‖) that quantifies the surface fluctuations.³⁰ 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=ξ ISM(x⃗)$ as a Voronoi tessellation³¹ that Fourier transformed gives the CW amplitudes $ξ̂ ISM(q⃗‖)$⁠. 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₀, i.e., the integral of the first peak in the IDP (inset in Fig. 1).

With early insight, Stillinger⁴⁴ 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₀ to describe the outermost liquid layer. Other physical properties (kinetics, dynamics,³¹ and even heat transport⁴⁵) 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₀ that depend on the temperature and the molecular interactions.³¹ 

Figure 2 presents γ(q_‖) at T* = 0.763, over the range 0.5 ≤ n₀ ≤ 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₀ ≈ 1. That n₀ is well above the ISM optimal choice, n₀ = 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₀ (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₀ 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), $ρ in ISM(z)$ and γ^ISM(q_‖).

The MD snapshot in Fig. 3 (bottom) gives a pictorial idea on how the choice of n₀ is reflected in the IS. The coarse shape of $z=ξ ISM(x⃗)$⁠, (i.e., the surface undulations with longer wavelengths) is quite robust, similar for the two choices of n₀ 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₀ becomes relevant. From a single MD snapshot, we could hardly determine which value of n₀ 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.75 (leading to K^ISM = 0.28 ± 0.05) than with n₀ = 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 ISM(q‖,qz)≡S T(q‖,qz)−S BW(q‖,qz)$ that is consistent with the molecular definition for $z=ξ ISM(x⃗)$⁠. This exercise had been done²⁸ 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).⁴⁶ 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,²⁸ 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 ISM(q‖,qz)$ with HD proposal $Sbg HD(q‖,qz)$⁠, and it lead to a molecular interpretation of the function γ(q_‖) that is extracted from surface diffraction experiments.

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‖′=[q‖2+qz2⁡sin2⁡ϕ]1/2≫q‖$⁠, then most of the CW contribution would be cancelled, $S BW(q‖′,qz′)≪S BW(q‖,qz)$⁠, while $qz′=qz⁡cos⁡ϕ$ 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‖′$ may still be important and $ΔS T(q‖,qz;ϕ)≡S T(q‖,qz)−S T(q‖′,qz′)$ 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 $ΔS BW(q‖,qz;ϕ)≡S BW(q‖,qz)−S BW(q‖′,qz′)$⁠.

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‖′=2q‖$⁠, reducing to a half the $1/q‖2$ CW factor in (11). As q_z increases, $q‖′$ 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.

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)|² 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σ⁻² (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 ISM(q‖,qz)≡S T(q‖,qz)−S BW(q‖,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 ISM>Sbg HD$ between our ISM background and HD proposal. Notice that the original HD surface term $δSbg HD(q‖,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 $δSbg HD(q‖,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 G$ in S_lv(q), keeping fix N_l + N_v = N_T, with the values of $δnl≡(Nl G−Nl)/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 ISM$ and $Sbg HD$⁠. 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 G$⁠, Δn_l = 0), and change the ISM parameter n₀ for the LJ surface, to achieve $Sbg ISM≈Sbg HD$⁠. However, with the insight of the N-particle structure in MD, we know that the mesoscopic view of the liquid surface given by $z=ξ ISM(x⃗)$ would be spoiled if we take n₀ = 1, instead of its optimal value n₀ ≈ 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 ISM$ and $Sbg HD$ 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 ISM$⁠, 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.²⁵ The accuracy of our proposed fit $Sbg ISM(q‖,qz)≈Slv(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 $δSbg HD(q‖,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²⁸) molecular interaction model.

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 ISM(q‖,qz)$ is not far from the proposal of Höfling and Dietrich²⁵ $Sbg HD(q‖,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 $Δnl≡(Nl G−Nl)/A$ to fit $Slv(q)≈Sbg ISM(q‖,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₀ 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.³¹ 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;²⁵ 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,²⁵ who compared the length $ℓ≡K/γ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.

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 G$ and $Nv G$ (taking into account the exponential decay from the X-rays adsorption²⁵) 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‖2$ order, leaving K as the only free parameter. Given γ_o and |Φ(q_z)|², 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 $γ(q‖)=γo+Kq‖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 $Σ(x)≡2πβγo(S(0)−S(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⁻¹ ≪ q_‖ ≪ σ⁻¹, a still simpler approximation to (13) may be obtained as

with the gamma function and the exponent $η(q‖)=qz2/(2πβγ(q‖))$⁠, which carries all the dependence on γ(q_‖). For low q_z, η(q_‖) ≪ 1 and $qz2Γ(η(q‖)/2)≈4πβγ(q‖)+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_‖)].⁴⁷ 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⁻¹ ≪ q_‖ ≪ σ⁻¹ 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⁻¹, 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 T(q‖,qz)≥Sbg0$⁠, since otherwise the assumed CW signal $H0≡S T−Sbg0$ 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 HD(q‖,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⁰. 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 T(q‖,0)=Sbg0$⁠, γ(q_‖) has to raise back [since H⁰(q_‖, 0) goes to zero], producing a shape qualitatively similar to that reported as enhanced capillary waves in some experimental studies.^14,15

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.⁵ 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),²⁵ 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 HD$⁠, 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.⁴⁴ Notice that HD setting $Nl=Nl G$ 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=ξ ISM(x⃗)$ that optimizes the intrinsic layering structure at the liquid surface, as it had been foreseen by Stillinger long time ago.⁴⁴ 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=ξ(x⃗)$ 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⁶ 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 ISM(q‖,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 $∼1/r6$ potential tail,⁴⁸ observed the dispersion force effects in γ^ISM(q_‖), as generically predicted by Napiórkowski and Dietrich,⁸ but its quantitative effect is small⁴⁹ 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.

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²⁸) molecular interaction model. In light of this result, we discuss enhanced CW [γ(q_‖) < γ_o].

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).

The authors have no conflicts to disclose.

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).

The data that support the findings of this study are available from the corresponding author upon reasonable request.