Simulations of water near extended hydrophobic spherical solutes have revealed the presence of a region of depleted density and accompanying enhanced density fluctuations. The physical origin of both phenomena has remained somewhat obscure. We investigate these effects employing a mesoscopic binding potential analysis, classical density functional theory (DFT) calculations for a simple Lennard-Jones solvent, and Grand Canonical Monte Carlo (GCMC) simulations of a monatomic water (mw) model. We argue that the density depletion and enhanced fluctuations are near-critical phenomena. Specifically, we show that they can be viewed as remnants of the critical drying surface phase transition that occurs at bulk liquid–vapor coexistence in the macroscopic planar limit, i.e., as the solute radius *R*_{s} → *∞*. Focusing on the radial density profile *ρ*(*r*) and a sensitive spatial measure of fluctuations, the local compressibility profile *χ*(*r*), our binding potential analysis provides explicit predictions for the manner in which the key features of *ρ*(*r*) and *χ*(*r*) scale with *R*_{s}, the strength of solute–water attraction *ɛ*_{sf}, and the deviation from liquid–vapor coexistence of the chemical potential, *δμ*. These scaling predictions are confirmed by our DFT calculations and GCMC simulations. As such, our theory provides a firm basis for understanding the physics of hydrophobic solvation.

## I. INTRODUCTION

A myriad of chemical and biological phenomena and processes involve a solute dissolved in a solvent such as water. Solutes may range greatly in size and form, from atoms or simple molecules to macromolecules such as proteins. Of particular interest is when the solute has an aversion to the solvent: It is then termed “solvophobic” (or “hydrophobic” in the particular case of water). Given the importance of water as a solvent, hydrophobic solvation is a topic of widespread relevance in physical chemistry and biochemistry.^{1–8} Typically, hydrophobicity is manifested by nonpolar solutes, but large solutes can exhibit both polar and nonpolar regions, leading to amphiphilic behavior. For example, it is argued hydrophobicity drives amphiphilic molecules to self-assemble into micelles^{9} and proteins to bind with ligands.^{10}

In order to understand the basics of hydrophobicity, and how some of the complex phenomena mentioned above might occur, a physical theory is required that explains how the solvent/water orders in the vicinity of a given solvophobe/hydrophobe. For water, it is known that length scales play an important role, and for solutes whose size is comparable to that of a water molecule, the solvent behavior will generally differ from that near a larger solute.^{11–13} Specifically, for the case of a very small hydrophobic solute, the hydrogen bond water network can flex to accommodate the solute, resulting in minimal disruption to the local water structure. However, for an extended hydrophobe, i.e., one whose size is substantially greater than that of a water molecule, the hydrogen bond network becomes disrupted. Computer simulations show that when the ratio of solute diameter to water molecule diameter *σ*_{s}/*σ*_{w} ≳ 3, a region of depleted water density develops around the solute together with an enhancement in the magnitude of density fluctuations, both measured relative to bulk water.^{2–4,6,8,11,14–17} Although the extent and magnitude of these effects has been argued to increase with the degree of hydrophobicity, to date no precise explanation has been offered for their physical origin and dependence on the solute’s size, its material properties and the deviation of the state point of the solvent from bulk liquid–vapor coexistence. There is also little insight into how the strength of local water density fluctuations depends on the distance from the solute. More generally, it is unclear as to whether the phenomena observed for extended hydrophobic solutes arise from the “special” hydrogen bonded network character of water, as seems to be implied by some authors, or is simply a particular case of a more universal solvophobic behavior.

Beyond its influence on nanoscale solvation processes, hydrophobicity also plays an important role on macroscopic length scales. A familiar phenomenon is the ability of certain materials such as the leaves of the lotus plant or fluorinated surfaces such as Teflon to cause liquid water to form sessile drops with large contact angles. For planar substrates and a fluid at bulk vapor–liquid coexistence, the angle *θ* at which the surface of the drop makes contact with the substrate is determined thermodynamically from the three interfacial tensions via Young’s equation. A contact angle *θ* > 90° indicates hydrophobic behavior, with the extreme hydrophobic limit corresponding to *θ* → 180°. Planar substrates for which *θ* > 140° are typically referred to as “superhydrophobic.” Experiments in which water is in contact with a planar superhydrophobic surface have established that a region of depleted water density forms adjacent to the surface,^{18–20} although its extent has been a matter of debate.^{21} Simulation studies indicate that in models of water close to weakly attractive planar substrates, density fluctuations are strongly enhanced relative to the bulk, e.g., Ref. 22 and references therein. However, quantitative measures and the physical origin of the fluctuations remain obscure.

Recent theoretical and simulation studies of hydrophobicity have adopted a new viewpoint that rationalizes the observed phenomenology of water at superhydrophobic planar substrates in terms of the physics of surface phase transitions, specifically the phenomenon of *critical drying*.^{23–26} For a solvent/fluid that is at vapor–liquid coexistence and in contact with a single infinite planar substrate, the contact angle grows continuously as the attractive strength of the substrate–solvent potential *ɛ*_{sf} is reduced. The drying point, $\epsilon sfd$, marks the extreme hydrophobic limit for which the contact angle attains *θ* = 180°.

Rather than considering the experimental scenario of a sessile liquid drop that arises when the number of solvent molecules *N* is fixed, theoretical approaches generally utilize a grand canonical ensemble description in which the chemical potential *μ* is prescribed. Here, one considers the adsorption of solvent/fluid at the planar substrate characterized by a density profile *ρ*(*z*) that measures the average local number density of solvent molecules as a function of the perpendicular distance *z* from the substrate. Grand Canonical Monte Carlo (GCMC) simulations of such a system find that as the drying point is approached on lowering *ɛ*_{sf} toward $\epsilon sfd$ from above, *ρ*(*z*) first displays a depleted low-density region near the substrate followed by the formation of a thin vapor layer, intruding between the bulk liquid and the substrate, which grows in thickness. Very close to the transition, the profile acquires a form similar to that of a free vapor–liquid interface. Simulations display strongly fluctuating vapor bubbles at the substrate, the sizes of which span many length scales.^{25} These bubbles are a manifestation of near-critical behavior; they are associated with the observed enhancement of density fluctuations^{23,25} in the vicinity of the substrate. The typical lateral extent of bubbles is measured by a correlation length *ξ*_{‖} that diverges in a power-law fashion as $\epsilon sf\u2192\epsilon sfd$ from above.^{25}^{,} *ξ*_{‖} is in turn linked to an important spatial measure of density fluctuations, the local compressibility profile,^{27} defined as $\chi (z)=\u2202\rho (z)/\u2202\mu T$. Specifically, one can show^{25} that for distances *z* in the interfacial region $\chi (z)\u223c\rho \u2032(z)\xi \Vert 2$, where *ρ*′ is the gradient of the density profile. This implies that the maximum in *χ*(*z*) also diverges as $\epsilon sf\u2192\epsilon sfd$. For very large but finite values of *ξ*_{‖}, the system can be regarded as comprising a liquid that is separated from the substrate by a vapor film of equilibrium thickness *ℓ*_{eq}, the magnitude of which also diverges as $\epsilon sf\u2192\epsilon sfd$. Extensive details concerning critical drying behavior, including the definitions and values of various critical exponents, are set out by Evans *et al.*^{25} That work used density functional theory (DFT) and GCMC simulation to show that in the vicinity of a critical drying point, the form of *χ*(*z*) provides a sensitive measure of the spatial variation of local density fluctuations near a hydrophobic substrate, while the magnitude of its maximum quantifies the degree of hydrophobicity.^{23–25,27} The local compressibility $\chi (r)=\u2202\rho (r)/\u2202\mu T$ is generally the correlator (covariance)^{23} of the total number operator and the local number density operator at position **r**. Recently, other measures of local density fluctuations have been defined in terms of other correlators.^{28} Investigations^{29} for a model liquid at a planar substrate near critical drying show that the local compressibility drives the form of other measures and provides the sharpest indicator of the strength of fluctuations; so, we choose to employ this particular measure in our present study.

All solvents, including water,^{23} are expected to exhibit critical drying when the solvent is (a) at liquid–vapor coexistence and (b) in contact with a single infinite hydrophobic planar substrate having $\epsilon sf=\epsilon sfd$. However, water at standard temperature and pressure (STP), i.e., ambient conditions, is not quite at coexistence.^{30} Specifically, ambient water has a very small but nonzero supersaturation: the chemical potential deviation *δμ* from coexistence *μ*_{co} is given by *βδμ* = *β*(*μ* − *μ*_{co}) ≈ 10^{−3}, *β* = 1/*k*_{B}*T*. Moreover, in experiments it is not currently possible to realize substrates that correspond to the extreme hydrophobic limit, i.e., which have $\delta \epsilon sf\u2261\epsilon sf\u2212\epsilon sfd=0$, although for materials that are the most superhydrophobic *δɛ*_{sf} is very small. Accordingly the true critical drying point is not attained in real water. This means that experimentally, the contact angle of a water drop *θ* < 180°, and in the grand canonical setting this implies a vapor layer of finite thickness *ℓ*_{eq} forms and the maximum of *χ*(*z*) remains finite. However, the small values of both *βδμ* and *δɛ*_{sf} for water at STP and in contact with a superhydrophobic surface imply that such a system is *near-critical* and hence its properties can be rationalized in terms of critical point scaling concepts.^{23}

As will prove pertinent to the properties of hydrophobic solvation, we note that an important feature of the phenomenology of surface phase transitions is the role played by the *range* of particle interactions.^{26,31–33} The range of both fluid–fluid (ff) and substrate–fluid (sf) interactions is relevant and one must distinguish between (i) long-ranged (LR) interactions, in which the full power-law tail of the dispersion interaction potential pertains, and (ii) short-ranged (SR) interactions. The latter arise when the interaction potential is explicitly truncated at a few molecular diameters as is commonly done in simulations of neutral fluids. The effective interactions in Coulombic systems with explicit charge decay exponentially due to screening. Exponentially decaying effective interactions also classify as SR. Generic surface phase diagrams for planar substrates have been presented in Ref. 26, which reveal dramatically the effects that the interaction range can have on the qualitative nature of the phase behavior for both drying and wetting transitions.

This brief discussion summarizes how hydrophobicity is manifest on length scales varying from the nanoscale to macroscopic planar substrates. Intriguingly, there are common features, namely, the depletion of water density near the hydrophobe and an enhancement of density fluctuations. However, it is fair to argue that there is a more detailed and systematic understanding of the physics of hydrophobicity on macroscopic length scales such as superhydrophobic substrates than there is for microscopic systems such as an extended hydrophobic solute. Moreover, to date, there is no comprehensive theory that unifies the observed phenomena across disparate length scales and that includes a proper account of the effects of the range of solvent–solvent interactions.

In the present work, we attempt to develop such a theory by considering the properties of solvents in the vicinity of a single spherical solvophobic solute of radius *R*_{s}. By treating the solute curvature $Rs\u22121$ as a scaling field that measures deviations from surface criticality, we construct a mesoscopic theory of critical drying, which generalizes earlier work on drying transitions at planar substrates by incorporating an incipient drying (vapor) film around the model solute. The theory assumes that for sufficiently small curvature (large *R*_{s}), the system is near-critical and provides specific predictions for the scaling properties of *ℓ*_{eq} and the maximum of $\chi (r)=\u2202\rho (r)/\u2202\mu T$ in terms of *R*_{s}, *βδμ*, *ɛ*_{sf} for (i) the experimentally relevant case of truly LR fluid–fluid (ff), or solvent–solvent, interactions, making contact with the general phenomenology for drying and wetting transitions mentioned above, and (ii) the case of SR ff interactions typical of simulations that utilize a truncated pair potential interaction between solvent particles. We test the predictions using two microscopic approaches: classical DFT for a simple model solute–solvent system and GCMC simulations of a monatomic water model in contact with a hydrophobic spherical solute. The results confirm the scaling predictions and thus corroborate our hypothesis that the density depletion and enhanced density fluctuations occurring near a solvophobic solute can be regarded as remnants of the critical drying transition. Our study reveals that for these properties, there are no qualitative differences between water and a Lennard-Jones (LJ) solvent, demonstrating that for extended solutes water behaves like any other solvent and hydrophobicity is merely a particular case of universal solvophobic behavior. In particular, by forging the link between critical drying and the properties of solvents at solvophobic solutes with finite radius, our theory provides a firm basis for rationalizing the observed enhancement of the density fluctuations that occur in water at a hydrophobic solute.

Our paper is arranged as follows: Sec. II describes the two model systems that we investigate. The first is a weakly attractive spherical solute particle immersed in a Lennard-Jones solvent, which we study via DFT, and the second is an equivalent solute immersed in a solvent of monatomic water, which we study via GCMC simulation. In Sec. III, we introduce our mesoscopic binding potential theory and present the key scaling predictions for how the thickness of the vapor film *ℓ*_{eq}, which determines the adsorption, and the maximum in the local compressibility depend on the solute radius, undersaturation, and solute–solvent attraction and how these quantities depend on the range of solvent–solvent interactions. Section IV describes our DFT calculations and presents a detailed comparison of the results with the scaling predictions. In Sec. V, we present results of our simulations of a monatomic water model, once again making contact with the scaling predictions. We conclude in Sec. VI with a discussion of our results and their repercussions. A short report on some of this work appeared previously.^{34}

## II. MICROSCOPIC MODELS FOR THE SOLUTE–SOLVENT SYSTEM

In this section, we introduce the two microscopic models that we study in this work. The first model employs a Lennard-Jones (LJ) description of solute–solvent and solute–solute interactions investigated via DFT. The second utilizes a popular monatomic water (mw) model that we investigate via GCMC.

### A. Lennard-Jones solvent

^{25}The setup is outlined below and a detailed description can be found in the thesis of Coe.

^{35}As shown in Fig. 1, a single spherical solute particle of radius

*R*

_{s}resides at the origin and is imagined to be composed of smaller “virtual” particles distributed homogeneously with density

*ρ*

_{s}. Within the mean-field DFT treatment of attraction, solvent particles interact with each other via a pairwise potential of the general form

^{25,36,37}

*r*= |

**r**−

**r**′| and where

**r**and

**r**′ denote the position vectors of a pair of particles.

*r*

_{min}= 2

^{1/6}

*σ*is the distance corresponding to the minimum,

*r*

_{c}is the cutoff radius,

*σ*is the LJ diameter, and

*ɛ*is the well depth of the LJ interaction.

^{35}by integrating over the angular degrees of freedom in which the fluid has homogeneous density to yield

^{35,38}

*σ*

_{s}and well depth

*ɛ*

_{s}, over the volume of the solute

^{35}gives rise to the net solute–solvent external potential

*R*

_{s}and

*r*are as shown in Fig. 1, and

*r*

_{+}=

*r*+

*r*

_{min,sf}, where

*r*

_{min,sf}is the minimum of the external potential. The external potential is shifted for numerical reasons such that the minimum occurs at the impenetrable surface of the solute. It is incorporated in the DFT calculations in the standard manner.

^{35,38}

As mentioned in Sec. I, solvent–solvent interactions are manifestly short ranged (SR) when the interparticle truncation distance is finite and short. In our DFT calculations, we set *r*_{c} = 2.5*σ* when we study the case of SR solvent–solvent interactions. In order to approximate the LR interaction, i.e., ultimate algebraic decay of dispersion forces, we set *r*_{c} = 200*σ*. Such a large truncation distance is sufficient to capture accurately the LR decay.

### B. Monatomic water model

^{39}This coarse-grained model represents a water molecule as a single particle and reproduces the tetrahedral network structure of liquid water using a parameterization of the Stillinger–Weber potential. Within the mw model, particles interact via the potential

^{39}

*ϕ*

_{mw,2}, and three-body,

*ϕ*

_{mw,3}, potentials are

*A*= 7.049 556 277,

*B*= 0.602 224 558 4,

*γ*= 1.2 are constants that determine the form and scale of the potential,

*λ*= 23.15 is the tetrahedrality parameter,

*θ*

_{0}= 109.47° is the angle favored between waters,

*a*= 1.8 sets the cutoff radius,

*σ*

_{mw}= 2.3925 Å is the diameter of a mw particle, and

*ɛ*

_{mw}= 6.189 kcal mol

^{−1}is the mw–mw (water–water) interaction strength. The mw solvent we employ in simulation is evidently SR: both the pair and three-body potentials decay exponentially.

## III. BINDING POTENTIAL ANALYSIS AND ITS PREDICTIONS

Binding potential, or effective interfacial potential, theory is a mesoscopic (coarse-grained) approach for understanding the dependence of the thermodynamic and certain structural properties of an inhomogeneous fluid on parameters, such as the temperature, chemical potential as well as the geometrical and material properties of a substrate such as the strength of substrate–fluid attraction *ɛ*_{sf}. The approach has a long history in the physics of interfacial phenomena. Most pertinent to the present investigation is its deployment in several recent studies of drying for fluids at planar surfaces^{24,25} and its use in ascertaining the surface phase diagrams for wetting and drying for different combinations of substrate–fluid and fluid–fluid interaction range.^{26} A binding potential analysis was also applied earlier to the study of drying around (very) large spherical particles/colloids.^{40–42} These studies focused on the singular behavior of the free energy of solvation and the adsorption in the limit $Rs\u22121\u21920$. Here, we build upon these studies, applying the binding potential approach to a smaller solute but one that is still much larger than the size of a solvent particle. We consider how the density profile, as characterized by the thickness of the vapor film *ℓ*_{eq}, and the local compressibility profile *χ*(*r*) depend on the chemical potential and temperature of the solvent as well as on the radius of the solute and the solute–solvent attractive strength. We shall assume that the solute–solvent interactions (which we refer to as solute–fluid sf) are LR as in Eq. (3), which is the experimentally relevant situation and one commonly adopted in simulations. For the solvent–solvent (ff) interactions we consider (a) the case of SR ff interactions typical of fluid simulations that utilize a truncated pair potential and (b) the experimentally relevant case of truly LR ff interactions.

*R*

_{s}in contact with a liquid that—like water—has a small supersaturation, i.e., the chemical potential

*μ*is slightly above the value for liquid–vapor coexistence

*μ*

_{co}for some prescribed subcritical temperature

*T*<

*T*

_{c}. Provided the solute–fluid attractive strength is sufficiently weak, a vapor film of width

*ℓ*will form around the solute, intruding between the solute surface and the liquid. Following standard treatments,

^{43,44}we employ a sharp-kink approximation to describe the density profile of the fluid around the solute as

*ρ*

_{v}and

*ρ*

_{l}are the coexisting densities of the vapor and liquid, respectively. Such a description assumes that the region of excluded volume between the surface particles of the solute and the fluid particles

*dw*is incorporated into the radius of the solute, such that

*R*

_{s}is the effective radius and $\rho (Rs+)$ is the first nonzero fluid density, as shown in Fig. 2(b); see also Fig. 1. The equivalent planar surface system, explored in previous work,

^{25}is given in Fig. 2(a).

_{ex}for such a system can be written as

^{41,45}

*γ*is the surface tension and

*A*the surface area of the solute–vapor (subscript

*sv*) and liquid–vapor (subscript

*lv*) interfaces,

*δμ*= (

*μ*−

*μ*

_{co}), Δ

*ρ*= (

*ρ*

_{l}−

*ρ*

_{v}),

*V*

_{v}is the volume of the vapor and

*ω*(

*ℓ*|

*R*

_{s}) is the binding potential—the contribution to the free energy that arises from the

*ℓ*-dependent interaction of the incipient vapor–liquid interface and the solute. Following earlier treatments,

^{41,45}we assume that the solute is sufficiently large that

*γ*

_{sv}(

*R*

_{s}) and

*γ*

_{lv}(

*R*

_{s}+

*ℓ*) can be approximated by their planar surface equivalents,

*γ*

_{sv}and

*γ*

_{lv}, respectively. Furthermore, we note that it has been shown previously

^{41}that $\omega (\u2113|Rs)=\omega (\u2113)(1+\u2113/Rs+O(\u21132ln(\u2113/2Rs)/Rs2))$. It follows that the binding potential for a curved surface/solute can be approximated in terms of that for a planar surface

*ω*(

*ℓ*) for large solutes. Adopting these approximations and substituting $Asv=4\pi Rs2$, $Alv=4\pi (Rs+\u2113)2$, and $Vv=4\pi (Rs+\u2113)3/3$ gives

As emphasized in our short report,^{34} $p\u0303$ combines the pressure of the intruding volume of supersaturated vapor with the Laplace pressure arising from the curvature of the (incipient) liquid–vapor interface. It is clear that within the large solute approximation, the two contributions play equivalent roles and the only dependence of the excess grand potential on *R*_{s} is via the Laplace pressure. Thus, in the same approximation, the physics of drying/wetting at the spherical solute is determined by the form of the *planar* binding potential *ω*(*ℓ*), which encapsulates the *l* dependence of the interactions between the liquid–vapor interface and the planar substrate. As we recall below, the specific form of *ω*(*ℓ*) depends on the functional form of fluid–fluid (ff) and solute–fluid (sf) interactions, primarily on whether these are LR or SR in character.

### A. SR ff, LR sf

*ω*(

*ℓ*) takes the form

^{24–26}

*a*(

*T*) is independent of

*l*, depends on temperature, and has dimensions of energy per unit area.

*ξ*

_{b}is the bulk correlation length of the (intruding) vapor phase and

*b*(

*T*) has dimensions of energy and is given by

*b*

_{o}=

*π*Δ

*ρ*/3, where, as previously,

*ρ*

_{s}is the density of (virtual) solute particles and

*ɛ*

_{s}and

*σ*

_{s}denote, respectively, the well depth and diameter of the (virtual) solute particle–fluid particle interaction. Hereafter, for simplicity, we neglect higher order terms (H.O.T).

*R*

_{s}dependence in Eq. (9) forges a potential link between hydro/solvophobicity phenomena on microscopic and macroscopic scales. In the case of a hard solute, where

*ɛ*

_{s}= 0, one finds that $\u2113eq\u223c\u2212lnp\u0303$. As shown by Evans

*et al.*,

^{42}this leads to the identification of two regimes of scaling, separated by the length scale of capillary evaporation

*R*

_{c}= 2

*γ*

_{lv}/

*δμ*Δ

*ρ*,

*δμ*= 0, two regimes of scaling can be identified,

*et al.*

^{25}These results, which pertain to different physical limiting cases, highlight that there are regions of the parameter space (

*δμ*,

*ɛ*

_{s},

*R*

_{s}) for which the behavior of

*ℓ*

_{eq}is dependent predominately on only one parameter. We note that as shown previously,

^{25}for the planar substrate case

*R*

_{s}=

*∞*, Eqs. (15) and (14) imply that critical drying for a system with SR ff, LR sf interactions occurs for

*δμ*= 0,

*ɛ*

_{s}= 0.

The magnitude of the local compressibility at *ℓ* = *ℓ*_{eq} provides a useful measure of the scale of density fluctuations in the neighborhood of a hydro/solvophobic surface. Within the binding potential analysis, this can be obtained by assuming that the density profile is a smooth function of the distance from the solute *ρ*(*r*) = *S*(*r* − (*R*_{s} + *ℓ*)).

*ℓ*

_{eq}can then be found using

^{25}

*ρ*′ is the spatial derivative of the profile. Substituting

*ℓ*

_{eq}given in Eq. (13) yields

*ℓ*

_{eq}, we can identify three regimes in which the behavior of

*χ*(

*ℓ*

_{eq}|

*R*

_{s}) is controlled predominately by

*R*

_{s},

*δμ*, or

*ɛ*

_{s}. For a hard solute, where

*ɛ*

_{s}= 0, we find

^{27}At bulk coexistence

*δμ*= 0, we find

*et al.*

^{25}

### B. LR ff, LR sf

^{26}and references therein,

*b*

_{o}=

*π*Δ

*ρ*/3, as defined previously, and

*z*

_{min}= (2/5)

^{1/6}

*σ*

_{s}is the location of the minimum in the (planar) substrate–fluid potential; see Ref. 35 for more details. Clearly,

*c*(

*T*) is positive for all temperatures. Assuming that

*ρ*

_{s}is constant, the sign of

*b*

_{LR}(

*T*) can change from negative to positive upon increasing

*T*as the vapor density

*ρ*

_{v}increases or, indeed, as the attraction strength

*ɛ*

_{s}decreases. The behavior of

*b*

_{LR}(

*T*) determines the location of the minimum in

*ω*

_{ex}(

*ℓ*) that corresponds to the equilibrium film width

*ℓ*

_{eq}. In particular, the drying temperature

*T*

_{d}, at which the thickness of the vapor film

*ℓ*

_{eq}→

*∞*, is determined by the condition

*b*

_{LR}(

*T*

_{d}) = 0. If one fixes the temperature, it follows that the critical drying point occurs for substrate–fluid attraction strength given by

^{26}Summarizing, critical drying for a system with LR ff, LR sf interactions occurs at bulk coexistence

*δμ*= 0 in the planar limit

*R*

_{s}=

*∞*when the attraction strength is given by Eq. (22).

*b*

_{LR}(

*T*) in the approach to the critical drying point from below allows us to define a dimensionless measure of the deviation from this point,

^{40}The convenience stems from the fact that

*b*

_{LR}(

*T*) ∝

*t*′.

*ℓ*gives

*R*

_{s},

*δμ*, and

*ɛ*

_{s}individually dominate the behavior of

*ℓ*

_{eq}.

*χ*(

*ℓ*

_{eq}|

*R*

_{s}),

^{40}that

*ℓ*

_{eq}can be expressed as a scaling function. An appropriate form is

*χ*(

*ℓ*

_{eq}|

*R*

_{s}) can be written as

*c*(

*T*). Note that while the result for

*b*

_{LR}(

*T*) remains valid beyond the sharp-kink approximation, this is not the case for

*c*(

*T*).

^{46}

## IV. RESULTS FROM CLASSICAL DFT CALCULATIONS FOR A LENNARD-JONES SOLVENT

Our DFT calculations^{47} are based on the familiar Rosenfeld functional for the HS functional and the standard mean-field treatment of attraction,^{36,38,48} i.e., they implement the same free energy functional as described in earlier papers for the case of a planar substrate.^{25,27} Our calculations adopt a system of the form shown in Fig. 1 with LR sf interactions described by Eq. (3). As mentioned earlier, for the ff interactions we consider (i) the SR ff case in which the LJ interparticle potential is truncated at *r*_{c} = 2.5*σ* and left unshifted and (ii) the LR ff case in which the true long-ranged potential is approximated by taking *r*_{c} = 200*σ*. The geometry required to treat a spherical solute gives rise to specific weight functions as described by Roth.^{38}

We perform our DFT calculations at the subcritical temperature *T* = 0.775*T*_{c} in accord with previous work.^{24,25} Further details of the bulk and coexistence state points that we have studied at this temperature are set out in the thesis of Coe,^{35} which also provides a guide to numerical implementation of DFT for a spherical solute. Measuring the density profile *ρ*(*r*) and the local compressibility profile *χ*(*r*) for various values of $\epsilon sf=2\pi \rho s\epsilon s\sigma s3/3$, the effective solute–fluid attraction strength, *δμ*, *T*, and *R*_{s}, provides insight into how the solvophobic response of the solvent is influenced by the solute properties and the proximity of the solvent to bulk coexistence.

### A. Profiles of density and local compressibility

Figures 3(a)–3(f) demonstrate the effect on *ρ*(*r*) and *χ*(*r*) of varying *R*_{s} (various colors) and *βδμ* (from left to right) for a (hard) solute with *ɛ*_{sf} = 0, which is the value for critical drying in this case of LR sf, SR ff interactions. As *R*_{s} → *∞*, *ρ*(*r*) , and *χ*(*r*) tend smoothly to the profiles for the planar substrate, thereby illustrating smooth connection between microscopic and macroscopic solvophobicity. The extent of the depleted density region and the magnitude of the local compressibility increase as *βδμ* is lowered toward zero and the deviation from critical drying is reduced.

These observations are in accord with the binding potential analysis of Sec. III A, which predicts, for the hard solute, a length scale *R*_{c} = 2*γ*_{lv}/*δμ*Δ*ρ* that separates behavior dependent on and independent of the curvature of the solute: For solutes having *R*_{s} ≫ *R*_{c}, the density and local compressibility profiles should be close to those of a planar substrate. Figures 3(a)–3(f) permit a test of this prediction. For *βδμ* = 10^{−3}, *R*_{c} ≈ 913*σ* and hence we would expect profiles for *R*_{s} > 1000*σ* to be indistinguishable from those of the planar substrate, as is indeed confirmed by the DFT results. Corresponding behavior is seen in Figs. 3(c)–3(f) for the cases of *βδμ* = 10^{−4} and *βδμ* = 10^{−5}, for which *R*_{c} = 9134*σ* and *R*_{c} = 91 340*σ*, respectively.

In Figs. 3(g)–3(l), *βδμ* is held constant at 10^{−3} while varying *δɛ*_{sf} (colors) and *R*_{s} (left to right). In this case, the binding potential analysis of Sec. III A again predicts two regimes in which the behavior is dependent and independent of curvature, though in this case the theory delivers no convenient expression for the crossover point. Nevertheless, separate regimes can be identified in Figs. 3(g)–3(l) by comparing the variation in *χ*(*r*) relating to *δɛ*_{sf}/*ɛ* ≥ 0.6 to those of *δɛ*_{sf}/*ɛ* < 0.6 when moving from left to right, increasing *R*_{s}. In the former regime, the scale and form of *χ*(*r*) vary little, while for the latter regime, the position and height of the maximum of *χ*(*r*) increase substantially. Note that *ρ*(*r*) exhibits a weaker evolution, confirming that *χ*(*r*) is by far the more sensitive indicator of solvophobicity.

This latter observation is pertinent when attempting to define a solvophobic substrate. Consider the case of *δɛ*_{sf}/*ɛ* = 1.2 in figures (g)–(k). In all cases, the corresponding density profiles show pronounced oscillations, which appear to originate around the bulk density and since these exhibit a weakly enhanced contact density, it is tempting to interpret such behavior as indicative of a solvophilic solute. However, the local compressibility profiles provide important new insight. While these profiles also exhibit oscillations, these are not centered on the bulk compressibility, as they would if the substrate were truly solvophilic.^{27} We note that for a planar substrate and the fluid at bulk coexistence, *βδμ* = 0.0, *δɛ*_{sf}/*ɛ* = 1.2 corresponds to a contact angle of $\u2248107.4\xb0$ and such a substrate would be designated solvophobic. It follows that local density fluctuations, as measured by *χ*(*r*), appear to be a far more reliable indicator of the degree of solvophobicity than the density profile alone.

Figures 3(m)–3(r) for LR ff, LR sf interactions demonstrate similar features as those in Figs. 3(a)–3(l) for a system with SR ff, LR sf interactions. Figures 3(m) and 3(n) demonstrate the influence of varying *R*_{s} for constant *βδμ* = 10^{−3} and *δɛ*_{sf} = 0. As in the SR ff, LR sf case, curvature dependent and independent regimes can be separated by the value of the parameter *R*_{c}, which for this case is *R*_{c} ≈ 1138*σ*. Figures 3(o)–3(r) compare the influence of varying *δɛ*_{sf}/*ɛ* and *R*_{s} and again the regime dependent on curvature occurs for values of *δɛ*_{sf}/*ɛ* < 0.6. While the general forms of the density and local compressibility profiles for LR ff, LR sf interactions differ little from those of SR ff, LR sf interactions, overall the magnitude of density fluctuations and the extent of density depletion are larger, in agreement with the predictions of Sec. III.

### B. Testing the scaling predictions

*ℓ*

_{eq}and thus

*χ*(

*ℓ*

_{eq};

*R*

_{s}). To do so, we define

*ℓ*

_{eq}in the standard way, e.g.,

^{49}

*ρ*

_{b}is the density of the bulk liquid.

*χ*(

*ℓ*

_{eq};

*R*

_{s}) can then be found by performing the derivative of the density profile w.r.t.

*μ*at

*r*=

*ℓ*

_{eq}.

We employ DFT measurements of *ℓ*_{eq} and *χ*(*ℓ*_{eq}; *R*_{s}) for a large range of values of *δμ*, *ɛ*_{sf}, *R*_{s} in order to perform detailed tests of the scaling predictions of the binding potential analysis.

#### 1. SR ff, LR sf interactions

Results for *ℓ*_{eq} with SR ff, LR sf interactions are shown in Fig. 4, for systems with parameters in the range (10^{−6} ≤ *βδμ* ≤ 10^{−3}, 0.0 ≤ *δɛ*_{sf}/*ɛ* ≤ 1.0, 10*σ* ≤ *R*_{s} ≤ 10^{8}*σ*). Results for planar substrates are also included. Within the figure, color is used to indicate the degree of solvophobicity by associating each value of *δɛ*_{sf}/*ɛ* with the corresponding value of the contact angle that would pertain for a planar system at liquid–vapor coexistence.

Within Fig. 4, values of *ℓ*_{eq} obtained from the DFT results for the case of SR ff, LR sf interactions are compared to Eq. (13). Excellent agreement with the predicted (linear) form is found for *δɛ*_{sf}/*ɛ* < 0.4, for a wide range of parameters—any deviation from the linear relationship is associated with solutes of radius *R*_{s} < 20*σ*. The latter indicates a limit in the size of solute for which the effects of the drying critical point can be felt, and one might consider whether such change of behavior could be related to the change in solvation behavior often predicted to occur for solutes of radius about 1 nm dissolved in water. For *δɛ*_{sf}/*ɛ* > 0.4, the agreement between the binding potential prediction and DFT results is not as good. Considering the density profiles in Figs. 3(h), 3(j), and 3(l), this is unsurprising; their form is far from what might be reasonably described by the sharp-kink approximation upon which the binding potential predictions are based. We note that such values of *δɛ*_{sf}/*ɛ* correspond to contact angles of $<150\xb0$, indicating that the effects of the drying critical point are most strongly felt for very weak *sf* attraction corresponding to the supersolvophobic regime.

Figure 5 compares the predictions of the binding potential analysis for *χ*(*ℓ*_{eq}; *R*_{s}) to DFT results for the same systems as in Fig. 4. Again, we see excellent agreement between (17) and the DFT results; however, the linear behavior is found over a far more limited range: *δɛ*_{sf}/*ɛ* ≤ 0.2. Any deviation within this range is for *R*_{s} < 20*σ*, as in Fig. 4. For *δɛ*_{sf}/*ɛ* = 0.3, there is a clear discrepancy between the prediction and the DFT results. Here, it is important to note limitations in the binding potential analysis. Consider the values of *ℓ*_{eq} for which a clear deviation begins—from Fig. 4, we see these correspond to *ℓ*_{eq} < 2*σ*. The bulk vapor correlation length at this temperature is *ξ*_{b} = 0.51*σ*, hence for the binding potential prediction to be physical, 1 − 3*ξ*_{b}/*ℓ*_{eq} > 0 and therefore *ℓ*_{eq} > 1.53*σ*. For smaller values of *ℓ*_{eq}, the binding potential prediction for *χ*(*ℓ*_{eq}; *R*_{s}) is no longer physical. One might attempt to include the neglected higher order terms; however, these depend on the shape of the liquid–vapor interface and are difficult to calculate. The crucial point is that we are pushing the binding potential treatment to its extremes.

#### 2. LR ff, LR sf interactions

The binding potential prediction for *ℓ*_{eq} for a system with LR ff, LR sf interaction is given in Eq. (26) and involves the scaling function $L$. We plot our DFT results for systems at fixed *βδμ* = 10^{−3}, with parameters in the range 0.0 ≤ *δɛ*_{sf}/*ɛ* ≤ 1.0, 10*σ* ≤ *R*_{s} ≤ 10^{8}*σ* and employing the arguments of this scaling form, in Fig. 6. We choose to make such a comparison, i.e., with the scaling form of *ℓ*_{eq}, because of the difficulty in determining *c*(*T*) accurately. Three temperatures, *T* = 0.7*T*_{c}, 0.775*T*_{c}, and 0.85*T*_{c}, are considered, with the arrow indicating the direction of increasing temperature. Color is used to indicate *δɛ*_{sf}/*ɛ* and the corresponding contact angle is given in the inset.

The temperature dependence of the scaling function is inherent in Eq. (24); the constants *b* and *c* are both temperature dependent. While the formula for *b*(*T*) is easily calculated and this coefficient in the binding potential expansion is expected to be valid beyond the sharp-kink approximation, the coefficient *c*(*T*) is dependent on the form of the liquid–vapor interface^{46} and leads to some ambiguity in determining the explicit scaling function. We do not attempt to ascertain $L$ but note this was attempted for the special case of *βδμ* = 0.0 and very large solutes by Stewart and Evans.^{40} Here, we focus simply on data collapse.

Overall, the data collapse predicted from the binding potential for *ℓ*_{eq} is confirmed by the DFT results in the regime *δɛ*_{sf}/*ɛ* < 0.6 for the three temperatures; any deviation corresponds typically to cases where *R*_{s} < 20*σ*. For *δɛ*_{sf}/*ɛ* > 0.6, there is clear deviation from the predicted functional form, which becomes more pronounced as *δɛ*_{sf}/*ɛ* is increased further. As in the case of SR ff, LR sf interactions, density profiles of such systems cannot be represented accurately using a sharp-kink approximation and hence this deviation is unsurprising. The values of *δɛ*_{sf}/*ɛ* for which data collapse is best obeyed correspond to large contact angles—see inset. Again, this suggests that the influence of the drying critical point is most pronounced for solute–fluid interaction strengths for which the contact angle is *θ* > 150°.

Turning finally to our DFT results for *χ*(*ℓ*_{eq}; *R*_{s}) for the LR ff, LR sf case, we tested the predicted scaling of Eq. (27) in the main part of Fig. 7. However, it is also possible to utilize the relationship in Eq. (25) and the inset compares this prediction to the DFT results. Again, the arrow indicates the direction of increasing temperature. Overall, there is clear consistency between the binding potential predictions and the DFT results for *δɛ*_{sf}/*ɛ* < 0.6, similar to the case of SR ff, LR sf interactions.

### C. Contour plots

As predicted in Sec. III, and observed in the profiles considered in Sec. IV A, there are regimes of parameter space for which the individual parameters *R*_{s}, *δμ*, and *ɛ*_{s} dominate the behavior of *ℓ*_{eq} and *χ*(*ℓ*_{eq}). Such behavior can be visualized more readily in the contour plots of Fig. 8. Figures 8(a) and 8(b) compare values of *ℓ*_{eq} obtained from DFT calculations for systems with SR ff, LR sf and with LR ff, LR sf interactions, respectively, for varying *R*_{s} and *δɛ*_{sf} at fixed *βδμ* = 10^{−3}. Figures 8(d) and 8(e) show the corresponding values of *χ*(*ℓ*_{eq}; *R*_{s}) normalized to the bulk values. The regions of parameter space in which individual parameters dominate are immediately apparent: When *R*_{s} is small, the contours are largely horizontal, indicating that changing *δɛ*_{sf}/*ɛ* has little influence on the behavior of *ℓ*_{eq} and *χ*(*ℓ*_{eq}; *R*_{s}). However, when *R*_{s} is large, the contours are almost vertical, indicating that *R*_{s} has little influence on the behavior of *ℓ*_{eq} and *χ*(*ℓ*_{eq}; *R*_{s}) in this region. From such plots, we can make numerical estimates of the crossover length scale for the change in behavior, say for a given choice of *δɛ*_{sf}, which was not possible from the binding potential analysis alone.

Figures 8(c) and 8(f) compare values of *ℓ*_{eq} and *χ*(*ℓ*_{eq}; *R*_{s}), respectively, for a system with SR ff, LR sf interactions at constant *δɛ*_{sf} = 0.0 and various *βδμ* and *R*_{s}. As was discussed in Secs. III A and IV A, for this case we expect the behavior of both quantities to depend almost solely on *R*_{s} when *R*_{s} < *R*_{c} and on *δμ* when *R*_{s} > *R*_{c}. The contours in Figs. 8(c) and 8(f) suggest this to be the case. As an example, we consider the case *βδμ* = 10^{−4} for which *R*_{c} ≈ 9134*σ*. Figures 8(c) and 8(f) confirm the crossover between different regimes indeed occurs around this value of *R*_{s}.

## V. RESULTS FROM GCMC SIMULATIONS OF THE mw SOLVENT

### A. Coexistence properties and simulation state points

mw is a relatively recent water model that has been shown to reproduce accurately many of the properties of water under ambient conditions while being much faster to implement than more established water models such as Extended Simple Point Charge (SPC/E).^{39} Recently, we have presented the first highly accurate liquid–vapor phase diagram for mw^{50} that we have measured via GCMC simulations.^{51} The temperature–density projection is reproduced in the upper panel of Fig. 9 and shows that mw has a critical temperature of 917.6 K, which exceeds that of water (647.1 K) by some 50%. Clearly, mw is not an accurate model for water at all temperatures. However, as shown in our previous work,^{50} mw appears to obey a law of corresponding states aligning with real water that other models such as TIP4P and SPC/E do not achieve to the same degree.^{50} Specifically, when the temperature–density phase diagram (coexistence curve) is scaled by the critical temperature and critical density, a data collapse is observed onto the similarly scaled phase diagram for real water, as shown in the lower panel of Fig. 9. This finding suggests that in seeking to study mw water under conditions equivalent to those of ambient real water, it is reasonable to employ the same scaled temperature as ambient water, i.e., *T*/*T*_{c} = 0.46. For mw, this corresponds to a simulation temperature of 426 K, which we adopt in our simulation studies below. We note that GCMC simulations of mw are substantially more computationally efficient at 426 K than at 300 K, although we have also considered the latter case as we mention below.

As mentioned in the Introduction, water at ambient conditions exhibits a small oversaturation^{30} of approximately *βδμ* ≈ 10^{−3}. When attempting to model accurately hydrophobic solvation, it is important to employ a realistic value of the oversaturation because this sets the deviation from liquid–vapor coexistence at which critical drying occurs. For near-critical planar systems, the magnitude of response functions depends strongly on the deviation from criticality and we expect the same to be true for solvophobic hydration at large spherical solutes. Our previously reported accurate measurements of the coexistence properties of mw in the *μ*-*T* plane^{50} allow us to control precisely the oversaturation for our model and thus impose a value appropriate to ambient water.

### B. Profiles of local density and compressibility

As for the solvophobic LJ case, understanding the influence of varying the parameters *R*_{s}, *δμ*, and *δɛ*_{sf} on the hydrophobic response of mw can be gained from measurements of *ρ*(*r*) and *χ*(*r*). Figures 10(a) and 10(b) show the variation with *R*_{s} and *ɛ*_{sf} at fixed *βδμ* = 10^{−3}, the value of the oversaturation pertaining to ambient water. The temperature is fixed at *T* = 426 K, which reproduces the value of *T*/*T*_{c} for ambient water as discussed in Sec. V A. Solutes of size *σ*_{mw} ≤ *R*_{s} ≤ 17*σ*_{mw} were studied corresponding to 0.239 25 nm ≤ *R*_{s} ≤ 4.067 25 nm in physical water units. This range was chosen to span the widely reported qualitative change in hydrophobic behavior that supposedly occurs around *R*_{s} ≈ 1 nm as discussed in Sec. I as well as to incorporate solutes whose size approaches that of small proteins. Owing to the large number of mw particles required, the computational effort required to study larger solutes becomes prohibitive.

The hydrophobic system considered here has interactions that correspond to SR ff, LR sf, see Sec. II B, and hence in the planar limit $Rs\u22121\u21920$, critical drying is expected to occur when *ɛ*_{sf} = 0.0. Profiles for values of *ɛ*_{sf} that correspond to realistic hydrophobic solutes^{53} are shown in Figs. 10(a) and 10(b). In all cases, for very small solutes, *R*_{s} < 1 nm, we find no region of depleted density and enhanced density fluctuations, in agreement with previous work.^{1} As *R*_{s} is increased beyond 1 nm, a region of depleted density begins to form and density fluctuations are enhanced on a similar length scale. The extent of the former and magnitude of the latter grow with *R*_{s}. Oscillations in the density profile, which are indicative of liquid packing effects, are dampened as the solute radius increases. Oscillations in the local compressibility profiles are no longer centered on *χ*_{b} as *R*_{s} increases. Note that similar behavior was observed in our DFT calculations for a solvophobic system, Sec. IV A. The apparent change in hydrophobic behavior occurring at around *R*_{s} = 1 nm is consistent with that reported in previous simulation studies of the density profiles of atomistic water models,^{1,9,54,55} strengthening our confidence that mw is a suitable water model for studies of hydrophobic solvation. While the limit of a planar substrate is not explored here, we note that a previous study of SPC/E water^{23} found very similar density and local compressibility profiles to those of Figs. 10(a) and 10(b). This reinforces further the connection between hydrophobicity on the macroscopic and microscopic length scales.

Although the majority of our result for mw were obtained at the fractional temperature *T*/*T*_{c} = 0.46 as water (which for mw corresponds to *T* = 426 K - see Sec. V A), we have also investigated other temperatures. Figures 10(c)–10(h) compare density and local compressibility profiles for three temperatures: 300 K, the ambient temperature for which mw was parameterized; 360 K, the temperature at which mw almost exactly reproduces the liquid–vapor surface tension of water; and 426 K. As the temperature is lowered, one observes that oscillations within both the density and local compressibility profiles become more prominent, and the magnitude of the maximum in *χ*(*r*) increases. For each temperature, the same general behavior is observed and mirrors that seen for a general solvophobic system as studied by DFT, Sec. IV A.

The binding potential analysis suggested three regimes within which individual parameters dominate the expected scaling behavior. The existence of these regimes was confirmed in DFT for a general solvophobic system by examining the density and local compressibility profiles calculated in Sec. IV A. Limitations on the radius of solute that can be studied via GCMC simulation prevent as full an exploration of such regimes for the hydrophobic (mw) case. However, it is still possible to confirm the existence of certain limiting behavior. Consider, for example, the local compressibility profiles *χ*(*r*) in Fig. 10(a). As *R*_{s} is reduced, the variation of *χ*(*r*) with *ɛ*_{sf} slows, which is similar to the angled horizontal contours of Fig. 8(d). For the hard solute, *ɛ*_{sf} = 0.0, we expect the behavior of the density and local compressibility to be almost solely dependent on *R*_{s} when *R*_{s} < *R*_{c}, as was shown to be the case for solvophobic systems in Figs. 3(a)–3(f). For mw at *T* = 426 K, *R*_{c} ≈ 0.583 *µ*m when *βδμ* = 10^{−3} and *R*_{c} ≈ 5.83 *µ*m when *βδμ* = 10^{−4}. Thus, for all values of *R*_{s} investigated in our simulations, we expect the density and local compressibility profiles to be almost independent of *βδμ*. Figure 11 confirms this to be the case: The profiles are near indistinguishable for *βδμ* = 10^{−3} and *βδμ* = 10^{−4}.

### C. Testing the binding potential predictions for *ℓ*_{eq} and *χ*(*ℓ*_{eq})

In contrast to DFT, our mw simulations cannot access the very large values of *R*_{s} that are required to verify fully the relationships predicted by our binding potential analysis: Eqs. (13) and (17). However, as discussed in Sec. III A, for the case *ɛ*_{sf} = 0.0, the behavior at small *R*_{s} is expected to be dominated by the Laplace $(Rs\u22121)$ term entering the effective pressure $p\u0303$, so we predict that *ℓ*_{eq} ∼ ln(*R*_{s}) and *χ*(*ℓ*_{eq}; *R*_{s}) ∼ *R*_{s} [see Eqs. (14) and (18)]. The binding potential analysis predicts that this scaling should also hold when *ɛ*_{sf} is sufficiently small. These predictions are tested and verified in Figs. 12 and 13 for several vales of *ɛ*_{sf}. The fact that the predicted scaling behavior is observed when *R*_{s} > 1 nm lends weight to our assertion that the critical drying point controls the properties of microscopic hydrophobicity on sufficiently large length scales. However, it is also striking that the predicted scaling of *ℓ*_{eq} and *χ*(*ℓ*_{eq}; *R*_{s}) appears to work quite well down to small solute radii where *ℓ*_{eq} is calculated to be a small fraction of the water diameter *σ*_{mw}. In this regime, there is no discernible vapor “film.” Rather, there are regions of density depletion, shown in Figs. 10 and 11, that extend across only short distances from the surface of the solute.

### D. Nature of the depleted density region

Experimental studies of hydrophobicity at a planar substrate have revealed the presence of a region of depleted water density at distances within a few molecular diameters of the substrate. However, the precise extent of this region remains controversial. X-ray reflectivity studies measure only the net depletion of electron density. Results and interpretation remain subject to debate.^{18–21} Atomic force microscopy studies^{56,57} appear to provide evidence that the depleted density region takes the form of “nanobubbles.” The formation of microscopic bubbles very close to the drying transition was observed in simulation studies of a Lennard-Jones fluid at a solvophobic planar substrate^{25} in which a rich fractal-like bubble structure was observed. As our present results suggest a common connection between features of hydrophobicity on macroscopic and microscopic length scales, it is interesting to consider the nature of the mw water structure at the surface of a hydrophobic solute.

Figure 14 presents simulation snapshots of cross sections through the center of the simulation box for six solutes of various sizes. In each case, the solute is “hard,” i.e., *ɛ*_{sf} = 0.0. mw solvent particles are shown in blue, with the shade of blue representing the depth of the particle from the foreground—lighter particles are farther away. In all cases, the size ratio of the mw particles to the solute is to scale.

As *R*_{s} increases, bubbles (low-density regions) form at the surface of the solute and hence our simulations are in line with the nanobubble picture of hydrophobicity.^{56,57} Bubbles occur over a large extent of the solute surface, particularly when *R*_{s} > 9*σ*_{mw}. Note that the bubbles are localized very close to the solute; they do not seem to extend far from the surface. This observation accords with the general expectation that correlations parallel to the substrate diverge faster than in the perpendicular direction on approaching the drying critical point.^{25} Although not presented here, the bubbles move across the surface of the solute frequently during the course of a simulation, changing in location and spatial extent. Increasing the solute–solvent attraction, *ɛ*_{sf}, has the effect of suppressing the bubbles and reducing their size.

We note that the complex bubble nanostructure that we observe in GCMC simulations is, of course, not accessible to either the binding potential or DFT calculations. The former coarse grains into a single variable *l*, the thickness of the intruding vapor film, and the latter averages over all the bubble configurations to yield an average density profile *ρ*(*r*) reflecting the average depletion arising from nanobubble formation. Nevertheless, one expects on general grounds^{25} that such theories will capture the correct large length scaling behavior for the thickness of the depleted region and the magnitude of the local compressibility. This appears to be the case, as borne out by the results of our simulations of mw.

## VI. DISCUSSION AND CONCLUSIONS

A detailed understanding at the molecular level of the behavior of water in the vicinity of a hydrophobic solute is an important goal in areas ranging from solution chemistry to biophysics. For the case of an extended strongly hydrophobic solute, previous work has reported a region of depleted density and enhanced density fluctuations in water near the solute’s surface. However, the physical origin of these effects has remained obscure. Progress in understanding the phenomena has been hindered by the use of imprecise nomenclature. Previous work points to the proximity of ambient water to liquid–vapor coexistence and implies that this leads to a “dewetting transition”^{4,5,9} or sometimes a “drying transition”^{1,14} of water near a hydrophobic solute or in the region between two such solutes.^{58} However, these terms seem to be merely shorthand for the appearance of a region of depleted water density; the precise definition of the “transitions” alluded to and their relevance to the phenomenology of hydrophobicity was not clarified. Here, we are careful to define precisely what is a drying transition and, in particular, what identifies a critical drying transition and why this is important for phenomena associated with hydrophobic and, more generally, solvophobic solutes.

In the present work, we have studied in detail the density depletion and enhanced compressibility close to a model spherical solute employing a combination of meso- and microscopic theoretical and computational methods. We hypothesized that these phenomena are attributable to the critical drying surface phase transition that occurs at liquid–vapor coexistence for a very weakly interacting solute in the limit of an infinite solute radius, i.e., a planar substrate. Quite generally on approaching a critical point, we expect the strength of fluctuations to grow with increase in the correlation length and to diverge precisely at criticality. It follows that enhancement of (density) fluctuations should be observed in a significant region of parameter space surrounding a critical point. For typical experimental systems in which water under standard temperature and pressure is in contact with an extended strongly hydrophobic solute, the oversaturation Δ*μ*, the solute curvature $Rs\u22121$, and the solute–solvent attraction *ɛ*_{sf} are all sufficiently small for the system to qualify as “near” to the critical drying point. Accordingly, we might expect the system’s behavior to be driven by enhanced local density fluctuations and, possibly, the emergence of a surface vapor region. Suitable measures should exhibit near-critical scaling behavior.

In order to investigate this proposal and make quantitative predictions, we constructed a mesoscopic binding potential theory that allows us to explore how the size of the solute, its interaction strength with water, and the degree of oversaturation determine the extent of density depletion, as measured by *ρ*(*r*), and the magnitude of local density fluctuations, as measured by *χ*(*r*). The resulting mean-field scaling predictions were tested via classical DFT calculations for a generalized solute–solvent system. From the evidence presented in the profiles of *ρ*(*r*) and *χ*(*r*) in Fig. 3, the comparison of binding potential predictions to DFT results in Figs. 4–7 and the contour plots of Fig. 8, it is clear that the predictions of the macroscopic binding potential analysis for the near-critical scaling of *ℓ*_{eq} and *χ*(*ℓ*_{eq}; *R*_{s}) are in agreement with DFT results for a wide variety of solvophobic systems that extend down to microscopic solute sizes. For microscopic solutes (*R*_{s} ≲ 10^{3}*σ*), our results suggest that the magnitude of local density fluctuations near a microscopic solute is most sensitive to changes in *R*_{s}, with small variations in *T*, *δμ*, and solvophobicity, as measured by *ɛ*_{sf}, having limited effects. This observation is pertinent with regard to proteins—the size ratio of a (small) protein to a water molecule is such that scaling behavior would be expected to fall in the curvature dominated regime. In turn, this implies that small variations in other parameters, such as temperature, would have little effect on the hydrophobic behavior, e.g., density depletion and fluctuations. As density fluctuations are sometimes conjectured to facilitate protein folding, this insight might potentially prove useful for understanding protein folding.

Our binding potential scaling predictions were tested further via GCMC simulation studies of a monatomic water model. Although the limited range of solute radii accessible to simulation permitted a less comprehensive test than for DFT, principal aspects of the scaling were nevertheless verified. The agreement between simulations and the binding potential and DFT predictions provides confirmation of the expectation that mean-field scaling should apply in such systems.^{25,35}

The simulations provide additional molecular-level insight into the nature of the local configurational structure of water near the solute surface and how this engenders the enhancement in local compressibility. We found that elongated vapor bubbles form at the solute surface, whose position and size fluctuate strongly during the course of a simulation run. The bubbles become larger and more distinct with increasing solute radius in a situation reminiscent of what is observed for a planar substrate.^{25} Thus, despite the impression conveyed by simulation measurements of the density profile *ρ*(*r*) (cf. Fig. 10), hydrophobicity does not immediately lead to the emergence of a smooth liquid–vapor-like interface around a large solute, at least not unless $Rs\u22121,\delta \mu ,\delta \epsilon sf$ are all sufficiently small that the largest bubbles encompass much of the solute’s surface area. Of course, here we shall be close to critical drying. Furthermore, our simulation results show no evidence of the strongly hydrogen bonded hydration shell that has been proposed to form around small hydrophobic solutes, e.g., Ref. 7. Should such a structure develop, we believe this can only occur for solutes smaller than those we studied here.

Taken together, we believe that our results shed new light on the nature and origin of hydrophobic solvation phenomena and provide a firm basis for rationalizing how properties on microscopic length scales depend on the solute size and the strength of solute–water attraction. As indicated above, in future work it would be interesting to investigate whether the insights gained here might facilitate an improved understanding of physical processes near hydrophobic entities that are believed to be mediated by local density fluctuations, such as those occurring in protein dynamics.

## ACKNOWLEDGMENTS

This work used the facilities of the Advanced Computing Research Centre, University of Bristol. We thank F. Turci for valuable discussions. R.E. acknowledges the support received under Leverhulme Trust Grant No. EM-2020-029\4.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Mary Kathryn Coe**: Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal). **Robert Evan**: Conceptualization (equal); Investigation (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). **Nigel B. Wilding**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

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

## REFERENCES

*in situ*x-ray study of the hydrophobic gap at the water–octadecyl-trichlorosilane interface

*Fundamentals of Inhomogeneous Fluids*

*Phase Transitions and Critical Phenomena*

*Les Houches 1988 Session XLVIII Liquids at Interfaces*

*Fluid Interfacial Phenomena*

*NIST Chemistry WebBook, NIST Standard Reference Database Number 69*

We note that the term “dewetting” is most commonly used in a rather different context of the rupturing of a thin film of liquid on a planar substrate.