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.
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 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, , 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 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 from above.25, ξ‖ is in turn linked to an important spatial measure of density fluctuations, the local compressibility profile,27 defined as . Specifically, one can show25 that for distances z in the interfacial region , where ρ′ is the gradient of the density profile. This implies that the maximum in χ(z) also diverges as . 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 . 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 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 . 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 , 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 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 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
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
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.
B. Monatomic water model
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 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 . 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.
As emphasized in our short report,34 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.
A. SR ff, LR sf
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 + ℓ)).
B. LR ff, LR sf
IV. RESULTS FROM CLASSICAL DFT CALCULATIONS FOR A LENNARD-JONES SOLVENT
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 , 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.
A. Profiles of density and local compressibility
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.
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 Rs ≫ Rc, 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 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.
B. Testing the scaling predictions
We employ DFT measurements of ℓeq and χ(ℓeq; Rs) 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.
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 , 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; Rs) 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 χ(ℓeq; Rs) 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 . 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.
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 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 χ(ℓeq; Rs) 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 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 χ(ℓeq; Rs) 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 χ(ℓeq; Rs). However, when Rs is large, the contours are almost vertical, indicating that Rs has little influence on the behavior of ℓeq and χ(ℓeq; Rs) 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; Rs), 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.
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 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.
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.
B. Profiles of local density and compressibility
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 σmw ≤ Rs ≤ 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.
The hydrophobic system considered here has interactions that correspond to SR ff, LR sf, see Sec. II B, and hence in the planar limit , 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.
C. Testing the binding potential predictions for ℓeq and χ(ℓeq)
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 term entering the effective pressure , so we predict that ℓeq ∼ ln(Rs) and χ(ℓeq; Rs) ∼ 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 χ(ℓeq; Rs) 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 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.
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.
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 , 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; Rs) 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 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
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.