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 Rs. 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 Rs, 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.

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 micelles9 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, ε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 ε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 fluctuations23,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 εsfεsfd from above.25, ξ is in turn linked to an important spatial measure of density fluctuations, the local compressibility profile,27 defined as χ(z)=ρ(z)/μT. Specifically, one can show25 that for distances z in the interfacial region χ(z)ρ(z)ξ2, where ρ′ is the gradient of the density profile. This implies that the maximum in χ(z) also diverges as εsfε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 εsfε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 χ(r)=ρ(r)/μ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 Investigations29 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 εsf=ε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/kBT. Moreover, in experiments it is not currently possible to realize substrates that correspond to the extreme hydrophobic limit, i.e., which have δεsfεsfε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 Rs. By treating the solute curvature Rs1 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 Rs), the system is near-critical and provides specific predictions for the scaling properties of eq and the maximum of χ(r)=ρ(r)/μT in terms of Rs, βδμ, ɛ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 

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.

This system is a direct extension of that used in a previous study of solvophobic planar substrates.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 Rs 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 form25,36,37
(1)
where r = |rr′| and where r and r′ denote the position vectors of a pair of particles. rmin = 21/6σ is the distance corresponding to the minimum, rc is the cutoff radius, σ is the LJ diameter, and ɛ is the well depth of the LJ interaction.
FIG. 1.

Sketch of the model system used in DFT and simulation. A smooth and impenetrable solute of radius Rs is centered on the origin and composed of smaller particles distributed homogeneously with density ρs. Solvent particles are represented by the blue circles of radius R. The smaller particles comprising the solute interact with the solvent particles via a 6–12 LJ potential of well depth ɛs and diameter σs. The net external (radial) potential experienced by solvent particles outside the impenetrable zone, Vext(r), is shown by the dashed black line. The radially symmetric density profile of the fluid is shown by the solid blue line. The first nonzero density or contact density, ρ(Rs+), occurs at a distance Rs as indicated.

FIG. 1.

Sketch of the model system used in DFT and simulation. A smooth and impenetrable solute of radius Rs is centered on the origin and composed of smaller particles distributed homogeneously with density ρs. Solvent particles are represented by the blue circles of radius R. The smaller particles comprising the solute interact with the solvent particles via a 6–12 LJ potential of well depth ɛs and diameter σs. The net external (radial) potential experienced by solvent particles outside the impenetrable zone, Vext(r), is shown by the dashed black line. The radially symmetric density profile of the fluid is shown by the solid blue line. The first nonzero density or contact density, ρ(Rs+), occurs at a distance Rs as indicated.

Close modal
These pairwise potentials give rise to an attractive potential at each point in the solvent whose form is derived35 by integrating over the angular degrees of freedom in which the fluid has homogeneous density to yield
(2)
This effective attractive potential is then incorporated in the DFT calculations (see Sec. IV) in the standard mean-field fashion.35,38
Integrating the (virtual) solute particle–fluid particle LJ pair potential, with diameter σs and well depth ɛs, over the volume of the solute35 gives rise to the net solute–solvent external potential
(3)
where εsf=2πρsεsσs3/3 is the effective solute–fluid attraction strength, Rs and r are as shown in Fig. 1, and r+ = r + rmin,sf, where rmin,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 rc = 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 rc = 200σ. Such a large truncation distance is sufficient to capture accurately the LR decay.

Our GCMC simulations utilize a popular monatomic water (mw) model proposed by Molinero and Moore.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 potential39 
(4)
where the two-body, ϕmw,2, and three-body, ϕmw,3, potentials are
(5)
(6)
and 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.

The present work considers a periodic cubic simulation box of mw particles that contains a spherical solute of radius Rs centered on the origin as in Fig. 1. The mw particles interact with the solute via a potential having the same form of Eq. (3) used in our DFT calculations.

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 surfaces24,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 Rs10. 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.

Our model system comprises a solute of radius Rs 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 < Tc. 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
(7)
where ρ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 Rs is the effective radius and ρ(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).
FIG. 2.

Illustration of the sharp-kink (SK) approximation for the density profile employed in the binding potential analysis. (a) For a planar substrate, the system is assumed to consist of three slabs, representing the substrate (gray), vapor (light blue), and liquid (darker blue). The densities of the vapor and liquid are assumed to take their coexistence values ρv and ρl, respectively. A region of excluded volume of width dw exists between the solute and vapor arising from the finite size of the particles. The density of the fluid is measured from the center of the fluid particle and the width of the vapor slab is . (b) The system is assumed to consist of a large solute (gray), a vapor shell (light blue) of width , and liquid (darker blue). Again, a region of excluded volume (white) of width dw exists between the solute and vapor, leading to an effective solute radius Rs.

FIG. 2.

Illustration of the sharp-kink (SK) approximation for the density profile employed in the binding potential analysis. (a) For a planar substrate, the system is assumed to consist of three slabs, representing the substrate (gray), vapor (light blue), and liquid (darker blue). The densities of the vapor and liquid are assumed to take their coexistence values ρv and ρl, respectively. A region of excluded volume of width dw exists between the solute and vapor arising from the finite size of the particles. The density of the fluid is measured from the center of the fluid particle and the width of the vapor slab is . (b) The system is assumed to consist of a large solute (gray), a vapor shell (light blue) of width , and liquid (darker blue). Again, a region of excluded volume (white) of width dw exists between the solute and vapor, leading to an effective solute radius Rs.

Close modal
The excess grand potential Ωex for such a system can be written as41,45
(8)
where γ is the surface tension and A the surface area of the solute–vapor (subscript sv) and liquid–vapor (subscript lv) interfaces, δμ = (μμco), Δρ = (ρlρv), Vv is the volume of the vapor and ω(|Rs) 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(Rs) and γlv(Rs + ) can be approximated by their planar surface equivalents, γsv and γlv, respectively. Furthermore, we note that it has been shown previously41 that ω(|Rs)=ω()(1+/Rs+O(2ln(/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πRs2, Alv=4π(Rs+)2, and Vv=4π(Rs+)3/3 gives
(9)
where under the large solute approximation, terms of order O(2/Rs2) and higher have been neglected. Here, p̃ is the effective pressure, defined as
(10)

As emphasized in our short report,34  p̃ 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 Rs 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.

For a system with SR ff, LR sf interactions, the binding potential ω() takes the form24–26 
(11)
where 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
(12)
with bo = πΔρ/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).
Substituting into Eq. (9) and minimizing with respect to yields an expression for the equilibrium vapor layer thickness, eq,
(13)
which reduces to the expression for a planar surface26 in the limit Rs10.
Including the Rs 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 eqlnp̃. 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 Rc = 2γlv/δμΔρ,
(14)
Similarly, at bulk coexistence δμ = 0, two regimes of scaling can be identified,
(15)
where the latter corresponds to the planar case mentioned in the work of Evans et al.25 These results, which pertain to different physical limiting cases, highlight that there are regions of the parameter space (δμ, ɛs, Rs) 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 Rs = , 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 − (Rs + )).

The local compressibility at eq can then be found using25 
(16)
where ρ′ is the spatial derivative of the profile. Substituting eq given in Eq. (13) yields
(17)
As expected, in the limit Rs10, this reduces to the result for a planar substrate.
Similarly to eq, we can identify three regimes in which the behavior of χ(eq|Rs) is controlled predominately by Rs, δμ, or ɛs. For a hard solute, where ɛs = 0, we find
(18)
where the latter is the result for drying at a hard planar wall, e.g., see the work of Evans and Stewart.27 At bulk coexistence δμ = 0, we find
(19)
where the latter was identified, for a planar wall, by Evans et al.25 
For the case of LR ff, LR sf interactions, the binding potential takes the form26 and references therein,
(20)
where
(21)
where bo = πΔρ/3, as defined previously, and zmin = (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 bLR(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 bLR(T) determines the location of the minimum in ωex() that corresponds to the equilibrium film width eq. In particular, the drying temperature Td, at which the thickness of the vapor film eq, is determined by the condition bLR(Td) = 0. If one fixes the temperature, it follows that the critical drying point occurs for substrate–fluid attraction strength given by
(22)
as derived previously.26 Summarizing, critical drying for a system with LR ff, LR sf interactions occurs at bulk coexistence δμ = 0 in the planar limit Rs = when the attraction strength is given by Eq. (22).
The variation of bLR(T) in the approach to the critical drying point from below allows us to define a dimensionless measure of the deviation from this point,
(23)
which vanishes at critical drying and which we term the effective reduced temperature. This definition is similar to that adopted by Stewart and Evans.40 The convenience stems from the fact that bLR(T) ∝ t′.
Turning to the case of a spherical solute, we substitute Eq. (20) into Eq. (9) and then minimizing w.r.t. gives
(24)
Once again, this equation reduces to that for a planar surface as Rs10. As in the SR ff, LR sf case, it is possible to identify three regimes in which Rs, δμ, and ɛs individually dominate the behavior of eq.
Differentiating Eq. (24) yields an expression for χ(eq|Rs),
(25)
For the case of LR ff, LR sf interactions, it has been shown previously40 that eq can be expressed as a scaling function. An appropriate form is
(26)
where L obeys Eq. (24). Adopting this result, χ(eq|Rs) can be written as
(27)
where L is the derivative of the scaling function. This result obeys Eq. (25). Scaling forms for this case of LR ff, LR sf interactions are particularly useful owing to the difficulty in calculating c(T). Note that while the result for bLR(T) remains valid beyond the sharp-kink approximation, this is not the case for c(T).46 

Our DFT calculations47 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 rc = 2.5σ and left unshifted and (ii) the LR ff case in which the true long-ranged potential is approximated by taking rc = 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.775Tc 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 εsf=2πρsεsσs3/3, the effective solute–fluid attraction strength, δμ, T, and Rs, 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.

Figures 3(a)3(f) demonstrate the effect on ρ(r) and χ(r) of varying Rs (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 Rs, ρ(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.

FIG. 3.

(a)–(l) DFT results for ρ(r) and χ(r), normalized to their bulk values, for the case of SR ff, LR sf; (m)–(r) the corresponding results for LR ff, LR sf. Top row: Varying Rs and βδμ at constant δɛsf = 0.0. Middle row: Varying δɛsf and Rs at constant βδμ = 10−3. Bottom row: Varying δɛsf and Rs at constant βδμ = 10−3 for LR ff, LR sf.

FIG. 3.

(a)–(l) DFT results for ρ(r) and χ(r), normalized to their bulk values, for the case of SR ff, LR sf; (m)–(r) the corresponding results for LR ff, LR sf. Top row: Varying Rs and βδμ at constant δɛsf = 0.0. Middle row: Varying δɛsf and Rs at constant βδμ = 10−3. Bottom row: Varying δɛsf and Rs at constant βδμ = 10−3 for LR ff, LR sf.

Close modal

These observations are in accord with the binding potential analysis of Sec. III A, which predicts, for the hard solute, a length scale Rc = 2γlv/δμΔρ that separates behavior dependent on and independent of the curvature of the solute: For solutes having RsRc, 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, Rc ≈ 913σ and hence we would expect profiles for Rs > 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 Rc = 9134σ and Rc = 91 340σ, respectively.

In Figs. 3(g)3(l), βδμ is held constant at 10−3 while varying δɛsf (colors) and Rs (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 Rs. 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 107.4° 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 Rs 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 Rc, which for this case is Rc ≈ 1138σ. Figures 3(o)3(r) compare the influence of varying δɛsf/ɛ and Rs 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.

From the density and local compressibility profiles obtained from DFT, it is possible to extract eq and thus χ(eqRs). To do so, we define eq in the standard way, e.g.,49 
(28)
where Γ is the Gibbs excess adsorption, obtained from the calculated density profiles using numerical integration as follows:
(29)
in the case of a curved substrate/solute and
(30)
in the case of a planar substrate. In each case, ρb is the density of the bulk liquid. χ(eqRs) 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 χ(eqRs) for a large range of values of δμ, ɛsf, Rs 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σRs ≤ 108σ). 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.

FIG. 4.

Comparison of DFT results for eq obtained using Eq. (28) and the (linear) scaling prediction of Eq. (13), for the case of SR ff, LR sf interactions. The temperature is fixed at T = 0.775Tc, while (10−6βδμ ≤ 10−3, 0.0 ≤ δɛsf/ɛ ≤ 1.0, 10σRs ≤ 108σ). The value of δɛsf and the associated contact angle θ for each result are indicated by the color bar.

FIG. 4.

Comparison of DFT results for eq obtained using Eq. (28) and the (linear) scaling prediction of Eq. (13), for the case of SR ff, LR sf interactions. The temperature is fixed at T = 0.775Tc, while (10−6βδμ ≤ 10−3, 0.0 ≤ δɛsf/ɛ ≤ 1.0, 10σRs ≤ 108σ). The value of δɛsf and the associated contact angle θ for each result are indicated by the color bar.

Close modal

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 Rs < 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°, 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 χ(eqRs) 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 Rs < 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 χ(eqRs) 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.

FIG. 5.

Comparison of DFT results for χ(eqRs) and the (linear) scaling prediction of Eq. (17), for the case of SR ff, LR sf interactions. The temperature is fixed at T = 0.775Tc, while (10−6βδμ ≤ 10−3, 0.0 ≤ δɛsf/ɛ ≤ 0.3, 10σRs ≤ 108σ). The value of δɛsf and the associated contact angle θ for each result are indicated by the color bar.

FIG. 5.

Comparison of DFT results for χ(eqRs) and the (linear) scaling prediction of Eq. (17), for the case of SR ff, LR sf interactions. The temperature is fixed at T = 0.775Tc, while (10−6βδμ ≤ 10−3, 0.0 ≤ δɛsf/ɛ ≤ 0.3, 10σRs ≤ 108σ). The value of δɛsf and the associated contact angle θ for each result are indicated by the color bar.

Close modal

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σRs ≤ 108σ 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.7Tc, 0.775Tc, and 0.85Tc, 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.

FIG. 6.

DFT results for eq plotted according to the predicted scaling form of Eq. (26) for LR ff, LR sf interactions. Results are shown for three temperatures, T = 0.7Tc, 0.775Tc, and 0.85Tc, with the red arrow indicating the direction of increasing temperature. In each case, βδμ = 10−3, while 0.0 ≤ δɛsf/ɛ ≤ 1.0, as indicated by the color bar and 10σRs ≤ 108σ. The inset shows the contact angle that would pertain in the planar limit for each δɛsf and T considered.

FIG. 6.

DFT results for eq plotted according to the predicted scaling form of Eq. (26) for LR ff, LR sf interactions. Results are shown for three temperatures, T = 0.7Tc, 0.775Tc, and 0.85Tc, with the red arrow indicating the direction of increasing temperature. In each case, βδμ = 10−3, while 0.0 ≤ δɛsf/ɛ ≤ 1.0, as indicated by the color bar and 10σRs ≤ 108σ. The inset shows the contact angle that would pertain in the planar limit for each δɛsf and T considered.

Close modal

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 interface46 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 Rs < 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 χ(eqRs) 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.

FIG. 7.

DFT results for χ(eqRs) plotted according to the predicted scaling form of Eq. (27) for LR ff, LR sf interactions. The inset makes the comparison for an alternative scaling form of Eq. (25). In each case, βδμ = 10−3, while 0.0 ≤ δɛsf/ɛ ≤ 1.0, as indicated by the color bar, and 10σRs ≤ 108σ. Results are shown for three temperatures, T = 0.7Tc, 0.775Tc, and 0.85Tc, with the red arrows indicating the direction of increasing temperature. Corresponding values of the contact angle can be found in the inset of Fig. 6.

FIG. 7.

DFT results for χ(eqRs) plotted according to the predicted scaling form of Eq. (27) for LR ff, LR sf interactions. The inset makes the comparison for an alternative scaling form of Eq. (25). In each case, βδμ = 10−3, while 0.0 ≤ δɛsf/ɛ ≤ 1.0, as indicated by the color bar, and 10σRs ≤ 108σ. Results are shown for three temperatures, T = 0.7Tc, 0.775Tc, and 0.85Tc, with the red arrows indicating the direction of increasing temperature. Corresponding values of the contact angle can be found in the inset of Fig. 6.

Close modal

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 Rs, δμ, 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 Rs and δɛsf at fixed βδμ = 10−3. Figures 8(d) and 8(e) show the corresponding values of χ(eqRs) normalized to the bulk values. The regions of parameter space in which individual parameters dominate are immediately apparent: When Rs is small, the contours are largely horizontal, indicating that changing δɛsf/ɛ has little influence on the behavior of eq and χ(eqRs). However, when Rs is large, the contours are almost vertical, indicating that Rs has little influence on the behavior of eq and χ(eqRs) 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.

FIG. 8.

Contour plots of DFT results showing the dependence of eq and χ(eq) on Rs1 and δɛsf. The temperature T = 0.775Tc is fixed. In plots (a), (b), (d), and (e), βδμ = 10−3 and the contact angle that pertains in the planar limit at bulk coexistence for a given δɛsf is included. (a) eq for SR ff interactions, (b) eq for LR ff, (d) χ(eq) for SR ff, and (e) χ(eq) for LR ff. (c) eq for SR ff with δɛsf = 0.0 and increasing βδμ. (f) χ(eq) for SR ff with δɛsf = 0.0 and increasing βδμ.

FIG. 8.

Contour plots of DFT results showing the dependence of eq and χ(eq) on Rs1 and δɛsf. The temperature T = 0.775Tc is fixed. In plots (a), (b), (d), and (e), βδμ = 10−3 and the contact angle that pertains in the planar limit at bulk coexistence for a given δɛsf is included. (a) eq for SR ff interactions, (b) eq for LR ff, (d) χ(eq) for SR ff, and (e) χ(eq) for LR ff. (c) eq for SR ff with δɛsf = 0.0 and increasing βδμ. (f) χ(eq) for SR ff with δɛsf = 0.0 and increasing βδμ.

Close modal

Figures 8(c) and 8(f) compare values of eq and χ(eqRs), respectively, for a system with SR ff, LR sf interactions at constant δɛsf = 0.0 and various βδμ and Rs. 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 Rs when Rs < Rc and on δμ when Rs > Rc. 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 Rc ≈ 9134σ. Figures 8(c) and 8(f) confirm the crossover between different regimes indeed occurs around this value of Rs.

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 mw50 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/Tc = 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.

FIG. 9.

(a) Upper panel: liquid–vapor phase diagram (binodal) of mw as measured by GCMC simulation. Statistical uncertainties are smaller than the data points. The critical point is marked by a cross and the dotted line denotes the coexistence diameter. Both physical and reduced units are given. In the latter, σmw = 2.3925 Å, while ɛmw = 6.189 kcal/mol. (b) Lower panel: comparison of the liquid–vapor phase diagrams of water (pink squares) and mw (blue circles) with temperature T and density ρ scaled by their bulk critical values. The dotted gray line indicates the value of T/Tc for water at ambient conditions. Data for water are taken from Ref. 52.

FIG. 9.

(a) Upper panel: liquid–vapor phase diagram (binodal) of mw as measured by GCMC simulation. Statistical uncertainties are smaller than the data points. The critical point is marked by a cross and the dotted line denotes the coexistence diameter. Both physical and reduced units are given. In the latter, σmw = 2.3925 Å, while ɛmw = 6.189 kcal/mol. (b) Lower panel: comparison of the liquid–vapor phase diagrams of water (pink squares) and mw (blue circles) with temperature T and density ρ scaled by their bulk critical values. The dotted gray line indicates the value of T/Tc for water at ambient conditions. Data for water are taken from Ref. 52.

Close modal

As mentioned in the Introduction, water at ambient conditions exhibits a small oversaturation30 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 plane50 allow us to control precisely the oversaturation for our model and thus impose a value appropriate to ambient water.

As for the solvophobic LJ case, understanding the influence of varying the parameters Rs, δμ, 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 Rs 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/Tc for ambient water as discussed in Sec. V A. Solutes of size σmwRs ≤ 17σmw were studied corresponding to 0.239 25 nm ≤ Rs ≤ 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 Rs ≈ 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.

FIG. 10.

GCMC results for mw. (a) local compressibility χ(r) normalized by its bulk value for various values of solute radius Rs, as indicated on the abscissa, and various values of solute–fluid attraction ɛsf/ɛmw as indicated in the key. In units of kcal/mol, these values are from lower to higher: ɛsf = 0.0, 0.030 95, 0.123 78, 0.247 56. The temperature T = 426 K and supersaturation βδμ = 10−3 are fixed. (b) The corresponding local density profile ρ(r) normalized by its bulk value. (c)–(h) Variation of χ(r) and ρ(r) with temperature T for ɛsf = 0.0 and various Rs as indicated in the key. The supersaturation is βδμ = 10−3.

FIG. 10.

GCMC results for mw. (a) local compressibility χ(r) normalized by its bulk value for various values of solute radius Rs, as indicated on the abscissa, and various values of solute–fluid attraction ɛsf/ɛmw as indicated in the key. In units of kcal/mol, these values are from lower to higher: ɛsf = 0.0, 0.030 95, 0.123 78, 0.247 56. The temperature T = 426 K and supersaturation βδμ = 10−3 are fixed. (b) The corresponding local density profile ρ(r) normalized by its bulk value. (c)–(h) Variation of χ(r) and ρ(r) with temperature T for ɛsf = 0.0 and various Rs as indicated in the key. The supersaturation is βδμ = 10−3.

Close modal

The hydrophobic system considered here has interactions that correspond to SR ff, LR sf, see Sec. II B, and hence in the planar limit Rs10, critical drying is expected to occur when ɛsf = 0.0. Profiles for values of ɛsf that correspond to realistic hydrophobic solutes53 are shown in Figs. 10(a) and 10(b). In all cases, for very small solutes, Rs < 1 nm, we find no region of depleted density and enhanced density fluctuations, in agreement with previous work.1 As Rs 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 Rs. 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 Rs 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 Rs = 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 water23 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/Tc = 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 Rs 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 Rs when Rs < Rc, as was shown to be the case for solvophobic systems in Figs. 3(a)3(f). For mw at T = 426 K, Rc ≈ 0.583 µm when βδμ = 10−3 and Rc ≈ 5.83 µm when βδμ = 10−4. Thus, for all values of Rs 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.

FIG. 11.

GCMC results for density (lower) and local compressibility (upper) profiles, scaled by their bulk values, for mw at T = 426 K and ɛsf = 0.0, and for various Rs and two values of βδμ. In the case of the density profile, uncertainties do not exceed linewidth.

FIG. 11.

GCMC results for density (lower) and local compressibility (upper) profiles, scaled by their bulk values, for mw at T = 426 K and ɛsf = 0.0, and for various Rs and two values of βδμ. In the case of the density profile, uncertainties do not exceed linewidth.

Close modal

In contrast to DFT, our mw simulations cannot access the very large values of Rs 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 Rs is expected to be dominated by the Laplace (Rs1) term entering the effective pressure p̃, so we predict that eq ∼ ln(Rs) and χ(eqRs) ∼ Rs [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 Rs > 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 χ(eqRs) 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.

FIG. 12.

Results for eq obtained from GCMC simulations of mw for fixed T = 426 K and βδμ = 10−3; ɛsf is given by the symbol. Each dotted line is a linear fit to the scaling prediction.

FIG. 12.

Results for eq obtained from GCMC simulations of mw for fixed T = 426 K and βδμ = 10−3; ɛsf is given by the symbol. Each dotted line is a linear fit to the scaling prediction.

Close modal
FIG. 13.

Results for χ(eqRs) obtained from GCMC simulations of mw for fixed T = 426 K and βδμ = 10−3; ɛsf is given by the symbol. Each dotted line is a linear fit to the scaling prediction.

FIG. 13.

Results for χ(eqRs) obtained from GCMC simulations of mw for fixed T = 426 K and βδμ = 10−3; ɛsf is given by the symbol. Each dotted line is a linear fit to the scaling prediction.

Close modal

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 studies56,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 substrate25 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.

FIG. 14.

Snapshots of mw at hard solutes, ɛsf = 0.0; T = 426 K, βδμ = 10−3. (a) Rs = σmw, (b) Rs = 5σmw, (c) Rs = 9σmw, (d) Rs = 13σmw, (e) Rs = 15σmw, (f) Rs = 17σmw. As the radius Rs increases, one observes the formation of more “nanobubbles” next to the solute.

FIG. 14.

Snapshots of mw at hard solutes, ɛsf = 0.0; T = 426 K, βδμ = 10−3. (a) Rs = σmw, (b) Rs = 5σmw, (c) Rs = 9σmw, (d) Rs = 13σmw, (e) Rs = 15σmw, (f) Rs = 17σmw. As the radius Rs increases, one observes the formation of more “nanobubbles” next to the solute.

Close modal

As Rs 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 Rs > 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 grounds25 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.

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 Rs1, 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. 47 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 χ(eqRs) are in agreement with DFT results for a wide variety of solvophobic systems that extend down to microscopic solute sizes. For microscopic solutes (Rs ≲ 103σ), our results suggest that the magnitude of local density fluctuations near a microscopic solute is most sensitive to changes in Rs, 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 Rs1,δμ,δε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.

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.

The authors have no conflicts to disclose.

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).

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

1.
D. M.
Huang
and
D.
Chandler
, “
The hydrophobic effect and the influence of solute-solvent attractions
,”
J. Phys. Chem. B
106
,
2047
2053
(
2002
).
2.
A.
Oleinikova
and
I.
Brovchenko
, “
Thermodynamic properties of hydration water around solutes: Effect of solute size and water–solute interaction
,”
J. Phys. Chem. B
116
,
14650
14659
(
2012
).
3.
J.
Mittal
and
G.
Hummer
, “
Static and dynamic correlations in water at hydrophobic interfaces
,”
Proc. Natl. Acad. Sci. U. S. A.
105
,
20130
20135
(
2008
).
4.
S.
Sarupria
and
S.
Garde
, “
Quantifying water density fluctuations and compressibility of hydration shells of hydrophobic solutes and proteins
,”
Phys. Rev. Lett.
103
,
037803
(
2009
).
5.
A. J.
Patel
,
P.
Varilly
,
S. N.
Jamadagni
,
M. F.
Hagan
,
D.
Chandler
, and
S.
Garde
, “
Sitting at the edge: How biomolecules use hydrophobicity to tune their interactions and function
,”
J. Phys. Chem. B
116
,
2498
2503
(
2012
).
6.
S.
Vaikuntanathan
,
G.
Rotskoff
,
A.
Hudson
, and
P. L.
Geissler
, “
Necessity of capillary modes in a minimal model of nanoscale hydrophobic solvation
,”
Proc. Natl. Acad. Sci. U. S. A.
113
,
E2224
E2230
(
2016
).
7.
I.
Bischofberger
,
D. C. E.
Calzolari
,
P.
De Los Rios
,
I.
Jelezarov
, and
V.
Trappe
, “
Hydrophobic hydration of poly-N-isopropyl acrylamide: A matter of the mean energetic state of water
,”
Sci. Rep.
4
,
4377
(
2014
).
8.
N. B.
Rego
and
A. J.
Patel
, “
Understanding hydrophobic effects: Insights from water density fluctuations
,”
Annu. Rev. Condens. Matter Phys.
13
,
303
324
(
2022
).
9.
D.
Chandler
, “
Interfaces and the driving force of hydrophobic assembly
,”
Nature
437
,
640
647
(
2005
).
10.
J.
Qvist
,
M.
Davidovic
,
D.
Hamelberg
, and
B.
Halle
, “
A dry ligand-binding cavity in a solvated protein
,”
Proc. Natl. Acad. Sci. U. S. A.
105
,
6296
6301
(
2008
).
11.
K.
Lum
,
D.
Chandler
, and
J. D.
Weeks
, “
Hydrophobicity at small and large length scales
,”
J. Phys. Chem. B
103
,
4570
4577
(
1999
).
12.
N. T.
Southall
and
K. A.
Dill
, “
The mechanism of hydrophobic solvation depends on solute radius
,”
J. Phys. Chem. B
104
,
1326
1331
(
2000
).
13.
F. H.
Stillinger
, “
Structure in aqueous solutions of nonpolar solutes from the standpoint of scaled-particle theory
,”
J. Solution Chem.
2
,
141
158
(
1973
).
14.
D. M.
Huang
and
D.
Chandler
, “
Cavity formation and the drying transition in the Lennard-Jones fluid
,”
Phys. Rev. E
61
,
1501
1506
(
2000
).
15.
H.
Acharya
,
S.
Vembanur
,
S. N.
Jamadagni
, and
S.
Garde
, “
Mapping hydrophobicity at the nanoscale: Applications to heterogeneous surface and proteins
,”
Faraday Discuss.
146
,
353
365
(
2010
).
16.
S. I.
Mamatkulov
,
P. K.
Khabibullaev
, and
R. R.
Netz
, “
Water at hydrophobic substrates: Curvature, pressure and temperature effects
,”
Langmuir
20
,
4756
4763
(
2004
).
17.
A. J.
Patel
,
P.
Varilly
, and
D.
Chandler
, “
Fluctuations of water near extended hydrophobic and hydrophilic surfaces
,”
J. Phys. Chem. B
114
,
1632
1637
(
2010
).
18.
M.
Mezger
,
H.
Reichert
,
S.
Schöder
,
J.
Okasinski
,
H.
Schröder
,
H.
Dosch
,
D.
Palms
,
J.
Ralston
, and
V.
Honkimäki
, “
High-resolution in situ x-ray study of the hydrophobic gap at the water–octadecyl-trichlorosilane interface
,”
Proc. Natl. Acad. Sci. U. S. A.
103
,
18401
18404
(
2006
).
19.
M.
Mezger
,
F.
Sedlmeier
,
D.
Horinek
,
H.
Reichert
,
D.
Pontoni
, and
H.
Dosch
, “
On the origin of the hydrophobic water gap: An X-ray reflectivity and md simulation study
,”
J. Am. Chem. Soc.
132
,
6735
6741
(
2010
).
20.
B. M.
Ocko
,
A.
Dhinojwala
, and
J.
Daillant
, “
Comment on “How water meets a hydrophobic surface”
,”
Phys. Rev. Lett.
101
,
039601
(
2008
).
21.
S.
Chattopadhyay
,
A.
Uysal
,
B.
Stripe
,
Y.-G.
Ha
,
T. J.
Marks
,
E. A.
Karapetrova
, and
P.
Dutta
, “
How water meets a very hydrophobic surface
,”
Phys. Rev. Lett.
105
,
037803
(
2010
).
22.
A. P.
Willard
and
D.
Chandler
, “
The molecular structure of the interface between water and a hydrophobic substrate is liquid-vapour like
,”
J. Chem. Phys.
141
,
18C519
(
2014
).
23.
R.
Evans
and
N. B.
Wilding
, “
Quantifying density fluctuations in water at a hydrophobic surface: Evidence for critical drying
,”
Phys. Rev. Lett.
115
,
016103
(
2015
).
24.
R.
Evans
,
M. C.
Stewart
, and
N. B.
Wilding
, “
Critical drying of liquids
,”
Phys. Rev. Lett.
117
,
176102
(
2016
).
25.
R.
Evans
,
M. C.
Stewart
, and
N. B.
Wilding
, “
Drying and wetting transitions of a Lennard-Jones fluid: Simulations and density functional theory
,”
J. Chem. Phys.
147
,
044701
(
2017
).
26.
R.
Evans
,
M. C.
Stewart
, and
N. B.
Wilding
, “
From hydrophilic to superhydrophobic surfaces: A unified picture of the wetting and drying of liquids
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
23901
23908
(
2019
).
27.
R.
Evans
and
M. C.
Stewart
, “
The local compressibility of liquids near non-adsorbing substrates: A useful measure of solvophobicity and hydrophobicity?
,”
J. Phys.: Condens. Matter
27
,
194111
(
2015
).
28.
T.
Eckert
,
N. C. X.
Stuhlmüller
,
F.
Sammüller
, and
M.
Schmidt
, “
Fluctuation profiles in inhomogeneous fluids
,”
Phys. Rev. Lett.
125
,
268004
(
2020
).
29.
M. K.
Coe
,
R.
Evans
, and
N. B.
Wilding
, “
Measures of fluctuations for a liquid near critical drying
,”
Phys. Rev. E
105
,
044801
(
2022
).
30.
C. A.
Cerdeiriña
,
P. G.
Debenedetti
,
P. J.
Rossky
, and
N.
Giovambattista
, “
Evaporation length scales of confined water and some common organic liquids
,”
J. Phys. Chem. Lett.
2
,
1000
1003
(
2011
).
31.
C.
Ebner
,
W. F.
Saam
, and
A. K.
Sen
, “
Critical and multicritical wetting phenomena in systems with long-range forces
,”
Phys. Rev. B
31
,
6134
6136
(
1985
).
32.
C.
Ebner
and
W. F.
Saam
, “
Effect of long-range forces on wetting near bulk critical temperatures: An Ising-model study
,”
Phys. Rev. B
35
,
1822
1834
(
1987
).
33.
M. P.
Nightingale
,
W. F.
Saam
, and
M.
Schick
, “
Wetting and growth behaviors in adsorbed systems with long-range forces
,”
Phys. Rev. B
30
,
3830
3840
(
1984
).
34.
M. K.
Coe
,
R.
Evans
, and
N. B.
Wilding
, “
Density depletion and enhanced fluctuations in water near hydrophobic solutes: Identifying the underlying physics
,”
Phys. Rev. Lett.
128
,
045501
(
2022
).
35.
M. K.
Coe
, “
Hydrophobicity across length scales: The role of surface criticality
,” Ph.D. thesis,
University of Bristol
,
2021
.
36.
R.
Evans
, “
Density functionals in the theory of nonuniform fluids
,” in
Fundamentals of Inhomogeneous Fluids
, edited by
D.
Henderson
(
Marcel Dekker
,
1992
), pp.
85
175
.
37.
M. C.
Stewart
, “
Effect of substrate geometry on interfacial phase transitions
,” Ph.D. thesis,
University of Bristol
,
2006
.
38.
R.
Roth
, “
Fundamental measure theory for hard-sphere mixtures: A review
,”
J. Phys.: Condens. Matter
22
,
063102
(
2010
).
39.
V.
Molinero
and
E. B.
Moore
, “
Water modelled as an intermediate element between carbon and silicon
,”
J. Phys. Chem. B
113
,
4008
4016
(
2009
).
40.
M. C.
Stewart
and
R.
Evans
, “
Critical drying at a spherical substrate
,”
J. Phys.: Condens. Matter
17
,
S3499
S3505
(
2005
).
41.
M. C.
Stewart
and
R.
Evans
, “
Wetting and drying at a curved substrate: Long-ranged forces
,”
Phys. Rev. E
71
,
011602
(
2005
).
42.
R.
Evans
,
J. R.
Henderson
, and
R.
Roth
, “
Nonanalytic curvature contributions to solvation free energies: Influence of drying
,”
J. Chem. Phys.
121
,
12074
12084
(
2004
).
43.
S.
Dietrich
, “
Wetting phenomena
,” in
Phase Transitions and Critical Phenomena
, edited by
C.
Domb
and
J. L.
Lebowitz
(
Academic Press
,
1988
), Vol. 12.
44.
M.
Schick
, “
Introduction to wetting phenomena
,” in
Les Houches 1988 Session XLVIII Liquids at Interfaces
, edited by
J.
Charvolin
,
J. F.
Joanny
, and
J.
Zinn-Justin
(
North-Holland
,
1990
).
45.
T.
Bieker
and
S.
Dietrich
, “
Wetting of curved surfaces
,”
Physica A
252
,
85
137
(
1998
).
46.
S.
Dietrich
and
M.
Napiórkowski
, “
Analytic results for wetting transitions in the presence of van der Waals tails
,”
Phys. Rev. A
43
,
1861
1885
(
1991
).
47.
M. K.
Coe
, “
cDFT Package
,” https://github.com/marykcoe/cDFT_Package,
2021
(Online).
48.
R.
Evans
, “
The nature of the liquid-vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids
,”
Adv. Phys.
28
,
143
200
(
1979
).
49.
D. E.
Sullivan
and
M. M. T.
da Gama
, “
Wetting transitions and multilayer adsorption at fluid interfaces
,” in
Fluid Interfacial Phenomena
, edited by
C. A.
Croxton
(
John Wiley & Sons
,
1986
).
50.
M. K.
Coe
,
R.
Evans
, and
N. B.
Wilding
, “
The coexistence curve and surface tension of a monatomic water model
,”
J. Chem. Phys.
156
,
154505
(
2022
).
51.
M. K.
Coe
, “
GCMC Package for mw water
,” https://github.com/marykcoe/Mont_Carlo,
2021
(Online).
52.
E. W.
Lemmon
,
M. O.
McLinden
, and
D. G.
Friend
, “
Thermophysical properties of fluid systems
,” in
NIST Chemistry WebBook, NIST Standard Reference Database Number 69
, edited by
P. J.
Linstrom
and
W. G.
Mallard
(
National Institute of Standards and Technology
Gaithersburg MD
,
2021
); accessed 29 June 2021.
53.
S. N.
Jamadagni
,
R.
Godawat
, and
S.
Garde
, “
Hydrophobicity of proteins and interfaces: Insights from density fluctuations
,”
Annu. Rev. Chem. Biomol. Eng.
2
,
147
171
(
2011
).
54.
D. M.
Huang
,
P. L.
Geissler
, and
D.
Chandler
, “
Scaling of hydrophobic solvation free energies
,”
J. Phys. Chem. B
105
,
6704
6709
(
2001
).
55.
U.
Schnupf
and
J. W.
Brady
, “
Water structuring above solutes with planar hydrophobic surfaces
,”
Phys. Chem. Chem. Phys.
19
,
11851
11863
(
2017
).
56.
J. W. G.
Tyrrell
and
P.
Attard
, “
Images of nanobubbles on hydrophobic surfaces and their interactions
,”
Phys. Rev. Lett.
87
,
176104
(
2001
).
57.
R.
Steitz
,
T.
Gutberlet
,
T.
Hauss
,
B.
Klösgen
,
R.
Krastev
,
S.
Schemmel
,
A. C.
Simonsen
, and
G. H.
Findenegg
, “
Nanobubbles and their precursor layer at the interface of water against a hydrophobic substrate
,”
Langmuir
19
,
2409
2418
(
2003
).
58.

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.

Published open access through an agreement with University of Bristol School of Physics