Local measures of fluctuations in inhomogeneous liquids: statistical mechanics and illustrative applications

We show in detail how three one-body fluctuation profiles, namely the local compressibility, the local thermal susceptibility, and the reduced density, can be obtained from a statistical mechanical many-body description of classical particle-based systems. We present several different and equivalent routes to the definition of each fluctuation profile, facilitating their explicit numerical calculation in inhomogeneous equilibrium systems. This underlying framework is used for the derivation of further properties such as hard wall contact theorems and novel types of inhomogeneous one-body Ornstein–Zernike equations. The practical accessibility of all three fluctuation profiles is exemplified by grand canonical Monte Carlo simulations that we present for hard sphere, Gaussian core and Lennard–Jones fluids in confinement.


I. INTRODUCTION
The description of inhomogeneous soft matter is an important and challenging task both from a theoretical perspective as well as in practical applications arising throughout physics, chemistry and biology.Many naturally occuring phenomena such as for instance waterrepelling plant surfaces [1], the organization of cells and membranes [2] and protein interactions [3] are attributed to be largely influenced by microscopic mechanical effects that occur upon immersion into a liquid, most commonly water.Some of the underlying principles, e.g. that of hydrophobicity, are even deemed a necessity for the emergence of living matter [4].Likewise, these phenomena can be transferred to technological purposes, as was done e.g. in the design of self-cleaning materials [5,6] or in the construction of specialized nanostructures [7].
It is apparent that a bulk description is insufficient in situations where the behavior of a fluid is determined by an inhomogeneous external environment, e.g.due to being in contact with a substrate or by confinement in cavities or pores [8][9][10].Instead, to properly understand many of the above mentioned phenomena occuring in nature, but also to be able to tailor the behavior of soft matter for technological applications, a framework that accounts for microscopic spatial correlations imposed by an external potential landscape is needed.
From the viewpoint of statistical mechanics, molecular computer simulations [11] as well as density functional theory (DFT) [12,13] are popular and powerful tools to reveal the structure of soft matter.DFT is based on a minimization principle of the grand canonical potential Ω[ρ], which is expressed universally and in principle exactly as a functional of the spatially resolved one-body density profile ρ(r).Functional minimization of Ω[ρ] with respect to ρ(r), carried out at fixed chemical potential µ, temperature T , and external potential, then yields the physical equilibrium density profile.Thus, ρ(r) is allevi-ated from a mere observable to a central quantity that contains, via the connection to the grand potential Ω[ρ], all relevant thermodynamic information.
Given this situation, it is perhaps all the more surprising that, in recent years, not only the density profile ρ(r) was considered in many of the fundamental phenomena mentioned above.Instead, when solvophobic effects dominate the behavior of fluids, the local compressibility χ µ (r) = ∂ρ(r)/∂µ has been shown to uncover much more information about structural correlations than the density profile alone [14][15][16][17][18][19][20][21][22][23].Solvophobicity, also called hydrophobicity if the solute is water, thereby refers to the tendency of a liquid phase to avoid contact with a substrate or solute.This effect can occur over a broad range of length scales, from the solvation of single molecules to density depletion at macroscopic interfaces [9].
Given the simple definition of the local compressibility as a parametric derivative of the density profile with respect to the chemical potential, it is clear that χ µ (r) is easily accessible in DFT calculations, and that it also can be obtained straightforwardly from molecular simulation.Recently, further developments that go beyond the density profile as the primary means of characterizing inhomogeneous fluids have emerged.Notable examples are the use of force density profiles within smart sampling techniques for molecular simulations [24][25][26][27], or reformulations of DFT to better account for two-body correlations [28][29][30][31].
Despite the demonstration of its practical relevance, the definition of χ µ (r) is somewhat ad hoc.The question arises whether the observation of χ µ (r) being a better indicator of solvophobicity than ρ(r) is merely relevant in specific situations, or whether it points to a more general theoretical framework.Indeed, it has been shown by Evans and Wilding that the local compressibility is a suitable and rigorous quantification of localized fluctuations of the particle number [15].Thus, when these local fluctuations govern the behavior of inhomogeneous fluids -as was shown to be the case for the multitude of processes related to hydrophobicity [32] -the local compressibility enables a thorough study of the resulting physical effects.

arXiv:2305.07297v3 [cond-mat.stat-mech] 25 Jul 2023
Based on the above findings and the apparent importance of χ µ (r), one could still ask whether there are different types of density fluctuations that may be captured by additional one-body fluctuation fields, and whether a more fundamental mathematical framework than the mere chemical potential derivative of the density profile can be found.In previous work, we proposed to complement the local compressibility by two additional observables such that a full set of fluctuation profiles is obtained, which reveals even more insight into the nature of local fluctuations [33].For the grand canonical ensemble, those additional one-body profiles are the local thermal susceptibility χ T (r) and the reduced density χ (r).
Coe et al. [21] have investigated the properties of these fluctuation profiles for critical drying based on DFT and on a mesoscopic treatment.The authors concluded that for state points in the critical drying regime, all three measures, when scaled by their respective bulk values, show very similar behavior and the ratio χ T (z)/χ µ (z) is constant when the distance z from the substrate is chosen to be in the gas-liquid interfacial region.Such behavior can be expected, based on general scaling arguments of surface thermodynamics [21].Furthermore, addressing the microscopic behavior of electrolytes, Cats and van Roij [34] argue that the differential capacitance is a probe for the electric double layer structure and the bulk composition of an electrolyte.Inspired by Ref. [33], they formulated an expression for the differential capacitance by correlating the local density with the global charge in the system.Given the significance of the fluctuation profiles in these previous works, we return to fundamentals and show in detail how the three fluctuation profiles can be obtained systematically and consistently from a manybody description of inhomogeneous equilibrium systems.
This paper is organized as follows.Starting from statistical mechanics, as laid out in Section II A, we derive functional generator expressions in Section II B, which give fluctuation measures as response functions to variation of the external potential V ext (r).Via a rewriting of derivatives, these expressions can be reformulated as in practice more easily accessible parametric derivatives with respect to the thermodynamic variables.In Section II C, we demonstrate that microscopic covariance forms exist for all three fluctuation profiles, which constitutes a further method for their numerical computation.In Section II D, one-body Orstein-Zernike equations for χ µ (r) and χ T (r) are derived.The fluctuation profiles are described for the hard sphere model in Section II E and hard wall contact theorems are given in Section II F. While our derivations are formulated in a grand canonical setting, we explicitly point out differences that occur when using the canonical ensemble in Section II G.
Based on grand canonical Monte Carlo (GCMC) simulations as described in Section III, we consider several standard situations using different interparticle interaction potentials, including hard spheres in Section III B, and Lennard-Jones particles and the Gaussian core model in Section III C. The behaviour of the fluctuation pro-files is studied primarily near substrates, where we observe markedly different behavior depending on both the particle and substrate type.The effects of changing the thermodynamic state point are considered, which reveals characteristic behaviour of the fluctuation profiles as a means to probe and predict emerging phase transitions.Lastly, we demonstrate the numerical calculation of fluctuation profiles for different ensembles in Section III D and conclude in Section IV.

A. Grand ensemble
We base our derivations of identities and properties of fluctuation profiles on the grand canonical ensemble (cf.Section II G for a generalization to the canonical ensemble).In the following, first our notation is introduced.
We consider systems of N particles described by the Hamiltonian with mass m, particle coordinates r N ≡ {r 1 , . . ., r N }, and momenta p N ≡ {p 1 , . . ., p N } in D spatial dimensions.The particles interact both via an interparticle potential u(r N ) and with an external potential V ext (r), where r is a generic position variable.In equilibrium, the grand canonical probability distribution function attains the Boltzmann form Here, Ξ denotes the grand partition sum given by Ξ = Tr e −β(H−µN ) where β = 1/(k B T ) with Boltzmann's constant k B , temperature T and chemical potential µ.
The grand canonical trace is defined as and it represents the sum of all particle numbers N and phase space integral over all particle positions and momenta; h denotes Planck's constant.The grand canonical average of a phase space function A r N , p N is defined as The appropriate thermodynamic potential is the grand potential Ω = −k B T ln Ξ, which is related to the internal energy U , entropy S and average number of particles N via Parametric differentiation of Ω with respect to its natural thermodynamic variables yields Note that due to Eq. ( 6), the decomposition (5) of the grand potential possesses the structure of a Legendre transform.The state of a grand canonical system is specified by the values of the variables T , µ and V .Since we consider the presence of a general external potential, we can express any finite system by including hard walls in V ext (r).Thus, we keep the volume V fixed, without loss of generality, and omit the variable V in the following.

B. Functional and parametric derivatives
We first show that one-body fluctuation profiles can be generated systematically as functional derivatives of appropriate global objects with respect to the external potential V ext (r).Additionally, it is demonstrated that the resulting fields can be obtained equivalently from parametric derivatives of the density profile with respect to thermodynamic variables.
We first recall the well-known relation [12] δΩ δV ext (r) T,µ = ρ(r), (7) which expresses that the density profile ρ(r) ≡ N i=1 δ(r i − r) is generated as a response function of the grand potential Ω upon changing V ext (r).Using the decomposition (5) of the grand potential into U , S and N , we take the functional derivative on the left hand side of Eq. ( 7) separately for each term, As before we have made explicit in the notation that the functional derivative is evaluated at constant thermodynamic state point T, µ.Equation ( 8) is a fundamental relation that we take as motivation to naturally consider the three terms on the right hand side separately.We first focus on the last two terms on the right hand side of Eq. ( 8) in order to reformulate the functional derivatives of S and N as parametric derivatives of the density profile.In the latter case, this is achieved by recalling N = −∂Ω/∂µ from relation (6), which is inserted into the functional derivative δ N /δV ext (r).Then, we perform an exchange of functional and parametric derivative, thereby paying attention to the variables that are kept fixed.After this exchange, the density profile can be identified due to Eq. ( 7) and one obtains the identity Here, the external potential is kept fixed when evaluating the parametric derivative.We can identify from Eq. ( 9) the local compressibility as a first fluctuation profile, which has been introduced by Evans and Stewart [14] in the parametric derivative form (11).The above derivation reveals that χ µ (r) can be defined equivalently as a functional generator expression of the average number of particles N with respect to changes in the external potential [15,33].
For the functional derivative of the entropic term in Eq. ( 8), analogous reasoning as in Eq. ( 9) leads to Therefore, a second fluctuation profile emerges from Eq. ( 12), which we define [33] as the local thermal susceptibility Formally, χ T (r) possesses an analogous structure to χ µ (r), as it can be acquired either from a functional generator expression with respect to variation of V ext (r), or alternatively as a parametric derivative of ρ(r).
We define the remaining first term of the right hand side of Eq. ( 8) as the reduced density profile: In contrast to the availability of S or N as a derivative of the grand potential, cf.Eq. ( 6), the internal energy U possesses no such form.Therefore, also no reformulation into a single parametric derivative of the density profile exists for χ (r).However, by rearranging Eq. ( 8) and inserting the functional generator definitions in Eqs.(10), (13) and (15), we can express the reduced density as We recall that χ µ (r) and χ T (r) can be formulated equivalently as parametric derivatives of ρ(r), cf.Eqs. ( 11) and ( 14); thus the reduced density profile can be obtained without the need of explicitly evaluating the functional derivative (15).Additionally, Eq. ( 16) possesses for fixed position r the structure of a local Legendre transform This is analogous to the global Legendre structure of the grand potential resulting from inserting Eq. ( 6) into Eq.( 5).As a summary, starting from Eq. ( 8), we have obtained the local compressibility χ µ (r), the local thermal susceptibility χ T (r) and the reduced density χ (r), which constitute a set of fluctuation profiles in the grand canonical ensemble.The local compressibility χ µ (r) as given via its parametric derivative form in Eq. ( 11) has been studied extensively in Refs.[15][16][17][18][19][20].The authors of these works conclude that χ µ (r) is a better measure of drying than is ρ(r).In the case of drying, the local density of a liquid adsorbed against a solvophobic wall falls below the fluid bulk density even all the way down to the vapour density.This vapour layer in front of the wall induces particularly large localized fluctuations in the number of particles, which are captured by χ µ (r) due to Eq. (11).The above derivation shows that additional fluctuation profiles χ T (r) and χ (r) exist, which correspond to entropic and energetic fluctuations, cf.Eq. (13) and Eq. ( 15) respectively.By providing means to capture those fluctuations, the local thermal susceptibility χ T (r) and the reduced density χ (r) therefore complement χ µ (r) in a way that promises to aid the analysis of fluctuations in inhomogeneous fluids, see Ref. [21,22] for an in-depth analysis for the case of critical drying.
For completeness, the bulk counterparts of χ µ (r) and χ T (r) are defined by replacing the one-body density ρ(r) with the bulk density ρ b = N /V in the parametric derivative definitions ( 11) and ( 14), i.e. χ b µ = ∂ρ b /∂µ and χ b T = ∂ρ b /∂T .Using Eq. ( 6), where the inequality follows due to Ω(µ) being a concave function.Hence, χ b µ is always positive even when interparticle interactions are present.In contrast to the bulk value of the chemical susceptibility being positive, this property does not hold locally for χ µ (r), as will be shown in the simulation results in Section III.
The bulk thermal susceptibility χ b T can be positive as well as negative since it is the off-diagonal element ∂ 2 Ω/(∂µ∂T ) of the Hessian matrix of Ω.Hence, χ T (r), the localized counterpart to χ b T , can also take on positive as well as negative values.

C. Correlator expressions
Expressing χ µ (r) and χ T (r) via Eqs.( 11) and ( 14) as parametric derivatives of ρ(r) already provides an intuitive interpretation of these functions as a local response of the density profile to changes in µ and T , respectively.Unlike the functional generator expressions, parametric derivatives are easily accessible in simulations, albeit still needing multiple runs carried out at neighboring state points with different values of µ or T or relying on histogram methods.Especially near phase coexistence, small shifts of the thermodynamic state point can induce large responses in ρ(r) [14,15], which makes the finite differences difficult to evaluate accurately.In the following, an alternative method for the derivation of fluctuation profiles is presented, which avoids the evaluation of parametric derivatives and reveals the accessibility of χ µ (r), χ T (r) and χ (r) at fixed thermodynamic parameters µ and T .
The basis of the reasoning is the analytical evaluation of the functional derivatives in Eq. ( 8).This is possible in equilibrium, as the grand canonical probability distribution function Ψ is known for a system described by a Hamiltonian H, cf. the familiar Boltzmann form (2) of Ψ.We show in the following that functional differentiation then yields exact microscopic expressions for χ µ (r) [15], χ T (r) and χ (r), which contain covariances of the density operator ρ(r) = N i=1 δ(r i − r) with N and H.The covariance of two phase space functions A r N , p N and B r N , p N is thereby defined as where we recall the definition of the grand ensemble average (4) as indicated by the angular brackets.
To evaluate the functional generator expressions given in Eqs. ( 10), ( 13) and ( 15), we need to differentiate N , S and U functionally with respect to V ext (r).Thereby, functional derivatives of H, Ξ and Ψ are required in turn.We present these in detail in Appendix A.
We start with the response of the average number of particles N = Tr N Ψ with respect to changes in the external potential according to Eq. ( 10) and obtain For the derivation of δΨ/δV ext (r) = −β [ρ(r) − ρ(r)] Ψ, see Eq. (A3) in Appendix A. The density profile ρ(r) in Eq. ( 22) can be taken out of the average, which yields Eq. ( 24) after identification of the covariance.As the left hand side of Eq. ( 21) is the negative functional generator expression (10) of the local compressibility χ µ (r) we obtain the identity [15,17] It is striking that the total number of particles N appears in the covariance, despite χ µ (r) being a localized (onebody) quantity.From inspecting Eq. ( 24), we observe a general mechanism: if an operator A is not explicitly dependent on V ext (r), e.g.A = N , then the functional derivative δ A /δV ext (r)| T,µ yields −β cov (A, ρ(r)).
We turn to the first term of the right hand side of Eq. ( 8) and hence consider the response of U = H = Tr HΨ to changes in V ext (r).Clearly, the Hamiltonian (1) has an explicit dependence on V ext (r).Therefore, unlike in Eq. ( 25), additional terms besides the covariance of H and −β ρ(r) are expected to occur.Using the product rule of differentiation, we evaluate For the derivation of δH/δV ext (r) = ρ(r), we refer the reader to Eq. (A1) in Appendix A. Besides the energy-density covariance, an additional ρ(r) contribution emerges in Eq. ( 28), which is the "direct response" of H upon changes in V ext (r).We recall the definition (15) of χ (r) and find from Eq. ( 28) the covariance form of the reduced density profile We next evaluate the response of the entropy S with respect to changes in V ext (r) to derive the covariance form of the local thermal susceptibility χ T (r).The entropy operator is defined [35] as and the total entropy then follows from S = Ŝ = −k B Tr Ψ ln Ψ.In contrast to the particle number or the Hamiltonian, Ŝ depends explicitly on the many-body distribution function Ψ as given in Eq. ( 2).Functional differentiation of S with respect to V ext (r) then leads to The product rule has been used to obtain Eq. ( 32) and the functional derivative δΨ/δV ext (r) = −β [ρ(r) − ρ(r)] Ψ has been evaluated in Eq. (A3).This yields Eq. (33), where one can successively identify the entropy operator (30) and the covariance (20) to obtain Eq. (35).Comparison with the functional generator expression (13) for χ T (r) establishes Due to Eq. ( 36), χ T (r) can be interpreted as a residual entropy, and it is related to the entropy densities introduced in Refs.[36,37].While Eq. ( 36) constitutes a successful rewriting of χ T (r) into a covariance form, its practical use is not fully revealed yet.This is due to the fact that Ŝ cannot be obtained straightforwardly from molecular simulations as it depends explicity on the partition function, cf.Eq. ( 30).However, one can express Eq. ( 36) via more easily accessible quantities as follows.
We insert the definition (30) of the entropy operator into Eq.( 36) which yields where the Boltzmann form (2) of Ψ has been used in Eq. (37).As ln Ξ = −βΩ is a constant with respect to the phase space variables, its covariance with ρ(r) has vanished in Eq. ( 39), and one arrives at Equation (40) only contains covariances of the density operator with H and N respectively.Contrary to Eq. ( 36), these covariances can be sampled in a standard manner in GCMC simulations without the need for methods such as thermodynamic integration [11,38].
Hence, practically useful covariance forms of χ µ (r), χ T (r) and χ (r) as given in Eqs. ( 25), ( 29) and ( 40) enable the sampling of all three fluctuation profiles directly within a GCMC simulation at fixed thermodynamic state point.This method is applied in Section III, where the fluctuation profiles of prototypical fluid models in confinement are investigated within GCMC simulations.

D. Fluctuation Ornstein-Zernike relations
The standard inhomogeneous two-body Ornstein-Zernike (OZ) relation is an integral equation that connects the pair correlation function h(r 1 , r 2 ) with the direct correlation function c 2 (r 1 , r 2 ) and the density profile as follows [13,39] (for a pedagogical introduction, see e.g.Ref. [40]) (41) We show that analogous integral equations can be derived for the fluctuation profiles χ µ (r) and χ T (r) and that these relations are of simpler, one-body nature.Due to their resemblance to the original OZ equation ( 41), we will also refer to the resulting integral equations for χ µ (r) and χ T (r) as fluctuation Ornstein-Zernike equations.
We proceed similar to the concepts of Refs.[41,42] and start with the fundamental minimization principle of the grand potential Ω[ρ], which is expressed as a functional of a (trial) density profile ρ(r) that attains its minimum at the physical equilibrium density profile, The grand potential is composed of ideal and excess intrinsic free energy functionals, F id [ρ] and F exc [ρ] respectively, as well as of external contributions [12], Explicit evaluation of the functional derivative in (42) after insertion of Eq. ( 43) yields (upon exchanging the left and right hand sides) which is the standard Euler-Lagrange equation of DFT.We have used the exact expression for the ideal gas free energy functional T and can identify the second term on the right hand side of Eq. ( 44) as a contribution proportional to the one-body direct correlation function c 1 (r) = −βδF exc [ρ]/δρ(r) [13].
We now differentiate both sides of Eq. ( 44) with respect to µ upon keeping T and V ext (r) fixed.This amounts to the physical situation of isothermally changing µ in the given system, The functional chain rule has been applied to obtain the second term in Eq. ( 46) and χ µ (r) is identified in Eq. ( 47) via its parametric derivative definition (11).Furthermore, c 2 (r, r ) = −βδ 2 F exc /δρ(r)δρ(r ) denotes the twobody direct correlation function [13].Multiplying Eq. ( 47) by βρ(r) then yields a one-body Ornstein-Zernike relation for the local compressibility To obtain a similar equation for χ T (r), the parametric derivative of the Euler-Lagrange Eq. ( 44) with respect to T is calculated, which yields The derivatives of V ext (r) and µ with respect to T vanish, since neither V ext (r) nor µ depend on T .The dependence of δF exc [ρ]/δρ(r) = −k B T c 1 (r) on T results in the second line of Eq. ( 50).The spatial integral over c 2 (r, r ) and χ T (r ) emerges again from the functional chain rule similar to Eq. ( 46).Additionally, the thermal wavelength Λ contains an explicit dependence on T , such that which is considered in the differentiation of ln(Λ D ρ(r)) in Eq. ( 49) to obtain the kinetic contributions (proportional to D/2) to χ T (r).For the convention of the thermal wavelength set to the particle diameter Λ = σ, all of the following equations still apply, but with D set to zero.
Taking the temperature dependence of Λ into account explicitly also leads to a transformation of the chemical potential µ, which is further illustrated in Appendix B.
Multiplying both sides of Eq. ( 50) with βρ(r) and rearranging terms results in an Ornstein-Zernike relation for the local thermal susceptibility The fluctuation OZ Eqs. ( 48) and ( 52) can be simplified further by performing a splitting of χ µ (r) and χ T (r) into ideal (superscript id) and excess (superscript exc) contributions.Details of an analysis of the local compressibility and local thermal susceptibility for the ideal gas are given in Appendix C. The results are For an interacting system, the excess parts of χ µ (r) and χ T (r) are defined by evaluating Eqs. ( 53) and ( 54) at the equilibrium density profile ρ(r) of the interacting system and considering the difference to the full fluctuation profiles, i.e.
One recognizes that this splitting appears naturally in the fluctuation OZ Eqs. ( 48) and (52), which can hence be rewritten succinctly as

E. Hard sphere fluctuation profiles
For model fluids consisting of hard spheres, further simplifications arise.We start with the fluctuation Ornstein-Zernike equation (52), which simplifies for hard spheres to since ∂c 1 (r)/∂T = 0 holds in this case.Furthermore, as the internal interaction potential vanishes for all allowed microstates in the case of hard sphere interactions, i.e. u(r N ) = 0 in Eq. ( 1), the correlator expression for the local compressibility (40) reduces to The equipartition theorem has been used for the evaluation of the kinetic energy part in Eq. (60).Despite the local fluctuations being one-body fields, Eq. (62) reveals explicitly that they also hold information about many-body correlation effects in interacting systems.The two-body density-density covariance H 2 (r, r ) = cov (ρ(r), ρ(r )) mediates non-local effects of the external potential, which are hence also captured in χ T (r).
For the case of a bulk system, i.e. for vanishing external potential V ext (r) = 0, or in systems with only hard walls, we find by rearrangement of Eq. ( 62), which makes the ratio of thermal and chemical fluctuations a constant, independent of position.Equation ( 63) can be considered to be a contact theorem, as it (trivially) relates the contact value of fluctuation profiles at the hard wall to bulk quantities of the fluid.In the following section, a more general contact theorem for χ T (r) is derived that holds for arbitrary particle interactions, and which reduces to Eq. ( 63) in the case of hard spheres.

F. Fluctuation hard wall contact theorems
Evans and Stewart [14] derived the following contact theorem for the local compressibility, for a planar geometry with distance x from the hard wall and x → ∞ for the bulk fluid.Following their derivation, we show that a contact theorem can be formulated for the local thermal susceptibility, which establishes connections to the internal energy of the system.To derive the contact value of χ T (r), we differentiate the well-known contact theorem of the density [30], with respect to temperature, which yields Here p(µ, T ) indicates the bulk pressure.Substituting the Gibbs-Duhem equation in the first term on the right hand side of Eq. ( 66) and Ω = −pV into the second term results in from which we obtain Equation ( 70) therefore connects at given T and µ the contact value of the thermal susceptibility to the mean bulk internal energy per volume, H /V , and to the bulk density ρ(∞).
Using a splitting of χ T (r), we can show that the theorem (70) holds separately for the relevant contributions.To perform the splitting of χ T (r), we recall its covariance form (40) and insert the Hamiltonian (1), which consists of internal, external and kinetic contributions.For the local thermal susceptibility, this decomposition yields where the kinetic term has been evaluated with the equipartition theorem, When V ext (r) is a hard wall potential, the second term on the right hand side of Eq. ( 71) vanishes since V ext (r) = 0 for all allowed particle positions.Then, for the contact location, Eq. ( 71) can be evaluated by inserting Eq. ( 64) and Eq.(70).A rearrangement yields χ T,int (0 which constitutes a separate contact theorem for the part of χ T (r) that is generated only by the internal interaction u(r N ) of the particles and that we hence indicate as χ T,int (r).
As the thermodynamic relations only hold in the thermodynamic limit, one expects that the above contact theorems are not satisfied in small systems, but converge in simulations for large box sizes, where effects of the wall become negligible for the bulk fluid.A numerical analysis and validation of the contact theorems is shown in Section III C.

G. Canonical ensemble
Throughout the previous sections, the grand canonical ensemble has been used to define and investigate fluctuation profiles.In the following, we show that suitable fluctuation profiles can be derived for the canonical ensemble, where we take as a starting point the Helmholtz free energy F N = U N − T S N as the relevant thermodynamic potential.All relevant canonical quantities are denoted by a subscript N .
We follow a route similar to that taken in Section II B and consider the functional derivative of F N (instead of the grand potential Ω considered in Section II B) with respect to the external potential V ext (r).This leads to the fundamental splitting analogous to Eq. ( 8), as the basis for the derivation of canonical fluctuation profiles.
The functional derivatives are evaluated as in Section II B, but using the canonical trace and probability distribution.This results in the canonical covariance expressions where ρ N (r) is the canonical density profile and cov N (•, •) denotes the covariance (20) but build with canonical averages.As the decomposition (74) of the free energy does not contain the term µ N , the local compressibility χ µ (r) cannot be defined as a corresponding functional derivative in the canonical ensemble.
The form of the canonical covariances Eqs.(75) and (76) formally coincides with the grand canonical covariance identities (29) and (40) for χ (r) and χ T (r), respectively.However, as different ensemble averages are used in the covariances in Eqs. ( 13) and (76) for χ T (r) and χ T,N (r), respectively, χ T (r) and χ T,N (r) are really distinct from each other.The same holds analogously for χ (r) and χ ,N (r), and they are therefore also distinct quantities.This is contrary to observables that are defined as simple averages, such as the density profile ρ(r) = i δ(r − r i ) , since they generally lose their ensemble-dependence in the thermodynamic limit.We will investigate the ensemble differences of the fluctuation profiles below in Section III D.
Similar to Section II F, one can proceed with a decomposition of χ T,N (r) into In the canonical ensemble, integrals over χ T,N (r) vanish due to N being fixed, which can be proven by using the expression in Eq. ( 76), This holds separately as well for the internal and external contribution of χ T (r), cf.Eq. ( 71).We note that in principle, one could still get the grand canonical χ T (r) from canonical data by taking the partial derivative route in Eq. ( 14), which states that ρ(r) must be derived with respect to T along constant µ.Practically, this is hindered by the fact that the particle number N rather than the chemical potential µ is controlled in the canonical ensemble, thus requiring some measurement of µ [43].

A. Numerical methods
As we have shown in Sections II B and II C on the basis of the grand canonical ensemble, several possible techniques exist to make the fluctuation profiles accessible in molecular simulations.We recall that one option is the explicit evaluation of the parametric derivatives of the density profile ρ(r) via finite difference versions of Eqs. ( 11), ( 14) and (16).For this, one must obtain ρ(r) for slightly different temperatures T ± ∆T and chemical potentials µ ± ∆µ to yield an estimate for the fluctuation profiles at a target thermodynamic state point T, µ.This requires carrying out multiple distinct simulations.Arguably more convenient are the correlator expressions Eqs. ( 25), ( 29) and (36), which enable the sampling of all fluctuation profiles in a single molecular simulation at fixed temperature T and fixed chemical potential µ.
Recall that the definition (13) of χ T (r) contains thermal quantities due to the appearance of the entropy operator Ŝ = −k B ln Ψ.We reemphasize that nevertheless χ T (r) can be reformulated in equilibrium via Eq.( 40) as an expression that avoids the explicit evaluation of the partition function or a thermodynamic potential.The local thermal susceptibility can thus be sampled in a straightforward manner similar to what is possible for χ µ (r) using only averages that are easily accessible in a standard simulation.
In the following, we present results of grand canonical Monte Carlo simulations [11] both for fluids in bulk and in planar confinement.We obtain the fluctuation profiles via the correlator sampling scheme.Apart from the practical superiority in this context, the correlator sampling method also leads to reduced error propagation of compound expressions such as χ T (r)/χ µ (r), which will become important below for hard sphere systems.Nevertheless, we have verified consistency with the partial derivative route.
We examine fluids confined by planar walls perpendicular to the x-direction and retain a planar geometry by imposing periodic boundary conditions along the y-and z-axis.The walls are modeled by an external potential V ext (x) for which two common forms will be considered.We only give the expressions of V ext (x) for the left wall at x = 0 in the following.
First, we consider the behavior at a planar hard wall as represented by the external potential Second, a truncated and shifted Lennard-Jones potential is imposed as a representative soft wall given by (83) for x ≤ x c and V ext (x) = 0 for x > x c with cutoff distance x c , energy scale ε w and interaction length σ x .
In the following, the influence of these external potentials on the spatial structure of the fluctuation profiles is examined for systems consisting of hard sphere, Gaussian and Lennard-Jones particles.We remind the reader that we keep the full temperature dependence of the thermal de Broglie wavelength Λ = Λ(T ).

B. Hard Particles
We first turn to the arguably simplest model of a confined fluid, namely hard spheres of diameter σ confined by two parallel planar hard walls.The distance between the walls is chosen sufficiently large such that a mutual influence is negligible and the simulation effectively serves as an approximation of a bulk fluid in contact with a single hard wall.Figure 1 shows the numerical results for all three fluctuation profiles obtained from this simulation.The quantities χ T (r) and χ µ (r) show strong oscillatory behaviour near the hard wall.They even become negative around x ≈ 0.5σ.
Although the strictly positive density profile ρ(r) is a linear combination of χ T (r), χ µ (r) and χ (r), the fluctuation profiles do not have to obey the constraint of positivity separately.For χ µ (r), this stands in contrast to the behaviour of the corresponding bulk value, which is necessarily positive, as shown in section II B. Negative values can also be seen in results by Evans and Stewart [14] for the local compressibility χ µ (r) in Lennard-Jones liquids near a solvophilic wall.Our simulation results show that this effect does not require attractive interactions -be they from internal or external contributions.Local maxima and minima coincide approximately for all fluctuation profiles, although small deviations can still be noticed.Interestingly, while ρ(x) > 0 is composed of the three profiles, χT (x), χµ(x) and χ (x), each one of them does not have to be strictly positive.Here, negative values of χT (x) and χµ(x) occur at x ≈ 0.5σ.Note that χµ(r) must be positive in bulk due to Eq. ( 19).The sign change rather also occurs in an athermal situation where only hard core repulsions exist.
For hard spheres, we have shown in Section II E that the fluctuation ratio χ T (r)/χ µ (r) is a natural quantity to consider, and that this ratio is expected to show markedly different behavior depending on the type of the external potential.In Fig. 2, this ratio is shown both for the case of a hard wall as well as for a Lennard-Jones wall given by Eq. ( 83) with σ x = σ, x c = 2σ and ε w = k B T .For the hard wall, the ratio χ T (r)/χ µ (r) remains spatially constant and attains a value consistent with Eq. ( 63).Conversely, in the case of the soft Lennard-Jones wall the fluctuation ratio deviates from its constant bulk value and shows oscillatory behavior.
In the vicinity of the wall, the value of the susceptibility ratio increases and resembles the shape of the external potential, reproducing its minimum at x ≈ 1σ.Especially for higher packing fractions, i.e. for larger values of µ, oscillations that reach a few hard sphere diameters into the bulk fluid can be observed.This behaviour can be traced back to Eq. ( 62), which shows that effects of the external potential are mediated non-locally via H 2 (r, r ) and contribute directly to the fluctuation ratio.We next turn to the investigation of fluctuation profiles for soft particles, where additional contributions to the Hamiltonian due to a finite interparticle pair potential arise.As two prototypical microscopic models, we choose i) the Gaussian core model as a purely repulsive particle type that smoothly interpolates from an ideal gas to strongly repulsive interactions and ii) the Lennard-Jones fluid as a standard model fluid with both repulsive and attractive contributions where the presence of the later generates a liquid-vapor phase transition.In the following, confinement of those fluid models by two parallel planar hard walls is considered.The local fluctuation profiles are means to probe the vicinity of the thermodynamic state point, cf.Section II B. Therefore, it is natural to investigate in detail the behavior upon changes in T, µ in a grand canonical setting.Depending on the prevailing bulk phase, the potential proximity to bulk and surface phase transitions, and the influence of external potentials, markedly different structure is thereby expected.

C. Lennard-Jones and Gaussian Particles
For simplicity, we keep the chemical potential µ fixed and vary the temperature T across a suitable range.For the Gaussian core model, this enables the analysis of the emergence of structure in the fluctuation fields when gradually decreasing T , as the influence of the interparticle potential barrier becomes increasingly dominant.The energy scale of the Gaussian interparticle interaction potential is ε, the standard deviation is σ and the cutoff radius is 3σ.
The results for temperatures ranging from k B T /ε = 0.1 to 2.5, and for fixed chemical potential µ = 0 are shown in Fig. 3.All three fluctuation profiles are nearly featureless for high temperature.Such behaviour is expected due to the small influence of the interparticle interactions at high temperature.For lower temperatures, the gradual formation of oscillations can be seen.The observed structure of fluctuation profiles for the Gaussian core model in this situation is -as anticipated -similar to the structure observed e.g. in Fig. 1 for hard spheres.
Next, we consider the Lennard-Jones fluid with energy scale ε, length scale σ and cutoff radius 2.5σ, and show results in Fig. 4. We proceed similarly to above and investigate a temperature sweep from k B T /ε = 0.82 to 1.2.Within this range, a liquid-vapor bulk phase transition occurs, the effects of which are reflected clearly in the fluctuation profiles.While the fluctuation profiles in the vapor phase show no significant structure, the sudden appearance of oscillations can be observed when the system transitions to a liquid.Note that in the density profile, apart from a global jump to larger values due to the increased mean density, very little structural information is revealed as ρ(r) remains rather smooth.
This discrepancy between an almost featureless density profile and a dominant signal in the fluctuation profiles is especially obvious when surface effects are accounted for.For temperatures just below the bulk liquid-vapor phase transition, a desorbing fluid is observed, for which a density depletion is visible near the wall.At even lower temperatures, the liquid becomes adsorbing, i.e. the density increases in the vicinity of the walls and shows oscillations.Especially for the desorbing liquid, the fluctuation profiles attain values that are orders of magnitude larger than in the middle of the box and display oscillating behavior, while the density profile still has no significant structural features.This phenomenon is known and well studied for the local compressibility [14], and in accordance with Ref. [21][22][23].We verify here that such an indicative feature is indeed universal to all fluctuation profiles.
A complete set of fluctuation profiles reveals information about general variations of the thermodynamic state point (i.e.not just for variations of µ as in [14]).It is therefore natural to consider -apart from the local compressibility -the full set of fluctuation profiles to study phase transitions and criticality in confined fluids.This might therefore be especially important in situations where the considered system is particularly sensitive to changes in a specific thermodynamic variable.The construction and measurement of the respective fluctuation profile then enables a much more detailed investigation than the density profile alone.We point the reader again to Ref. [21] for their investigation of critical drying.
We note that we have identified the bulk behavior by considering the middle of the box.This assumption is based on the local density approximation (LDA), i.e. the premise that locally, for sufficiently small spatial variation of V ext (r), the value of a one-body field attains the value of the corresponding bulk quantity.That this approximation holds for fluctuation profiles is shown numerically in Appendix D.
In Section II F, we have proven that contact theorems cannot only be derived for χ µ (r) as already done in [14], but also for the local thermal susceptibility χ T (r).As a result, the contact values of the constituent parts of χ T (r) are intricately linked to contributions to the global internal energy of the system.
In the following, we demonstrate the validity of the contact theorem (70) numerically for the Lennard-Jones vapor and liquid as well as for the Gaussian core model fluid.To eliminate finite size effects, we perform simulations for a range of box sizes with fixed cross-section 6σ × 6σ but variable length L and extrapolate σ/L → 0 by fitting a linear function to the contact values thus obtained.Additionally, bulk simulations are conducted for comparison to the respective bulk value that the theorem connects to the contact values.Due to the splitting of χ T (r) according to Eq. ( 71) and the validity of the contact theorem (66) for χ µ (r), we focus here on the contact theorem (73) for the only nontrivial contribution χ T,int (0 + ).The results are shown in Fig. 5, and it is evident that the contact theorem (73) indeed holds in all considered cases.Within numerical accuracy the contact values χ T,int (0 + ) do not change with the size of the simulation box, which means that even for very small boxes one can estimate the total internal energy of a bulk fluid very accurately.

D. Comparing Canonical and Grand Canonical Systems
As laid out in Section II G, fluctuation profiles can be defined canonically, analogous to the procedure based on the grand canonical ensemble, by considering a splitting of the appropriate thermodynamic potential.We demonstrated that, when taking the covariance route, this yields distinct quantities that differ from the grand canonical counterpart.Still, the canonical profiles may share the same qualitative behaviour in characterizing density variations as a response to changing the thermodynamic state point and thus, e.g., indicate proximity to phase transi- Comparison of ρ(x) (black lines) and the thermal susceptibility (blue lines) in canonical (dashed lines) and grand canonical (solid lines) systems of Lennard-Jones particles with liquid-vapor coexistence stabilized by a sinusoidal external potential with a small amplitude of 0.2ε.We show the full canonical and grand canonical profiles (x) and χT (x) as well as the contribution χT (x) in the grand canonical system that arises solely due to potential energy contributions, see Eq. ( 71).Both systems are simulated in a box of size 25 × 6 × 6σ 3 at kBT = 0.85ε.The chemical potential in the grand canonical system is set to µ = −11.8ε,which we have taken as an empirical estimate for the coexistence value at the specified temperature and which leads to a stable liquid-vapor interface in the considered case.In the canonical system, N = 401 is chosen to match the mean density of the grand canonical system.tions and critical phenomena.
Recall that it is still possible to obtain ensembleinvariant fluctuation profiles by resorting to the partial derivative route, cf.Eqs.(11) and (14), where one must then pay attention to the variables that shall be kept constant (which is cumbersome, but technically possible, in ensembles where those quantities are dependent variables).For the canonical ensemble, apart from the reduced density χ ,N (r) the thermal fluctuation profile χ T,N (r) emerges from the splitting of the Helmholtz free energy F N = U N − T S N .Here we focus on the Lennard-Jones system near liquid-vapor coexistence.The spatial coexistence is stabilized by a weak sinusoidal external potential in both the canonical and the grand canonical case, separating the system into a liquid and a vapor region.
In order to be able to compare fluctuation profiles of a canonical and grand canonical system, we use the same interaction potential, impose the same temperature and match the mean density of the grand canonical system to the canonical density.In the case of small particle numbers where ensemble differences are significant, it is generally not possible to obtain the same density profile in both systems if the same external potential is employed.However, in the thermodynamic limit, ensemble differences in ρ(r) vanish, such that in the finite simulated systems we are able to achieve nearly the same density profile in sufficiently large canonical and grand canonical simulation domains.For approaches that relate canonical and grand canonical DFT to each other, see Refs.[44][45][46][47].
We note that χ T (r) in the grand canonical ensemble consists of multiple contributions as shown in Eq. (71).A similar splitting can be performed for the canonical χ T,N (r).To facilitate a concise comparison, we show both the full thermal susceptibility profiles χ T,N (r) and χ T (r) in the canonical and grand canonical ensemble as well as only the constituent of χ T (r) that arises from contributions of the potential energy, which we denote by χT (r) in the following.Note that formally, the correlator expressions for χ T,N (r) and χT (r) contain the same terms and only differ in the form of the ensemble average, cf.Eqs. ( 71) and (77).Figure 6 shows these profiles as well as the density in the phase-separated Lennard-Jones systems.
As expected, the densities of the canonical and grand canonical systems almost coincide.The local thermal susceptibilities χ T,N (r) and χ T (r) on the other hand are very different, which manifests the anticipated ensemble difference.Nevertheless, for both systems, we expect the extrema of the local thermal susceptibility to occur in regions of substantial particle fluctuation.This can be verified in Fig. 6, and showcases that observables defined by the formalism of Section II are indeed capable of identifying such regions.We observe a larger magnitude of the grand canonical thermal susceptibility χ T (r) as compared to the canonical χ T,N (r), which can be attributed to the variability of the particle number N in the grand ensemble.Note that this also leads to an indirect effect which results in larger fluctuations of the internal interaction potential u(r N ).To further interpret the results in the grand canonical case, it is instructive to consider, say, a growth of the liquid region, thereby shifting the location of the interface and leading to an increase of the total particle number of the system.As a result, more particles within the fluid phase reside in the well of the pair potentials of their neighbors and decrease the internal energy u(r N ).Therefore, fluctuations in the particle number are anti-correlated with the value of the total internal energy and a significant negative peak forms in χT (r) in the vicinity of the interface.In the total profile χ T (r), this effect is counteracted by the additional terms proportional to χ µ (r), see Eq. (71).In the canonical case, fluctuations in the particle number cannot occur, thus leading to a reduced variation of u(r N ).Nevertheless, the dominant scenario for particles within the interface is the escape to the gas phase, thereby increasing the potential energy of the system.Hence, one observes a positive peak of χ T,N (r) at the interface, which is slightly shifted towards the gas phase.In total, while χ T (r) and χ T,N (r) differ quantitatively, they are equally suitable for localizing fluctuations that are consistent with the considered ensemble.This information cannot be easily determined by analyzing the density profile alone.

IV. CONCLUSION AND OUTLOOK
We have described different statistical mechanical approaches to the derivation of a complete set of fluctuation profiles.For the grand canonical ensemble, this set comprises the local compressibility χ µ (r), the local thermal susceptibility χ T (r) and the reduced density χ (r).
The first way of defining these fields is based on the Legendre structure of the thermodynamic potential as laid out in Section II B. A seperate evaluation of functional derivatives that emerge from the splitting Eq. ( 8) of Ω reveals generator expressions for fluctuation profiles.We get those definitions in the grand canonical ensemble by utilizing the fundamental relation Ω = U − T S − µ N of the grand potential, which yields the generator Eqs. ( 10), ( 13) and ( 15) for χ µ (r), χ T (r) and χ (r), respectively.This constitutes a formalism for the derivation of suitable fluctuation fields also for the canonical ensemble in Section II G, and it shows that those fields -similar to the density profile ρ(r) -attain a fundamental status due to their connection to thermodynamics.
By an exchange of the order of functional and thermodynamic parametric derivatives, we have constructed an equivalent definition of fluctuation profiles via the parametric derivatives (11) and ( 14) of the density profile.This reveals an intuitive interpretation of the fluctuation profiles as indicators of local susceptibilities of the density upon changes in µ, T , thereby probing the vicinity of the thermodynamic state point.This route also provides a straightforward way of accessing these fields from simulation data: by the evaluation of finite differences of the density profile, obtained from multiple simulations close to the target thermodynamic state point, one is able to obtain numerical results for χ µ (r), χ T (r) and χ (r).
A simpler numerical alternative for the evaluation of fluctuation profiles opens up when considering the third possible definition based on covariances as in Section II C. Explicit evaluation of the functional generators is possible (as the underlying distribution function is known in equilibrium), and yields Eqs. ( 25), ( 29) and (40) as covariance expressions for χ µ (r), χ (r) and χ T (r), respectively.As all fields obtained from such functional generator expression can be written as covariances, their role as spatially localized measures of fluctuations is again manifest.Also, the covariance route provides a convenient way of numerically determining fluctuation profiles, as only a single simulation run at fixed thermodynamic state point is needed.
Having multiple definitions of the fluctuation profiles at hand, one can analytically evaluate them for the ideal gas, as done in Appendix C. From there, for arbitrarily interacting systems, inhomogeneous Ornstein-Zernike equations were derived in Section II D for the excess (i.e. over ideal) parts of χ µ (r) and χ T (r), which are of simpler one-body type as compared to the standard inhomogeneous OZ equation.If appropriate closure relations could be found in future work, one might get a feasible scheme for obtaining the inhomogeneous fluctuation profilesakin to liquid integral equation theory [13] -from these Ornstein-Zernike equations.A more immediate application arises by considering the fluctuation OZ equations as easily accessible sum rules, which could prove to be useful in the development and verification of novel numerical methods and machine learning techniques [48].
We specialized to the hard sphere model fluid in Section II E and derived contact theorems for the fluctuation profiles of fluids in contact with a hard wall in Section II F. The derivation of a contact theorem for the local compressibility χ µ (r) has been done in [14]; we complement this derivation with one for the local chemical susceptibility χ T (r).Thereby, a natural splitting of χ T (r) into additive contributions inherited from the Hamiltonian is found, and we show that hard wall contact theorems can be derived.
Next we turned to numerical investigation via grand canonical Monte Carlo simulations.Our focus was on showcasing fluctuation profiles for a range of different fluids in confinement, and to verify some of the predicted properties such as the sensitivity to phase transitions and the contact theorems.To prove the universal accessibility, we have considered hard core, Gaussian core and Lennard-Jones model fluids as three fundamentally different particle types.We discriminated between the behavior found at hard walls and at soft walls of Lennard-Jones type.
Especially for fluids undergoing phase transitions or showing surface effects, we saw that fluctuation profiles provide additional structural information that can hardly be revealed on the basis of the density profile alone.The fluctuation profiles are therefore suitable indicators of such phenomena.While this has also been described in other works for the local compressibility χ µ (r), we emphasize that a complete set of fluctuation profiles gives additional information of the vicinity of the thermodynamic state point, and is therefore natural to consider in its entirety.This might be especially important if the system under consideration is predominantly affected by changes in one specific thermodynamic variable.
Lastly, a numerical comparison of canonical and grand canonical systems was given for a spatially phaseseparated Lennard-Jones fluid in Section III D. While the relevant fluctuation profiles do not coincide quantitatively as they carry an intrinsic ensemble difference (cf.Section II G), their qualitative behavior matches and can be deemed universal in that it localizes regions of large fluctuations (in the presented case the liquid-vapor interface).
As an outlook on possible future work, we mention the investigation of fluctuation profiles for more sophisticated particle models and especially for water due to its crucial role in nature.For the latter, a variety of atomistic [49] and coarse-grained [50] interaction potentials exist, and one could assess the resulting fluctuation profiles as an additional criterion of their accuracy.We note that density fluctuations near a curved hydrophobic solute have recently been investigated for a monatomic water model [51] via the local compressibility by Coe et al. [21].The generalization to nonequilibrium interfacial thermodynamics (see e.g.Refs.[52,53] for recent work) is an open question for future work as is the investigation of thermal Noether invariance [54][55][56][57] on local fluctuations.
Additionally, in this work, we focused on the liquid gas regime and refrained from considering the freezing transition, which occurs e.g. for the Gaussian core model both into a bcc as well as an fcc crystal [58].Still, the fluctuation profiles could show intricate behavior in this case as well, and may be considered as precursors of crystallization in the future.Local density approximation for χT (x) (blue squares, solid line), and χµ(x) (green pluses, dashed line) in a bulk system with external potential Eq. (D1).Lines depict inhomogeneous simulation data while symbols are values obtained via LDA.We show (a) a Lennard-Jones liquid in a system of size 12 × 6 × 6σ 3 at µ = −11ε and kBT = 0.9ε with α = 1ε, (b) a Lennard-Jones gas in a system of size 25 × 6 × 6σ 3 at µ = −12ε and kBT = 0.85ε with α = 0.2ε, and (c) a Gaussian fluid in a system of size 25 × 6 × 6σ 3 at µ = −4ε and kBT = 0.5ε with α = 0.2ε.
One can qualitatively suppose that LDA is also a viable approximation for fluctuation profiles by comparison of the behavior of bulk and confined systems.This is not clear a priori, as the covariance expressions contain thermodynamic quantities (e.g.S, N , U ).
To carry out a quantitative comparison, we impose a modulated external potential in a bulk system of size L × 6σ × 6σ where α scales the amplitude of the modulation.We then compare the obtained profiles at position x with bulk systems that are simulated at the same effective chemical potential.This procedure recovers a bulk density equal to the value of the density profile in the inhomogeneous system at position x.As we show, LDA is a viable approximation for the fluctuation profiles as well.
In a very dilute Lennard-Jones gas (ρ b ≈ 0.025 σ −3 ) with α = 0.2 ε and L = 25 σ shown in Fig. 8(b), values of bulk systems coincide with the local values of the fluctuation profiles in the perturbed system.This can also be verified for a Lennard-Jones liquid (ρ b ≈ 0.78 σ −3 ) in Fig. 8(a) and for a dilute fluid of Gaussian particles (ρ b ≈ 0.14 σ −3 ) in Fig. 8(c).

Figure 1 .
Figure1.The three fluctuation profiles and the density profile for a hard sphere bulk fluid in contact with a hard wall at µ = −5.0kBT .Layering near the wall is observable in all three fluctuation profiles.Local maxima and minima coincide approximately for all fluctuation profiles, although small deviations can still be noticed.Interestingly, while ρ(x) > 0 is composed of the three profiles, χT (x), χµ(x) and χ (x), each one of them does not have to be strictly positive.Here, negative values of χT (x) and χµ(x) occur at x ≈ 0.5σ.Note that χµ(r) must be positive in bulk due to Eq. (19).

Figure 2 .
Figure 2. The ratio χT (x)/χµ(x) for the hard sphere fluid in contact with (a) a hard wall and (b) a Lennard-Jones wall is shown for varying chemical potential µ/(kBT ) = −8, −9, −10.While the fluctuation ratio remains spatially constant in (a) and attains values according to Eq. (63), it shows oscillatory behavior in (b).This is due to Eq. (62), and shows that nonlocal effects of a soft external potential can be observed in the fluctuation ratio.
core model fluid confined between hard walls at x = 0 and 12σ.The density profile ρ(x) (a), reduced density χ (x) (b), local compressibility χµ(x) (c), and local thermal susceptibility χT (x) (d) are shown.The chemical potential µ = 0 is fixed and the temperature range kBT /ε ∈ {0.1, 0.3, 0.6, 1.0, 1.5, 2.5} is considered.Arrows and the color gradient from blue to red thereby indicate increasing temperature.The inset in (a) shows the density profile for low temperatures in more detail, whereas the insets in (b), (c) and (d) depict the respective fluctuation profiles divided by the normalized density to illustrate structural variations larger than the density.

Figure 4 .
Figure 4. Planar system of Lennard-Jones particles confined between hard walls at x = 0 and 12σ.The density profile ρ(x) (a), reduced density χ (x) (b), local compressibility χµ(x) (c), and local thermal susceptibility χT (x) (d) are shown for fixed chemical potential µ = −11.5εand varying temperature kBT /ε ∈ {0.82, 0.83, 0.84, 0.86, 1.2}.Arrows and color gradient from blue to red thereby represent increasing temperature.The transition from vapor (dotted) to desorbing liquid (dash-dotted) to adsorbing liquid (solid) can be observed, see Appendix B for a detailed account of the consequences that arise due to including the full temperature dependence of the thermal wavelength Λ. Insets show the respective fields for the vapor state.

2 )Figure 5 .
Figure 5. Contact value χT,int(0 + ) (red crosses) and mean internal interaction energy per volume (orange pluses) in a Lennard-Jones liquid (bottom), Lennard-Jones gas (middle) and Gaussian fluid (top) as a function of system size.Linear fits to the data are included as a guide to the eye.Results for the Lennard-Jones gas are scaled up by a factor of 2000.Arrows on the left side indicate the mean internal energy per volume in a completely periodic bulk system with L = 25σ.

Figure 6 .
Figure 6.Comparison of ρ(x) (black lines) and the thermal susceptibility (blue lines) in canonical (dashed lines) and grand canonical (solid lines) systems of Lennard-Jones particles with liquid-vapor coexistence stabilized by a sinusoidal external potential with a small amplitude of 0.2ε.We show the full canonical and grand canonical profiles (x) and χT (x) as well as the contribution χT (x) in the grand canonical system that arises solely due to potential energy contributions, see Eq. (71).Both systems are simulated in a box of size 25 × 6 × 6σ 3 at kBT = 0.85ε.The chemical potential in the grand canonical system is set to µ = −11.8ε,which we have taken as an empirical estimate for the coexistence value at the specified temperature and which leads to a stable liquid-vapor interface in the considered case.In the canonical system, N = 401 is chosen to match the mean density of the grand canonical system.