# Understanding the physics of hydrophobic solvation

Mary K. Coe  Robert Evans  Nigel B. Wilding  H.H. Wills Physics Laboratory,
Royal Fort, University of Bristol, Bristol BS8 1TL, U.K.
<EMAIL_ADDRESS>
###### Abstract

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 (LJ) solvent and Grand Canonical Monte Carlo (GCMC) simulations
of a monatomic water (mw) model. We argue that the density depletion and
enhanced fluctuations are near-critical phenomena. Specifically, we show that
they can be viewed as remnants of the critical drying surface phase transition
that occurs at bulk liquid-vapor coexistence in the macroscopic planar limit,
i.e. as the solute radius $R_{s}\to\infty$. Focusing on the radial density
profile $\rho(r)$ and a sensitive spatial measure of fluctuations, the local
compressibility profile $\chi(r)$, our binding potential analysis provides
explicit predictions for the manner in which the key features of $\rho(r)$ and
$\chi(r)$ scale with $R_{s}$, the strength of solute-water attraction
$\varepsilon_{sf}$, and the deviation from liquid-vapor coexistence of the
chemical potential, $\delta\mu$. 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 Huang and
Chandler (2002); Oleinikova and Brovchenko (2012); Mittal and Hummer (2008);
Sarurpria and Garde (2009); Patel  _et al._ (2012); Vaikuntanathan  _et al._
(2016); Bischofberger  _et al._ (2014); Rego and Patel (2022). Typically
hydrophobicity is manifested by non-polar solutes, but large solutes can
exhibit both polar and non-polar regions leading to amphiphilic behaviour. For
example, it is argued hydrophobicity drives amphiphilic molecules to self
assemble into micelles Chandler (2005a) and proteins to bind with ligands
Qvist  _et al._ (2008).

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 behaviour will generally differ from that near a larger
solute Lum, Chandler, and Weeks (1999); Southall and Dill (2000); Stillinger
(1973). 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
$\sigma_{s}/\sigma_{w}\gtrsim 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 Lum, Chandler, and Weeks
(1999); Huang and Chandler (2000a); Sarupria and Garde (2009); Acharya  _et
al._ (2010); Mamatkulov and Khabibullaev (2004); Patel, Varilly, and Chandler
(2010); Mittal and Hummer (2008); Oleinikova and Brovchenko (2012);
Vaikuntanathan  _et al._ (2016); Rego and Patel (2022). 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 behaviour.

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 $\theta$ 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
$\theta>90^{\circ}$ indicates hydrophobic behaviour, with the extreme
hydrophobic limit corresponding to $\theta\to 180^{\circ}$. Planar substrates
for which $\theta>140^{\circ}$ 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 surfaceMezger  _et al._ (2006, 2010); Ocko,
Dhinojwala, and Daillant (2008), although its extent has been a matter of
debateChattopadhyay  _et al._ (2010). 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 Willard and
Chandler (2014) 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 rationalises the observed phenomenology of water at
superhydrophobic planar substrates in terms of the physics of surface phase
transitions, specifically the phenomenon of critical drying Evans and Wilding
(2015); Evans, Stewart, and Wilding (2016, 2017, 2019). 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 $\epsilon_{sf}$ is reduced. The
drying point, $\epsilon_{sf}^{d}$ marks the extreme hydrophobic limit for
which the contact angle attains $\theta=180^{\circ}$.

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 $\mu$ is prescribed. Here one considers the adsorption
of solvent/fluid at the planar substrate characterised by a density profile
$\rho(z)$ that measures the average local number density of solvent molecules
as a function of the perpendicular distance $z$ from the substrate. GCMC
simulations of such a system find that as the drying point is approached on
lowering $\epsilon_{sf}$ towards $\epsilon_{sf}^{d}$ from above, $\rho(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 lengthscales Evans, Stewart, and Wilding (2017). These bubbles
are a manifestation of near-critical behaviour; they are associated with the
observed enhancement of density fluctuations Evans and Wilding (2015); Evans,
Stewart, and Wilding (2017) in the vicinity of the substrate. The typical
lateral extent of bubbles is measured by a correlation length
$\xi_{\parallel}$ that diverges in a power law fashion as
$\epsilon_{sf}\to\epsilon_{sf}^{d}$ from above Evans, Stewart, and Wilding
(2017). $\xi_{\parallel}$ is in turn linked to an important spatial measure of
density fluctuations, the local compressibility profile Evans and Stewart
(2015), defined as $\chi(z)=\left(\partial\rho(z)/\partial\mu\right)_{T}$.
Specifically, one can show Evans, Stewart, and Wilding (2017) that for
distances $z$ in the interfacial region
$\chi(z)\sim\rho^{\prime}(z)\xi_{\parallel}^{2}$, where $\rho^{\prime}$ is the
gradient of the density profile. This implies that the maximum in $\chi(z)$
also diverges as $\epsilon_{sf}\to\epsilon_{sf}^{d}$. For very large but
finite values of $\xi_{\parallel}$, the system can be regarded as comprising a
liquid that is separated from the substrate by a vapor film of equilibrium
thickness $\ell_{eq}$, the magnitude of which also diverges as
$\epsilon_{sf}\to\epsilon_{sf}^{d}$. Extensive details concerning critical
drying behaviour, including the definitions and values of various critical
exponents, are set out by Evans et al. Evans, Stewart, and Wilding (2017).
That work used DFT and GCMC simulation to show that in the vicinity of a
critical drying point, the form of $\chi(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 Evans and Stewart (2015); Evans and Wilding (2015); Evans,
Stewart, and Wilding (2016, 2017). The local compressibility
$\chi(\textbf{r})=\left(\partial\rho(\textbf{r})/\partial\mu\right)_{T}$ is
generally the correlator (covariance) Evans and Wilding (2015) 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 Eckert  _et al._ (2020). InvestigationsCoe, Evans, and
Wilding (2022a) for a model liquid at a planar substrate near critical drying
show that the local compressibility drives the form of the 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 Evans and Wilding (2015), are expected to
exhibit critical drying when the solvent is (a) at liquid-vapor coexistence
and (b) in contact with a single infinite hydrophobic planar substrate having
$\epsilon_{sf}=\epsilon_{sf}^{d}$. However, water at standard temperature and
pressure (STP) i.e. ambient conditions, is not quite at coexistence Cerdeiriña
_et al._ (2011). Specifically ambient water has a very small but non-zero
supersaturation: the chemical potential deviation $\delta\mu$ from coexistence
$\mu_{co}$ is given by $\beta\delta\mu=\beta(\mu-\mu_{co})\approx
10^{-3},\;\beta=1/k_{B}T$. Moreover, in experiments it is not currently
possible to realize substrates that correspond to the extreme hydrophobic
limit, i.e. which have
$\delta\epsilon_{sf}\equiv\epsilon_{sf}-\epsilon_{sf}^{d}=0$, although for
materials that are the most superhydrophobic $\delta\epsilon_{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
$\theta<180^{\circ}$, and in the grand canonical setting this implies a vapor
layer of finite thickness $\ell_{eq}$ forms and the maximum of $\chi(z)$
remains finite. However, the small values of both $\beta\delta\mu$ and
$\delta\epsilon_{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 rationalised in terms of critical point scaling concepts Evans and Wilding
(2015).

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 Ebner, Saam, and Sen
(1985); Ebner and Saam (1987); Nightingale, Saam, and Schick (1984); Evans,
Stewart, and Wilding (2019). The range of both fluid-fluid (ff) and substrate-
fluid (sf) interactions are 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. Evans, Stewart, and Wilding (2019)
which reveal dramatically the effects that the interaction range can have on
the qualitative nature of the phase behaviour for both drying and wetting
transitions.

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

In the present work we attempt to develop such a theory by considering the
properties of solvents in the vicinity of a single spherical solvophobic
solute of radius $R_{s}$. By treating the solute curvature $R_{s}^{-1}$ as a
scaling field that measures deviations from surface criticality we construct a
mesoscopic theory of critical drying which generalizes earlier work on drying
transitions at planar substrates by incorporating an incipient drying (vapor)
film around the model solute. The theory assumes that for sufficiently small
curvature (large $R_{s}$) the system is near-critical and provides specific
predictions for the scaling properties of $\ell_{eq}$ and the maximum of
$\chi(r)=\left(\partial\rho(r)/\partial\mu\right)_{T}$ in terms of
$R_{s},\beta\delta\mu,\epsilon_{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 namely 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 solvent, demonstrating that for extended
solutes water behaves like any other solvent and hydrophobicity is merely a
particular case of universal solvophobic behaviour. 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 rationalising
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
vapour film $\ell_{eq}$ , that 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. Sec. 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 Coe, Evans, and Wilding (2022b).

##  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 which is investigated via DFT. The
second utilizes a popular monatomic water (mw) model which we investigate via
GCMC.

###  II.1 Lennard-Jones solvent

Figure 1: Sketch of the model system used in DFT and simulation. A smooth and
impenetrable solute of radius $R_{s}$ is centred on the origin and composed of
smaller particles distributed homogeneously with density $\rho_{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 $\varepsilon_{s}$ and diameter $\sigma_{s}$. The
net external (radial) potential experienced by solvent particles outside the
impenetrable zone, $V_{ext}(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 non-zero density or contact density, $\rho(R_{s}^{+})$, occurs
at a distance $R_{s}$ as indicated.

This system is a direct extension of that used in a previous study of
solvophobic planar substrates Evans, Stewart, and Wilding (2017). The setup is
outlined below and a detailed description can be found in the thesis of Coe
Coe (2021a). As shown in Fig. 1, a single spherical solute particle of radius
$R_{s}$ resides at the origin and is imagined to be composed of smaller
’virtual’ particles distributed homogeneously with density $\rho_{s}$. Within
the mean-field DFT treatment of attraction solvent particles interact with
each other via a pairwise potential of the general form Evans (1992); Stewart
(2006); Evans, Stewart, and Wilding (2017)

$\phi_{att}(r)=\begin{cases}-\varepsilon&r<r_{\rm min}\\\
4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]&r_{\rm
min}<r<r_{c}\\\ 0&r>r_{c}\end{cases}$ (1)

where $r=|\mathbf{r}-\mathbf{r}^{\prime}|$ and where $\mathbf{r}$ and
$\mathbf{r}^{\prime}$ denote the position vectors of a pair of particles.
$r_{\rm min}=2^{1/6}\sigma$ is the distance corresponding to the minimum,
$r_{c}$ is the cut-off radius, $\sigma$ is the LJ diameter, and $\varepsilon$
is the well-depth of the LJ interaction.

These pairwise potentials give rise to an attractive potential at each point
in the solvent whose form is derived Coe (2021a) by integrating over the
angular degrees of freedom in which the fluid has homogeneous density to
yield:

$\overline{\phi}_{att}(|r-r^{\prime}|)=\begin{cases}\frac{\pi\varepsilon}{rr^{\prime}}\left[(r-r^{\prime})^{2}-r_{\rm
min}^{2}+\frac{4}{5}\sigma^{12}\left(\frac{1}{r_{\rm
min}^{10}}-\frac{1}{r_{c}^{10}}\right)-2\sigma^{6}\left(\frac{1}{r_{\rm
min}^{4}}-\frac{1}{r_{c}^{4}}\right)\right]&|r-r^{\prime}|<r_{\rm min}\\\
\frac{\pi\varepsilon}{rr^{\prime}}\left[\frac{4}{5}\sigma^{12}\left(\frac{1}{(r-r^{\prime})^{10}}-\frac{1}{r_{c}^{10}}\right)-2\sigma^{6}\left(\frac{1}{(r-r^{\prime})^{4}}-\frac{1}{r_{c}^{4}}\right)\right]&r_{\rm
min}<|r-r^{\prime}|<r_{c}\\\ 0&|r-r^{\prime}|>r_{c}\end{cases}$ (2)

This effective attractive potential is then incorporated in the DFT
calculations (see Sec. IV) in the standard mean field fashion Roth (2010); Coe
(2021a).

Integrating the (virtual) solute particle -fluid particle LJ pair potential ,
with diameter $\sigma_{s}$ and well-depth $\varepsilon_{s}$, over the volume
of the solute Coe (2021a) gives rise to the net solute-solvent external
potential

$V_{ext}(r)=\begin{cases}\infty&r<R_{s}\\\
\varepsilon_{sf}\left[\frac{2\sigma_{s}^{9}}{15}\left(\frac{1}{(r_{+}-R_{s})^{9}}-\frac{1}{(r_{+}+R_{s})^{9}}\right)\right.\\\
+\frac{3\sigma_{s}^{9}}{20r_{+}}\left(\frac{1}{(r_{+}+R_{s})^{8}}-\frac{1}{(r_{+}-R_{s})^{8}}\right)\\\
+\sigma_{s}^{3}\left(\frac{1}{(r_{+}+R_{s})^{3}}-\frac{1}{(r_{+}-R_{s})^{3}}\right)\\\
+\left.\frac{3\sigma_{s}^{3}}{2r_{+}}\left(\frac{1}{(r_{+}-R_{s})^{2}}-\frac{1}{(r_{+}+R_{s})^{2}}\right)\right]&r>R_{s}\end{cases}$
(3)

where $\varepsilon_{sf}=2\pi\rho_{s}\varepsilon_{s}\sigma_{s}^{3}/3$ is the
effective solute-fluid attraction strength, $R_{s}$, $r$ are as shown in Fig
1, and $r_{+}=r+r_{\rm min,sf}$, where $r_{\rm min,sf}$ is the minimum of the
external potential. The external potential is shifted for numerical reasons
such that the minimum occurs at the impenetrable surface of the solute. It is
incorporated in the DFT calculations in the standard mannerRoth (2010); Coe
(2021a).

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

###  II.2 Monatomic water model

Our GCMC simulations utilise a popular monatomic water (mw) model proposed by
Molinero and Moore Molinero and Moore (2009). This coarse-grained model
represents a water molecule as a single particle and reproduces the
tetrahedral network structure of liquid water using a parameterization of the
Stillinger-Weber potential. Within the mw model, particles interact via the
potential Molinero and Moore (2009)

$\phi_{mw}(\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{r}_{k},\theta_{ijk})=\sum_{i}\sum_{j>i}\phi_{mw,2}(\mathbf{r}_{i},\mathbf{r}_{j})\\\
+\sum_{i}\sum_{j\neq
i}\sum_{k>j}\phi_{mw,3}(\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{r}_{k},\theta_{ijk})$
(4)

where the two-body, $\phi_{mw,2}$ and three-body, $\phi_{mw,3}$, potentials
are

$\phi_{mw,2}(\mathbf{r}_{i},\mathbf{r}_{j})=\\\
A\varepsilon_{mw}\left[B\left(\frac{\sigma_{mw}}{r}\right)^{4}-1\right]\exp\left(\frac{\sigma_{mw}}{|\mathbf{r}_{i}-\mathbf{r}_{j}|-a\sigma_{mw}}\right)\;,$
(5)
$\phi_{mw,3}(\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{r}_{k},\theta_{ijk})=\\\
\lambda\varepsilon_{mw}\left[\cos\theta_{ijk}-\cos\theta_{0}\right]^{2}\exp\left(\frac{\gamma\sigma_{mw}}{|\mathbf{r}_{i}-\mathbf{r}_{j}|-a\sigma_{mw}}\right)\\\
\times\exp\left(\frac{\gamma\sigma_{mw}}{|\mathbf{r}_{i}-\mathbf{r}_{k}|-a\sigma_{mw}}\right)\;,$
(6)

and $A=7.049556277$, $B=0.6022245584$, $\gamma=1.2$ are constants which
determine the form and scale of the potential, $\lambda=23.15$ is the
tetrahedrality parameter, $\theta_{0}=109.47^{\circ}$ is the angle favoured
between waters, $a=1.8$ sets the cut-off radius, $\sigma_{mw}=2.3925$Å is the
diameter of a mw particle, and
$\varepsilon_{mw}=6.189\;\mathrm{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 $R_{s}$ which is centered on the
origin as in Fig. 1. The mw particles interact with the solute via a potential
having the same form Eq. (3) used in our DFT calculations.

##  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 $\epsilon_{sf}$. The approach has a long history in
the physics of interfacial phenomena. Most pertinent to the present
investigation is its deployment in several recent studies of drying for fluids
at planar surfaces Evans, Stewart, and Wilding (2016, 2017) and its use in
ascertaining the surface phase diagrams for wetting and drying for different
combinations of substrate-fluid and fluid-fluid interaction range Evans,
Stewart, and Wilding (2019). A binding potential analysis was also applied
earlier to the study of drying around (very) large spherical
particles/colloids Stewart and Evans (2005a, b); Evans, Henderson, and Roth
(2004). These studies focused on the singular behaviour of the free energy of
solvation and the adsorption in the limit $R_{s}^{-1}\rightarrow 0$. 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 $\ell_{eq}$, and the local compressibility profile $\chi(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 $R_{s}$ in contact with a liquid
which –like water– has a small supersaturation, i.e. the chemical potential
$\mu$ is slightly above the value for liquid-vapour coexistence $\mu_{co}$ for
some prescribed subcritical temperature $T<T_{c}$. Provided the solute-fluid
attractive strength is sufficiently weak, a vapor film of width $\ell$ will
form around the solute, intruding between the solute surface and the liquid.
Following standard treatments Dietrich (1988); Schick (1990), we employ a
sharp-kink approximation to describe the density profile of the fluid around
the solute as

$\rho(r)=\begin{cases}\rho_{v}&R_{s}<r<R_{s}+\ell\\\
\rho_{l}&R_{s}+\ell<r\end{cases}$ (7)

where $\rho_{v}$ and $\rho_{l}$ are the coexisting densities of the vapor and
liquid, respectively. Such a description assumes that the region of excluded
volume between the surface particles of the solute and the fluid particles
$dw$ is incorporated into the radius of the solute, such that $R_{s}$ is the
effective radius and $\rho(R_{s}^{+})$ is the first non-zero fluid density, as
shown in fig. 2(b); see also Fig.1. The equivalent planar surface system,
explored in previous work Evans, Stewart, and Wilding (2017), is given in fig.
2(a).

Figure 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
(grey), vapour (light blue) and liquid (darker blue). The densities of the
vapor and liquid are assumed to take their coexistence values $\rho_{v}$ and
$\rho_{l}$, respectively. A region of excluded volume of width $dw$ exists
between the solute and vapour arising from the finite size of the particles.
The density of the fluid is measured from the centre of the fluid particle and
the width of the vapour slab is $\ell$. (b) The system is assumed to consist
of a large solute (grey), a vapour shell (light blue) of width $\ell$, and
liquid (darker blue). Again a region of excluded volume (white) of width $dw$
exists between the solute and vapour leading to an effective solute radius
$R_{s}$.

The excess grand potential $\Omega_{ex}$ for such a system can be written as
Bieker and Dietrich (1998); Stewart and Evans (2005b)

$\Omega_{ex}=\gamma_{sv}(R_{s})A_{sv}+\gamma_{lv}(R_{s}+\ell)A_{lv}+\omega(\ell|R_{s})A_{sv}+\delta\mu\Delta\rho
V_{v}$ (8)

where $\gamma$ is the surface tension and $A$ the surface area of the solute-
vapour (subscript $sv$) and liquid-vapour (subscript $lv$) interfaces,
$\delta\mu=(\mu-\mu_{co})$, $\Delta\rho=(\rho_{l}-\rho_{v})$, $V_{v}$ is the
volume of the vapour and $\omega(\ell|R_{s})$ is the binding potential – the
contribution to the free energy that arises from the $\ell$-dependent
interaction of the incipient vapor-liquid interface and the solute. Following
earlier treatments Stewart and Evans (2005b); Bieker and Dietrich (1998), we
assume that the solute is sufficiently large that $\gamma_{sv}(R_{s})$ and
$\gamma_{lv}(R_{s}+\ell)$ can be approximated by their planar surface
equivalents, $\gamma_{sv}$ and $\gamma_{lv}$ respectively. Furthermore, we
note that it has been shown previously Stewart and Evans (2005b) that
$\omega(\ell|R_{s})=\omega(\ell)(1+\ell/R_{s}+\mathcal{O}(\ell^{2}ln(\ell/2R_{s})/R_{s}^{2}))$.
It follows that the binding potential for a curved surface/solute can be
approximated in terms of that for a planar surface $\omega(\ell)$ for large
solutes. Adopting these approximations and substituting $A_{sv}=4\pi
R_{s}^{2}$, $A_{lv}=4\pi(R_{s}+\ell)^{2}$ and $V_{v}=4\pi(R_{s}+\ell)^{3}/3$
gives

$\omega_{ex}(\ell|R_{s})\equiv\frac{\Omega_{ex}}{A_{sv}}\approx\gamma_{sv}+\gamma_{lv}+\omega(\ell)+\tilde{p}\ell$
(9)

where under the large solute approximation, terms of order
$\mathcal{O}(\ell^{2}/R_{s}^{2})$ and higher have been neglected. Here
$\tilde{p}$ is the effective pressure, defined as

$\tilde{p}=\frac{2\gamma_{lv}}{R_{s}}+\delta\mu\Delta\rho\>.$ (10)

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

###  III.1 SR ff, LR sf

For a system with SR ff, LR sf interactions the binding potential
$\omega(\ell)$ takes the form Evans, Stewart, and Wilding (2019, 2017, 2016)

$\omega_{SR}(\ell)=a(T)e^{-\ell/\xi_{b}}+\frac{b(T)}{\ell^{2}}+H.O.T$ (11)

where $a(T)$ is independent of $l$, depends on temperature and has dimensions
of energy per unit area. $\xi_{b}$ is the bulk correlation length of the
(intruding) vapor phase and $b(T)$ has dimensions of energy and is given by

$b(T)=-b_{o}\rho_{s}\varepsilon_{s}\sigma_{s}^{6}$ (12)

with $b_{o}=\pi\Delta\rho/3$,where, as previously, $\rho_{s}$ is the density
of (virtual) solute particles and $\varepsilon_{s}$ and $\sigma_{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 eqn. (9) and minimising with respect to $\ell$ yields an
expression for the equilibrium vapor layer thickness, $\ell_{eq}$:

$-\frac{\ell_{eq}}{\xi_{b}}=\ln\left(\frac{\xi_{b}}{a}\right)+\ln\left(\tilde{p}-\frac{2b}{\ell_{eq}^{3}}\right)\>.$
(13)

which reduces to the expression for a planar surfaceEvans, Stewart, and
Wilding (2019) in the limit $R_{s}^{-1}\rightarrow 0$.

Including the $R_{s}$ dependence in Eq. 9 forges a potential link between
hydro/solvophobicity phenomena on microscopic and macroscopic scales. In the
case of a hard solute, where $\varepsilon_{s}=0$, one finds that
$\ell_{eq}\sim-\ln\tilde{p}$. As shown by Evans et. al. Evans, Henderson, and
Roth (2004), this leads to the identification of two regimes of scaling,
separated by the length-scale of capillary evaporation
$R_{c}=2\gamma_{lv}/\delta\mu\Delta\rho$

$\ell_{eq}(\varepsilon_{s}=0|R_{s})\sim\begin{cases}\ln R_{s}&R_{s}\ll
R_{c}\\\ -\ln\delta\mu&R_{c}\ll R_{s}\end{cases}$ (14)

Similarly, at bulk coexistence $\delta\mu=0$, two regimes of scaling can be
identified :

$\ell_{eq}(\delta\mu=0|R_{s})\sim\begin{cases}\ln
R_{s}&R_{s}\ll\frac{\gamma_{lv}\ell_{eq}^{3}}{|b|}\\\
-\ln\varepsilon_{s}+3\ln\ell_{eq}&\frac{\gamma_{lv}\ell_{eq}^{3}}{|b|}\ll
R_{s}\end{cases}$ (15)

where the latter corresponds to the planar case Evans et. al. Evans, Stewart,
and Wilding (2017). These results, which pertain to different physical
limiting cases, highlight that there are regions of the parameter space
$(\delta\mu,\varepsilon_{s},R_{s})$ for which the behaviour of $\ell_{eq}$ is
dependent predominately on only one parameter. We note that as shown
previously Evans, Stewart, and Wilding (2017), for the planar substrate case
$R_{s}=\infty$, Eqs. 15 and 14 imply that critical drying for a system with SR
ff, LR sf interactions occurs for $\delta\mu=0,\varepsilon_{s}=0$.

The magnitude of the local compressibility at $\ell=\ell_{eq}$ provides a
useful measure of the scale of density fluctuations in the neighbourhood 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 $\rho(r)=S(r-(R_{s}+\ell))$.

The local compressibility at $\ell_{eq}$ can then be found using Evans,
Stewart, and Wilding (2017)

$\chi(\ell_{eq}|R_{s})=-\rho^{\prime}(R_{s}+\ell_{eq})\left.\frac{\partial\ell_{eq}}{\partial\mu}\right|_{T}$
(16)

where $\rho^{\prime}$ is the spatial derivative of the profile. Substituting
$\ell_{eq}$ given in eqn. (13) yields

$\chi(\ell_{eq}|R_{s})=\xi_{b}\Delta\rho\rho^{\prime}(R_{s}+\ell_{eq})\left(\tilde{p}-\frac{2b}{\ell_{eq}^{3}}\left(1-\frac{3\xi}{\ell_{eq}}\right)\right)^{-1}$
(17)

As expected, in the limit $R_{s}^{-1}\rightarrow 0$, this reduces to the
result for a planar substrate.

Similarly to $\ell_{eq}$, we can identify three regimes in which the behaviour
of $\chi(\ell_{eq}|R_{s})$ is controlled predominately by $R_{s}$, $\delta\mu$
or $\epsilon_{s}$. For a hard solute, where $\varepsilon_{s}=0$, we find

$\chi(\ell_{eq}(\varepsilon_{s}=0)|R_{s})\sim\begin{cases}R_{s}&R_{s}\ll
R_{c}\\\ \delta\mu^{-1}&R_{c}\ll R_{s}\end{cases}$ (18)

where the latter is the result for drying at a hard planar wall, e.g. Evans
and Stewart Evans and Stewart (2015). At bulk coexistence $\delta\mu=0$, we
find

$\chi(\ell_{eq}(\delta\mu=0)|R_{s})\sim\begin{cases}R_{s}&R_{s}\ll\frac{\gamma_{lv}\ell_{eq}^{3}}{|b|}\\\
\varepsilon_{s}^{-1}&\frac{\gamma_{lv}\ell_{eq}^{3}}{|b|}\ll R_{s}\end{cases}$
(19)

where the latter was identified, for a planar wall, by Evans et. al. Evans,
Stewart, and Wilding (2017).

###  III.2 LR ff, LR sf

For the case of LR ff LR sf interactions, the binding potential takes the form
Evans, Stewart, and Wilding (2019) and references therein

$\omega_{LR}(\ell)=\frac{b_{LR}(T)}{\ell^{2}}+\frac{c(T)}{\ell^{3}}+H.O.T$
(20)

where

$\displaystyle b_{LR}(T)$ $\displaystyle=$ $\displaystyle
b_{o}(\rho_{v}\varepsilon\sigma^{6}-\rho_{s}\varepsilon_{s}\sigma_{s}^{6}),$
$\displaystyle c(T)$ $\displaystyle=$ $\displaystyle
2(dw+z_{min})\rho_{s}\varepsilon_{s}\sigma_{s}^{6}b_{o}\>,$ (21)

where $b_{o}=\pi\Delta\rho/3$, as defined previously and
$z_{min}=(2/5)^{1/6}\sigma_{s}$ is the location of the minimum in the (planar)
substrate -fluid potential; see Ref. 39 for more details. Clearly $c(T)$ is
positive for all temperatures. Assuming that $\rho_{s}$ is constant, the sign
of $b_{LR}(T)$ can change from negative to positive upon increasing $T$ as the
vapour density $\rho_{v}$ increases or, indeed, as the attraction strength
$\epsilon_{s}$ decreases. The behaviour of $b_{LR}(T)$ determines the location
of the minimum in $\omega_{ex}(\ell)$ that corresponds to the equilibrium film
width $\ell_{eq}$. In particular, the drying temperature $T_{d}$, at which the
thickness of the vapour film $\ell_{eq}\to\infty$, is determined by the
condition $b_{LR}(T_{d})=0$. If one fixes the temperature, it follows that the
critical drying point occurs for substrate-fluid attraction strength

$\varepsilon_{sf}^{d}=\frac{2\pi\rho_{v}\sigma^{6}}{3\sigma_{s}^{3}}\varepsilon$
(22)

as derived previously Evans, Stewart, and Wilding (2019). Summarizing,
critical drying for a system with LR ff, LR sf interactions occurs at bulk
coexistence $\delta\mu=0$ in the planar limit $R_{s}=\infty$ when the
attraction strength is given by Eq.(22).

The variation of $b_{LR}(T)$ in the approach to the critical drying point from
below allows us to define a dimensionless measure of the deviation from this
point :

$t^{\prime}=(\rho_{v}\varepsilon\sigma^{6}-\rho_{s}\varepsilon_{s}\sigma_{s}^{6})\varepsilon^{-1}\sigma^{-3}$
(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
Stewart and Evans (2005a). The convenience stems from the fact that
$b_{LR}(T)\propto t^{\prime}$.

Turning to the case of a spherical solute, we substitute eqn. (20) into eqn.
(9) and then minimising w.r.t. $\ell$ gives

$\frac{2b}{\ell_{eq}^{3}}+\frac{3c}{\ell_{eq}^{4}}=\tilde{p}$ (24)

Once again, this equation reduces to that for a planar surface as
$R_{s}^{-1}\rightarrow 0$. As in the SR ff, LR sf case, it is possible to
identify three regimes in which $R_{s}$, $\delta\mu$ and $\epsilon_{s}$
individually dominate the behaviour of $\ell_{eq}$.

Differentiating eqn. (24) yields an expression for $\chi(\ell_{eq}|R_{s})$:

$\chi(\ell_{eq}|R_{s})=\Delta\rho\;\rho^{\prime}(R_{s}+\ell_{eq})\frac{\ell_{eq}}{2}\left(2\tilde{p}-\frac{b}{\ell_{eq}^{3}}\right)^{-1}$
(25)

For the case of LR ff, LR sf interactions, it has been shown previouslyStewart
and Evans (2005a) that $\ell_{eq}$ can be expressed as a scaling function. An
appropriate form is

$\ell_{eq}=\frac{\sigma}{|t^{\prime}|}\mathcal{L}\left(\frac{\tilde{p}}{\varepsilon\Delta\rho|t^{\prime}|^{4}}\right)$
(26)

where $\mathcal{L}$ obeys equation (24). Adopting this result,
$\chi(\ell_{eq}|R_{s})$ can be written as

$\chi(\ell_{eq}|R_{s})=-\rho^{\prime}(R_{s}+\ell_{eq})\frac{\sigma}{\varepsilon|t^{\prime}|^{5}}\mathcal{L}^{\prime}\left(\frac{\tilde{p}}{\varepsilon\Delta\rho|t^{\prime}|^{4}}\right)$
(27)

where $\mathcal{L}^{\prime}$ is the derivative of the scaling function. This
result obeys eqn. (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 whilst the result for $b_{LR}(T)$ remains valid beyond the
sharp kink approximation this is not the case for $c(T)$ Dietrich and
Napiórkowski (1991).

##  IV Results from Classical DFT calculations for a Lennard-Jones solvent

Our DFT calculations Coe (2021b) are based on the familiar Rosenfeld
functional for the HS functional and the standard mean-field treatment of
attraction Evans (1979); Roth (2010); Evans (1992), i.e. they implement the
same free energy functional as described in earlier papers for the case of a
planar substrateEvans and Stewart (2015); Evans, Stewart, and Wilding (2017).
Our calculations adopt a system of the form shown in fig. 1 with LR sf
interactions described by Eq. (3). As mentioned earlier, for the ff
interactions we consider (i) the SR ff case in which the LJ interparticle
potential is truncated at $r_{c}=2.5\sigma$ and left unshifted, and (ii) the
LR ff case in which the true long ranged potential is approximated by taking
$r_{c}=200\sigma$. The geometry required to treat a spherical solute gives
rise to specific weight functions as described by Roth Roth (2010).

We perform our DFT calculations at the subcritical temperature $T=0.775T_{c}$
in accord with previous work Evans, Stewart, and Wilding (2016, 2017). Further
details of the bulk and coexistence state points that we have studied at this
temperature are set out in the thesis of Coe Coe (2021a) which also provides
guides to numerical implementation of DFT for a spherical solute. Measuring
the density profile $\rho(r)$ and the local compressibility profile $\chi(r)$
for various values of
$\varepsilon_{sf}=2\pi\rho_{s}\varepsilon_{s}\sigma_{s}^{3}/3$ , the effective
solute-fluid attraction strength, $\delta\mu,T$ and $R_{s}$, provides insight
into how the solvophobic response of the solvent is influenced by the solute
properties and the proximity of the solvent to bulk coexistence.

###  IV.1 Profiles of density and local compressibility

Figures 3(a)-(f) demonstrate the effect on $\rho(r)$ and $\chi(r)$ of varying
$R_{s}$ (various colours) and $\beta\delta\mu$ (from left to right) for a
(hard) solute with $\varepsilon_{sf}=0$ which is the value for critical drying
in this case of LR sf, SR ff interactions. As $R_{s}\rightarrow\infty$,
$\rho(r)$ and $\chi(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 $\beta\delta\mu$ is lowered
towards zero and the deviation from critical drying is reduced.

Figure 3: (a)-(l) display DFT results for $\rho(r)$ and $\chi(r)$, normalized
to their bulk values, for the case of SR ff, LR sf; (m)-(r) display the
corresponding results for LR ff, LR sf. Top row: Varying $R_{s}$ and
$\beta\delta\mu$ at constant $\delta\varepsilon_{sf}=0.0$. Middle Row: Varying
$\delta\varepsilon_{sf}$ and $R_{s}$ at constant $\beta\delta\mu=10^{-3}$.
Bottom Row: Varying $\delta\varepsilon_{sf}$ and $R_{s}$ at constant
$\beta\delta\mu=10^{-3}$ for LR ff, LR sf.

These observations accord with the binding potential analysis of Sec. III.1
which predicts, for the hard solute, a length scale
$R_{c}=2\gamma_{lv}/\delta\mu\Delta\rho$ that separates behaviour dependent
and independent of the curvature of the solute: for solutes having $R_{s}\gg
R_{c}$, the density and local compressibility profiles should be close to
those of a planar substrate. Figs. 3(a)-(f) permit a test of this prediction.
For $\beta\delta\mu=10^{-3}$, $R_{c}\approx 913\sigma$ and hence we would
expect profiles for $R_{s}>1000\sigma$ to be indistinguishable from those of
the planar substrate, as is indeed confirmed by the DFT results. Corresponding
behaviour is seen in figures 3(c),(d) and 3(e)-(f) for the cases of
$\beta\delta\mu=10^{-4}$ and $\beta\delta\mu=10^{-5}$, for which
$R_{c}=9134\sigma$ and $R_{c}=91340\sigma$, respectively.

In figures 3(g)-(l), $\beta\delta\mu$ is held constant at $10^{-3}$ whilst
varying $~{}\delta\varepsilon_{sf}$ (colours) and $R_{s}$ (left to right). In
this case, the binding potential analysis of section III.1 again predicts two
regimes in which the behaviour 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 figures
3(g)-(l) by comparing the variation in $\chi(r)$ relating to
$\delta\varepsilon_{sf}/\varepsilon\geq 0.6$ to those of
$\delta\varepsilon_{sf}/\varepsilon<0.6$ when moving from left to right,
increasing $R_{s}$ . In the former regime, the scale and form of $\chi(r)$
varies little, whilst for the latter regime, the position and height of the
maximum of $\chi(r)$ increase substantially. Note that $\rho(r)$ exhibits a
weaker evolution, confirming that $\chi(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 $\delta\varepsilon_{sf}/\varepsilon=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 behaviour as indicative of a solvophilic solute. However, the
local compressibility profiles provide important new insight. Whilst these
profiles also exhibit oscillations, these are not centred on the bulk
compressibility, as they would if the substrate were truly solvophilicEvans
and Stewart (2015). We note that for a planar substrate and the fluid at bulk
coexistence, $\beta\delta\mu=0.0$ , $\delta\varepsilon_{sf}/\varepsilon=1.2$
corresponds to a contact angle of $\approx 107.4^{\circ}$ and such a substrate
would be designated solvophobic. It follows that local density fluctuations,
as measured by $\chi(r)$, appear to be a far more reliable indicator of the
degree of solvophobicity than the density profile alone.

Figures 3(m)-(r) for LR ff, LR sf interactions demonstrate similar features as
those in figures 3(a)-(l) for a system with SR ff, LR sf interactions. Figures
3(m) and (n) demonstrate the influence of varying $R_{s}$ for constant
$\beta\delta\mu=10^{-3}$ and $\delta\varepsilon_{sf}=0$. As in the SR ff, LR
sf case, curvature dependent and independent regimes can be separated by the
value of the parameter $R_{c}$, which for this case is $R_{c}\approx
1138\sigma$. Figures 3(o)-(r) compare the influence of varying
$\delta\varepsilon_{sf}/\varepsilon$ and $R_{s}$ and again the regime
dependent on curvature occurs for values of
$\delta\varepsilon_{sf}/\varepsilon<0.6$. Whilst 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 the density fluctuations and extent of the density depletion are larger, in
agreement with the predictions of section III.

###  IV.2 Testing the scaling predictions

From the density and local compressibility profiles obtained from DFT, it is
possible to extract $\ell_{eq}$ and thus $\chi(\ell_{eq};R_{s})$. To do so, we
define $\ell_{eq}$ in the standard way, e.g. Sullivan and da Gama (1986)

$\ell_{eq}\equiv-\frac{\Gamma}{A_{sv}\Delta\rho}$ (28)

where $\Gamma$, is the Gibbs excess adsorption, obtained from the calculated
density profiles using numerical integration:

$\frac{\Gamma}{A_{sv}}=\frac{1}{R_{s}^{2}}\int_{R_{s}}^{\infty}\mathrm{d}r\;r^{2}(\rho(r)-\rho_{b})$
(29)

in the case of a curved substrate/solute and

$\frac{\Gamma}{A_{sv}}=\int_{0}^{\infty}\mathrm{d}z\;(\rho(z)-\rho_{b})$ (30)

in the case of a planar substrate. In each case, $\rho_{b}$ is the density of
the bulk liquid. $\chi(\ell_{eq};R_{s})$ can then be found by performing the
derivative of the density profile w.r.t. $\mu$ at $r=\ell_{eq}$.

We employ DFT measurements of $\ell_{eq}$ and $\chi(\ell_{eq};R_{s})$ for a
large range of values of $\delta\mu,\epsilon_{sf},R_{s}$ in order to perform
detailed tests of the scaling predictions of the binding potential analysis.

####  IV.2.1 SR ff, LR sf interactions

Results for $\ell_{eq}$ with SR ff LR sf interactions are shown in figure 4,
for systems with parameters in the range $(10^{-6}\leq\beta\delta\mu\leq
10^{-3},0.0\leq\delta\varepsilon_{sf}/\varepsilon\leq 1.0,10\sigma\leq
R_{s}\leq 10^{8}\sigma)$. Results for planar substrates are also included.
Within the figure, colour is used to indicate the degree of solvophobicity by
associating each value of $\delta\varepsilon_{sf}/\varepsilon$ with the
corresponding value of the contact angle which would pertain for a planar
system at liquid-vapor coexistence.

Within Fig. 4, values of $\ell_{eq}$ obtained from the DFT results for the
case of SR ff LR sf interactions are compared to equation (13). Excellent
agreement with the predicted (linear) form is found for
$\delta\varepsilon_{sf}/\varepsilon<0.4$, for a wide range of parameters - any
deviation from the linear relationship is associated with solutes of radius
$R_{s}<20\sigma$. The latter indicates a limit in size of solute for which the
effects of the drying critical point can be felt, and one might consider
whether such change of behaviour could be related to the change in solvation
behaviour often predicted to occur for solutes of radius about 1$\mathrm{nm}$
dissolved in water. For $\delta\varepsilon_{sf}/\varepsilon>0.4$ the agreement
between the binding potential prediction and DFT results is not as good.
Considering the density profiles in figures 3(h),(j) and (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 $\delta\varepsilon_{sf}/\varepsilon$
correspond to contact angles of $<150^{\circ}$ indicating that the effects of
the drying critical point are most strongly felt for very weak $sf$ attraction
corresponding to the supersolvophobic regime.

Figure 4: Comparison of DFT results for $\ell_{eq}$ obtained using eq. (28)
and the (linear) scaling prediction of equation (13), for the case of SR ff,
LR sf interactions. The temperature is fixed at $T=0.775T_{c}$, whilst
$(10^{-6}\leq\beta\delta\mu\leq
10^{-3},0.0\leq\delta\varepsilon_{sf}/\varepsilon\leq 1.0,10\sigma\leq
R_{s}\leq 10^{8}\sigma)$. The value of $\delta\varepsilon_{sf}$ and the
associated contact angle $\theta$ for each result are indicated by the colour
bar.

Figure 5 compares the predictions of the binding potential analysis for
$\chi(\ell_{eq};R_{s})$ to DFT results for the same systems as in figure 4.
Again we see excellent agreement between (17) and the DFT results however the
linear behaviour is found over a far more limited range :
$\delta\varepsilon_{sf}/\varepsilon\leq 0.2$. Any deviation within this range
is for $R_{s}<20\sigma$, as in figure 4. For
$\delta\varepsilon_{sf}/\varepsilon=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 $\ell_{eq}$ for
which the clear deviation begins - from figure 4 we see these correspond to
$\ell_{eq}<2\sigma$. The bulk vapor correlation length at this temperature is
$\xi_{b}=0.51\sigma$, hence for the binding potential prediction to be
physical, $1-3\xi_{b}/\ell_{eq}>0$ and therefore $\ell_{eq}>1.53\sigma$. For
smaller values of $\ell_{eq}$, the binding potential prediction for
$\chi(\ell_{eq};R_{s})$ is no longer physical. One might attempt to include
the neglected higher order terms, however these depend on the shape of the
liquid-vapor interface and are difficult to calculate. The crucial point is
that we are pushing the binding potential treatment to its extremes.

Figure 5: Comparison of DFT results for $\chi(\ell_{eq};R_{s})$ and the
(linear) scaling prediction of Eq. (17), for the case of SR ff, LR sf
interactions. The temperature is fixed at $T=0.775T_{c}$, whilst
$(10^{-6}\leq\beta\delta\mu\leq
10^{-3},0.0\leq\delta\varepsilon_{sf}/\varepsilon\leq 0.3,10\sigma\leq
R_{s}\leq 10^{8}\sigma)$. The value of $\delta\varepsilon_{sf}$ and the
associated contact angle $\theta$ for each result are indicated by the colour
bar.

####  IV.2.2 LR ff, LR sf interactions

The binding potential prediction for $\ell_{eq}$ for a system with LR ff LR sf
interaction is given in equation (26) and involves the scaling function
$\mathcal{L}$. We plot our DFT results for systems at fixed
$\beta\delta\mu=10^{-3}$ ,with parameters in the range
$0.0\leq\delta\varepsilon_{sf}/\varepsilon\leq 1.0,10\sigma\leq R_{s}\leq
10^{8}\sigma$ and employing the arguments of this scaling form, in figure 6.
We choose to make such a comparison, i.e. with the scaling form of
$\ell_{eq}$, because of the difficulty in determining accurately $c(T)$. Three
temperatures, $T=0.7T_{c},0.775T_{c}$ and $0.85T_{c}$ are considered, with the
arrow indicating the direction of increasing temperature. Colour is used to
indicate $\delta\varepsilon_{sf}/\varepsilon$ 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. Whilst the formula
for $b(T)$ is easily calculated and this coefficient in the binding potential
expansion is expected to be valid beyond the sharp-kink approximation the
coefficient $c(T)$ is dependent on the form of the liquid-vapor interface
Dietrich and Napiórkowski (1991) and leads to some ambiguity in determining
the explicit scaling function. We do not attempt to ascertain $\mathcal{L}$
but note this was attempted for the special case of $\beta\delta\mu=0.0$ and
very large solutes by Stewart and Evans Stewart and Evans (2005a). Here we
focus simply on data collapse.

Overall, the data collapse predicted from the binding potential for
$\ell_{eq}$ is confirmed by the DFT results in the regime
$\delta\varepsilon_{sf}/\varepsilon<0.6$ for the three temperatures; any
deviations correspond typically to cases where $R_{s}<20\sigma$. For
$\delta\varepsilon_{sf}/\varepsilon>0.6$, there is clear deviation from the
predicted functional form which becomes more pronounced as
$\delta\varepsilon_{sf}/\varepsilon$ 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 $\delta\varepsilon_{sf}/\varepsilon$
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 interactions strengths for which the contact
angle is $\theta>150^{\circ}$.

Figure 6: DFT results for $\ell_{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.7T_{c},0.775T_{c}$ and $0.85T_{c}$, with the red
arrow indicating the direction of increasing temperature. In each case
$\beta\delta\mu=10^{-3}$,whilst $0.0\leq\delta\varepsilon_{sf}/\varepsilon\leq
1.0$,as indicated by the colour bar and $10\sigma\leq R_{s}\leq 10^{8}\sigma$.
The inset shows the contact angle that would pertain in the planar limit for
each $\delta\varepsilon_{sf}$ and $T$ considered.

Turning finally to our DFT results for $\chi(\ell_{eq};R_{s})$ for the LR ff,
LR sf case, we tested the predicted scaling of eq. (27) in the main part of
figure 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
$\delta\varepsilon_{sf}/\varepsilon<0.6$, similar to the case of SR ff, LR sf
interactions.

Figure 7: DFT results for $\chi(\ell_{eq};R_{s})$ 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 eq. (25). In each case
$\beta\delta\mu=10^{-3}$, whilst
$0.0\leq\delta\varepsilon_{sf}/\varepsilon\leq 1.0$ , as indicated by the
colour bar, and $10\sigma\leq R_{s}\leq 10^{8}\sigma$. Results are shown for
three temperatures, $T=0.7T_{c},0.775T_{c}$ and $0.85T_{c}$, 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.

###  IV.3 Contour Plots

As predicted in section III, and observed in the profiles considered in
section IV.1, there are regimes of parameter space for which the individual
parameters $R_{s}$, $\delta\mu$ and $\epsilon_{s}$ dominate the behaviour of
$\ell_{eq}$ and $\chi(\ell_{eq})$. Such behaviour can be visualised more
readily in the contour plots of figure 8. Figures 8(a) and (b) compare values
of $\ell_{eq}$ obtained from DFT calculations for systems with SR ff, LR sf
and with LR ff, LR sf interactions, respectively, for varying $R_{s}$ and
$\delta\varepsilon_{sf}$ at fixed $\beta\delta\mu=10^{-3}$. Figures 8(d) and
(e) show the corresponding values of $\chi(\ell_{eq};R_{s})$ normalized to the
bulk values. The regions of parameter space in which individual parameters
dominate is immediately apparent: when $R_{s}$ is small, the contours are
largely horizontal, indicating that changing
$\delta\varepsilon_{sf}/\varepsilon$ has little influence on the behaviour of
$\ell_{eq}$ and $\chi(\ell_{eq};R_{s})$. However, when $R_{s}$ is large, the
contours are almost vertical, indicating that $R_{s}$ has little influence on
the behaviour of $\ell_{eq}$ and $\chi(\ell_{eq};R_{s})$ in this region. From
such plots we can make numerical estimates of the crossover length scale for
the change in behaviour, say for a given choice of $\delta\varepsilon_{sf}$,
which was not possible from the binding potential analysis alone.

Figures 8(c) and (f) compare values of $\ell_{eq}$ and
$\chi(\ell_{eq};R_{s})$, respectively, for a system with SR ff, LR sf
interactions at constant $\delta\varepsilon_{sf}=0.0$ and various
$\beta\delta\mu$ and $R_{s}$. As was discussed in sections III.1 and IV.1, for
this case we expect the behaviour of both quantities to depend almost solely
on $R_{s}$ when $R_{s}<R_{c}$ and on $\delta\mu$ when $R_{s}>R_{c}$. The
contours in figures 8(c) and (f) suggest this to be the case. As an example,
we consider the case $\beta\delta\mu=10^{-4}$ for which $R_{c}\approx
9134\sigma$. Figures 8(c) and (f), confirm the crossover between different
regimes indeed occurs around this value of $R_{s}$.

Figure 8: Contour plots of DFT results showing the dependence of $\ell_{eq}$
and $\chi(\ell_{eq})$ on $R_{s}^{-1}$ and $\delta\varepsilon_{sf}$.The
temperature $T=0.775T_{c}$ is fixed. In plots (a,b,d,e))
$\beta\delta\mu=10^{-3}$ and include the contact angle that pertains in the
planar limit at bulk coexistence for a given $\delta\varepsilon_{sf}$. (a)
$\ell_{eq}$ for SR ff interactions, (b) $\ell_{eq}$ for LR ff, (d)
$\chi(\ell_{eq})$ for SR ff and (e) $\chi(\ell_{eq})$ for LR ff. (c)
$\ell_{eq}$ for SR ff with $\delta\varepsilon_{sf}=0.0$ and increasing
$\beta\delta\mu$. (f) $\chi(\ell_{eq})$ for SR ff with
$\delta\varepsilon_{sf}=0.0$ and increasing $\beta\delta\mu$.

##  V Results from GCMC simulations of the mw solvent

###  V.1 Coexistence properties and simulation state points

mw is a relatively recent water model which has been shown to reproduce
accurately many of the properties of water under ambient conditions whilst
being far faster to implement than more established water models such as SPC/E
Molinero and Moore (2009). Recently, we have presented the first highly
accurate liquid-vapor phase diagram for mw Coe, Evans, and Wilding (2022c)
which we have measured via GCMC simulations Coe (2021c). 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 workCoe, Evans, and
Wilding (2022c), 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 Coe, Evans, and Wilding (2022c). Specifically when the
temperature-density phase diagram (coexistence curve) is scaled by the
critical temperature and critical density, a data collapse is observed onto
the similarly scaled phase diagram for real water, as shown in the lower panel
of Fig. 9. This finding suggests that in seeking to study mw water under
conditions equivalent to those of ambient real water, it is reasonable to
employ the same scaled temperature as ambient water, i.e. $T/T_{c}=0.46$. For
mw, this corresponds to a simulation temperature of $426$K, which we adopt in
our simulation studies below. We note that GCMC simulations of mw are
substantially more computationally efficient at $426$K than at $300$K,
although we have also considered the latter case as we mention below.

As mentioned in the Introduction, water at ambient conditions exhibits a small
oversaturation Cerdeiriña  _et al._ (2011) of approximately
$\beta\delta\mu\approx 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 $\mu$-$T$ planeCoe, Evans, and Wilding
(2022c) allow us to control precisely the oversaturation for our model and
thus impose a value appropriate to ambient water.

Figure 9: Upper plot: 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, $\sigma_{mw}=2.3925$Å, whilst
$\varepsilon_{mw}=6.189\;\mathrm{kcal/mol}$. Lower plot: Comparison of the
liquid-vapor phase diagrams of water (pink squares) and mW (blue circles) with
temperature $T$ and density $\rho$ scaled by their bulk critical values. The
dotted grey line indicates the value of $T/T_{c}$ for water at ambient
conditions. Data for water are taken from ref. Lemmon, McLinden, and Friend .

###  V.2 Profiles of local density and compressibility

As for the solvophobic LJ case, understanding the influence of varying the
parameters $R_{s}$,$\delta\mu$ and $\delta\epsilon_{sf}$ on the hydrophobic
response of mw can be gained from measurements of $\rho(r)$ and $\chi(r)$.
Figures 10(a) and (b) show the variation with $R_{s}$ and $\epsilon_{sf}$ at
fixed $\beta\delta\mu=10^{-3}$ , the value of the oversaturation pertaining to
ambient water. The temperature is fixed at $T=426\mathrm{K}$ which reproduces
the value of $T/T_{c}$ for ambient water as discussed in Sec. V.1. Solutes of
size $\sigma_{mw}\leq R_{s}\leq 17\sigma_{mw}$, were studied corresponding in
physical water units to $0.23925\mathrm{nm}\leq R_{s}\leq 4.06725\mathrm{nm}$.
This range was chosen to span the widely reported qualitative change in
hydrophobic behaviour that supposedly occurs around $R_{s}\approx
1\mathrm{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 $R_{s}^{-1}\to 0$
critical drying is expected to occur when $\varepsilon_{sf}=0.0$. Profiles for
values of $\varepsilon_{sf}$ that correspond to realistic hydrophobic solutes
Jamadagni, Godawat, and Garde (2011) are shown in Fig. 10(a,b). In all cases,
for very small solutes, $R_{s}<1\mathrm{nm}$, we find no region of depleted
density and enhanced density fluctuations, in agreement with previous work
Huang and Chandler (2002). As $R_{s}$ is increased beyond $1\mathrm{nm}$, a
region of depleted density begins to form, and density fluctuations are
enhanced on a similar length scale. The extent of the former and magnitude of
the latter grow with $R_{s}$. Oscillations in the density profile, which are
indicative of liquid packing effects, are dampened as the solute radius
increases. Oscillations in the local compressibility profiles are no longer
centred on $\chi_{b}$ as $R_{s}$ increases. Note that similar behaviour was
observed in our DFT calculations for a solvophobic system, sec. IV.1. The
apparent change in hydrophobic behaviour occurring at around
$R_{s}=1\mathrm{nm}$ is consistent with that reported in previous simulation
studies of the density profiles of atomistic water models Huang, Geissler, and
Chandler (2001); Chandler (2005a); Huang and Chandler (2002); Schnupf and
Brady (2017), strengthening our confidence that mw is a suitable water model
for studies of hydrophobic solvation. Whilst the limit of a planar substrate
is not explored here, we note that a previous study of SPC/E water Evans and
Wilding (2015) found very similar density and local compressibility profiles
to those of figures 10(a),(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 same
fractional temperature $T/T_{c}=0.46$ as water (which for mw corresponds to
$T=426$ K - see Sec. V.1), we have also investigated other temperatures.
Figures 10(c) to (h) compare density and local compressibility profiles for
three temperatures : $300\mathrm{K}$, the ambient temperature for which mw was
parameterised; $360\mathrm{K}$, the temperature at which mw almost exactly
reproduces the liquid-vapor surface tension of water; and $426\mathrm{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 $\chi(r)$ increases. For each temperature, the
same general behaviour is observed and mirrors that seen for a general
solvophobic system as studied by DFT, section IV.1.

The binding potential analysis suggested three regimes within which individual
parameters dominate the expected scaling behaviour. The existence of these
regimes was confirmed in DFT for a general solvophobic system by examining the
density and local compressibility profiles calculated in section IV.1.
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
behaviour. Consider, for example, the local compressibility profiles $\chi(r)$
in figure 10(a). As $R_{s}$ is reduced, the variation of $\chi(r)$ with
$\varepsilon_{sf}$ slows which is similar to the angled horizontal contours of
figure 8(d). For the hard solute , $\varepsilon_{sf}=0.0$, we expect the
behaviour of the density and local compressibility to be almost solely
dependent on $R_{s}$ when $R_{s}<R_{c}$, as was shown to be the case for
solvophobic systems in figures 3(a)-(f). For mw at $T=426\mathrm{K}$,
$R_{c}\approx 0.583\mu m$ when $\beta\delta\mu=10^{-3}$ and $R_{c}\approx
5.83\mu m$ when $\beta\delta\mu=10^{-4}$. Thus, for all values of $R_{s}$
investigated in our simulations we expect the density and local
compressibility profiles to be almost independent of $\beta\delta\mu$. Figure
11 confirms this to be the case: the profiles are near indistinguishable for
$\beta\delta\mu=10^{-3}$ and $\beta\delta\mu=10^{-4}$.

Figure 10: GCMC results for mw. (a) local compressibility $\chi(r)$ normalised
by its bulk value for various values of solute radius $R_{s}$, as indicated on
the abcissa, and various values of solute-fluid attraction
$\epsilon_{sf}/\epsilon_{mw}$ as indicated in the key. In units of kcal/mol
these values are from lower to higher:
$\varepsilon_{sf}=0.0,0.03095,0.12378,0.24756$. The temperature $T=426$K and
supersaturation $\beta\delta\mu=10^{-3}$ are fixed. (b) The corresponding
local density profile $\rho(r)$ normalised by its bulk value. (c)-(h)
Variation of $\chi(r)$and $\rho(r)$ with temperature $T$ for
$\varepsilon_{sf}=0.0$ and various $R_{s}$ as indicated in the key. The
supersaturation is $\beta\delta\mu=10^{-3}$.  Figure 11: GCMC results for
density (lower) and local compressibility (upper) profiles, scaled by their
bulk values, for mW at $T=426\mathrm{K}$ and $\varepsilon_{sf}=0.0$, and for
various $R_{s}$ and two values of $\beta\delta\mu$. In the case of the density
profile, uncertainties do not exceed line width.

###  V.3 Testing the Binding Potential Predictions for $\ell_{eq}$ and
$\chi(\ell_{eq})$

In contrast to DFT, our mw simulations cannot access the very large values of
$R_{s}$ that are required to verify fully the relationships predicted by our
binding potential analysis: namely eqs. (13) and (17). However, as discussed
in Sec. III.1, for the case $\varepsilon_{sf}=0.0$, the behaviour at small
$R_{s}$ is expected to be dominated by the Laplace ($R_{s}^{-1}$) term
entering the effective pressure $\tilde{p}$ so we predict that
$\ell_{eq}\sim\ln(R_{s})$ and $\chi(\ell_{eq};R_{s})\sim R_{s}$ (see eqs. (14)
and (18)). The binding potential analysis predicts that this scaling should
also hold when $\varepsilon_{sf}$ is sufficiently small. These predictions are
tested and verified in figures 12 and 13 for several values of
$\varepsilon_{sf}$. The fact that the predicted scaling behaviour is observed
when $R_{s}>1\mathrm{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 $\ell_{eq}$ and $\chi(\ell_{eq};R_{s})$ appears to work
quite well down to small solute radii where $\ell_{eq}$ is calculated to be a
small fraction of the water diameter $\sigma_{mw}$. In this regime there is no
discernible vapor ’film’. Rather there are regions of density depletion, shown
in Figures 10,11, that extend across only short distances from the surface of
the solute.

Figure 12: Results for $\ell_{eq}$ obtained from GCMC simulations of mw for
fixed $T=426\mathrm{K}$ and $\beta\delta\mu=10^{-3}$; $\varepsilon_{sf}$ is
given by the symbol. Each dotted line is a linear fit to the scaling
prediction. Figure 13: Results for $\chi(\ell_{eq};R_{s})$ obtained from GCMC
simulations of mw for fixed $T=426\mathrm{K}$ and $\beta\delta\mu=10^{-3}$;
$\varepsilon_{sf}$ is given by the symbol. Each dotted line is a linear fit to
the scaling prediction.

###  V.4 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 Mezger  _et al._ (2006, 2010); Ocko, Dhinojwala, and Daillant (2008);
Chattopadhyay  _et al._ (2010). Atomic force microscopy studies Tyrrell and
Attard (2001); Steitz  _et al._ (2003) 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 interface
Evans, Stewart, and Wilding (2017) 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 centre
of the simulation box for six solutes of various sizes. In each case the
solute is ’hard’, i.e. $\varepsilon_{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 further away. In all cases, the size
ratio of the mw particles to the solute is to scale.

As $R_{s}$ increases, bubbles (low density regions) form at the surface of the
solute and hence our simulations are in line with the nanobubble picture of
hydrophobicity Tyrrell and Attard (2001); Steitz  _et al._ (2003). Bubbles
occur over a large extent of the solute surface, particularly when
$R_{s}>9\sigma_{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 the approach
to the drying critical point Evans, Stewart, and Wilding (2017). 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 $\varepsilon_{sf}$ has the effect of
suppressing the bubbles and reducing their size.

Figure 14: Snapshots of mw at hard solutes, $\varepsilon_{sf}=0.0$ ; $T=426$K,
$\beta\delta\mu=10^{-3}$. (a) $R_{s}=\sigma_{mw}$, (b) $R_{s}=5\sigma_{mw}$,
(c) $R_{s}=9\sigma_{mw}$, (d) $R_{s}=13\sigma_{mw}$, (e)
$R_{s}=15\sigma_{mw}$, (f) $R_{s}=17\sigma_{mw}$. As the radius $R_{s}$
increases one observes the formation of more ’nanobubbles’ next to the solute.

We note that the complex bubble nanostructure that we observe in GCMC
simulations is, of course, not accessible to either 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 $\rho(r)$ reflecting
the average depletion arising from nanobubble formation. Nevertheless, one
expects on general grounds Evans, Stewart, and Wilding (2017) that such
theories will capture the correct large length scaling behaviour 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 behaviour 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. One barrier to
progress in explaining the phenomenology appears to be confusion in the
literature regarding the appropriate nomenclature for describing it. Previous
work points to the proximity of ambient water to liquid-vapor coexistence and
implies that this leads to a ‘dewetting transition’ Chandler (2005b); Patel
_et al._ (2012); Sarupria and Garde (2009) or sometimes a ‘drying transition’
Huang and Chandler (2000b, 2002) of water near a hydrophobic solute, or in the
region between two such solutes 111We note that the term ‘dewetting’ is most
commonly used in the quite different context of the rupturing of a thin film
of liquid on a planar substrate.. 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
hypothesised 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 the increasing 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 $\Delta\mu$, the solute
curvature $R_{s}^{-1}$ and the solute-solvent attraction $\epsilon_{sf}$ are
all sufficiently small for the system to qualify as ‘near’ to the critical
drying point. Accordingly we might expect behaviour to be driven by enhanced
local density fluctuations and, possibly, the emergence of a surface vapor
region. Suitable measures should exhibit near-critical scaling behaviour.

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 the density depletion, as
measured by $\rho(r)$, and the magnitude of local density fluctuations, as
measured by $\chi(r)$. The resulting mean field scaling predictions were
tested via classical DFT calculations for a generalised solute-solvent system.
From the evidence presented in the profiles of $\rho(r)$ and $\chi(r)$ in
figure 3, the comparison of binding potential predictions to DFT results in
figures 4, 5, 6 and 7 and the contour plots of figure 8, it is clear that the
predictions of the macroscopic binding potential analysis for the near-
critical scaling of $\ell_{eq}$ and $\chi(\ell_{eq};R_{s})$ are in agreement
with DFT results for a wide variety of solvophobic systems that extend down to
microscopic solute sizes. For microscopic solutes ($R_{s}\lesssim
10^{3}\sigma$), our results suggest that the magnitude of local density
fluctuations near a microscopic solute is most sensitive to changes in
$R_{s}$, with small variations in $T$, $\delta\mu$ and solvophobicity, as
measured by $\epsilon_{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 behaviour 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 behaviour, e.g. the density depletion and fluctuations. As density
fluctuations are sometimes conjectured to facilitate protein folding, this
insight might potentially provide 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 is expected
to apply in such systemsEvans, Stewart, and Wilding (2017); Coe (2021a).

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 Evans, Stewart, and
Wilding (2017). Thus despite the impression conveyed by simulation
measurements of the density profile $\rho(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
$R_{s}^{-1},\delta\mu,\delta\epsilon_{sf}$ are all sufficiently small that the
largest bubbles encompass much of the solute’s surface area. Of course, here
we shall be close to critical drying. We further further that our simulation
results show no evidence of the strongly hydrogen bonded hydration shell which
has been proposed to form around small hydrophobic solutes, e.g.
ref.Bischofberger  _et al._ (2014). 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
rationalising 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 occur in protein dynamics.

###### Acknowledgements.

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

## References

  * Huang and Chandler (2002) D. M. Huang and D. Chandler, “The hydrophobic effect and the influence of solute-solvent attractions,” J. Phys. Chem. B 106, 2047–53 (2002).
  * Oleinikova and Brovchenko (2012) 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).
  * Mittal and Hummer (2008) J. Mittal and G. Hummer, “Static and dynamic correlations in water at hydrophobic interfaces,” PNAS 105, 20130–20135 (2008).
  * Sarurpria and Garde (2009) S. Sarurpria and S. Garde, “Quantifying water density fluctuations and compressibility of hydration shells of hydrophobic solutes and proteins,” PRL 103, 037803 (2009).
  * Patel  _et al._ (2012) 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–503 (2012).
  * Vaikuntanathan  _et al._ (2016) S. Vaikuntanathan, G. Rotskoff, A. Hudson, and P. L. Geissler, “Necessity of capillary modes in a minimal model of nanoscale hydrophobic solvation,” PNAS 113, E2224–E2230 (2016).
  * Bischofberger  _et al._ (2014) 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,” Scientific Reports 4, 4377 (2014).
  * Rego and Patel (2022) N. B. Rego and A. J. Patel, “Understanding hydrophobic effects: Insights from water density fluctuations,” Annual Review of Condensed Matter Physics 13, 303–324 (2022).
  * Chandler (2005a) D. Chandler, “Interfaces and the driving force of hydrophobic assembly,” Nature 437, 640–647 (2005a).
  * Qvist  _et al._ (2008) J. Qvist, M. Davidovic, D. Hamelberg, and B. Halle, “A dry ligand-binding cavity in a solvated protein,” PNAS 105, 6296–301 (2008).
  * Lum, Chandler, and Weeks (1999) K. Lum, D. Chandler, and J. D. Weeks, “Hydrophobicity at small and large length scales,” J. Phys. Chem. B 103, 4570–7 (1999).
  * Southall and Dill (2000) N. T. Southall and K. A. Dill, “The mechanism of hydrophobic solvation depends on solute radius,” J. Phys. Chem. B 104, 1326–31 (2000).
  * Stillinger (1973) F. H. Stillinger, “Structure in aqueous solutions of nonpolar solutes from the standpoint of scaled-particle theory,” J. Solution Chem. 2, 141–58 (1973).
  * Huang and Chandler (2000a) D. M. Huang and D. Chandler, “Cavity formation and the drying transition in the Lennard-Jones fluid,” Phys. Rev. E 61, 1501–1506 (2000a).
  * Sarupria and Garde (2009) 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).
  * Acharya  _et al._ (2010) 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–65 (2010).
  * Mamatkulov and Khabibullaev (2004) S. I. Mamatkulov and P. K. Khabibullaev, “Water at hydrophobic substrates: Curvature, pressure and temperature effects,” Langmuir 20, 4756–63 (2004).
  * Patel, Varilly, and Chandler (2010) A. J. Patel, P. Varilly, and D. Chandler, “Fluctuations of water near extended hydrophobic and hydrophilic surfaces,” J. Phys. Chem. B 114, 1632–7 (2010).
  * Mezger  _et al._ (2006) 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,” PNAS 103, 18401–18404 (2006).
  * Mezger  _et al._ (2010) 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).
  * Ocko, Dhinojwala, and Daillant (2008) B. M. Ocko, A. Dhinojwala, and J. Daillant, “Comment on “how water meets a hydrophobic surface”,” Phys. Rev. Lett. 101, 039601 (2008).
  * Chattopadhyay  _et al._ (2010) 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).
  * Willard and Chandler (2014) 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).
  * Evans and Wilding (2015) R. Evans and N. B. Wilding, “Quantifying density fluctuations in water at a hydrophobic surface: Evidence for critical drying,” PRL 115, 016103 (2015).
  * Evans, Stewart, and Wilding (2016) R. Evans, M. C. Stewart, and N. B. Wilding, “Critical drying of liquids,” PRL 117, 176102 (2016).
  * Evans, Stewart, and Wilding (2017) 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).
  * Evans, Stewart, and Wilding (2019) R. Evans, M. C. Stewart, and N. B. Wilding, “From hydrophilic to superhydrophobic surfaces: a unified picture of the wetting and drying of liquids,” PNAS 116, 23901–8 (2019).
  * Evans and Stewart (2015) 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).
  * Eckert  _et al._ (2020) T. Eckert, N. C. X. Stuhlmüller, F. Sammüller, and M. Schmidt, “Fluctuation profiles in inhomogeneous fluids,” PRL 125, 268004 (2020).
  * Coe, Evans, and Wilding (2022a) M. K. Coe, R. Evans, and N. B. Wilding, “Measures of fluctuations for a liquid near critical drying,” Phys. Rev. E 105, 044801 (2022a).
  * Cerdeiriña  _et al._ (2011) 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–3 (2011).
  * Ebner, Saam, and Sen (1985) 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).
  * Ebner and Saam (1987) C. Ebner and W. Saam, “Effect of long-range forces on wetting near bulk critical temperatures: An ising-model study,” Phys. Rev. B 35, 1822–1834 (1987).
  * Nightingale, Saam, and Schick (1984) 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).
  * Coe, Evans, and Wilding (2022b) 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 (2022b).
  * Coe (2021a) M. K. Coe,  _Hydrophobicity Across Length Scales: The Role of Surface Criticality_ , Ph.D. thesis, University of Bristol (2021a).
  * Evans (1992) R. Evans, “Density functionals in the theory of nonuniform fluids,” in  _Fundamentals of Inhomogeneous Fluids_ , edited by D. Henderson (Marcel Dekker Inc., 1992) pp. 85–175.
  * Stewart (2006) M. C. Stewart,  _Effect of Substrate Geometry on Interfacial Phase Transitions_ , Ph.D. thesis, University of Bristol (2006).
  * Roth (2010) R. Roth, “Fundamental measure theory for hard-sphere mixtures: a review,” J. Phys.: Condens. Matter 22, 063102 (2010).
  * Molinero and Moore (2009) V. Molinero and E. B. Moore, “Water modelled as an intermediate element between carbon and silicon,” J. Phys. Chem. B 113, 4008–16 (2009).
  * Stewart and Evans (2005a) M. C. Stewart and R. Evans, “Critical drying at a spherical substrate,” J. Phys.: Condens. Matter 17, S3499–505 (2005a).
  * Stewart and Evans (2005b) M. C. Stewart and R. Evans, “Wetting and drying at a curved substrate: Long-ranged forces,” Phys. Rev. E 71, 011602 (2005b).
  * Evans, Henderson, and Roth (2004) R. Evans, J. R. Henderson, and R. Roth, “Nonanalytic curvature contributions to solvation free energies: Influence of drying,” J. Chem. Phys. 121, 12074–84 (2004).
  * Dietrich (1988) S. Dietrich, “Wetting phenomena,” in  _Phase Transitions and Critical Phenomena_ , Vol. 12, edited by C. Domb and J. L. Lebowitz (Academic Press, 1988).
  * Schick (1990) 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).
  * Bieker and Dietrich (1998) T. Bieker and S. Dietrich, “Wetting of curved surfaces,” Physica A 252, 85–137 (1998).
  * Dietrich and Napiórkowski (1991) S. Dietrich and M. Napiórkowski, “Analytic results for wetting transitions in the presence of van der waals tails,” Phys. Rev. A 43, 1861–85 (1991).
  * Coe (2021b) M. K. Coe, “cDFT Package,” https://github.com/marykcoe/cDFT_Package (2021b), online.
  * Evans (1979) R. Evans, “The nature of the liquid-vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids,” Advances in Physics 28, 143–200 (1979).
  * Sullivan and da Gama (1986) 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).
  * Coe, Evans, and Wilding (2022c) M. K. Coe, R. Evans, and N. B. Wilding, “The coexistence curve and surface tension of a monatomic water model,” The Journal of Chemical Physics 156, 154505 (2022c).
  * Coe (2021c) M. K. Coe, “GCMC Package for mw water,” https://github.com/marykcoe/Mont_Carlo (2021c), online.
  * (53) 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, 20899) accessed 29 June 2021.
  * Jamadagni, Godawat, and Garde (2011) S. N. Jamadagni, R. Godawat, and S. Garde, “Hydrophobicity of proteins and interfaces: Insights from density fluctuations,” Annu. Rev. Chem. Biomol. Eng. 2, 147–71 (2011).
  * Huang, Geissler, and Chandler (2001) D. M. Huang, P. L. Geissler, and D. Chandler, “Scaling of hydrophobic solvation free energies,” J. Phys. Chem. B 105, 6704–9 (2001).
  * Schnupf and Brady (2017) U. Schnupf and J. W. Brady, “Water structuring above solutes with planar hydrophobic surfaces,” Phys. Chem. Chem. Phys. 19, 11851–63 (2017).
  * Tyrrell and Attard (2001) J. W. G. Tyrrell and P. Attard, “Images of nanobubbles on hydrophobic surfaces and their interactions,” PRL 87, 176104 (2001).
  * Steitz  _et al._ (2003) 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 subtrate,” Langmuir 19, 2409–18 (2003).
  * Chandler (2005b) D. Chandler, “Interfaces and the driving force of hydrophobic assembly,” Nature 437, 640–647 (2005b).
  * Huang and Chandler (2000b) D. M. Huang and D. Chandler, “Cavity formation and the drying transition in the Lennard-Jones fluid,” Phys. Rev. E 61, 1501–6 (2000b).
  * Note (1) We note that the term ‘dewetting’ is most commonly used in the quite different context of the rupturing of a thin film of liquid on a planar substrate.