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 the measurements of the mean and variance of solvent particle number fluctuations in a set of contiguous subvolumes 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 as follows: (i) the local compressibility profile *χ*(**r**) of Evans and Stewart is considerably more sensitive to variations in the strength of local density fluctuations than the spatial fluctuation profile $F(r)$ and can resolve much more detailed structure; and (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 a 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 best-known examples arise in the context of 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 and the thickness of the film—and hence the Gibbs adsorption—diverge 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, for example, the reviews by Dietrich^{1} and Bonn *et al.*^{2} Drying is the analog 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 as follows: (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 (iii) 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 a 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 *μ*, i.e., *χ*(**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 $\xi \Vert 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,12} introduced 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 an 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 analog 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 simulations; 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. This work was important in identifying how measures other than the density profile itself might provide information about the degree of hydrophobicity of a given substrate. However, with this measure, 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 *χ*_{fl}(*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. Section 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.

## II. BACKGROUND AND METHODS

### A. Spatial fluctuation profile

*V*=

*L*

^{d}(with

*d*being the dimensionality of the system) defined by $\kappa T=\u22121V\u2202V\u2202pT$, with

*p*being the system pressure. Within the grand canonical (constant

*μVT*) 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 following well-known expression:

^{24}

^{,}

*ρ*

_{b}≡⟨

*ρ*⟩ = ⟨

*N*⟩/

*V*is the fixed bulk number density. Here, we note that since ⟨

*N*⟩ scales linearly with

*V*, 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 $\u27e8\rho 2\u27e9\u2212\u27e8\rho \u27e92=kBT\rho b2\kappa 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-

*NVT*and constant-

*NPT*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.

*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,

*κ*

_{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.

*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,

*V*(

*z*) apply now with regard to

*V*(

*r*).

*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, for example, 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-*NVT* ensemble yields an accurate estimate for *κ*_{T}. It was found that it failed to do so for subvolumes less than about ten 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–32}

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

^{8}is defined within the GCE and takes the following form:

*ρ*(

**r**)⟩ is the ensemble average of the instantaneous density profile $\rho (r)=\u2211i=1N\delta (r\u2212ri)$, where

**r**is the d-dimensional position vector and

*μ*is the system chemical potential. It is straightforward to show (see the 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 $\chi b=\rho b2\kappa T$.

The definition of Eq. (5) is formally for a continuous *d*-dimensional 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.

*ρ*(

**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

^{33}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,

*μ*in the chemical potential to yield an estimate for the average profile corresponding to

*μ*+ Δ

*μ*,

*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

*χ*(

**r**) can also be expressed in terms of a correlator,

^{3,10}which takes the following form:

*χ*(

**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. (9) 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:

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.A Lennard-Jones solvent in contact with a solvophobic spherical solute; the radial profiles

*χ*(*r*) and $F(r)$ are studied using spherical shell subvolumes.SPC/E water

^{34}in contact with a hydrophobic spherical solute; the profiles*χ*(**r**) and $F(r)$ are studied using cubic subvolumes.

### 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.446 10^{37} at a (reduced) temperature *T** = 1.0 = 0.842*T*_{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 long-range 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*. The histograms of the instantaneous density profile are formed as *ρ*(*z*) = *N*(*z*)/*V*(*z*). Figure 1 shows our results for the ensemble average ⟨*ρ*(*z*)⟩ [normalized by the independently determined average bulk liquid density *ρ*_{b} = 0.654(1)] for various choices of Δ*z*. The principal features of these 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 are described in Sec. II. The results for *χ*(*z*) are shown in Fig. 2(a) from which one sees (as already established in Ref. 3) that *χ*(*z*) exhibits a 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 $\chi b=\rho b2\kappa T$ as denoted by the red dashed 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 $(kBT\rho b)\u22121$ 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 $(kBT\rho b)\u22121F(z)$ fails to decay to its bulk value *κ*_{T} far from the walls. The 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 and state point, 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.125*σ*. 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 horizontal dashed 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. While 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 a comparison of profiles for *χ*(*r*) and $F(r)$ for three different values of Δ*r*. Figure 5(a) shows that *χ*(*r*) is insensitive to the subvolume size, and the detailed non-trivial structure of *χ*(*r*) is reproduced consistently for each Δ*r*. By contrast, the fluctuation profile $F(r)$ shown in 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*). While 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.^{39,40} The temperature was set to *T* = 300 K, and the chemical potential to its corresponding coexistence value *βμ* = −15.24.^{38} A solute of radius *R*_{s} = 14 Å was fixed at the origin of a simulation box of size $V=(40A\u030a)3$, and the solute–solvent interactions were assigned similarly to Sec. 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.5A\u030a)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 one-dimensional 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 $\u223c30%$ 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 Figs. 6(b) and 6(c), respectively. In the former case, a strong enhancement of local compressibility is seen, of which the spatial range 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 $(kBT\rho b)\u22121F(r)$ rises from zero to a constant value that is ∼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.5A\u030a)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 equation (4) evaluates as $F(r)=1\u2212\u27e8N(r)\u27e9$. 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 $(kBT\rho b)\u22121F(r)=(kBT\rho b)\u22121(1\u2212\u27e8N(r)\u27e9)$ shown in Fig. 6(c) given the measured value of the bulk molecular number density *ρ*_{b} = 3.348 × 10^{28} m^{−3}.

## 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 sub-particle 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 biomolecule, 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.^{41}

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 Secs. 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. 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,42,43} 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 were made to relate such a behavior to the degree of hydrophobicity of the substrate/solute.^{42} 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} In addition, 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^{39,40} and LAMMPS,^{44} 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.

## ACKNOWLEDGMENTS

N.B.W. is grateful to Tom Underwood for guidance in the use of DL_MONTE. The computer simulations were carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol, as well as the Isambard 2 UK National Tier-2 HPC Service (http://gw4.ac.uk/isambard/) operated by GW4 and the UK Met Office, and funded by EPSRC (EP/T022078/1). R.E. acknowledges the support of the Leverhulme Trust, Grant No. EM-2020-029\4.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Nigel B. Wilding**: Conceptualization (lead); Software (lead); Validation (lead); Writing – original draft (lead). **Robert Evans**: Conceptualization (supporting); Writing – original draft (supporting). **Francesco Turci**: Conceptualization (supporting); Investigation (supporting); Software (supporting); Writing – review & editing (supporting).

## DATA AVAILABILITY

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

### APPENDIX: CORRELATOR FOR THE LOCAL COMPRESSIBILITY

For completeness, we include a derivation of Eq. (9). In the grand canonical (constant-*μ*, *V*, *T*) ensemble, the particle number *N* fluctuates under the control of the chemical potential *μ*. Denote a configuration of *N* particles as ${ri}N$, where *i* = 0 … *N*. For an inhomogeneous fluid under the influence of an external potential *V*_{ext}(**r**) that is a function of the position vector **r**, the Hamiltonian can be written as $H({ri}N|\mu ,V,T)=E({ri}N)+\u222bdr\rho (r)Vext(r)\u2212\mu N$, where $E({ri}N)$ is the configurational energy. The grand partition function follows as $\Xi (\mu ,V,T)=\u2211N=0\u221e\u222b(\u220fi=0Ndri)e\u2212\beta H$, with $\beta =(kBT)\u22121$, where we have suppressed the phase space and combinatorial factors.

*β*

^{−1}ln Ξ w.r.t. the external potential,

*ρ*

_{b}, for which

*ρ*(

**r**) =

*ρ*

_{b}= ⟨

*N*⟩/

*V*and

*χ*(

**r**) =

*χ*

_{b}. Integrating Eq. (A7) w.r.t. the differential volume

*dV*, over the fixed system volume and noting that $\u222b0VdV=V$ and $\u222b0V\rho (r)dV=N$, one obtains

## REFERENCES

*Phase Transitions and Critical Phenomena*

*Theory of Simple Liquids*

*Understanding Molecular Simulation*

The magnetic analogs of *χ*(*z*) and $F(z)$ have been defined by Binder and co-workers to describe surface phase transitions in lattice-based magnetic models and termed *χ*_{n} and *χ*_{nn}, respectively. In accordance with our study, a much stronger magnetic fluctuation signal occurs for *χ*_{n} than for *χ*_{nn}.^{14} However, in contrast to fluids, the choice of subvolume size is uniquely predefined in the lattice context.