Structure of liquid–vapor interfaces: perspectives from liquid state theory, large-scale simulations, and potential grazing-incidence X-ray diffraction

Grazing-incidence X-ray diffraction (GIXRD) is a scattering technique which allows one to characterize the structure of fluid interfaces down to the molecular scale, including the measurement of the surface tension and of the interface roughness. However, the corresponding standard data analysis at non-zero wave numbers has been criticized as to be inconclusive because the scattering intensity is polluted by the unavoidable scattering from the bulk. Here we overcome this ambiguity by proposing a physically consistent model of the bulk contribution which is based on a minimal set of assumptions of experimental relevance. To this end, we derive an explicit integral expression for the background scattering, which can be determined numerically from the static structure factors of the coexisting bulk phases as independent input. Concerning the interpretation of GIXRD data inferred from computer simulations, we extend the model to account also for the finite sizes of the bulk phases, which are unavoidable in simulations. The corresponding leading-order correction beyond the dominant contribution to the scattered intensity is revealed by asymptotic analysis, which is characterized by the competition between the linear system size and the X-ray penetration depth in the case of simulations. Specifically, we have calculated the expected GIXRD intensity for scattering at the planar liquid–vapor interface of Lennard-Jones (LJ) fluids with truncated pair interactions via extensive, high-precision computer simulations. The reported data cover interfacial and bulk properties of fluid states along the whole liquid–vapor coexistence line. A sensitivity analysis shows that our findings are robust

Experimentally, scattering techniques such as X-ray reflectometry and grazing-incidence X-ray diffraction (GIXRD) a) Electronic mail: f.hoefling@fu-berlin.de provide the most detailed view of fluid interfaces.Whereas reflectometry allows one to infer interfacial density profiles, GIXRD probes density fluctuations [56][57][58] and surface structures 29,30 .This kind of scattering experiments under grazing incidence have been carried out for liquid surfaces of water 36 , molecular fluids 37 , and liquid metals 38 .The analysis of the scattered intensities requires accounting for the inevitable scattering from the fluid bulk phases, which depends on a suitable background model [59][60][61] , nowadays available within GIXRD analysis software 62 ; yet, such a modeling choice can potentially introduce ambiguities into the interpretation of the experimental data.
Focusing on thermal equilibria of coexisting liquid and vapor phases, the increase of the interfacial area by a macroscopic amount ∆A incurs a free energy cost of γ 0 ∆A, which defines the (macroscopic) surface tension γ 0 .At the molecular scale, the interface is roughened by thermal fluctuations, which are tamed by an analogous cost in free energy.This mesoscopic picture leads to the capillary wave (CW) theory [63][64][65] , which assumes the interface to be a smooth, membrane-like surface with a tension modulus which is governed by an interface Hamiltonian.On the other hand, by adopting the particle perspective, the theory of inhomogeneous fluids 64,66,67 considers (two-point) density fluctuations in order to characterize the interfacial region.For a free, planar interface, the classical result by Wertheim 68 and Weeks 69 states that fluctuations with wavevector q parallel to the surface lead to a scattering intensity proportional to k B T/(γ 0 q 2 ), which applies for small Figure 1.Top: cross-section of a snapshot of a simulated liquid-vapor interface configuration for a truncated Lennard-Jones (LJ) fluid at liquid-vapor bulk coexistence for the temperature T ≈ 0.8T c ; T c is the critical temperature of this fluid.Bottom: interface-related fluctuations are switched off in a Gedankenexperiment (see the main text).This snapshot depicts the initial configuration of the simulation, composed of two independent bulk configurations, before further equilibration.q = |q|, i.e., for macroscopic wavelengths 2π/q [see, cf., Fig. 8 and Eq. ( 66)]; the thermal energy scale is denoted by k B T as usual.This asymptotic result (q → 0) matches with the free energy cost of the excitation of corresponding CWs.It has stimulated elaborate derivations of effective interfacial Hamiltonians [33][34][35][70][71][72] in order to capture also the (anticipated) structure of the two-point correlation function at higher orders in q. These xtended CW theories include also a bending modulus κ of the interface 43,44 or, more generally, a wave number-dependent surface tension γ(q) which exhibits the relation γ(q → 0) = γ 0 .The leading correction to this macroscopic limit is either non-analytic (of the form q 2 log(q)) or quadratic (κq 2 ), depending on whether dispersion forces are present or not.Molecular dynamics (MD) simulations testing the behavior of γ(q) have frequently been analyzed in line with the CW picture by evaluating the fluctuation spectrum of the local interface position 20,[40][41][42][43][44][45][46][73][74][75] .This requires one to adopt a suitable definition of this position down to the molecular scale.Common choices are either a local Gibbs dividing surface or an elaborate algorithm to identify an intrinsic interface which separates CWs from bulk-like density fluctuations.However, in view of the various possible choices the published data for γ(q) disagree on even the sign of the bending coefficient κ (in the absence of dispersion forces).Moreover, it was shown within density functional theory (DFT) 50 that any effective interface Hamiltonian fails to reproduce the behavior of interfacial density correlations, which are obtained from simulations.
As an alternative which bypasses the ambiguities associated with the definition of a local interface position, we previously introduced an effective surface tension γ(q) as a proxy for interfacial density fluctuations which is entirely based on quantities amenable to scattering experiments 49 [see below, Eq. ( 8)].Similarly, as for the interpretation of scattering data, this approach hinges on a consistent model for the background scattering, which is the subject of the present study.It is based on a hypothetical liquid-vapor interface with all interface-related correlations switched off (Fig. 1), for which the scattered intensity can be worked out analytically.Somewhat unexpectedly, these open boundary conditions at the interface lead to corrections in the small-wave number density correlations due to bulk fluctuations 76 , which, as it will turn out, interfere with the small-q behavior of γ(q); a related boundary effect has been observed for the two-dimensional strip geometry 77,78 .Accounting for these corrections removes an inconsistency between the results for γ(q) obtained from the direct calculation of the density correlations and obtained from simulated scattering intensities.Moreover, it renders the bending coefficient in γ(q) positive for all temperatures at liquid-vapor coexistence, between the triple point and the critical point of the fluid.
The paper is organized as follows.In Sect.II we provide a number of relations which will be useful for the theoretical analysis; in particular, we discuss the separation of the GIXRD intensity into interfacial and background parts, provide the definition of γ(q), and connect with Ref. 76.The background scattering for the above mentioned reference system is analyzed in Sect.III with special emphasis on the non-commuting limits of infinite X-ray penetration depth and infinite sample width.In Sect.IV, these results are applied to the analysis of simulation data for scattering intensities from liquid-vapor interfaces.

A. GIXRD master formula
For diffraction experiments on planar interfaces under grazing incidence, the scattering intensity is proportional to the two-point correlation function of the atomic number density.The dependence of the scattering intensity on the lateral scattering vector q = (q x , q y ) follows the master formula [Eq.(2.68)  in Ref. 56] with q = |q| and the z-axis coinciding with the interface normal.
Here, we have dropped the atomic form factors for reasons of simplicity (i.e., implicitly considering point-like particles) and have omitted the reflection and transmission coefficients, which are constant amplitudes with respect to the lateral momentum transfer. 79The two-point number density correlation function G(q, z, z ′ ) characterizes the sample and is independent of the experimental setup probing it.Within the theory of fluids, this density-density correlation function is defined as 64,67 where ρ(r) is the microscopic number density and ϱ(r) = ⟨ ρ(r)⟩ is the mean density at point r = (R, z).It is convenient to split off the singular self-part of G(r, r ′ ) and to introduce the pair correlation function g(r, r ′ ) by In planar geometry, translational invariance parallel to the mean interface entails ϱ(r) = ϱ(z) and G(r, r ′ ) = G(R − R ′ , z, z ′ ).This suggests to consider the lateral Fourier transform.Accordingly, for the two-component vectors q and ∆R = R − R ′ , one introduces which enters into Eq.(1).G(q, z, z ′ ) depends only on the modulus q = |q| of the wavevector due to the isotropy of the sample in the planes parallel to the interface.The function f (z) is determined by the scattering geometry (i.e., the angles α i and α f of incoming and outgoing beams, respectively; see Fig. 2 in Ref. 56) and the mean density profile ϱ(z).It describes the decay of the evanescent wave on the liquid side of the interface, f (z < 0) ∼ exp(−κ|z|), with the penetration depth 1/κ(α i , α f ) 56,80 .For α i and α f approaching the critical angle of total reflection, the penetration depth diverges, and, in the thermodynamic limit, I(q) is dominated by scattering from the liquid bulk.Thus, experimental setups aim at minimizing the penetration depth 1/κ, which in practice can reach down to a few nanometres 36,37 or even below 81,82 .On the other hand, one has to ensure that the penetration depth 1/κ is considerably larger than the interfacial width ζ (cf.Eq. ( 50)); otherwise, the observation of interfacial fluctuations would be incomplete.Therefore, the condition to access interfacial properties fully is κζ ≪ 1, which is met by the typical experimental setups; in analytical and simulation work, it implies that the undamped limit, κ → 0, is considered.

B. Capillary wave divergence
Concerning the q-dependence of the scattering intensity, classical capillary wave theory predicts, for the density correlations at small wave numbers, that asymptotically 64,68,69 where γ 0 denotes the macroscopic surface tension and ℓ c = √︀ γ 0 /(mg∆ϱ) is the capillary length; ∆ϱ = ϱ ℓ − ϱ v is the number density contrast between the coexisting phases, and mg is the gravitational weight of a molecule of mass m.Clearly, G(q, z, z ′ ) diverges as q −2 for large capillary lengths, ℓ c ≫ q −1 , which is considered in the following.In combination with Eq. ( 1), the divergence is passed on to the scattered intensity: Since the bulk scattering remains finite as q → 0, the interfacial contribution dominates the signal for small q.In order to obtain information about the interface which goes beyond this divergence, one needs to unambiguously separate interfacerelated and bulk contributions to the scattering intensity: This leads one to introduce the interfacial structure factor as H(q) := lim κ→0 [I(q; κ) − I b (q; κ)] = I int (q; κ = 0) (7)   and an effective, wave number-dependent surface tension γ(q) by generalizing Eq. ( 5): The macroscopic surface tension is recovered as γ 0 = γ(q → 0) and the q-dependence of γ(q) quantifies the deviations from the divergence H(q → 0) ∼ q −2 , which can be attributed to interface-related density fluctuations.The division in Eq. ( 6) hinges on a consistent model for the bulk scattering I b (q; κ) which is based on clearly formulated assumptions and on experimentally accessible quantities only.A simple and commonly used model is based on the bulk structure factor of the liquid phase 36,37 : where χ T is the isothermal compressibility.We shall show below that this approximation, even in the small-q regime, creates a constant shift in the interfacial structure factor H(q) and thus a bias of O(q 2 ) in γ(q).Moreover, it leads to inconsistencies between experiments and simulations.

C. Structure factor of an open slab of liquid
In Ref. 76, we investigated the effect of open boundaries of a liquid sample of finite width L < ∞.It revealed a finite-size correction to the structure factor and anticipates the solution strategy followed in Sect.III.This study is based on virtually excavating a planar slab of width L from a homogeneous liquid, thereby imposing free boundary conditions on the newly created surfaces (see Fig. 1 of Ref. 76; the corresponding situation with equilibrated interfaces is depicted below in Fig. 4).The observable of interest is the slab structure factor, which, up to the prefactor, resembles Eq. ( 1) for a step function f (z) = 1 [0,L] (z), chosen as the indicator function of the interval [0, L]; the latter emerges in the limit κ → 0 and for a finite system.Using the two-point correlation function of a homogeneous liquid, G(q, z, z ′ ) = G ℓ (q, ∆z = z − z ′ ), yields ϱLS (q; L) = I ℓ (q; κ → 0, L) [see Eqs. ( 1) and ( 6)].The crucial step in the derivation is the observation that for a homogeneous fluid, the pair correlation function is fully determined by the bulk structure factor.This is expressed by the relation which is obtained by a partial Fourier back transform with k = (−q, k z ) denoting a three-component wave vector.The subscripts ℓ indicate that ϱ ℓ and S ℓ refer to the number density and the structure factor of the bulk liquid, respectively.Inserting G ℓ (|q|, ∆z) into Eq.( 10) and carrying out the integrals over z and z ′ yields 76 For wide slabs, one finds the asymptotic expansion whereJ 0 (q) is the leading-order correction integral Most importantly, the small-wave number limit J 0 (q → 0) is non-zero and can have either sign, depending on, e.g., the temperature (see below, Fig. 7).We stress that this kind of finite-size corrections does not appear for periodic boundary conditions at the surfaces z = 0 and z = L (as commonly used in simulations).In the latter case, the equation S per (q; L) = S ℓ (q) holds exactly.

A. Density-density correlation function
The interface-related contributions to the scattering intensity are unambiguously identified via the comparison with a background reference system in which all interface-induced perturbations are switched off.This idealized situation can be created by constructing an ensemble of particle configurations, which contain a planar interface and which are characterized solely by bulk correlations.Such configurations are obtained in a Gedankenexperiment by virtually cutting a homogeneous liquid sample of macroscopic size along a plane, denoted as z = 0, and by removing all molecules above the plane (i.e., with positions z > 0).The empty half space is then filled by a correspondingly treated sample of the coexisting vapor phase.This renders a flat liquid-vapor interface without any structural distortions in the vicinity of the plane z = 0 (Fig. 1).By construction, the local densities on opposite sides of the interface are independent of each other and, in particular, uncorrelated.On the other hand, for z and z ′ both being on the same side, the two-point correlation function of the reference system equals that of the respective bulk phase, no matter how close to the interface z and z ′ are chosen to be.Thus, we define the twopoint correlation function of this "naked interface" ensemble, which serves as background reference, as It describes an inhomogeneous system of two unperturbed and uncorrelated, coexisting bulk phases which share a flat, open boundary; due to the spatial homogeneity of the bulk phases, G ℓ and G v depend only on z − z ′ instead on z and z ′ separately.It is clear that the dynamic evolution of these reference configurations would necessarily lead to the usual equilibrium with fully developed interfacial fluctuations; however, statistical averages of the background configurations are understood as ensemble averages, not time averages.The bulk scattering follows by inserting this definition of a reference system [i.e., Eq. ( 15)] into the master formula, [Eq.(1)] which is carried out in Sect.III B. Concerning simulations of GIXRD experiments, the same model for the bulk correlations can be used with the only modification that the nonzero width of the liquid film must be taken into account.Accordingly, the scattering depends on both the penetration depth 1/κ and the sample width L with the relevant, non-commuting limits κ → 0 and L → ∞ (see Sect.III C).
We note that G b (q, z, z ′ ) as defined in Eq. ( 15) is discontinuous at the interface as the reference system changes abruptly between the two phases.On the other hand, the physically observed, full two-point correlation G(q, z, z ′ ) is a continuous function of z and z ′ for all wave numbers q.Any smooth interpolation of the background between the coexisting phases would necessarily make assumptions on the correlations in the interfacial region, e.g., it would introduce an unknown interpolation length. 83A discontinuity of G b (q, z, z ′ ) implies that the interfacial part G int (q, z, z ′ ) := G(q, z, z ′ ) − G b (q, z, z ′ ) is also discontinuous.We emphasize that this conceptually unavoidable discontinuity of G b and thus G int is not in conflict with theoretical constraints 51 ; most importantly, it does not contradict the asymptotically rigorous result for the CW divergence [Eq.( 4) with l c = 0]: which is continuous in z and z ′ since the mean density profile ϱ(z) is a smooth function.The same property carries over to G int (q, z, z ′ ), because G b (q, z, z ′ ) is a bounded function of q: It turns out that the continuity of G int (q, z, z ′ ) is not needed for this to hold; it is sufficient that the discontinuity is uniformly bounded in q (as a function of z and z ′ ).The line of arguments is as follows: Let us separate G int = G c int + G ∆ int into its continuous part G c int (q, z, z ′ ) and its discontinuous part G ∆ int (q, z, z ′ ), which is piecewise constant in z and z ′ .The latter is also the discontinuous part of G b , and, because G ℓ and G v are bounded, G ∆ int (q, z, z ′ ) is thus uniformly bounded in q.With this, in the limit q → 0, the discontinuous part drops out: B. Scattering from a macroscopic half-space of bulk liquid In this section, we turn to the specific setup of GIXRD experiments, i.e., a macroscopic liquid sample (L → ∞) and an evanescent wave on the liquid side (κ > 0).According to the master formula in Eq. ( 1), the GIXRD intensity for the bulk reference model [Eq.(15)] is the sum of two independent contributions stemming from the integrals over G b (q, z, z ′ ) in the domains z, z ′ < 0 (bulk liquid) and z, z ′ > 0 (bulk vapor).At low vapor pressure, the scattering from the latter side of the interface is negligible; alternatively, the result of Ref. 76 can be used directly because there is no damping of the propagating beam, i.e., f (z > 0) = 1.
The calculation of the bulk scattering from the liquid side (z < 0) is conceptually similar, albeit the algebraic expressions differ.Due to the exponential decay of the evanescent wave, i.e., f (z) = exp(−κ|z|) for z < 0, the integrals in Eq. ( 1) can be formulated explicitly: upon a change of variables u := z − z ′ = ∆z and v := −(z + z ′ ) and by taking into account the absolute value of the Jacobian The integral over v is simple, and replacing the bulk correlation function G ℓ (q, u) by Eq. ( 11) yields Finally, interchanging the integrations over u and k z leads to For a given static structure factor S ℓ (k) and wave numbers k < k max , the integral can be evaluated numerically as described in Ref. 76.To this end, we use the identity which can be truncated safely at large k z , recalling that S ℓ (k → ∞) = 1.As a by-product, one finds I ℓ (q → ∞) = ϱ ℓ /(2κ).
In the limit κ → 0, i.e., for scattering angles close to the angle of total reflection, one can identify a part of the integrand in Eq. ( 21) as a representation of Dirac's δ-distribution: which holds for a continuous and bounded test function φ(x).I ℓ (|q|) resembles the bulk structure factor evaluated at wave vectors q parallel to the interface: The next-to-leading order term O (︀ κ 0 )︀ within the asymptotic expansion of I ℓ (q) around κ = 0 is given by lim where we have taken the limit κ → 0 inside of the integral, as permitted by the theorem on monotone convergence, and we have made use of the definition ofJ 0 (q) [Eq.( 14)].Combining Eqs. ( 24) and (25), we arrive at one of our main results: The behavior of the bulk contribution I ℓ (q) is presented exemplarily in Fig. 2 for a Lennard-Jones (LJ) liquid along the liquid-vapor coexistence line at two temperatures: T * := k B T/ε = 0.70 slightly above the triple point temperature and T * = 1.15 in proximity of the critical temperature (T * c ≈ 1.22).The pair potential was truncated at r c = 3.5σ, and ε and σ refer to the interaction strength and range of the LJ potential, respectively; the phase diagram is shown in Fig. 3 and further details of the simulations are given in Appendix A. We calculated I ℓ (q) according to Eq. ( 22) from the bulk structure factors S ℓ (k) as input, using the same simulation data for S ℓ (k) as the ones in Ref. 76.The integration over k z in Eq. ( 22) was restricted to 0 ⩽ k z ⩽ k max = 50/σ, and S ℓ (k) was extrapolated to such large wave numbers as described in Ref. 76.The curves obtained for a range of penetration depths (0.01 ⩽ κσ ⩽ 1) convincingly corroborate the convergence of (2κ/ϱ ℓ ) I ℓ (q) to S ℓ (q) as κ → 0. The corrections, however, are significant for κσ ≳ 0.1: from Fig. 2 one infers that the value of I ℓ (q → 0) is increased over the value of S ℓ (q → 0) by a factor of up to ≈ 3.5 for T * = 0.70, which is close to the triple point, whereas it is suppressed by a factor of about 1.8 at the higher temperature (T * = 1.15).Moreover, the minimum in I ℓ (q) near qσ ≈ 2 becomes more shallow for increasing κ and seems to disappear at low temperatures.

C. GIXRD intensity from a liquid slab of finite width
For MD simulations of GIXRD intensities, the finite size of the system and, in particular, the finite width L of the bulk liquid must be accounted for properly (Fig. 4).One anticipates that the appropriate expression for the bulk scattering combines aspects both of Sects.II C and III B. We start the derivation of the expression for the bulk scattering by reiterating that the reference system is such that liquid and vapor regions are The normalization is taken such that in the limit κ → 0 the curves stay finite and approach S ℓ (q) [Eq.( 24 independent of each other and that there are no distortions of the microscopic density close to z = −L ℓ , 0, L v , where free boundary conditions are imposed.The finite widths of liquid and vapor slabs, respectively, are implemented as cutoffs in the weight function: Concerning the derivation of a formula for the background scattering, we again specialize to the liquid side (z < 0).The bulk contribution of the vapor side follows in the limit κ → 0 as discussed previously 76 (see Sect.II C).
Following the same route as before, we combine the GIXRD master formula [Eq.( 1)] for the two-point correlation function G ℓ (q, ∆z) of the homogeneous bulk with the truncated form of f (z) [Eq.(27)] which yields the finite-slab (L ℓ < ∞) version of Eq. ( 19): after carrying out the elementary integral over v in the first line.In the next step, we use Eq.(11) in order to substitute G ℓ (q, u) by the bulk structure factor S ℓ (k) and interchange the integrations over k z and u.Based on the integral we arrive at our central result for a finite system: the last line contains the contribution for q = 0, and its prefactor stems from the integral In the final step, we make use of the identity and recast Eq. ( 30) in a form similar to Eq. ( 22): where the 1 in the numerator of the integrand is balanced by the last term.With this, the remaining integration over k z is approximated well by a finite integration domain |k z | < k max , because S (k) ≃ 1 for k sufficiently large.Thus, Eq. ( 33) is suitable for a numerical evaluation; in particular, increasing k max decreases the truncation error.We have followed the procedure described in Ref. 76, which has already been used to integrate Eq. ( 22), choosing k max = 50/σ as before. 84t is elucidating to discuss certain limiting cases of Eq. ( 30).For large wave number q, one reads off from Eq. ( 33), using that S ℓ (q → ∞) = 1: In Eqs. ( 30) and ( 33), care is needed when taking the limits κ → 0 and L ℓ → ∞ at finite values of q.For κ fixed, the limit L ℓ → ∞ and the integration over k z may be interchanged due to the dominated convergence theorem (using that the numerator is bounded), which allows one to reproduce the previous results in Eqs. ( 21) and ( 24): The limit κ → 0 for L ℓ fixed is more intricate.The integrand in Eq. ( 30) is dominated by the integrable function ] for all κ ⩾ 0, where M is the maximum of S ℓ (•).[Proof: Put y = e −κL ℓ and a = cos(k z L ℓ ) and use the fact that the expression 1 + y 2 − 2ya is monotonically increasing for y ⩾ a, i.e., it is maximal at y = 1.]Therefore, the limit κ → 0 may be taken inside the integral so that where x = k z L ℓ .Except for the prefactor ϱ ℓ L ℓ , this is precisely the structure factor of a liquid slab [see Eq. ( 12)], and only for large L ℓ it approaches the bulk structure factor: This result is particularly relevant for simulations, because for them finite-size effects cannot be avoided, and because a finite value of L ℓ permits to probe the physically meaningful limit κ → 0 directly.In Fig. 5, we test the expression in Eq. ( 30) for the GIXRD intensity scattered from finite-sized liquid samples at low (T * = 0.70) and high (T * = 1.15) temperatures.As for Fig. 2, we have employed the bulk structure factors S ℓ (q) obtained from MD simulations.Alternatively, we have calculated I ℓ (q) within the simulations and via Eq.( 1), with the integration domain restricted to 0 ⩽ z, z ′ ⩽ L ℓ .These simulation data are in excellent agreement with the results from Eq. ( 30) for all parameter sets.The top row of the panels [(a),(b)] shows the convergence of I ℓ (q) upon increasing the sample width L → ∞ at a fixed, exemplary penetration depth κσ = 0.01, whereas in the bottom row [(c),(d)] the limit κ → 0 is taken at the fixed sample width L ℓ = 20σ.We note that in the latter case (κ → 0, L ℓ fixed), the scattered intensity does not converge to the bulk structure factor S ℓ (q) (black lines), but rather to the slab structure factor S (q; L ℓ ) (not shown), as expected from Eq. ( 37); the difference is O(L −1 ℓ ) and vanishes for macroscopic samples [see Eq. ( 13)].
It is noteworthy that the bulk contribution I ℓ (q) to the scattered intensity, including the case κ = 0, is non-additive: On the other hand, the expression ϱ ℓ L ℓ S ℓ (q) is trivially L ℓadditive, implying that it does not contain correlations in the transversal direction, i.e., between two particles at positions z z ′ .However, such correlations are contained in the scattered intensity.One anticipates an appreciable error in the interfacial structure factor H(q) [Eq.(7)] at small yet non-zero wave numbers q if the bulk contribution I ℓ (q) is approximated by I ℓ (q) ≈ ϱ ℓ L ℓ S ℓ (q) -which would correspond to assuming periodic boundary conditions at the interface (see the last paragraph of Sect.II C).A quantification of this error in H(q) follows directly from the expansion of S (q; L ℓ ) for large L ℓ [Eq.( 13)], which yields lim An analogous error in H(q) arises in the analysis of experimental data if only the leading order of the bulk scattering is subtracted from the scattered intensity [see Eq. ( 25)].However, the order of the limits κ → 0 and L ℓ → ∞ is different in the two cases, which results in the prefactor 2 on the r.h.s. of Eq. ( 39) relative to Eq. ( 25).The latter is understood by noting that κL ℓ in Eq. ( 30) is sent either to 0 (here) or to ∞ [Eq.( 25)].

D. Uniform asymptotic behavior
For the interpretation of both experiments and simulations, the asymptotic behavior of the bulk scattering for small (albeit non-zero) κ and large (finite or infinite) L ℓ is relevant.Accordingly, the issue arises whether Eq. ( 30) can be recast such that the leading asymptotic behavior is apparent from a single expression for both orders of the limits κ → 0 and L → ∞.Inspired by Eq. ( 34) and by the asymptotes of I ℓ (q) [Eqs.(35)  and ( 37)], we single out a term proportional to S ℓ (q) by using the identity in Eq. (32) and by rearranging the remaining integral: With the definition ofJ 0 (q) in Eq. ( 14), this reduces to The first term on the r.h.s. is the leading order for κ → 0 and L ℓ → ∞.One recovers both I ℓ (q) ≃ S ℓ (q)/(2κ) and I ℓ (q) ≃ ϱ ℓ L ℓ S ℓ (q), depending on the order of the limits.The second term is of O(1) w.r.t.κ → 0, L ℓ → ∞; it containsJ 0 (q), which is independent of κ and L ℓ and is given solely by S ℓ (q).The remaining two terms are higher-order corrections.For the third term, in the limit κ → 0 one finds Here, we have made use of Eq. ( 23) and of the expansion We note that S ′ ℓ (q)/q < ∞ for all q, including q → 0. Thus, the second line of Eq. ( 40) vanishes for κ → 0 at L ℓ fixed.
The term in the last line vanishes exponentially if the limit L ℓ → ∞ is taken first at fixed κ.For finite L ℓ , the remaining integral is bounded for all κ ⩾ 0 [e.g., by J 0 (q)] and it  is the Fourier cosine transform of an integrable and continuous function and thus decays as L −1 ℓ for L ℓ → ∞ [see also Eq. ( 13)].The limit κ → 0 can be interchanged with the integration over k z .Therefore, the last line of Eq. ( 40) behaves as e In summary, the asymptotic behavior of I ℓ (q) in the joint limit L ℓ → ∞ and κ → 0 reads: Taking L ℓ → ∞ first, one obtains for large penetration depths (i.e., κ → 0) Conversely, for L ℓ large but fixed and for κ → 0 one has The next-to-leading order terms of both expansions depend neither on κ nor on L ℓ , but differ by a factor of 2. Approximating the background intensity by the leading order only, as widely done in the analysis of experimental or simulation data (and also in theoretical treatments), leaves a contribution proportional to ϱ ℓ J 0 (q), but with a different prefactor, in the expression for the interface correlations.This introduces a systematic inconsistency when comparing results from experiments and simulations.(In the former results, unavoidably one has κ > 0 and an extrapolation κ → 0 is needed; in the latter results, one may put κ = 0 from the very beginning.)In particular, we shall find below [Eq.(55)] that this contribution pollutes the effective surface tension γ(q) at the order q 2 for small q.

IV. GIXRD FROM LIQUID-VAPOR INTERFACES
For GIXRD scattering from liquid-vapor interfaces, we consider the inhomogeneous reference system as described in Sect.III A, which is composed of two independent bulk phases at the coexisting densities and which exhibits the absence of correlations across the interface and unperturbed bulk correlations even close to the interface [see Eq. ( 15)].The evanescent wave on the liquid side is combined with a propagating wave on the vapor side so that in the master formula for the GIXRD intensity [Eq.( 1)].(Here, we do not consider the weak absorption of X-rays.).Evaluating the master formula with G(q, z, z ′ ) and f (z) as in Eqs. ( 15) and ( 46), respectively, yields the background intensity due to bulk scattering, with the liquid and the vapor contributions given by Eq. ( 21) and Eq. ( 36), respectively.The vapor contribution becomes particularly relevant at elevated temperatures, i.e., upon approaching the critical point.We also note that the limit κ → 0 of I v (q) is meaningful only for a finite width L v < ∞ of the vapor phase, as naturally encountered in simulations.For plotting purposes, we quote the large-q limit, which follows from Eqs. ( 34) and ( 47): Also for the CW divergence of G(q, z, z ′ ) in the classical Wertheim-Weeks theory, the presence of the function f (z) in the integrand of the GIXRD master formula implies a κdependent shift of the scattering amplitude.This is seen by combining Eqs. ( 1) and ( 4) and then inserting Eq. ( 46): where the integral can be expressed in terms of the digamma function ψ 0 (x) = Γ ′ (x)/Γ(x): Except in the close vicinity of the liquid-vapor critical point, it holds 1/κ ≫ ζ and one may approximate the second line by 1 − 1.39 κζ + O(κζ) 2 , which can be a noticeable correction to Eq. ( 5) at elevated temperatures; for the last result, we used ψ 0 (1) − ψ 0 (1/2) ≈ 1.39.
Let us examine the implications for γ(q) upon replacing the background intensity by the leading term, i.e., if only the divergent part of the bulk scattering is subtracted as it has been done in previous experimental studies [36][37][38] .To this end, we introduce an approximation ̃︀ H for the interface structure factor: and the corresponding surface tension, Equation (26) implies that ̃︀ H(q) = H(q) + ϱ ℓ J 0 (q).We emphasize that this latter q-dependent shift persists after taking the limit κ → 0. Thus, including the full expression for the background scattering results in a relative change of γ(q) given by The magnitude of this shift of the contribution O(q 2 ) to γ(q) is determined by two lengths: ϱ ℓ γ 0 /[(∆ϱ) 2 k B T ] and J 0 (q → 0), both of which are amenable to experimental investigations.As an example, we have estimated these quantities for liquid water at T = 27 ∘ C, where ∆ϱ ≈ ϱ ℓ , and have obtained γ 0 /(ϱ ℓ k B T ) ≈ 5.2 Å and J 0 (q → 0) ≈ 0.22 Å from evaluating Eq. ( 14) for available data of the bulk structure factor 85,86 .
In the following, we shall elucidate this shift further for simple theoretical models, which support the simulation results for LJ fluids reported in Ref. 49 and below in Sect.IV D.

B. Simple DFT models: square-gradient approximation
The Ornstein-Zernike form Table I.Interfacial properties of truncated LJ fluids (r c = 3.5σ) along the liquid-vapor coexistence line.The surface tension γ 0 was calculated from the anisotropy of the microscopic stress tensor.The length ℓ quantifies the contribution O(q 2 ) to γ(q) and was obtained from fits of Eq. ( 66) to the data in Fig. 10(a).The interfacial width ζ was determined from fits of Eq. (B1) to the density profile of inhomogeneous systems with a mean interfacial area of (100σ) 2 .The numbers in parentheses give the statistical uncertainty in the last digit.
of the structure factor is a common feature of densityfunctional theories based on the square-gradient approximation; here S 0 = ϱ ℓ k B T χ T , in terms of the isothermal compressibility χ T , and ξ is the OZ correlation length.S OZ (k) is a useful approximation for the structure factor of liquids within the range kξ ≪ 1 and at elevated temperatures.(This is valid for liquids with an appreciable compressibility, but not too close to their critical point.)For the integral corresponding to the leading-order correction [Eq.( 14)] one finds 76 Inserting this into Eq.( 26) yields, for the small-angle bulk scattering, This decreases the uncorrected, q-dependent surface tension γ(q) [Eq.( 55)]: The prefactor of the contribution O(q 2 ) has the dimension of a length squared and, upon approaching the critical point for d bulk dimensions and by using the exponent relation 87 2β + γ = νd.Thus, the correction to γ(q) scales as (qξ) 2 in square-gradient models.Noting that the macroscopic surface tension γ 0 = γ(q → 0) vanishes as 88 γ 0 ∼ |t| 2ν (see Fig. 6 and Table I), one finds that for T ↑ T c at O(q 2 ) the difference between γ(q) and γ(q) becomes independent of temperature.The relevance of the correction in Eq. ( 55) can be assessed by the comparison with the small-q behavior of γ(q).To this end, we consider a simplified DFT treatment of the liquidvapor interface based on the square-gradient approximation  I).The solid line depicts the critical scaling law γ 0 ≃ A γ |t| 2ν for t := (T − T c )/T c upon t ↑ 0 with the Ising exponent ν = 0.630.The inset shows the same data in a rectification plot of the critical law, yielding the critical temperature T * c = 1.220 ± 0.001 and the amplitude A γ = (2.554± 0.008)εσ −2 from a linear regression to the three data points for T * ⩾ 1.0.
with the symmetric double-parabola potential 50 .By defining γ(q) via the scattered intensity, for κ → 0 Parry et al. 50obtained γ(q → 0)/γ 0 ≃ 1+ 5 4 (qξ) 2 .Accounting for the correction of the bulk scattering due to the open boundary [Eq.( 59)] subtracts the contribution 1 4 (qξ) 2 , i.e., it decreases the contribution O(q 2 ) by 20%, leading to γ(q → 0)/γ 0 ≃ 1 + (qξ) 2 .We note that adopting a different definition of the q-dependent surface tension which is solely based on the two-point density correlation function G(q, z, z ′ ) would yield a different prefactor at order O(q 2 ), and such a definition would not require as a prerequisite a model for the bulk scattering as input 50 .However, up to date, experimentally G(q, z, z ′ ) is not directly accessible.

C. Hard-sphere approximation for the bulk liquid
At temperatures close to the triple point, the liquid phase is characterized by a low compressibility and a small correlation length; in particular, the bulk structure is dominated by the interparticle repulsion.In order to estimate the value ofJ 0 (q → 0) in this regime, we approximate the liquid bulk structure factor by the one of the Ashcroft-Lekner (AL) model for hard spheres 89 ; for simplicity, we ignore the contribution due to the attractive part of the pair potential 90 .We recall that hard spheres do not form liquid or vapor phases; accordingly there is no liquid-vapor interface.Nonetheless, one can expect that the AL model for S ℓ (k) renders a useful estimate of the bulk correctionJ 0 (q) for dense and nearly incompressible liquids.
In terms of the volume packing fraction η and the hard sphere diameter σ, the AL model reads  26)] as obtained from quadratures of Eq. ( 14) with bulk structure factors S ℓ (k) taken from (i) the square-gradient DFT [Eq.( 57 the direct correlation function and the coefficients Numerical integration of Eq. ( 14) for q = 0 yields the values of J 0 (0), which vary smoothly as function of the packing fraction, attaining their maximum ≈ 0.087 σ near η ≈ 0.2; at high packing fractionJ 0 (0) ≈ 0.070 σ for η = 0.45.The latter result should be compared withJ 0 (0) ≈ 0.063σ obtained from MD simulations for a LJ liquid at T * = 0.70 (see below for details).Moreover, using S AL (k) as input to Eq. ( 14) and numerically computing the full q-dependence ofJ 0 (q) yields a remarkably accurate approximation ofJ 0 (q) with S ℓ (k) obtained from the simulations (Fig. 7); for this comparison we used η = 0.45 and identified σ with σ.In particular, our analysis for the OZ and the AL model has revealed that J 0 (0) changes sign due to a subtle competition of excluded volume and long-ranged correlations.

D. Molecular dynamics simulations of Lennard-Jones fluids
Within MD simulations, we have numerically determined GIXRD intensities due to scattering off the liquid-vapor interface, based on Eq. ( 1) and a microscopic expression for G(q, z, z ′ ); details are given in Appendix A. For the calculation of γ(q), one has to account for the finite width of the liquid and vapor regions; close to T c , also the vapor contribution to the bulk scattering must not be neglected.On the other hand, the finite sizes of the bulk phases ensure that both I(q; κ) and I b (q; κ) remain finite as κ → 0 so that one can put κ = 0 already when calculating I(q; κ).In this case, the dimensionless quantity is given by the standard microscopic expression for the static structure factor [cf. Eq. (A2) for f (z) = 1], and the bulk contribution follows from the simpler expression in Eq. ( 12), which was derived in Ref. 76.In Eq. ( 63), A is the area of the planar mean interface and N is the number of particles in the simulation.We have followed both routes in order to test their consistency: (i) determine I(q; κ) and thus γ(q; κ) for a decreasing sequence of values of κ and take γ(q) = γ(q; κ → 0), and (ii) calculate S tot (q) and thus γ(q) directly.In Fig. 8, the decomposition of the scattered intensity into interfacial and background contributions [Eq.( 6)] is illustrated (a) for the temperature T * = 0.70 close to the triple point and (b) for the temperature T * = 1.15 in proximity of the liquid-vapor critical point (T * c ≈ 1.22), both for a penetration depth 1/κ = 10σ.(The analogous decomposition of S tot (q) for T * = 1.15 is provided by Fig. 2 in Ref. 49.)For T * = 0.70, a liquid slab of width L s = 25σ was simulated, but only particles with positions 0 ⩽ z ⩽ L ℓ = 20σ were admitted for the calculation of I(q; κ); here, z = 0 denotes the mean position of the interface.For the higher temperature, these values have been doubled in order to accommodate the much lower (macroscopic) surface tension and thus larger fluctuations of the local interface position.On the vapor side, only particles within slabs of L v = 50σ (T * = 0.70) and L v = 75σ (T * = 1.15) were taken into account.The background scattering was calculated according to Eq. ( 33) using the bulk structure factors S ℓ (q) and S v (q) as input, which were obtained from separate simulations of homogeneous fluids (see Fig. 5).Close to the triple point [Fig.8(a)], the CW divergence is clearly visible in the scattered intensity, I(q → 0; κ) ∼ q −2 , without adjusting any parameters.To this end, the prefactor of the divergence was taken from Eq. ( 51) using the macroscopic surface tension γ 0 (Fig. 6 and Table I) as determined independently from an integral over the stress tensor profile across the interface 64,91 ; the interfacial width ζ was obtained from the simulated mean density profiles ϱ(z) [see Eq. (B1)].For the interfacial part of the scattering I int (q; κ) ≈ H(q), this behavior of the CW divergence extends to a wide range of wave numbers qσ ≲ 2. At high temperature [Fig.8(b)], the CW divergence is barely visible in I(q; κ) itself, but it can clearly be recognized in I int (q; κ) for qσ ≲ 0.2.The slight mismatch between the predicted and the actual prefactors of the CW divergence (compare the thin and the thick black lines for qσ ≲ 0.2) disappears for larger penetration of the liquid side, e.g., κσ = 0.01.The mismatch is presumably due to higher order terms in Eq. ( 4); the apparently obvious cause, that large-amplitude CWs are not probed properly for insufficiently small κ, is already accounted for in Eq. ( 51).The sizable deviations of I int (q; κ) from the asymptotic behavior

GIXRD intensity
T D 1:15 Ä D 0:1 L `D 40 I.q/ I int .q/I b .q/S v .q/S `.q/(b) Figure 8. Simulated I(q; κ) = I int (q; κ) + I b (q; κ) (red line) from scattering off the corresponding liquid-vapor interface and its decomposition into bulk and interface contributions I b (q; κ) and I int (q; κ), respectively (gray dashed and thick black lines); all quantities shown are normalized by I b (q → ∞) [Eq.(48)].The two panels show results at the reduced temperatures (a) T * = 0.70 and (b) T * = 1.15, both for the penetration depth 1/κ = 10σ.The thin black line indicates the CW divergence, I(q → 0; κ) ∼ 1/(γ 0 q 2 ) [Eq. ( 51)], with the macroscopic surface tension γ 0 obtained independently from the simulated stress tensor [Fig.6 and Table I]; the deviation of I int (q; κ) from the CW divergence (thick vs. thin black lines), which is well developed in panel (b), gives rise to the wave number-dependent surface tension γ(q).The bulk contribution I b (q; κ) (gray dashed line) was calculated according to Eqs. ( 33) and ( 47), using the simulated bulk structure factors of the coexisting liquid (green line) and vapor (blue) phases S ℓ (q) and S v (q), respectively; the data are extrapolated to small q assuming the Ornstein-Zernike form (green and blue dotted lines).For the calculation of I(q; κ) only particles within a slab of width (a) L ℓ = 20σ and (b) L ℓ = 40σ, respectively, were considered on the liquid side.For the quantities I(q), S ℓ (q), and S v (q), lines connect actual simulation data (symbols, only shown for the five smallest wave numbers).In panel (b), the tiny wiggles in I int (q) at qσ ≳ 1 reflect the statistical uncertainty of the simulation data.at large wave numbers give rise to the q-dependent surface tension.
Figure 9 exhibits, for each wave number q, the convergence of γ(q; κ) to the physically meaningful limit γ(q) upon systematically decreasing the inverse penetration depth κ → 0. In simulations, γ(q) can be obtained also directly from S tot (q) (thick solid lines).However, if only the leading order of the background scattering is subtracted, the resulting γ(q) deviates from γ(q) = γ(q; κ → 0) as one would infer from the GIXRD data.In particular, from Eq. ( 13) one concludes that ̂︀ H(q) = H(q) + 2ϱ ℓ J 0 (q), which also differs from ̃︀ H(q) [see Eq. ( 53) and the following text].Thus, the discrepancy between γ(q) and γ(q) is larger for larger wave numbers [cf.Eq. ( 55)]; in the simulations, e.g., for T * = 0.70 and qσ ≈ 2.2, we find that γ(q) is almost 40% smaller than γ(q).
It was this inconsistency between the analysis of scattering data (as obtained from experiments) and the total structure factor (as obtained within MD simulations) which gave rise to the refined treatment of the bulk background as elaborated here.The inconsistency would not be lifted by considering only the divergent part of the background in the analysis of GIXRD intensities I(q; κ), which upon κ → 0 would yield γ(q) as introduced after Eq. ( 53).In retrospective, the issue is well understood by comparing Eqs. ( 44) and (45), which imply that γ(q) γ(q).This clarifies that the O(1)-term ∝ ϱ ℓ J 0 (q) must be included in the background contribution for a consistent ana- .Convergence of the effective q-dependent surface tension γ(q; κ → 0) as function of the inverse penetration depth κ for three, fixed wave numbers q at the temperature T * = 0.70.The data for γ(q; κ) were obtained via Eq.( 52) from MD simulations for the GIXRD intensity, using only a a fraction L ℓ = 20σ of the width L s of the bulk liquid in order to calculate I(q) and carrying out the decomposition shown in Fig. 8(a).Thick solid lines indicate the corresponding limits γ(q) calculated from S tot (q), with setting κ = 0 directly within the simulations.Dashed lines show γ(q) as obtained from S tot (q) when accounting only for the divergent part of the bulk scattering [see the main text and Eq. ( 64)].The thin blue line (qσ/2π = 0.1) represents the κ-dependence of γ(q → 0; κ) as implied by Eq. ( 51).
lysis of GIXRD data.This yields γ(q) irrespective of whether it was calculated along route (i) or route (ii), described at the end of the first paragraph in Sect.IV D. The q-dependent surface tension γ(q) of LJ fluids is shown in Fig. 10(a) for temperatures ranging from T * = 0.70 to 1.15.The data are taken from Ref. 49 and were obtained from the MD simulation results for S tot (q) [route (ii)] and from the full background contribution given in Eq. ( 37) and based on Eq. ( 12): The most notable effect is the enhancement of the effective surface tension at non-zero wave numbers upon increasing the temperature; this was discussed in Ref. 49.It can be rationalized by writing with a temperature-dependent length ℓ(T ), which increases monotonically as T is increased (Table I).This form of γ(q) is corroborated within a recent DFT treatment of liquid-vapor interfaces [52][53][54] , with ℓ(T ) ≃ ξ(T ) at temperatures close to the critical one.
A related observation was made for the curvaturedependence of the macroscopic surface tension, upon replacing 1/q by the radius of a spherical droplet 92 .The simulation results for ℓ(T ) exhibit a temperature dependence similar to the one of the OZ bulk correlation lengths of the coexisting liquid and vapor phases (Fig. 11).Moreover, anticipating the same critical scaling exponent as for the correlation length, ℓ(T ↑ T c ) ∼ |T − T c | −ν , the product γ 0 ℓ 2 is expected to converge to a constant.Interestingly, our data suggest that γ 0 ℓ 2 /(k B T ) → (4πω) −1 ≈ 0.09, where ω = lim T ↑T c k B T c /(4πγ 0 ξ 2 ℓ,v ) is a universal amplitude ratio [92][93][94] ; its most reliable estimate stems from Monte Carlo simulations of the three-dimensional Ising model 94 : ω ≈ 0.87.This would imply that indeed ℓ(T )/ξ ℓ,v (T ) → 1 as T ↑ T c .
The right panel of Fig. 10(b) shows the corresponding results for γ(q), which have been obtained from the same input data, but taking into account only the divergent part of the bulk scattering [Eq.(64)].(We recall that J 0 (0) changes sign as function of temperature, which is seen in Fig. 7).Whereas the data points shift slightly upwards at higher temperatures, as expected from Eq. ( 55) due toJ 0 (0) < 0 at these temperatures [Eq.( 57)], the repercussions are more significant at low temperatures (T * ≲ 0.80): opposed to the almost constant behavior of γ(q) up to qσ ≲ 2, γ(q) bends downwards, which results from J 0 (0) ≈ 0.063σ > 0 in this case and due to the small value of ℓ.Empirically, the data are described by γ(q) ≈ γ 0 (1 + Kq α ) with K < 0 and exponents α = 4 for T * = 0.80 and α = 2.5 for T * = 0.70.In particular, the data for γ(q) suggest that there is a distinguished temperature T 0 with 0.80 ≲ T * 0 ≲ 0.90 such that ℓ(T ) = 0 for T < T 0 , which appears to be implausible on physical grounds.This issue is removed by considering the full background scattering [Eq.( 65)], which includes the correction given byJ 0 (q) and which leads to γ(q) as shown in Fig. 10(a).

E. Sensitivity to the mean interface position
So far, for the interpretation of the simulation data, we have anticipated that the widths L ℓ and L v of the coexisting phases [Eq.( 27)] are known.Our protocol to construct the equilibrated inhomogeneous samples (see Appendix A) suggests a fixed ratio L v : L ℓ (e.g., 3:1 or 4:1) of the bulk phases, where we set L ℓ = L s and have chosen L s = 25σ or 50σ for the width of the pre-equilibrated slabs, depending on temperature.However, due to the broadening of the interface by capillary waves, these nominal values of L ℓ and L v are in general (slightly) different from the values that are deduced from the inhomogeneous sample.For the latter step, various definitions of the mean interface position were proposed and are used in the literature 40,42,51,64,65,74 .
A common choice is based on Gibbs' dividing surface (GDS), which in integral form is equivalent to is fixed to the extent of the simulation domain along the interface normal (i.e., the z-axis).Combining these two relations yields L ℓ for given N, A, ϱ ℓ , ϱ v , and L z .On the other hand, counting the particles in the inhomogeneous system, which was assembled from pre-equilibrated bulk slabs (Appendix A), yields the same expression for N/A as the GDS criterion.Hence, L ℓ according to the GDS definition agrees with the nominal values for L ℓ .
In the present analysis, we followed a different approach and exploited the fact that there are two well-separated interfaces in the simulation setup: we have determined L ℓ from fits to the simulated mean density profile ϱ sim (z) using an inflected sigmoidal function [Eq.(B1)].The obtained values for L ℓ are very close to the nominal values according to the GDS definition; the absolute deviations are less than 0.1σ at all considered temperatures, i.e., a difference of less than 4‰.In addition to the values for L ℓ , the fits produced precise values of the coexisting densities, ϱ ℓ and ϱ v , and the interfacial width ζ (see Tables I  and II), which allowed us to obtain an accurate estimate of the liquid-vapor critical point (see Fig. 3 and Appendix B), although it is not in our focus here.
The small ambiguity in the position of the mean interface has consequences for the q-dependent surface tension γ(q): A variation of the interface position by δL implies changing the widths of the bulk slabs from L ℓ and L v to L ℓ + δL and L v − δL, respectively.(δL can have either sign.)This modifies the background contribution, given by Eq. ( 65) in the context of the simulations, and thus the interfacial structure factor H(q), which, by virtue of Eq. ( 13), receives an additive contribution provided that L ℓ ≫ ξ ℓ .This means that γ(q), given by Eq. ( 8), is replaced by the adjusted expression For small wave numbers, one has δH(q) ≪ H(q) [Fig.8] so Figure 10.MD simulation results (symbols) for the wave number-dependent surface tension γ(q) of the truncated LJ liquid as obtained from S tot (q) after subtracting (a) the full background contribution [Eq.(65)] and (b) only its divergent part [Eq.( 64)]; the data points in panel (a) are reproduced from Ref. 49.The differences are most apparent at the two lowest temperatures.In order to facilitate the comparison, the solid lines are the same in both panels: for T * ⩾ 0.9, the lines show fits of the small-q behavior as in Eq. ( 66) to the data for γ(q) (panel (a)); for T * = 0.7 and 0.8, they represent empirical power law fits to the data for γ(q) (panel (b)).  2 ] in comparison to the bulk correlation lengths, ξ ℓ and ξ v , of the coexisting liquid and vapor phases (Table I).The inset tests the convergence γ 0 ℓ 2 /(k B T ) → (4πω) −1 ≈ 0.09 as T → T c , where ω ≈ 0.87 is a universal amplitude ratio [92][93][94] .that where T ] is a certain length.Inserting the values for these coefficients, as obtained in the simulations (Table II), yields Λ ≈ 0.074σ for T * = 0.70 and Λ ≈ 0.046σ for T * = 1.15.With the corresponding values of ℓ [Table I and Fig. 11] and assuming a physically meaningful range for δL, we conclude that the precise definition of the mean interface position has only a minor effect on the behavior of γ(q) for small q.
For large wave number, however, we have δH(q → ∞) = ∆ϱ δL.Thus, changing H(q) by such an amount has the potential to qualitatively modify the behavior of γ(q) for large q.An adjustment of L ℓ implies that a contribution ∝ 1/q 2 is added reciprocally to γ(q): At large q, taking δL < 0 leads to an increase of γ adj (q) relative to γ(q) and, for δL sufficiently large in magnitude, this can induce a bending upwards of γ adj (q) such that γ adj (q → ∞) → ∞, which is desirable in the context of CW theory employing effective surface Hamiltonians 34,35,46,71,74 .Recently, Hernández-Muñoz, Tarazona, and Chacón 75 argued that there are no interfacial density correlations at large wave number and they proposed to use the condition H(q → ∞) → 0 in order to tune certain parameters of the data analysis.In our case, this amounts to adjusting the length L ℓ by δL = −H(q → ∞)/∆ϱ, which would remove a putative 1/q 2 -contribution from γ(q) and one would indeed obtain that γ adj (q → ∞) → ∞.
For γ(q), we have tested this procedure for the simulation results for γ(q) shown in Fig. 10(a).The data appear to follow a decay γ(q) ≃ γ 0 (bq) −2 for large q, albeit only in a small wave number window; such a decay for large q would correspond to a non-zero limit of the interfacial structure factor, H(q → ∞) > 0 [see Eq. ( 8) and Fig. 8].From fits to the data for γ(q) in the range 2.6 ≲ qσ ≲ 3.8, we have obtained b = 0.41σ and 0.21σ at T * = 0.70 and 1.15, respectively.The assumed asymptotic form of γ(q) for large q is equivalent to H(q → ∞) ≃ k B T (∆ϱ) 2 b 2 /γ 0 =: H ∞ , and setting δH(q → ∞) = −H ∞ will remove such a spurious large-q contribution from H(q).Robustness of the q-dependent surface tension with respect to a variation of the width L ℓ of the liquid slab in the background contribution S b (q) by an amount δL.Results for this, a posteriori adjusted quantity γ adj (q) [Eqs.(67) and ( 68)], are shown for seven distinct values of δL and for the temperatures T * = 0.70 (disks) and T * = 1.15 (diamonds).These results are based on the data provided for γ(q) in Fig. 10(a).
The corresponding shift δL = H ∞ /∆ϱ for the two temperatures renders δL = 0.11σ and δL = 0.30σ, respectively.However, shifting the mean interface position by such an amount, has only a marginal influence on γ(q) and does not yield the desired qualitative change, namely that γ(q) bends upwards [Fig.12].Instead, γ(q), over the entire accessible range of wave numbers, is found to be rather robust against variations of L ℓ within a physically meaningful range.
The reason that the above choice for δL does not remove the apparent large-q decay of γ(q) can be understood by noting that the wave number window, within which one has γ(q) ≈ γ 0 (bq) −2 , is yet to the left of the first peak of the bulk structure (which is near qσ ≈ 6.8, see Fig. 8).There, S ℓ (q) ≈ S v (q) ≈ 1 does not hold in this regime, which was used to deduce the relation δL = H ∞ /∆ϱ.Therefore, in order to remove an apparent plateau in H(q) around a certain intermediate wave number q * , one has to consider the full qdependence of δH(q) given in Eq. ( 67), which suggests to set δL = −H(q * )/[ϱ ℓ S ℓ (q * ) − ϱ v S v (q * )].At T * = 0.70, reasonable estimates of S v (q * ) and S ℓ (q * ) are given by their values for q → 0. Using H(q * ) = H ∞ one finds δL = −b 2 /Λ in terms of the length Λ [introduced after Eq. ( 69)].This expression leads to δL ≈ 2.3σ, which is more than four times the interfacial width ζ and thus physically unplausible.
Based on recent insight into the resonance structure of the interfacial two-point correlations, Parry and Rascón 53,54,55 have proposed that the full wave number dependence of γ 0 /γ(q) is well approximated by a linear combination of the bulk structure factors S ℓ (q)/S ℓ (q → 0) and S v (q)/S v (q → 0) with suitable, weakly q-dependent weights to account for the liquidvapor asymmetry.Moreover, these DFT studies reveal that H(q → ∞) ∼ q −2 (with the exception of the overly simplified square-gradient models, for which S (q → ∞) ∼ q −2 and thus H(q → ∞) ∼ q −4 ).Concerning the present MD simulation data, we conclude that a finding of H(q → ∞) = H ∞ > 0 would indeed be in conflict with the above prediction.However, the observed decrease of γ(q) corresponds well with the increase of S ℓ (q) in the rising flank of its first peak (near qσ ≈ 4, see, e.g., Fig. 8).Moreover, the actual behavior of γ(q) for large wave number, i.e., qσ ≳ 5 cannot be obtained from the data due to the unavoidable statistical noise.Thus, from our data one cannot rule out that the actual γ(q) has a small, positive limit γ ∞ := γ(q → ∞) > 0 or, equivalently, that H(q → ∞) ∼ q −2 -which would be consistent with the DFT calculations.

V. SUMMARY AND CONCLUSIONS
In sum, we have discussed the wave number-dependence of the GIXRD intensity I(q; κ) due to scattering off liquid-vapor interfaces.We have proposed an unambiguous separation I(q; κ) = I int (q; κ) + I b (q; κ) into an interface-related contribution I int (q; κ) and the bulk background I b (q; κ), as illustrated in Fig. 8; κ is the inverse penetration length.The separation is based on a simple reference system for the coexisting bulk phases which avoids any assumption concerning the interfacial region.The essential ingredients are free boundary conditions for the bulk phases on both sides of the interface.This means that the reference system is composed of independent liquid and vapor phases and that their structures are identical to the respective bulk structures and are unperturbed by the presence of the interface.(Necessarily, such an idealized situation cannot occur in thermal equilibrium, but only on paper, becauseç apillary waves and other interfacial fluctuations would render any physical quantity to vary smoothly across the interface.)Accepting this simple reference model, it turns out that the background scattering I ℓ (q; κ) from, e.g., the liquid phase is not simply proportional to the structure factor S ℓ (q) of the bulk liquid [Eq.( 26)]; rather it is given as an integral over this function [Eq.(21) and Fig. 2].This is a consequence of the free boundary conditions and appears likewise in the static structure factor of a liquid slab of finite width 76 [Eqs.(12) and ( 13)].We note that any "continuous" model for the background scattering, i.e., one which imposes a continuous interpolation of G b (q, z, z ′ ) across the interface, would require knowledge of the microscopic structure of the interfacial region, at least on the length scale over which the interpolation takes place.Already at the level of the local mean density, the question how to switch between the different correlation lengths on the liquid and the vapor sides has no obvious answer without providing microscopic details.As was shown in Sect.III A, the discontinuity of the background contribution and thus of the interfacial part of the two-point correlation function G(q, z, z ′ ) is compatible with the asymptotically rigorous Wertheim-Weeks result for the CW divergence [Eq.( 4)].
The interfacial part of the scattering yields the interfacial structural factor H(q) = I int (q; κ → 0) for sufficiently deep sample penetration on the liquid side (κ −1 ≫ ζ, which is idealized as κ → 0).This expression of H(q) defines an effective, wave number-dependent surface tension γ(q) which is entirely based on density pair correlations [Eq.(8)].Only for small wave numbers, qℓ(T ) ≪ 1, the classical CW divergence, i.e., H(q) ∼ q −2 , is observed in the scattered intensity because in this regime γ(q) ≃ γ 0 [see Eq. ( 66) for the definition of ℓ(T )].Here, we have shown that considering merely the singular part of the background contribution, I ℓ (q; κ) ≈ (ϱ ℓ /2κ)S ℓ (q), as usually done in the analysis of GIXRD data, results in a different interfacial structure factor ̃︀ H(q) and, correspondingly, in a different surface tension γ(q).In particular, the neglected background terms do not drop out in the limit κ → 0 but modify the surface tension at order O(q 2 ) [Eq. (55)].The magnitude of this difference is controlled by the correction integralJ 0 (q), which is determined by the bulk structure factors [Eq.(14)] and which has the dimension of a length.Depending on the temperature,J 0 (q → 0) can have a positive (close to the triple point temperature T t ) or a negative sign (close to the critical temperature T c ).It turns out that at low temperatures, γ(q) and γ(q) exhibit qualitatively different q-dependences (Fig. 10).At higher temperatures, the relative difference is diminished due to the emergence of a contribution of O(q 2 ) in γ(q), which is characterized by a another length ℓ(T ) that grows as T is increased [Eq.(66)].
Based on MD simulations for the truncated LJ fluid with cutoff distance r c = 3.5σ, we have presented evidence that ℓ(T ) in fact diverges upon T ↑ T c and that ℓ(T ) approaches the bulk correlation lengths near criticality (Fig. 11).We observed further that the macroscopic surface tension γ 0 (T ) of this truncated LJ fluid happens to be described very well by the critical scaling law along the whole coexistence line, from the triple point to the critical point (Fig. 6); the same observation was made earlier 95,96 for a cutoff of r c = 2.5σ.Whereas the present study relies on LJ fluids as a generic test bed, analogous largescale MD simulations could be performed for other substances, which would permit the direct comparison between existing GIXRD data [36][37][38] and the simulation results; a similar program was carried out successfully for the bulk structure of water 97 .
At large wave number, γ(q) is only mildly affected by the change of temperature, which together with the increase of ℓ(T ) leads to a maximum in γ(q) at a certain wave number.This phenomenon was observed first in MD simulations 49 and has since been put on firm theoretical ground by Parry et al. [52][53][54][55] .These theoretical studies use DFT calculations for exactly solvable models, which give insight into the structure of the two-point correlation function of the inhomogeneous fluid and which suggests that G(q, z, z ′ ) can reliably be approximated using solely the bulk structure factors and related bulk properties 53,54 .The related expressions can, in principle, be translated into quantitative predictions for the scattered intensity I(q; κ), which could be tested against its small-κ counterpart S tot (q) in the simulations [Eq.(63)].The corresponding expressions, however, are complicated by the liquid-vapor asymmetry 55 and are not yet readily available in an explicit form.Nevertheless, our data for γ(q), covering the full range of the wave number, are qualitatively consistent with the corresponding expectations based on the DFT approximations for G(q, z, z ′ ); this includes the possibility that γ ∞ := γ(q → ∞) > 0, which cannot be resolved from the available data.A complementary, first principles route to G(q, z, z ′ ) has come into reach within a novel Barker-Henderson-like DFT treatment of inhomogeneous fluids 98 , albeit such an endeavor may be technically challenging.
Experiments with phase-separated colloidal suspensions can, in principle, render the knowledge of the three-dimensional positions of all colloids, given the tremendous advances in confocal microscopy during the past two decades, and thus provide experimental data for G(q, z, z ′ ).In previous experiments on polymer-colloid dispersions 5 , single scans of the focal plane perpendicular to the interface were used to obtain slices of the microscopic local density ρ(r) (compare Fig. 1).On this basis, capillary wave theory was then tested by assigning local interface positions and by calculating height-height correlation functions, closely resembling the traditional analysis of simulation data.Yet, the reconstruction of all three-dimensional particle positions from sequences of such focal scans appears to be an ambitious task.
Here, differential dynamic microscopy (DDM) [99][100][101] offers an alternative: it is based on the correlation of intensity images and can yield similar information as contained in the interfacial structure factor H(q) discussed in the present study.DDM is also applicable to dense suspensions which scatter multiple times if the confocal mode of the microscope is used 102 .In this case, the observation volume along the optical axis is restricted by the confocal depth, which introduces corrections in the obtained correlation functions which are analogous to the finite-κ and finite width effects discussed here and for bulk liquids 76 .We expect that a refined interpretation of confocal DDM data, accounting for such corrections, can be developed along the lines presented here.
DDM is also a suitable tool for the characterization of motile suspensions, with micro-organisms or synthetic microswimmers as constituents [102][103][104][105] .Despite being inherently out of equilibrium, such suspensions exhibit a motility-induced phase separation which shares certain universal features of the liquid-vapor transition 106,107 .A surface tension and a surface stiffness have been associated with simulation data for such phenomena 108,109 , although a debate about even the sign of the surface tension shows that active flows and mechanical contributions must be distinguished carefully in order to arrive at a consistent physical interpretation (see Ref. 110 and references therein).Similarly as for the equilibrium situation, a microscopic theory for the two-point density correlations in inhomogeneous active matter would be desirable in order to overcome the ambiguities which are associated with the notion of a fluctuating surface dividing the coexisting phases.
In a recent contribution, Hernández-Muñoz, Tarazona, and Chacón 75 discuss the predictions of extended CW theory for the surface diffraction at liquid-vapor interfaces with the fluctuating surface obtained from the intrinsic sampling method (ISM) 41,73 .This latter approach considers the many-particle structure in the interfacial region, which is accessible within simulations, in order to define a local interface position and, in this sense, goes beyond the mere use of pair correlations as considered here.Within both approaches to γ(q) (i.e., via multi-particle and via pair correlations), there is consensus that the "bending" contribution O(q 2 ) to the q-dependent surface tension should be positive; in particular, this should also hold for almost incompressible liquids at temperatures close to the triple point.Yet, the ISM values for the corresponding length ℓ are considerably larger than what we have found here.We have tested whether this difference can be diminished by tuning the slab width L ℓ of the bulk liquid; L ℓ enters the expression for the background contribution in the simulations [Eq.( 65)].However, for the investigated LJ fluid we find that, for all accessible wave numbers, γ(q) responds only marginally to changes of L ℓ within a physically plausible range (Fig. 12).Thus, we can conclude that our findings for γ(q) -as obtained from density pair correlations -are robust with respect to the details of the definition of the mean interface position.Moreover, the resulting shape of γ(q) is consistent with theoretical predictions [52][53][54][55] for G(q, z, z ′ ).We note that a different, ISM-based definition of γ(q) renders 46,74,75 γ(q) to diverge for large q.We conjecture that this apparent controversy on the large-q behavior of γ(q) is a consequence of whether the definition of γ(q) contains implicit information about three-and many-body correlations or not.This claim is motivated, first, by noting that the ISM approach relies on this additional information whereas the present analysis of the simulation data is restricted to the use of two-point density correlations.Second, we recall the good agreement between the data for γ(q) as obtained along this route and the above-mentioned DFT calculations for γ(q), both using essentially the same definition of γ(q) in terms of G(q, z, z ′ ).Finally, we note that GIXRD experiments on fluid interfaces merely probe two-point correlations, although confining the fluid in a disordered host lattice provides scope for GIXRD studies of higher-order correlations 111 .
In the context of extended CW theory, the physical interpretation of γ(q) is broader than serving just as a proxy for the interfacial two-point correlations; rather, γ(q) provides a mesoscopic characterization of liquid interfaces.In this picture [33][34][35]71 , the interface is thought of as a sharp surface which is locally dressed in an intrinsic density profile perpendicular to the surface (interpolating between the coexisting bulk phases as if there are no CWs). Th fluctuations of the local surface position are then governed by a corresponding surface Hamiltonian such that γ(q)q 2 is the free energy cost associated with surface corrugations of wave number q ("capillary waves").Naturally, such a mesoscopic description must break down at short distances (large wave numbers), for which the molecular discreteness becomes relevant.It has been demonstrated for solvable toy models 50 that one cannot unambiguously single out the naked CW contribution to the local density fluctuations at O(q 2 ) due to a non-local entanglement of the two; however, one can try to push the frontier as far as possible 72,74 .Overall, we can state that the long and slow-burning controversy on the concept of the wave numberdependent surface tension has been resolved, but care must be taken to respect its limitations and to not compare apples with pears.
The non-analytic contribution O(q 2 log(q)) to γ(q) due to dispersion forces is within the scope of mesoscopic surface Hamiltonians and unambiguously identifiable 33,35 .This contribution results in a minimum in γ(q) at mesoscopic wave numbers, corresponding to an enhancement of CW fluctuations.But the magnitude of this effect is not well understood yet: the minimum was found to be surprisingly shallow in simulations of untruncated LJ fluids 46 , but sizable in experimental data from GIXRD on various liquid surfaces [36][37][38] .For the latter, the correction discussed here has the potential to reduce the depth of the minimum to some extent, but we do not expect that it would qualitatively change the conclusions drawn from the experiments.For a direct comparison, simulation data for GIXRD on liquid surfaces of other substances than LJ fluids, e.g., water, would be of great value.It would also be of interest to highlight the role of dispersion forces in the interfacial density correlations, exploiting recent insight into their analytic structure [52][53][54][55] .
It is straightforward to extend the concepts developed here to fluid interfaces in phase-separating binary liquid mixtures 70,[112][113][114][115] , which would lay the basis for probing local changes of the composition (and its fluctuations) in the interfacial region.Within the Gaussian theory, the wave numberdependent surface tension γ(q) not only determines the fluctuations of the local interface height, but also the fluctuations of the local interface normal 116 .Therefore, in addition to GIXRD, the present results for γ(q) are relevant to a variety of further surface-specific experimental techniques such as fluorescence spectroscopy, infrared spectroscopy and linear dichroism, generation of Maxwell displacement current (MDC), second-harmonic generation (SHG), direct measurement of the tilt angle distribution, and laser scanning confocal miscroscopy 116 .Furthermore, it may prove fruitful to investigate the local density correlations under non-equilibrium conditions such as liquid-vapor interfaces in a temperature gradient [117][118][119] .Eventually, the wave number-dependent relaxation dynamics of capillary waves and interfacial fluctuations 16 may be probed within the framework put forward here.To this end, one merely needs to replace the static structure factors of the bulk by their corresponding intermediate scattering functions and to introduce a time lag between the factors of the two-point density correlation function [see Eq. (A1)].
the liquid-vapor critical point (T * c ≈ 1.22), and T * = 0.70, slightly above the triple point temperature.Periodic boundary conditions were applied at all faces of the cuboid simulation box with its edge lengths chosen as L x = L y < L z , so that stable, planar interfaces occur perpendicularly to the z-axis.We used L x = 100σ in order to obtain a large area A = L 2 x = 10 4 σ 2 of the mean interface and in order to access small wave numbers q, which must be integer multiples of 2π/L x .The MD simulations were carried out in the canonical ensemble with the software "HAL's MD package" 120 , which exploits the massively parallel architecture of high-end graphics processors and which is well suited for the study of long-wavelength and low-frequency phenomena in liquids 49,76,114,121 .(Concerning the relationship between the canonical and the grand canonical description of finite-size systems see Ref. 122.) Initial and final particle configurations, time series of observables, as well as correlation functions such as I(q; κ) were stored efficiently in the H5MD file format 123 .Further details on the simulations can be found in Ref. 49.
For the simulation of free liquid-vapor interfaces, a sufficiently thick film of bulk liquid is placed within the simulation domain, and the remaining space is filled with the coexisting vapor phase as to form two parallel, planar liquid-vapor interfaces (Fig. 4).The detailed protocol was as follows: (i) determine the coexisting liquid and vapor densities at the prescribed temperature; (ii) equilibrate the bulk phases of liquid and vapor independently, using slab-like, periodic boxes of width L s = 25σ (T * ⩽ 1.0) or 50σ (T * ⩾ 1.1); (iii) assemble the two phases, after squeezing the different configurations slightly in order to avoid particle overlaps at the boundaries between the phases (e.g., by an amount of 0.5σ along the z-direction via rescaling of the positions); (iv) equip the assembled system with periodic boundaries and let it relax to form an inhomogeneous fluid in equilibrium.
For steps (i) and (ii) only relatively short simulation runs are needed, whereas much longer simulations are required for step (iv) in order to ensure the equilibration of the capillary waves, especially at small wave numbers.In step (iii) we combined several replicas of the vapor phase with one slab of liquid so that the overall box size was L z = 5 × 25σ at low temperatures and L z = 4 × 50σ at high temperatures, yielding a total number of particles of N = 209 300 and 447 000, respectively.
k B T ε Table II.Bulk properties of truncated LJ fluids (r c = 3.5σ) along the liquid-vapor coexistence line.The densities ϱ ℓ and ϱ v of the coexisting liquid and vapor phases were determined from the density profiles obtained in simulations of the inhomogeneous system.The corresponding pressures P, the isothermal compressibilities χ (ℓ) T and χ (v) T , and the correlation lengths ξ ℓ and ξ v stem from separate simulations of the bulk phases; the last four quantities were calculated from OZ fits to the bulk structure factors [Eq.(56)].The numbers in parentheses give the measurement uncertainty in the last digit.
the true scaling behavior sets in only asymptotically for T ↑ T c .In particular, the temperature dependence of c V is found to be non-monotonic along the liquid branch of the coexistence line and has its minimum near T * ≈ 1.05 (data not shown; this calls for future research).
Along the transition line, we have also computed the pressures of the coexisting liquid and vapor phases from the bulk simulations [Table II and Fig. 3(c)], which served as a consistency check.Eventually, the critical pressure was obtained from a separate simulation of the bulk fluid at the quoted critical point (T c , ϱ c ), yielding P c = (0.11074 ± 0.00003)εσ −3 .

Figure 2 .
Figure2.Contribution of the bulk liquid to the GIXRD intensity for five penetration depths κ −1 (in units of the LJ diameter σ) as calculated from Eq. (22).As input serve the simulated bulk structure factors S ℓ (k) of the (truncated) LJ liquids at temperatures T * = k B T/ε = 0.70 (a) and T * = 1.15(b) in reduced units with the LJ energy scale ε.The normalization is taken such that in the limit κ → 0 the curves stay finite and approach S ℓ (q) [Eq.(24)].

Figure 3 .
Figure2.Contribution of the bulk liquid to the GIXRD intensity for five penetration depths κ −1 (in units of the LJ diameter σ) as calculated from Eq. (22).As input serve the simulated bulk structure factors S ℓ (k) of the (truncated) LJ liquids at temperatures T * = k B T/ε = 0.70 (a) and T * = 1.15(b) in reduced units with the LJ energy scale ε.The normalization is taken such that in the limit κ → 0 the curves stay finite and approach S ℓ (q) [Eq.(24)].

Figure 4 .
Figure 4. Snapshot of the simulated LJ liquid-vapor coexistence at T ≈ 0.8T c , comprising 210 000 particles in total; in the vapor phase, only every 10th particle has been drawn for clarity.The yellow frames indicate the mean positions of the two planar interfaces, which delimit the liquid slab of width L s = 25.0σ[see Eq. (B1)]; the area of each of the mean interfaces is (100σ) 2 .The lateral scattering vector q lies in the plane parallel to the interfaces (xy-plane).Along the direction normal to the interfaces (z-axis), the length of the simulation box is L z = 125σ, but only a part of it is shown.

Figure 5 .
Figure 5. Contribution to the GIXRD intensities collected from liquid slabs of finite width L ℓ and with open boundary conditions.Colored solid lines have been obtained from Eq. (33) using the bulk structure factor S ℓ (q) (black line) as input; symbols are results which follow from the direct evaluation of Eq. (1) within MD simulations for bulk liquids.The top panels [(a),(b)] show results for liquid slabs of four widths at fixed inverse penetration depth κσ = 0.01.The bottom panels [(c),(d)] show results for a liquid slab of fixed width L ℓ = 20σ with five penetration depths.The left panels [(a),(c)] refer to a liquid at the temperature T * = 0.70, i.e., close to the triple point, and the right panels [(b),(d)] correspond to T * = 1.15 near the liquid-vapor critical point.The thin horizontal lines at large q indicate the limits given in Eq. (34).

Figure 6 .
Figure 6.Temperature dependence of the macroscopic surface tension γ 0 along the liquid-vapor coexistence line from the triple point temperature T * t ≈ 0.65 . . .0.70 to the critical temperature T * c .The data points stem from MD simulations for truncated LJ fluids (r c = 3.5σ, TableI).The solid line depicts the critical scaling law γ 0 ≃ A γ |t| 2ν for t := (T − T c )/T c upon t ↑ 0 with the Ising exponent ν = 0.630.The inset shows the same data in a rectification plot of the critical law, yielding the critical temperature T * c = 1.220 ± 0.001 and the amplitude A γ = (2.554± 0.008)εσ −2 from a linear regression to the three data points for T * ⩾ 1.0.

Figure 7 .
Figure 7. Leading-order correction J 0 to the bulk background [Eq.(26)] as obtained from quadratures of Eq. (14) with bulk structure factors S ℓ (k) taken from (i) the square-gradient DFT [Eq.(57)], (ii) the Ashcroft-Lekner model (Sect.IV C), and MD simulations of a LJ liquid at the temperatures (iii) T * = 1.15 (close to the critical one, T * c ) and (iv) T * = 0.70 (close to the triple point T * t ).The parameters of the theoretical models were chosen to correspond to the simulated liquids: (i) S 0 = 1.18 and ξ = 1.25σ for the DFT model and (ii) σ = σ at packing fraction η = 0.45 for the hard sphere model.Parts of the figure are reproduced from Ref. 76.
Figure 7. Leading-order correction J 0 to the bulk background [Eq.(26)] as obtained from quadratures of Eq. (14) with bulk structure factors S ℓ (k) taken from (i) the square-gradient DFT [Eq.(57)], (ii) the Ashcroft-Lekner model (Sect.IV C), and MD simulations of a LJ liquid at the temperatures (iii) T * = 1.15 (close to the critical one, T * c ) and (iv) T * = 0.70 (close to the triple point T * t ).The parameters of the theoretical models were chosen to correspond to the simulated liquids: (i) S 0 = 1.18 and ξ = 1.25σ for the DFT model and (ii) σ = σ at packing fraction η = 0.45 for the hard sphere model.Parts of the figure are reproduced from Ref. 76.

Figure 9
Figure 9. Convergence of the effective q-dependent surface tension γ(q; κ → 0) as function of the inverse penetration depth κ for three, fixed wave numbers q at the temperature T * = 0.70.The data for γ(q; κ) were obtained via Eq.(52) from MD simulations for the GIXRD intensity, using only a a fraction L ℓ = 20σ of the width L s of the bulk liquid in order to calculate I(q) and carrying out the decomposition shown in Fig.8(a).Thick solid lines indicate the corresponding limits γ(q) calculated from S tot (q), with setting κ = 0 directly within the simulations.Dashed lines show γ(q) as obtained from S tot (q) when accounting only for the divergent part of the bulk scattering [see the main text and Eq.(64)].The thin blue line (qσ/2π = 0.1) represents the κ-dependence of γ(q → 0; κ) as implied by Eq. (51).

Figure 12 .
Figure 12.Robustness of the q-dependent surface tension with respect to a variation of the width L ℓ of the liquid slab in the background contribution S b (q) by an amount δL.Results for this, a posteriori adjusted quantity γ adj (q) [Eqs.(67) and (68)], are shown for seven distinct values of δL and for the temperatures T * = 0.70 (disks) and T * = 1.15 (diamonds).These results are based on the data provided for γ(q) in Fig.10(a).