What is the best simulation approach for measuring local density fluctuations near solvo/hydrophobes?

Measurements of local density fluctuations are crucial to characterizing the interfacial properties of equilibrium fluids. A specific case that has been well-explored involves the heightened compressibility of water near hydrophobic entities. Commonly, a spatial profile of local fluctuation strength is constructed from measurements of the mean and variance of solvent particle number fluctuations in a set of contiguous sub-volumes of the system adjacent to the solvo/hydrophobe. An alternative measure proposed by Evans and Stewart (J. Phys.: Condens. Matter 27, 194111 (2015)) defines a local compressibility profile in terms of the chemical potential derivative of the spatial number density profile. Using Grand Canonical Monte Carlo simulation, we compare and contrast the efficacy of these two approaches for a Lennard-Jones solvent at spherical and planar solvophobic interfaces, and SPC/E water at a hydrophobic spherical solute. Our principal findings are that: (i) the local compressibility profile $\chi({\bf r})$ of Evans and Stewart is considerably more sensitive to variations in the strength of local density fluctuations than the spatial fluctuation profile $F({\bf r})$ and can resolve much more detailed structure; (ii) while the local compressibility profile is essentially independent of the choice of spatial discretization used to construct the profile, the spatial fluctuation profile exhibits strong systematic dependence on the size of the subvolumes on which the profile is defined. We clarify the origin and nature of this finite-size effect.


I. INTRODUCTION
Fluids adsorbed at substrates or extended solutes typically display distinctive local density profiles together with enhanced density fluctuations.Perhaps the bestknown examples arise in the context of the continuous wetting and drying surface phase transitions that occur for a planar substrate of macroscopic interfacial area.A continuous wetting transition corresponds to adsorption from the bulk vapor whereby a film of liquid intrudes between the vapor and a (sufficiently attractive) substrate whose thickness -and hence the Gibbs excess adsorption -diverges continuously as a thermodynamic control parameter is varied.The approach to the transition is accompanied by growing density-density correlations parallel to the substrate, arising from capillary wave-like fluctuations, whose correlation length ξ ∥ diverges, see e.g. the reviews by Dietrich 1 and Bonn et.al 2 .Drying is the analogue of wetting, but now the bulk is liquid, and a film of vapor intrudes between it and the (repulsive or very weakly attractive) substrate.Recent theoretical and numerical studies by Evans et.al. 3,4 have shown that for realistic choices of fluid-fluid and substrate-fluid potentials a continuous drying transition is likely to occur at bulk coexistence by decreasing the strength of the attractive interaction between the substrate and the fluid.At a continuous drying transition, the thickness of the film of vapor and the negative of the Gibbs excess adsorption diverge, accompanied by growing density-density correlations characterized by a diverging correlation length ξ ∥ .Other situations where a diverging parallel correlation length occurs are: i) near a pre-wetting critical point, ii) a fluid confined between two identical parallel planar substrates on approaching a capillary critical point 5 and ii) a fluid in an asymmetric planar slit where one substrate prefers to wet and the other prefers to dry, e.g.Stewart and Evans 6 and references therein.For all these cases there is ample understanding of the nature of the surface criticality to enable proper analysis of simulation results.
The situation is different for adsorption at solutes of finite extent; the optimal way of measuring solvent density fluctuations is not immediately obvious.A suitable measure should treat the full range of solute size and shape; solutes might range from simple near spherical molecules to much larger objects such as colloidal particles and include more complex entities such as proteins or biological molecules.Often, the solvent considered is liquid water, and then one is concerned with hydrophobic solutes.Various attempts have been made to connect enhanced density fluctuations in water with the degree of hydrophobicity and also with drying phenomena; the recent review of Rego and Patel 7 provides many references.Of course, one cannot define a parallel correlation length for density-density correlations near a general solute.Rather, one requires a robust measure that describes how the strength of density fluctuations varies with the distance from the solute, the size of the solute, and the proximity of the thermodynamic state point of the solvent to bulk coexistence.
Experience with planar substrates indicates that there is a natural measure.This is the local compressibility χ(r) defined as the derivative of the equilibrium density profile ρ(r) with respect to (w.r.t) the chemical potential µ, ie.χ(r) = ∂ρ(r)/∂µ| T , which in bulk, where the density is constant, is proportional to the usual thermodynamic quantity, the isothermal compressibility.In planar substrate geometry, the maximum of χ(z) diverges (essentially) in the same way as ξ 2 ∥ on approaching a continuous surface transition.Employing classical Density Functional Theory (DFT) calculations for Lennard-Jones (LJ) fluids, Evans and Stewart 8 demonstrated that χ(z) provides a valuable measure for characterizing hydrophobicity, or more generally solvophobicity, at planar substrates.Specifically, they showed that for substrates where Young's contact angle θ is very large, but θ < π, the maximum in χ(z) is at least an order of magnitude larger than the bulk value χ b and occurs at distances z within one or two atomic diameters of the substrate.They argued that χ(z) is a much sharper indicator of the degree of solvophobicity of a substrate, as measured by the contact angle, than is the extent of density depletion, as extracted from the density profile ρ(z).Subsequent papers reinforced these ideas using Grand Canonical Monte Carlo (GCMC) simulations and DFT for LJ fluids 3,9 and simulations of SPC/E water 10 .Later studies by Eckert et.al. 11,12introduced a (closely related) local thermal susceptibility χ T (z) = ∂ρ(z)/∂T | µ that is the temperature derivative of the density profile.Coe et.al. 13 showed that on approaching critical drying, the singular behavior of χ(z) drives identical singular behavior in χ T (z).
Note that in an early paper 6 , and references therein, χ(z) was termed the local susceptibility since it is the direct analogue for continuum fluids of the layer magnetic susceptibility χ n in an Ising lattice subject to a surface magnetic field.Specifically, χ n = ∂m n /∂h is the derivative of the (average) magnetization m n , in the n-th layer away from the surface, w.r.t. the external magnetic field h; recall that h is equivalent to the chemical potential µ.The layer susceptibility provides a powerful measure of the strength of local magnetization fluctuations in layer n and measurements of χ n , using Monte Carlo simulations, have played a key role in elucidating the fundamental physics of wetting transitions and phase transitions arising from confinement 14,15 .
Arguments that χ(r) provides an effective measure of the degree of solvo/hydrophobicity at finite spherical solutes were presented in two recent papers 16,17 which draw upon ideas from earlier studies of critical drying at (very large) spherical particles 18 .Similarly to planar systems, a pronounced peak develops in χ(r) close to the solute whose height increases with the proximity to bulk coexistence and now also with the radius of the solute; the latter acts as a further variable in a comprehensive scaling analysis of solvent thermodynamics and density fluctuations at a solvophobe.Importantly, these studies, together with those for the planar cases, emphasize the lack of water-specific mechanisms for the behavior near an extended hydrophobe: the physics of the density fluctuations around a generic solvophobe should be the same.
Of course, there are other measures of density fluctuations near hydrophobes and Evans and Stewart 8 provide a summary.For models of water, or any fluid, confined in a slit pore, the mean square fluctuation (variance) of the total particle number provides a valuable measure of the overall compressibility of the system and is often investigated in GCE simulation, see the commentary by Bratko 19 and references therein.Here we focus on local measures.Acharya et.al. 20 (see also Sarupria and Garde 21 and the review by Jamadagni et.al. 22 ) attempted a definition of a local compressibility that involves a derivative of the density profile ρ(z) w.r.t.pressure.However, it is not clear precisely what the pressure is and how the derivative is performed.
In their paper, Acharya et.al. 20 also describe another quantity, denoted as χ f l (z), which measures the variance of particle number fluctuations in a slab of certain thickness ∆z located at a distance z from a substrate.This quantity provides a simple means of defining a spatial fluctuation profile and is essentially the quantity F(z) that we shall define in Sec.II A. Using molecular dynamics simulations, Willard and Chandler 23 measured F(z) for SPC/E water near a planar Lennard-Jones 12-6 wall and found that this quantity increased, for z close to the wall, with decreasing wall-fluid attraction, i.e. with increasing contact angle.They argued that F(z) provides a quantitative measure of the degree of hydrophobicity.Evans and Stewart 8 noted that their DFT results for χ(z) in the LJ fluid exhibited similar trends to those displayed in the SPC/E results for F(z), but could not perform direct comparisons.
In the present work, we perform a systematic comparison of the utility of the spatial fluctuation profile F(r) (i.e. the generalization to arbitrary geometry of F(z)), and the local compressibility profile χ(r), for probing local density fluctuations.Our paper is arranged as follows.Sec.II describes some pertinent background to fluctuations in particle number.We consider the definitions of F(r) and χ(r), emphasizing the distinction between these; χ(r) is proportional to the correlator (covariance) of the local particle number N (r) and the total particle number N whereas F(r) measures the correlator of N (r) with itself.In Sec.III we present GCMC results for three different physical systems: i) a LJ liquid confined between two planar solvophobic walls, ii) a LJ liquid at a spherical solvophobe, and iii) SPC/E water at a spherical hydrophobe.In each case, we compare the relative merits of the two profile measures in quantifying the strength and range of density fluctuations.We conclude in Sec.IV with a summary and discussion.

A. Spatial fluctuation profile
We start by considering the isothermal compressibility of a bulk (uniform) fluid of volume V = L d (with d the dimensionality of the system) defined by κ T = − 1 V ∂V ∂p T , with p the system pressure.Within the grand canonical (constant µV T ) ensemble (GCE), κ T can be related to the configurational average ⟨N ⟩ and variance ⟨N 2 ⟩−⟨N ⟩ 2 of the fluctuating total number of particles N by the wellknown expression 24 : where ρ b ≡ ⟨ρ⟩ = ⟨N ⟩/V is the fixed bulk number density.Here we note that since ⟨N ⟩ scales linearly with V , then for κ T to be intensive requires that ⟨N 2 ⟩ − ⟨N ⟩ 2 scales similarly.This is equivalent to requiring that the variance of the probability density function (pdf) of the number density ⟨ρ 2 ⟩ − ⟨ρ⟩ 2 = k B T ρ 2 b κ T /V , which is the standard result 25 for the finite-size scaling of Gaussian density fluctuations that arise as a consequence of the central limit theorem.
For statistical ensembles in which the total particle number N is fixed, such as the constant-N V T and constant-N P T ensembles, it is clear that Eq.( 1) cannot be applied at the level of the total system.However, an estimate of κ T can potentially be obtained by evaluating Eq. ( 1) over a domain or 'subvolume' 25 .In the limit in which the size of the subvolume v is large on the scale of the particle size, but small compared to the total size of the system so that v/V → 0, fluctuations in the number of subvolume particles N v will effectively be grand canonical in form.This observation has motivated several authors to attempt to evaluate via constant-N molecular dynamics simulations the spatial dependence of the strength of density fluctuations in inhomogeneous fluids at substrates or solutes by applying Eq. ( 1) to each of a set of contiguous subvolumes that cover the spatial region of interest 20,21,23,26 .Below we consider how this strategy is implemented in practice for various solute/solvent geometries of interest.
As a first example, consider the case of a threedimensional (d = 3) slit geometry in which a solvent occupies the space between a pair of planar substrates located at z = 0 and z = L, with the system assumed periodic in the x-y plane so that the density profile varies only in the z direction.It is natural to discretize the system in the z direction into identical thin parallel slabs of thickness ∆z, each having equal subvolume V (z) = L 2 ∆z, where we have chosen the parallel area=L 2 .Evaluating the mean and variance of the solvent particle number fluctuations in each slab yields a histogram of the fluctuation profile: Here we note in analogy to the discussion for κ T above, that since ⟨N (z)⟩ ∼ V (z), the numerator of Eq. 2 should similarly scale linearly with V (z) for F(z) to provide an intensive measure of the local fluctuations.Another common scenario considers an extended spherical solute particle fixed at the origin and immersed in a solvent.For this geometry, it is natural to discretize the space surrounding the solute into contiguous spherical shells (concentric with the origin), each of which encompasses the space between some radius r and r + ∆r and has a subvolume V (r) = 4πr 2 ∆r.Evaluating the mean and variance of the solvent particle number fluctuations within each shell yields a histogram of the radial fluctuation profile: for which similar scaling considerations as for V (z) apply now with regard to V (r).
The approach can readily be extended to the case of an irregular solute, such as a protein molecule, that lacks the symmetry of the above examples.Here, the spatial variation of the fluctuation of the solvent density in d = 3 can be mapped using a fluctuation profile F(r) calculated w.r.t. a position vector r which we assume can range over the discrete lattice vectors of a space-filling structure of eg.cubic subcells of volume V (r) = (∆l) 3 .Evaluating ⟨N (r)⟩ and ⟨N 2 (r)⟩ for fluctuations in the solvent particle number in each subcell yields a histogram F(r): Each of the above approaches to defining a discretized fluctuation profile entails a choice for the subvolume size and shape.The shape may be suggested by the geometry of the problem, but the size needs to be sufficiently small to resolve the pertinent features of the local density fluctuations.Furthermore, when operating in a fixed N ensemble, the subvolume should be much smaller than the system size; otherwise, density fluctuations will be suppressed, leading to biased results.However, it transpires that if one chooses a subvolume that is smaller than some multiple of the diameter of the solvent fluid particles, then the results can also be heavily biased.Aspects of the issues of measuring fluctuations via subvolumes have previously been highlighted by Villamaina and Trizac 27 and Román et al 28 who considered the extent to which Eq. ( 1) when applied to a square subvolume of a two-dimensional uniform fluid in the constant-N V T ensemble yields an accurate estimate for κ T .It was found that it failed to do so for subvolumes less than about 10 particle diameters in linear size even in the limit when this size was much smaller than that of the system as a whole 27 .Related studies of uniform three-dimensional fluids have considered whether one can correct measurements of subvolume compressibility to yield estimates of the bulk compressibility [29][30][31] .
The finding that subvolume estimates of κ T can exhibit serious finite-size effects is a consequence of the central limit theorem: when the subvolume size is insufficiently large, the local density fluctuations deviate from the Gaussian form that pertains to very large subvolumes.This is true even for situations that are far removed from a bulk or surface critical point so that the correlation length for density-density fluctuations remains of order the particle size.The consequence for fluctuation profiles is that the numerator in each of Eqs.(2-4) scales nonlinearly with the subvolume, thus engendering a subvolume size dependence of the fluctuation profile.However, to date, the detailed consequences of this subvolume finite-size effect for the sensitivity, resolution, and accuracy of measurements of local density fluctuations in inhomogeneous fluids have not been investigated.We address some of these here.

B. Local compressibility profile
The local compressibility profile as introduced by Evans and Stewart 8 is defined within the GCE and takes the form Here ⟨ρ(r)⟩ is the ensemble average of the instantaneous density profile ρ(r) where r is the ddimensional position vector and µ is the system chemical potential.It is straightforward to show (see appendix) that for a bulk system for which χ(r) = χ b is constant by translational invariance, the local compressibility profile is related to the bulk isothermal compressibility via χ b = ρ 2 b κ T .The definition eq. 5 is formally for a continuous ddimensional density profile.However, in practice χ(r) is accumulated as a histogram by discretizing the space of interest into a set of subvolumes of interest.Then r belongs to a discrete set of vectors running over the subvolumes, which may, for example, be cubic cells, spherical shells, or planar slabs as outlined in Sec.II A. In the latter two cases, one obtains one-dimensional profiles χ(r) and χ(z) respectively.
In simulations ⟨ρ(r)⟩ is calculated as a configurational average of histograms ρ(r) = N (r)/V (r) obtained by binning the solvent particle positions into the set of subvolumes 3,10 .A convenient approach for obtaining χ(r) is to explicitly perform the numerical derivative as per Eq. ( 5), by taking a finite difference.The apparent need to perform two simulations at different µ to achieve this can be neatly avoided by employing histogram reweighting 32 of a single simulation.In short, a GCE simulation is performed at some chemical potential µ of interest from which one accumulates an uncorrelated sequence of M simultaneous measurements of ρ(r) and the total number of particles N .From this sequence, one first forms an estimate for the ensemble-averaged density profile: Then this same sequence of measurements is reweighted w.r.t. a small notional change ∆µ in the chemical potential to yield an estimate for the average profile corresponding to µ + ∆µ: where N i is the total number of particles for measurement i.Since ∆µ is notional, it can be chosen to be arbitrarily small (we use ∆µ = 10 −4 ).The local compressibility profile follows simply as A straightforward derivation starting from the grand partition function (see the appendix) shows that χ(r) can also be expressed in terms of a correlator 3,10 which takes the form This expression provides an alternative route to that just described for measuring χ(r).Furthermore, it serves to expose the similarities and differences between the two approaches for probing local fluctuations represented by the quantities χ(r) and F(r).Comparing Eqs.(10) and (4) shows that χ(r) and F(r) are fundamentally distinct in character.Specifically, χ(r) correlates the instantaneous density profile ρ(r), or subvolume particle number N (r), with the instantaneous total particle number N , while F(r) correlates the subvolume particle number with itself.Thus, the numerical value of the product ⟨N (r)N ⟩ far exceeds that of ⟨N 2 (r)⟩.Consequently the central limit theorem (and, by extension, the linear scaling with V (r) at fixed V of the numerator in Eq. 10) is satisfied to a much greater degree for χ(r) than is the case for F(r), eq. 4. As we shall see, this means that χ(r) produces the bulk compressibility value even for a minimal choice of the subvolume size on which it is defined, in sharp contrast to F(r).
We end this Section by noting that the integral of the difference (χ(r) − χ b ) over the volume available to the fluid measures the surface excess compressibility χ e , which is the derivative of the Gibbs excess adsorption w.r.t.chemical potential.χ e is proportional to the difference between the variance of the total number of particles in the inhomogeneous fluid and the corresponding variance in the bulk fluid at the same chemical potential.This quantity provides a powerful measure of the integrated strength of density fluctuations as demonstrated explicitly for planar systems 8 .

III. MONTE CARLO SIMULATION RESULTS
We have investigated the relative merits of the spatial fluctuation profile F and the local compressibility χ via Monte Carlo simulations of three distinct physical setups: (i) A Lennard-Jones solvent confined by a pair of solvophobic planar walls; the profiles χ(z) and F(z) normal to the walls are studied using planar slab subvolumes.
(ii) A Lennard-Jones solvent in contact with a solvophobic spherical solute; the radial profiles χ(r) and F(r) are studied using spherical shell subvolumes.
(iii) SPC/E water 33 in contact with a hydrophobic spherical solute; the profiles χ(r) and F(r) are studied using cubic subvolumes.
In each instance, we work within the GCE 34 , employing the venerable algorithm of Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller 35 , which is celebrated in this special issue.

A. Fluctuations in a slit with solvophobic planar walls investigated with planar slab subvolumes
We construct a slit geometry from a cubic simulation box of volume L 3 by placing planar walls at z = 0 and z = L, with periodic boundary conditions in the x and y directions.For the solvent, we employ a 12-6 Lennard-Jones (LJ) fluid that is truncated at r c = 2.5σ (where σ is the LJ diameter) and left unshifted.The chemical potential µ is tuned to the conditions of bulk liquid-vapor coexistence: the (reduced) µ * = −3.44610 36at a (reduced) temperature T * = 1.0 = 0.842T c .The fluid interacts with both walls through a wall-fluid potential that is infinitely repulsive at z = 0 and z = L and has a longrange attraction of the modified 9-3 LJ form previously studied by Evans et al. 3 .By making the dimensionless wall-fluid attraction very weak (we chose ϵ w = 0.01, in the notation of eq. 5 of ref. 3 ) we render the wall strongly solvophobic, giving rise to a state that is very close to critical drying 3 .For such a state one expects substantial depletion of the density at each wall and greatly enhanced density fluctuations.Note that in our slit setup, the liquid is metastable w.r.t.capillary evaporation, but this does not prevent us from studying the near-drying region at the walls (see the DFT study in Fig. 8 of ref. 8 and Fig. 6 of ref. 3 ).
We discretize the space between the walls into L/∆z contiguous planar slab subvolumes of prescribed equal thickness ∆z and volume V (z) = L 2 ∆z.Histograms of the instantaneous density profile are formed as ρ(z) = N (z)/V (z).profiles are a considerable depletion in density near the solvophobic walls reflecting the incipient drying regions, and a relaxation to the bulk liquid density far from the walls.It should also be noted that within statistical uncertainties the profiles are independent of ∆z.
We now present our measurements of the local compressibility profile χ(z) and the spatial fluctuation profile F(z).The definitions and methodologies for calculating these quantities have been described in Sec.II.Results for χ(z) are shown in Fig. 2(a) from which one sees (as already established in ref. 3 ) that χ(z) exhibits strong enhancement -in this case by a factor of 75 compared to the bulk-in the region of the incipient vapor-liquid interface that forms near the slit walls.Far from the walls, χ(z) decays precisely to its bulk value χ b = ρ 2 b κ T as denoted by the dashed red line (κ T = 0.432(3) being determined by an independent simulation of the liquid in a large periodic system).Similarly to ⟨ρ(z)⟩, we see that the profiles for χ(z) are essentially independent of the choice of ∆z.
Matters are quite different for the spatial fluctuation profile F(z) shown in Fig. 2(b) and scaled by a factor of (k B T ρ b ) −1 to allow comparison with the bulk compressibility κ T whose value is indicated by the red horizontal dashed line.F(z) was calculated using the same sequence of configurations and slab subvolume sizes as for χ(z).Nonetheless, the profile characteristics are quite distinct.Principally there is a strong dependence on subvolume thickness ∆z.Although all profiles exhibit peaks near the wall, these are narrower than those of χ(z) and grow with increasing ∆z.The peak-to-trough enhancement factor for F(z) is considerably less than for χ(z), showing that the former is much less sensitive to variations in local density fluctuations than the latter.It is also less accurate: the profile for (k B T ρ b ) −1 F(z) fails to decay to its bulk value κ T far from the walls.This latter observation accords with the bulk studies of Villamaina and Trizac 27 .While further increasing the slab subvolume thickness ∆z would help to ameliorate these issues with F(z), this would come at the cost of a reduction in the resolution of the peak in this measure of density fluctuations.

B. Fluctuations at a spherical solvophobic solute investigated with spherical shell subvolumes
Next, we consider the case of a solvophobic spherical solute immersed in the same LJ solvent used in sec.III A. The center of the solute particle is fixed at the origin and for the solute-solvent interaction we employ a potential of the form given in Eq. ( 3) of Coe et al 17 , with the attractive well depth set to the very small (reduced) value ϵ sf = 0.01 to render the solute strongly solvophobic.We partition the volume around the solute into concentric spherical shell subvolumes (centered on the solute) extending from radius r to r + ∆r.The instantaneous radial density profile is measured as ρ(r) = N (r)/V (r) with V (r) = 4πr 2 ∆r.
In Fig.
(3) we show ⟨ρ(r)⟩ normalized by its bulk value ρ b for two different solute particles of radii R s = 3σ and R s = 5σ.The profiles were accumulated using spherical shells of equal thickness ∆r = 0.0125σ.One sees that in both cases close to the solute particle there is a pronounced depletion in solvent density reflecting the strong solvophobicity.The density profile in this region also exhibits 'kinks', particularly for the smaller solute.These occur on length scales of about σ and are remnants of the packing structure that occurs in bulk fluids.Away from the solute, the density relaxes to its bulk value ρ b as indicated by the red dashed horizontal line.
The corresponding local compressibility profiles χ(r) for the two solute radii are shown in Fig. (4)a.Close to the solute there is a great enhancement of the local compressibility compared to its bulk value.The peak in χ(r) for the larger solute particle extends over a greater range of radii than that for the smaller one, in accord with previous findings 17 .Evident also is considerable structure in the profile.Whilst this is inherited from the density profile 3 the structure is much richer and more pronounced, with subsidiary peaks in evidence.Far from the solute particle, the local compressibility decays smoothly and accurately to its bulk value χ b .
The corresponding scaled spatial fluctuation profiles are shown in Fig. ( 4)b, and were accumulated using the same subvolumes as for χ(r).Here one sees that F(r) exhibits a much weaker response to the enhanced fluctuations than does χ(r), and as in the slit case, the scaled profile fails to decay to the value of the bulk compressibility.Here the failure is dramatically clear.
To investigate the effects of varying the subvolume shell thickness on the results, we present in fig. 5  contrast the fluctuation profile F(r), Fig. 5(b) proves to be very much less sensitive to the enhanced fluctuations near the solute surface than χ(r).Specifically, F(r) has a much smaller peak compared to its limiting value, fails to reflect accurately the range of the enhanced density fluctuations, and does not resolve the structural features picked up by χ(r).Whilst the sensitivity of F(r) increases with the subvolume size, this is at the cost of spatial resolution and does not approach the level of detail that is provided by χ(r).
C. Fluctuations at a hydrophobic spherical solute investigated with cubic subvolumes for SPC/E water In this third example, we compare the local compressibility profile and the spatial fluctuation profile in cubic subvolumes.This is a scenario that one would likely adopt when seeking to map enhanced density fluctuations near an irregularly shaped entity that exhibits hydrophobic regions, as is the case for some large biomolecules.
The model that we have studied to illustrate this case is SPC/E water at a spherical solute, which we have simulated in the GCE using the open-source multipurpose Monte Carlo simulation engine DL_MONTE 38,39 .The temperature was set to T = 300K and the chemical potential to its corresponding coexistence value βµ = −15.24 37.A solute of radius R s = 14Å was fixed at the origin of a simulation box of size V = (40Å) 3 and the solute-solvent interactions were assigned similarly to Sect.III B by employing a potential of the form given in eq.( 3) of Coe et al 17 between the solute particle and the oxygen atom of the water molecules, with the attractive well depth fixed, as in the previous example, by the very small (reduced) value ϵ sf = 0.01.The system was discretized into equal cubic subcells of size V (r) = (0.5Å) 3 and both the local compressibility profile χ(r) and the spatial fluctuation profile F(r) were accumulated w.r.t.these.Although the subcells are the fundamental units on which measurements were performed, we can improve statistical sampling and generate a onedimensional profile for comparison purposes by exploiting the symmetry of our system to average both χ(r) and F(r) spherically, yielding profiles for χ(r) and F(r).We stress, however, that the latter are not to be confused with the profiles that would have resulted from choosing spherical shell subvolumes in the manner of Sec.III B.
The density profile for this system is shown in Fig. 6(a) and displays a ∼ 30% depletion in the water density close to the solute and a weak oscillatory decay to the bulk.The comparison of the spherically averaged forms of χ(r) and F(r) is shown in Fig. 6(b) and (c) respectively.In the former case, a strong enhancement of local compressibility is seen, the spatial range of which correlates with the depletion in the water density close to the solute, and which decays in a similar oscillatory fashion to the bulk value far from the surface.By contrast F(r), while calculated in the same simulation and with the same set of subcells as used to calculate χ(r), shows no signal of enhanced fluctuations at all -the spherically averaged profile (k B T ρ b ) −1 F(r) rises from zero to a constant value which is approximately 12 times that of the bulk compressibility κ T = 5.92 × 10 −10 m 3 J −1 .This complete insensitivity of F(r) to the density fluctuations is of course traceable to the small cell size V (r) = (0.5Å) 3 used in this case.While Fig. 6(b) shows that this choice of cell size is necessary and warranted to resolve the pertinent features of the local compressibility profile χ(r), it results in a complete lack of signal in F(r).
One can better understand this finding by considering the statistics of very small subvolumes 27 .When V (r) is smaller than the particle size, cell occupancy is limited to N (r) = 0 or 1. Accordingly, the probability distribution function for occupation is binomial with some mean occupancy ⟨N (r)⟩.Given this, one readily finds that the fluctuation profile eq. 4 evaluates as F(r) = 1 − ⟨N (r)⟩.When the subvolume cells are considerably smaller than the particle size or the solvent density is low, as can be the case in proximity to an extended solvo/hydrophobe, then ⟨N (r)⟩ ≪ 1 and the statistics approach the Poissonian limit for which the variance and mean of the fluctuations in N (r) are equal.Accordingly, F(r) approaches unity, which explains the absence of its response to density fluctuations as reflected in Fig. 6(c).For position vectors r corresponding to bulk liquid SPC/E water, we find ⟨N (r)⟩ = 0.0041(1) which is consistent with the limiting (large r ) value of

IV. SUMMARY AND DISCUSSION
In this paper, we have compared the utility of two distinct approaches for measuring the strength of local density fluctuations in fluids near solvo/hydrophobes.The local compressibility χ(r) introduced by Evans and Stewart 8 in this context has been shown to provide high levels of sensitivity, resolution, and accuracy regardless of the geometry or length scales (subvolume size) for which it is applied.In particular as evidenced by Figs.4(a) and 6(b), χ(r) can readily resolve the detailed features of the local compressibility even when they occur on subparticle length scales.As mentioned in the Introduction, this measure is well rooted in the statistical physics of interfacial phenomena, and its behavior in the vicinity of surface phase transitions is well established.
By contrast, the spatial fluctuation profile F(r), previously employed by several authors to study density fluctuations in simulations of water near hydrophobes, appears to be a considerably inferior measure in all respects.Specifically, it generally exhibits a much weaker response to the magnitude of local density fluctuations than χ(r).This can only be mitigated by increasing the subvolume size, which comes at the cost of a loss of spatial resolution.The severity of this trade-off depends somewhat on the geometry.For a planar solvo/hydrophobic substrate the subvolumes that one can reasonably employ may not be so small as to render inadequate the sensitivity and resolution of the spatial fluctuation profile.However, for spherical solvo/hydrophobic solutes where the subvolumes in proximity to the surface are necessarily small, the sensitivity and resolution are considerably worse than those provided by the local compressibility.If one attempts to map the spatial fluctuations around an irregular hydrophobe such as a bio-molecule with acceptable resolution, then the spatial fluctuation profile will likely fail to provide an adequate signal of the local density fluctuations, while the local compressibility retains its full utility 40 The differences between the two methods are traceable -in part at least-to the degree to which they satisfy the central limit theorem in the limiting case of a bulk fluid, as discussed in Sections II and III.On the smallest length scales the occupation statistics of subvolumes are binomial, tending to the Poissonian limit.Gaussian density fluctuations emerge in the bulk only for sufficiently large subvolumes whose linear extent greatly exceeds the local correlation length, and which contain correspondingly large numbers of particles.Comparing the correlators that define the two approaches, eqs.( 10) and (4), one sees that the typical magnitude of the product ⟨N (r)N ⟩ in the correlator for χ(r) far exceeds that of ⟨N 2 (r)⟩ appearing in F(r).It follows that the linear scaling, with subvolume size, of the numerator in the correlator that is required to yield an intensive measure of fluctuation strength is readily achieved for χ(r) but not for F(r).This is why χ(r) yields accurate estimates of the bulk compressibility irrespective of the choice of subvolume size.And these benefits of accuracy and high resolution extend beyond the case of purely bulk fluctuations to yield greater sensitivity to the non-Gaussian near-critical fluctuations that are the root cause of enhanced density fluctuations near extended solvo/hydrophobes 17 .
While the spatial fluctuation profile F(r) has been employed in several papers to study enhanced density fluctuation near extended hydrophobes 20,21,23 , we note that others 7,22,41 refer instead to the form of the pdf of the subvolume particle number P (N v ).The appearance of enhanced fluctuations in the form of a non-Gaussian tail at small N v has been reported and attempts made to relate such behavior to the degree of hydrophobicity of the substrate/ solute 41 .However, this tail is apparently visible only for quite large subvolumes 7 , which seems to rule out the use of P (N v ) to create a high-resolution fluctuation profile.
Finally, we note that a key feature of the local compressibility χ(r) is that it is defined within the GCE.This open ensemble is doubtlessly optimal for studying density fluctuations in inhomogeneous systems because fluctuations can occur on all length scales up to and including the system size.Furthermore, the GCE lends itself to accurate positioning of the thermodynamic state of interest relative to (bulk) phase coexistence.Recall that the deviation from bulk liquid-vapor coexistence, measured by the chemical potential deviation, is a crucial ingredient in ascertaining the origin, range and strength of solvo/hydrophobicity-induced density fluctuations 3 .Additionally, as we have seen, the GCE offers practical efficiencies such as the ability to calculate χ(r) via histogram reweighting (see Sec. II B).Grand Canonical Monte Carlo is implemented in several general-purpose molecular simulation engines such as DL_MONTE 38,39 and LAMMPS 42 and is thus readily accessible for a wide range of applications including complex molecules.Nevertheless, it must be acknowledged that many simulation studies of hydrophobic phenomena are not performed grand canonically but rather using molecular dynamics in a closed (constant-N ) ensemble.While we expect the deficiencies of the spatial fluctuation profile that we have identified in the GCE to be at least as pronounced in closed ensembles as they are in the GCE, future work could usefully consider the benefits offered by analogs of the local compressibility to closed ensembles.Significant steps in this direction have recently been reported by Eckert et al 12 who consider their local thermal susceptibility in both the canonical and grand ensembles.

Fig. 1 FIG. 1 .
FIG.1.Monte Carlo results for the normalized density profile ⟨ρ(z)⟩/ρ b calculated for a truncated LJ liquid at bulk liquidvapor coexistence and reduced temperature T * = 1.0 confined by a pair of solvophobic planar walls.Results are shown for three choices of subvolume size V (z) = L 2 ∆z with L = 25σ and ∆z given in the legend.The overall system volume is V = 25σ 3 and the dashed red line indicates ⟨ρ(z)⟩ = ρ b where the bulk value is determined from an independent simulation of a large fully periodic system.Statistical errors are smaller than the symbol sizes.

FIG. 2 .
FIG. 2. (a)Monte Carlo results for the local compressibility profile χ(z) normalized by its bulk value χ b .These correspond to the density profile shown in Figure1.Results are shown for three choices of subvolume size V (z) = L 2 ∆z with L = 25σ and ∆z given in the legend.The overall system volume is V = 25σ 3 and the dashed red horizontal line corresponds to χ(z)/χ b = 1.(b) The scaled fluctuation profile (kBT ρ b ) −1 F(z) evaluated for the same subvolumes considered in (a).Note the strong dependence on the subvolume thickness ∆z and that the profile does not decay in the slit middle to the correct value of the bulk compressibility κT = 0.432(3)  in reduced LJ units, as indicated by the dashed red horizontal line.In all cases, statistical errors are smaller than the symbol sizes.

FIG. 3 .
FIG.3.Monte Carlo results for the radial density profiles for the truncated LJ fluid at liquid-vapor coexistence for reduced temperature T * = 1.0.The LJ solvent is in contact with a very weakly attractive spherical solute and results are shown for solute radii Rs = 3.0σ and Rs = 5.0σ.The radial profiles were accumulated using spherical shell subvolumes of thickness ∆r = 0.125σ.The overall system volume is V = 25σ 3 and the dashed red horizontal line indicates ρ b , the bulk value.Note the structure ('kinks') in the density profiles which stems from packing effects.Statistical errors are smaller than the symbol sizes.

FIG. 4 .
FIG. 4. (a) Monte Carlo results for the radial local compressibility profile χ(r) for the truncated LJ fluid at liquid-vapor coexistence at T * = 1.0 and in contact with a very weakly attractive spherical solute of radius Rs.Results are shown for Rs = 3σ and Rs = 5σ.The subvolume shell thickness is ∆r = 0.125σ and the dashed red horizontal line corresponds to the bulk value.Note the structure in χ(r) which is related to the kinks in ρ(r) seen in fig. 3. (b) Corresponding results for the scaled radial fluctuation profile F(r) for the same ∆r = 0.125σ.This profile resolves neither the extended range of enhanced density fluctuations seen in χ(r) nor the detailed structure.Moreover, it does not decay to the correct value of the bulk compressibility κT = 0.432(3) in reduced LJ units, as indicated by the dashed red horizontal line.In all cases, statistical errors are comparable with the symbol sizes.

FIG. 5 .
FIG. 5. (a) Monte Carlo results for the radial local compressibility profile χ(r) for the truncated LJ fluid at liquid-vapor coexistence and T * = 1.0 in contact with a very weakly attractive spherical solute of radius Rs = 5.0σ.The dashed red horizontal line corresponds to the bulk value.Data are shown for three values of the spherical shell subvolume thickness as indicated in the legend.(b)The scaled radial fluctuation profile (kBT ρ b ) −1 F(r) evaluated for the same subvolumes considered in (a).Note the strong dependence on the subvolume shell thickness ∆r and that far from the solute the profile does not decay to the correct value of the bulk compressibility κT = 0.432(3) (in reduced LJ units), as indicated by the dashed red horizontal line.In all cases, statistical errors are comparable with the symbol sizes.

FIG. 6 .
FIG.6.Monte Carlo results for SPC/E water in contact with a very weakly attractive spherical solute of radius Rs = 14Å in a periodic box of volume V = (40Å)3 .The temperature T = 300K and µ is set to its bulk coexistence value 37 .(a) The radial density profile ⟨ρ(r)⟩/ρ b ; the inset shows a configurational snapshot.(b) The local compressibility profile χ(r).(c) The scaled radial fluctuation profile F(r) in units of m 3 J −1 .Each of ⟨ρ(r)⟩, χ(r) and F(r) was calculated for a 3D lattice of cubic subvolumes of size ∆V = (0.5Å)3 , and then averaged spherically to produce the radial profiles shown.Statistical errors are comparable with the symbol sizes.The dashed red horizontal line in (a) and (b) corresponds to the bulk value; that in (c) corresponds to the bulk compressibility κT = 5.92 × 10 −10 m 3 J −1 .Note the lack of response of F(r) to enhanced density fluctuations near the hydrophobic solute.