Pathways to dewetting in hydrophobic confinement
Richard C. Remsinga,1, Erte Xia,1, Srivathsan Vembanurb, Sumit Sharmac, Pablo G. Debenedettic, Shekhar Gardeb,
and Amish J. Patela,2

aDepartment of Chemical & Biomolecular Engineering, University of Pennsylvania, Philadelphia, PA 19104; bHoward P. Isermann Department of Chemical &
Biological Engineering, and Center for Biotechnology and Interdisciplinary Studies, Rensselaer Polytechnic Institute, Troy, NY 12180; and cDepartment of
Chemical & Biological Engineering, Princeton University, Princeton, NJ 08544

Edited by Ken A. Dill, Stony Brook University, Stony Brook, NY, and approved May 19, 2015 (received for review February 16, 2015).

Liquid water can become metastable with respect to its vapor in
hydrophobic confinement. The resulting dewetting transitions are
often impeded by large kinetic barriers. According to macroscopic
theory, such barriers arise from the free energy required to
nucleate a critical vapor tube that spans the region between two
hydrophobic surfaces—tubes with smaller radii collapse, whereas
larger ones grow to dry the entire confined region. Using exten-
sive molecular simulations of water between two nanoscopic hydro-
phobic surfaces, in conjunction with advanced sampling techniques,
here we show that for intersurface separations that thermodynami-
cally favor dewetting, the barrier to dewetting does not correspond
to the formation of a (classical) critical vapor tube. Instead, it corre-
sponds to an abrupt transition from an isolated cavity adjacent to one
of the confining surfaces to a gap-spanning vapor tube that is already
larger than the critical vapor tube anticipated by macroscopic theory.
Correspondingly, the barrier to dewetting is also smaller than the
classical expectation. We show that the peculiar nature of water
density fluctuations adjacent to extended hydrophobic surfaces—
namely, the enhanced likelihood of observing low-density fluctua-
tions relative to Gaussian statistics—facilitates this nonclassical
behavior. By stabilizing isolated cavities relative to vapor tubes,
enhanced water density fluctuations thus stabilize novel path-
ways, which circumvent the classical barriers and offer diminished
resistance to dewetting. Our results thus suggest a key role for
fluctuations in speeding up the kinetics of numerous phenomena
ranging from Cassie–Wenzel transitions on superhydrophobic
surfaces, to hydrophobically driven biomolecular folding and assembly.

capillary evaporation | fluctuations | kinetic barriers | assembly

The favorable interactions between two extended hydrophobicsurfaces drive numerous biomolecular and colloidal assem-
blies (1–5), and have been the subject of several theoretical,
computational, and experimental inquiries (6–22). Examples include
the association of small proteins to form multimeric protein com-
plexes, of amphiphlic block copolymers, dendrimers, or proteins to
form vesicular suprastructures, and of patchy colloidal particles into
complex crystalline lattices (23–27). When two such hydrophobic
surfaces approach each other, water between them becomes meta-
stable with respect to its vapor at a critical separation, dc, that can be
quite large (8, 9, 28–30). For nanometer-sized surfaces at ambient
conditions, dc is proportional to the characteristic size of the hy-
drophobic object, whereas for micron-sized and larger surfaces,
dc ∼ 1  μm (29, 30). However, due to the presence of large kinetic
barriers separating the metastable wet and the stable dry states, the
system persists in the wet state, and a dewetting transition is trig-
gered only at much smaller separations (∼ 1 nm) (13, 22, 28, 30).
To uncover the mechanism of dewetting, a number of theoretical

and simulation studies have focused on the thermodynamics as well
as the kinetics of dewetting in the volume between two parallel
hydrophobic surfaces that are separated by a fixed distance, d < dc
(8, 10–16, 18–21). These studies have highlighted that the bottle-
neck to dewetting is the formation of a roughly cylindrical, critical
vapor tube spanning the region between the surfaces (11, 14, 15). A
barrier in the free energetics of vapor tube formation as a function
of tube radius is also supported by macroscopic interfacial ther-
modynamics, wherein the barrier arises primarily from a competition

between the favorable solid–vapor and unfavorable liquid–vapor
surface energies (Eq. 1 and Fig. 1). Thus, the classical mechanism for
the dewetting transition prescribes that a vapor tube that spans the
volume between the two surfaces must first be nucleated, and if
the vapor tube is larger than a certain critical size, it will grow until
the entire confined volume is dry (9).
Although it has been recognized that water density fluctua-

tions must play a crucial role in nucleating vapor tubes (14, 15),
the precise mechanism by which these tubes are formed is not
clear. To understand how vapor tubes are formed and to investigate
their role in the dewetting process, here we use molecular simula-
tions in conjunction with enhanced sampling methods (31, 32) to
characterize the free energetics of water density fluctuations in
the region between two nanoscopic hydrophobic surfaces. Such a
characterization of water density fluctuations in bulk water and
at interfaces has already provided much insight into the physics
of hydrophobic hydration and interactions (5, 13, 31–44). In
particular, both simulations and theory have shown that the
likelihood of observing low-density fluctuations adjacent to ex-
tended hydrophobic surfaces is enhanced relative to Gaussian
statistics (13, 31, 36–38, 42). Further, the intricate coupling be-
tween enhanced solvent fluctuations and dewetting kinetics has
been highlighted by both coarse-grained (45–47) and atomistic
simulations (48–51).
Here we show that such enhanced water density fluctuations

influence the pathways to dewetting in hydrophobic confinement
by stabilizing isolated cavities adjacent to one of the confining
surfaces with respect to vapor tubes. As the density in the con-
fined region is decreased, the stability of isolated cavities relative
to vapor tubes also decreases, and at a particular density, isolated

Significance

Dewetting in hydrophobic confinement plays an important role
in diverse phenomena, ranging from protein folding and as-
sembly, to the heterogeneous nucleation of vapor bubbles and
superhydrophobicity. Using molecular simulations, we find
that dewetting proceeds through the formation of isolated
cavities adjacent to one of the confining surfaces. These iso-
lated cavities are stabilized by enhanced water density fluctu-
ations, and their growth is uphill in free energy. Upon growing
to a certain size, the isolated cavities transition abruptly into
supercritical vapor tubes that span the confined region, and
grow spontaneously. Consequently, this nonclassical pathway
results in lower free energy barriers than anticipated by mac-
roscopic theory, with important implications for the kinetics of
dewetting and hence for water-mediated self-assembly.

Author contributions: S.V., S.S., P.G.D., S.G., and A.J.P. designed research; R.C.R., E.X., and
A.J.P. performed research; R.C.R., E.X., and A.J.P. analyzed data; and R.C.R., E.X., P.G.D., S.G.,
and A.J.P. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.
1R.C.R. and E.X. contributed equally to this work.
2To whom correspondence should be addressed. Email: amish.patel@seas.upenn.edu.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.
1073/pnas.1503302112/-/DCSupplemental.

www.pnas.org/cgi/doi/10.1073/pnas.1503302112 PNAS | July 7, 2015 | vol. 112 | no. 27 | 8181–8186

A
P
P
LI
ED

P
H
Y
SI
C
A
L

SC
IE
N
C
ES

http://crossmark.crossref.org/dialog/?doi=10.1073/pnas.1503302112&domain=pdf
mailto:amish.patel@seas.upenn.edu
http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental
http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental
www.pnas.org/cgi/doi/10.1073/pnas.1503302112


cavities abruptly transition to vapor tubes. Surprisingly, for d K dc,
that is, separations for which dewetting is thermodynamically
favorable, we find that the nascent vapor tubes formed from the
isolated cavities are already larger than the corresponding criti-
cal vapor tubes predicted by classical theory. Because the newly
formed vapor tube is supercritical, it grows spontaneously. Im-
portantly, because the formation of this supercritical vapor tube
involves a nonclassical pathway that circumvents the critical vapor
tube altogether, the process entails a smaller free energetic cost.
Our results thus point to smaller kinetic barriers to dewetting than
predicted by macroscopic theory.

Macroscopic Theory
According to classical interfacial thermodynamics, the free energy
for creating a cylindrical vapor tube of radius r, which spans the
volume between two surfaces separated by a distance d, is given by:

ΔGthðr; dÞ = π
�
r2dΔP + 2rdγ + 2r2γ cos θ + 4rλ

�
, [1]

where ΔP is the difference between the system pressure and the
saturation pressure, γ is the liquid–vapor surface tension, θ is the
contact angle, and λ is the line tension. For nanoscopic surfaces,
the pressure-volume contribution is negligible at ambient condi-
tions (29, 52, 53), whereas the line tension contribution can be
important (20, 54). The term containing cos θ is negative for
hydrophobic surfaces and favors dewetting, whereas the term
corresponding to formation of the vapor–liquid area is unfavor-
able. The functional form of ΔGthðr; dÞ given in Eq. 1 is illus-
trated in Fig. 1D for three d values. In each case, a barrier
separates the liquid ðr = 0Þ and vapor (large r) basins, supporting
the notion of dewetting mediated by the nucleation and growth
of a vapor tube; both the critical vapor tube radius and the
barrier height increase with increasing d.

Density-Dependent Free Energy Profiles Feature a Kink
To investigate how water density fluctuations influence the
mechanism of dewetting in hydrophobic confinement, here we
perform molecular dynamics simulations of water confined be-
tween two roughly square hydrophobic surfaces of size L = 4 nm,
separated by a distance, d, as shown in Fig. 1A. Water in the
confined region is in equilibrium with a reservoir of water, which
in turn is in coexistence with its vapor ðΔP = 0Þ (31, 48). We

characterize the statistics of water density fluctuations in the
confined volume using indirect umbrella sampling (INDUS) (31,
32); that is, we estimate the free energy, ΔG, of observing N
water molecules in that volume. The free energy, ΔGðN; dÞ, thus
estimated is shown in Fig. 2A for a range of separations, d, with
the free energy of the liquid basin, N = Nliq, being set to zero in
each case. Over the entire range of separations considered, the
free energy profile displays distinctive liquid (high N) and vapor
(low N) basins with barriers separating them. Interestingly, the
free energy profiles also feature a kink, that is, an abrupt change
in the slope of ΔGðN; dÞ is observed at a particular value of
N between the liquid and vapor basins; we refer to this value of N
as Nkink. This discontinuity in slope is seen more clearly in the
derivatives of the free energy, shown in Fig. 2B. Small errors in ΔG
are amplified if simple finite differences are used to evaluate the
derivatives; we therefore smooth the free energy profiles before
evaluating the derivatives. Details of the smoothing procedure as
well as the unsmoothed derivatives are shown in the SI Appendix.

Kink Separates the Vapor Tube and Isolated Cavity
Ensembles
To investigate the significance of the kink in the free energy, we
characterize configurations corresponding to N on either side of
Nkink. We do so by building upon the instantaneous interface
method of Willard and Chandler (55) to identify isosurfaces

A C

B

D

Fig. 1. (A–C) Simulation snapshots of water (shown in red/white) in confine-
ment between two square hydrophobic surfaces (shown in cyan) of size L = 4 nm
that are separated by a distance of d = 20 Å; configurations highlighting (A) the
liquid basin, (B) a cylindrical vapor tube of radius, r, that spans the confined
region, and (C) the vapor basin are shown. In the front views, only one of the
confining surfaces is shown. (D) Macroscopic theory predicts a free energetic
barrier to vapor tube formation (Eq. 1), suggesting that a vapor tube larger than
a critical size must be nucleated before dewetting can proceed.

A

B

Fig. 2. (A) The simulated free energy profiles, βΔGðN; dÞ, as a function of the
number of water molecules between the surfaces, N, display marked kinks
(highlighted by circles). Here, β = 1=kBT, with kB being the Boltzmann constant
and T being the temperature. The size of the largest error bar is also shown for
d = 23 Å and N = 400. (B) The kinks are also apparent in the smoothed derivatives
of the free energy profiles, which display a sharp decrease in the vicinity of Nkink.

8182 | www.pnas.org/cgi/doi/10.1073/pnas.1503302112 Remsing et al.

http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental/pnas.1503302112.sapp.pdf
www.pnas.org/cgi/doi/10.1073/pnas.1503302112


that encompass the dewetted regions, that is, the regions from
which water is absent. The details of the method are included in
the SI Appendix. As illustrated in Fig. 3A for d = 20 Å and
N = Nkink − 12, characteristic configurations with N K Nkink con-
tain a vapor tube, that is, the dewetted region (in purple) clearly
spans the confined volume between the two surfaces (side view,
left). Water molecules (not shown for clarity) occupy the entire re-
gion between the surfaces not shown in purple and are also pre-
sent outside the confinement region. In contrast, for configurations
with N J Nkink, isolated cavities are observed adjacent to one sur-
face or the other; however, as seen in Fig. 3B for N = Nkink + 3, the
cavities do not span the region between the surfaces to form
vapor tubes. Movies corresponding to the configurations shown
in Fig. 3 A and B are included in the Supporting Information as
Movies S1–S4.

The configurations shown in Fig. 3 A and B suggest that
N = Nkink marks the boundary between the vapor tube and iso-
lated cavity ensembles. To put this notion on a quantitative
footing, we define indicator functions, htube and hcav, which are 1
if a given configuration has a vapor tube or an isolated cavity,
respectively, and 0 otherwise. The average value of the indicator
function hhtubeiN, subject to the constraint that the number of
water molecules in confinement is N, is shown in Fig. 3C for the
entire range of N values, and for several separations, d. Because
the configurations were generated in the presence of biasing
potentials, care must be exercised in evaluating the averages
shown in Fig. 3C; details of the averaging procedure as well as
the criteria used in the definition of the indicator functions can
be found in the SI Appendix. For a given separation, hhtubeiN,
which is the probability of observing a vapor tube conditional on
the number of waters in confinement being N, is a sigmoidal
function of N, decreasing sharply from 1 at low N to 0 for high N.
We define Ntube to be the value of N at which hhtubeiN undergoes
a sharp transition; in particular, where hhtubeiN crosses 0.5. As
shown in Fig. 3D, for the entire range of d values studied here
ð11  Å ≤ d ≤ 25  ÅÞ, Ntube is equal to Nkink, formalizing the notion
that the kink in ΔGðN; dÞ demarcates configurations that display
vapor tubes and those that do not.
Analogous to hhtubeiN, the conditional average hhcaviN quan-

tifies the fraction of configurations with N confined waters,
which feature isolated cavities. For all separations, hhcaviN dis-
plays a sharp increase, followed by a gradual decrease; see Fig.
3E. The sharp increase occurs in the vicinity of N = Nkink as vapor
tubes give way to isolated cavities, whereas the gradual decrease
corresponds to a crossover from isolated cavities to uniform
configurations. To specify the location of this crossover, we de-
fine Ncav to be the value of N where hhcaviN crosses 0.5 (with a
negative slope). Despite these common features in the functional
form of hhcaviN, there are subtle differences in hhcaviN at small
and large separations, which nevertheless have interesting con-
sequences. For the largest separations, a well-defined plateau at
hhcaviN =   1 separates the sharp increase in hhcaviN and its gradual
decrease; this plateau demarcates the range of N values, which
reliably feature isolated cavities. In contrast, the plateau at
hhcaviN =   1 is absent for the smaller separations. Instead, the
sharp increase in hhcaviN and its gradual decrease overlap in their
range of N values. This overlap suggests that configurations
corresponding to N = Nkink are not limited to those with vapor
tubes or isolated cavities, but could also be uniform. Given that
hhtubeiN = 0.5 at N = Ntube by definition, a value of hhcaviNtube <   0.5
would correspond to a nonzero likelihood of observing uniform
configurations with N = Ntube = Nkink. As shown in Fig. 3F, that
is indeed the case for d ≤ 16 Å, suggesting that although the
nucleation of a vapor tube must proceed through the formation
of isolated cavities for d > 16 Å, vapor tubes may be nucleated
directly from uniform configurations for smaller separations.

Free Energetics of Vapor Tubes and Isolated Cavities
Given the abrupt change in the dewetted morphologies at
N   =   Nkink, we expect the functional form of the free energetics
for N > Nkink and N < Nkink to be different. Fig. 3 C and D
collectively show that configurations with N K Nkink feature a va-
por tube spanning the confined region, consistent with classical
arguments. Although the vapor tube undergoes extensive shape
fluctuations, we find its average shape to be roughly cylindrical.
Coarse-grained density maps of select configurations reflecting the
average vapor tube shape are included in the SI Appendix. To
compare the free energetics of the vapor tubes obtained from our
simulations to macroscopic theory, we first transform the number
of waters in the confined region, N, to an approximate vapor tube
radius, r, using the simple relation, πr2=L2 = ðNliq − NÞ=Nliq. The
values of vapor tube radii thus obtained are consistent with the
average radii of the vapor tubes observed in our simulations for

A

C D

E F

B

Fig. 3. Instantaneous interfaces encompassing dewetted regions (shown in
purple) between the hydrophobic surfaces (shown in cyan) separated by
d = 20 Å highlight the presence of (A) a vapor tube for N = Nkink − 12, and
(B) an isolated cavity for N = Nkink + 3. Water molecules not shown for clarity.
(C) Average of the binary vapor tube indicator function, htube, conditioned
on the number of waters in confinement being N, displays a sharp transition
from 1 to 0 as N is increased. The color scheme is the same as that in Fig. 2.
The value of N corresponding to hhtubeiN = 0.5 (dashed line) is defined as
Ntube. (D) Ntube is identical to the location of the kink in the free energy
profiles, Nkink, as shown by the agreement between the simulation data and
a straight line (dashed). This agreement confirms that the kink demarcates
conformations with and without vapor tubes. (E) Conditional average of the
isolated cavity indicator function, hhcaviN, shows a sharp increase in the vi-
cinity of Nkink (the square symbols correspond to Ntube), followed by a
gradual decrease at larger N values, and eventually vanishes around N = Nliq.
(F) For the larger d values, hhtubeiNtube = hhcaviNtube = 0.5. However, for the
smaller d values, hhcaviNtube < 0.5, suggesting the possibility of direct vapor
tube nucleation without isolated cavities as intermediates.

Remsing et al. PNAS | July 7, 2015 | vol. 112 | no. 27 | 8183

A
P
P
LI
ED

P
H
Y
SI
C
A
L

SC
IE
N
C
ES

http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental/pnas.1503302112.sapp.pdf
http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental/pnas.201503302SI.pdf?targetid=nameddest=STXT
http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental
http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental
http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental
http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental/pnas.1503302112.sapp.pdf
http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental/pnas.1503302112.sapp.pdf


N < Nkink, as shown in the SI Appendix. Having a one-to-one re-
lation between N and r allows us to transform the simulated free
energies, ΔGðN; dÞ, into r-dependent free energies, ΔGðr; dÞ in the
region rkink < r < L=2; the corresponding free energies are shown as
symbols in Fig. 4A. The lines are fits to the ΔGðrÞ data using the
macroscopic expectation, ΔGfitðrÞ = ΔGthðrÞ + 2kBT lnð1 − 2r=LÞ,
where the logarithmic term corresponds to the translational en-
tropy of the vapor tube. ΔGðrÞ is fit separately for each d value and
yields values of γ and λ that are reasonable. Values of γ are in the
range of 12.2 − 15.8  kBT=nm2, comparable to the reported value of
14.5  kBT=nm2 for the water model that we use (56). Our fits yield
−λ=γ in the 6.5 − 7.5  Å range, in accord with a recently reported
experimental value of λ = −30 pN (54), which yields −λ=γ = 4.2  Å.
These agreements are remarkable considering the simplicity of the
model that we use as well as the assumptions that we make (cylin-
drical vapor tube shape, constant surface tension independent of
vapor tube curvature, etc.), and suggest that the energetics of the
vapor tube are well-described by classical macroscopic theory. Fur-
ther details of our fitting procedure, the values of the fit parameters

for each d, as well as our attempts to fit the simulation data to other
reasonable expressions of GthðrÞ can be found in the SI Appendix.
To investigate the free energetics of the isolated cavity en-

semble, in Fig. 4B, we focus on ΔGðN > Nkink; dÞ. In the liquid
basin, that is, in the vicinity of N = Nliq, the free energy (symbols)
is parabolic (solid lines), indicating that the underlying density
fluctuations are Gaussian. Although ΔG remains harmonic for
N > Nliq, it crosses over to being roughly linear (dashed lines) for
N < Nliq. Such a crossover from parabolic to linear has also been
observed adjacent to single extended hydrophobic surfaces (31,
41, 42), and corresponds to the interfacial water undergoing a
collective dewetting transition to expel water from a nanometer-
sized cavity (13, 38, 42). Indeed, as shown in Fig. 4B, the location
of the crossover agrees well with the corresponding Ncav values
(squares) discussed in the previous section. Interestingly, the
slopes of the linear fat tails are similar for all d values (in the
range of Nkink + 20 ≤ N ≤ Nliq − 50); the dashed lines shown in
Fig. 4B are linear fits with a slope of −0.396    kBT per water.
Additional details of the fitting procedure and the values of the
parameters obtained are provided in the SI Appendix. The dif-
ference between the values of Nliq and the x intercepts of the
linear fits is also approximately the same for all d values, and is
equal to 30  ±   5 waters. Thus, the free energy for forming an
isolated cavity of a given size (as quantified by the number of
waters displaced from the confined region, Nliq − N), is independent
of the separation between the surfaces; that is, the free energetics of
isolated cavity formation adjacent to one hydrophobic surface are
largely unaffected by the presence of the other confining surface. In
contrast, the free energetics of vapor tube formation clearly depend
on the intersurface separation, d. As a result, the location of the
kink, where isolated cavities become metastable with respect to
vapor tubes, also depends on d.

Nonclassical Dewetting Mechanism Reduces Barriers
Fig. 5A summarizes our findings for the dewetting mechanism
presented thus far; in addition to the simulated ΔGðN; dÞ for
d = 15  Å, it highlights the metastable branches of the vapor tube
and isolated cavity ensemble free energies, anticipated from the
fits shown in Fig. 4. It is clear that the system minimizes its free
energy at all times by staying on the branch with the lower free
energy; at N = Nkink, where the two free energy profiles in-
tersect, the system jumps from the isolated cavity to the vapor
tube ensemble. Importantly, the nonclassical path leading up
to the formation of nascent vapor tubes ðN > NkinkÞ can result
in smaller barriers to dewetting than anticipated by classical
theory, as shown in Fig. 5B. For d = 14  Å, the classical barrier
(critical vapor tube) appears in the metastable segment of the
free energy profile. The system thus circumvents the classical
barrier, and instead adopts the path involving isolated cavities,
which give way to vapor tubes only at N = Nkink; these nascent
vapor tubes are larger than the critical vapor tube, so their
subsequent growth is downhill in energy. Thus, the barrier to
dewetting is the free energetic cost for forming these nascent,
supercritical vapor tubes, which is clearly smaller than the
classical barrier. In Fig. 5C, we illustrate that the nascent vapor
tubes are not supercritical for all separations; for d = 23  Å, the
classical barrier appears in the stable segment of the free en-
ergy profile. Thus, although the kink in ΔGðN; dÞ again marks
the formation of a vapor tube for d = 23  Å, the vapor tube
formed is smaller than the critical vapor tube, and must grow
further in a process that is uphill in energy, before dewetting
can proceed. As a result, the barrier to dewetting coincides
with that predicted by macroscopic theory.
To uncover the separation at which the system transitions

from a supercritical to a subcritical nascent vapor tube, in Fig.
5D, we plot Nkink and Nmax as a function of d. Here, Nmax cor-
responds to the value of N between the liquid and vapor basins
where ΔGðN; dÞ is the highest. For small values of d, Nmax = Nkink,

A

B

Fig. 4. (A) βΔGðN ≤ Nkink; dÞ is recast as βΔGðr; dÞ, the free energy to form a
vapor tube of radius, r. The points were obtained from the simulated free
energy profiles by using the relation πr2=L2 = ðNliq − NÞ=Nliq in the region
rkink < r < L=2, and the lines are fits to macroscopic theory. (B) The portion of
the free energy corresponding to the liquid basin ðN ≥ NkinkÞ is parabolic at
high N (Gaussian fluctuations), but linear at low N (fat tails in water number
distributions). The crossover is gradual and occurs in the vicinity of Ncav (square
symbols), that is, the value of N for which hhcaviN is 0.5 with a negative slope. The
linear regions have roughly the same slope for all separations.

8184 | www.pnas.org/cgi/doi/10.1073/pnas.1503302112 Remsing et al.

http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental/pnas.1503302112.sapp.pdf
http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental/pnas.1503302112.sapp.pdf
http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental/pnas.1503302112.sapp.pdf
www.pnas.org/cgi/doi/10.1073/pnas.1503302112


indicating that ΔGðN; dÞ is the highest at the kink, consis-
tent with the formation of supercritical nascent vapor tubes. For
larger values of d, there is an additional maximum in ΔG at
Nmax < Nkink, suggesting that the newly formed vapor tubes are
smaller than the corresponding critical vapor tubes. Interestingly,
the separation at which the system transitions from nonclassical
to classical dewetting barriers is close to the separation, dc, at
which there is coexistence between the liquid and vapor. Thus,
for the separations with a thermodynamically favorable driving
force for dewetting, the mechanism for dewetting is manifestly
nonclassical, corresponding to the formation of supercritical vapor
tubes from isolated cavities, and requiring a smaller free energetic
barrier than anticipated by macroscopic theory.

Discussion and Outlook
Our results highlight that water density fluctuations can play a
central role in determining the pathways to dewetting in hydro-
phobic confinement. Enhanced water density fluctuations in the
vicinity of hydrophobic surfaces stabilize isolated cavities relative
to vapor tubes for N > Nkink; the system resides in the classical
vapor tube ensemble only for N < Nkink. Our results thus high-
light a breakdown of classical nucleation theory that stems from
the stabilization of a nonclassical species (isolated cavities) and
leads to an alternative pathway with a lower barrier. Similar
findings have been reported for crystal nucleation, wherein a
nonclassical pathway is stabilized by liquidlike solute clusters,
which results in a reduced barrier (57, 58).
Although the free energetics of both the isolated cavity and vapor

tube ensembles are well described by the number of waters in
confinement, N, this simple order parameter may not be sufficient to
describe the transition from one ensemble to the other. Indeed,
ΔGðN; dÞ represents a projection of a complex free energy land-
scape onto the single parameter, N, and a kink in ΔGðN; dÞ strongly
suggests that the transition from an isolated cavity to a vapor tube
involves order parameter(s) that are orthogonal to N. Although
beyond the scope of the present work, it would be interesting to
uncover these order parameters that define the transition state en-
semble along with N. Determining these additional parameters will
require path sampling techniques in conjunction with a character-
ization of commitment probabilities (14); it is conceivable that ad-
ditional barriers may present themselves in these order parameters.
Dewetting in nanoscopic hydrophobic confinement plays an im-

portant role in biology; ranging from the assembly of multimeric
proteins and the collapse of the hydrophobic protein core, to the
vapor-lock gating of ion channels and the specific binding of ligands
to hydrophobic grooves on their binding partners. In particular,

recent work has highlighted the importance of including the solvent
coordinate, N, in describing the kinetics of hydrophobically driven
collapse and assembly (49–51, 59). Our results show that water
density fluctuations stabilize nonclassical pathways, which reduce
the barriers along the N coordinate, and should therefore enhance
the kinetics of dewetting-mediated biophysical phenomena.
Dewetting in hydrophobic confinement can also be important

in a host of nonbiological phenomena, ranging from heteroge-
neous nucleation of vapor bubbles and contact line pinning, to
the Cassie–Wenzel transitions on textured surfaces (60). These
phenomena involve intricate confinement geometries, which
could result in complex pathways involving one or more transi-
tions between various dewetted morphologies, akin to the tran-
sition between isolated cavities and vapor tubes observed here.
Fluctuation-mediated pathways ought to also reduce dewetting
barriers associated with these diverse phenomena.

Materials and Methods
We simulate the SPC/E (extended simple point charge) model of water in con-
finement between two square hydrophobic surfaces of size L = 4 nm, for a range
of separations, d (Fig. 1A), ranging from 11 to 25 Å, chosen to span the entire
range of d values with both liquid and vapor basins. The surfaces are composed
of 1,008 atoms each, which are arranged on a hexagonal lattice with a spacing
of 1.4 Å. The surface atoms interact with the water oxygens through the Len-
nard–Jones potential with the parameters, σ = 3.283 Å and e = 0.121 kJ/mol (see
refs. 20, 21 for further details). As shown in the SI Appendix, this well depth was
chosen so that a water droplet on the hydrophobic surface makes a contact
angle, θ ≈ 120°, which is characteristic of alkyl-terminated self-assembled
monolayer surfaces (37). For water, we have chosen the SPC/E model (61) be-
cause it adequately captures the experimentally known features of water such as
surface tension, compressibility, and vapor–liquid equation of state near ambient
conditions, which are important in the study of hydrophobic effects (5, 62). Our
simulations contain roughly 10,000–15,000 water molecules and were performed
in the canonical ensemble, thermostated at T = 300 K using the canonical ve-
locity rescaling thermostat (63). We use a periodic simulation box with the hy-
drophobic surfaces of interest fixed at the center of the box, and a buffering
liquid–vapor interface nucleated at the top of the box with the help of a wall of
purely repulsive particles. The buffering interface ensures that the system is at
the saturation pressure of SPC/E water at 300 K (31, 48); free energies obtained
with such a construct have been shown to be nearly indistinguishable from
those obtained in the NPT ensemble at a pressure of 1 bar (32). Short-ranged
interactions were truncated at 1 nm; whereas, long-ranged electrostatic
interactions were computed using the particle mesh Ewald method (64). The
bonds in water were constrained using the SHAKE algorithm (65). To study
the free energetics of dewetting, we select the cuboid shaped ðL × L × dÞ
observation volume between the hydrophobic surfaces, and estimate the free
energies, ΔGðN; dÞ, using the INDUS method (31, 32). Each biased simulation was
run for 6 ns and the first 1 ns was discarded for equilibration.

A B C D

Fig. 5. (A) The simulated βΔGðN; dÞ for d = 15  Å (solid) is shown along with the expected metastable branches of the free energies corresponding to the vapor tube
(dot-dashed) and isolated cavity (dashed) ensembles. The metastable branches are anticipated from the fits shown in Fig. 4. The system minimizes its free energy by
localizing to the ensemble with the lower free energy. (B) For d = 14 Å, the nascent vapor tube formed at the kink is larger than the critical vapor tube anticipated by
macroscopic theory, and therefore grows spontaneously. As a result, the corresponding barrier to dewetting is smaller than that predicted by macroscopic theory.
(C) For larger separations, here d = 23 Å, the newly formed vapor tube is subcritical, and has to grow larger for dewetting to proceed. (D) Comparison of the location
of the kink, Nkink, with the location of the barrier between the liquid and vapor basins, Nmax. For small d, the barrier (point of highest ΔG) occurs at the kink, so that
Nmax ≈ Nkink. In contrast, for larger d values, the barrier occurs in the vapor tube segment of the simulated free energy profile and corresponds to the classical critical
vapor tube, so that Nmax < Nkink. The transition from nonclassical to classical behavior occurs in the vicinity of the coexistence separation, dc.

Remsing et al. PNAS | July 7, 2015 | vol. 112 | no. 27 | 8185

A
P
P
LI
ED

P
H
Y
SI
C
A
L

SC
IE
N
C
ES

http://www.pnas.org/lookup/suppl/doi:10.1073/pnas.1503302112/-/DCSupplemental/pnas.1503302112.sapp.pdf


ACKNOWLEDGMENTS. R.C.R. and A.J.P. were supported in part by financial
support from the National Science Foundation (NSF) through a Seed Grant
from the University of Pennsylvania Materials Research Science and Engineering

Center (NSF UPENN MRSEC DMR 11-20901). S.G. was supported by the NSF
(CBET-1159990). P.G.D. gratefully acknowledges support from the NSF (CBET-
1263565 and CHE-1213343).

1. Tanford C (1973) The Hydrophobic Effect - Formation of Micelles and Biological
Membranes (Wiley Interscience, New York).

2. Kauzmann W (1959) Some factors in the interpretation of protein denaturation. Adv
Protein Chem 14:1–63.

3. Stillinger FH (1973) Structure in aqueous solutions of nonpolar solutes from the
standpoint of scaled-particle theory. J Solution Chem 2:141–158.

4. Israelachvili JN (2011) Intermolecular and surface forces (Academic, Waltham, MA),
revised third edition.

5. Chandler D (2005) Interfaces and the driving force of hydrophobic assembly. Nature
437(7059):640–647.

6. Israelachvili J, Pashley R (1982) The hydrophobic interaction is long range, decaying
exponentially with distance. Nature 300(5890):341–342.

7. Christenson HK, Claesson PM (1988) Cavitation and the interaction between macro-
scopic hydrophobic surfaces. Science 239(4838):390–392.

8. Bérard DR, Attard P, Patey GN (1993) Cavitation of a Lennard-Jones fluid between hard
walls, and the possible relevance to the attraction measured between hydrophobic sur-
faces. J Chem Phys 98:7236–7244.

9. Parker JL, Claesson PM, Attard P (1994) Bubbles, cavities, and the long-ranged at-
traction between hydrophobic surfaces. J Phys Chem 98:8468–8480.

10. Wallqvist A, Berne BJ (1995) Computer simulation of hydrophobic hydration forces on
stacked plates at short range. J Phys Chem 99:2893–2899.

11. Lum K, Luzar A (1997) Pathway to surface-induced phase transition of a confined
fluid. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics 56:R6283–R6286.

12. Lum K, Chandler D (1998) Phase diagram and free energies of vapor films and tubes
for a confined fluid. Int J Thermophys 19:845–855.

13. Lum K, Chandler D, Weeks JD (1999) Hydrophobicity at small and large length scales.
J Phys Chem B 103:4570–4577.

14. Bolhuis PG, Chandler D (2000) Transition path sampling of cavitation between mo-
lecular scale solvophobic surfaces. J Chem Phys 113:8154–8160.

15. Leung K, Luzar A, Bratko D (2003) Dynamics of capillary drying in water. Phys Rev Lett
90(6):065502.

16. Huang X, Margulis CJ, Berne BJ (2003) Dewetting-induced collapse of hydrophobic
particles. Proc Natl Acad Sci USA 100(21):11953–11958.

17. Urbic T, Vlachy V, Dill KA (2006) Confined water: A Mercedes-Benz model study.
J Phys Chem B 110(10):4963–4970.

18. Choudhury N, Pettitt BM (2007) The dewetting transition and the hydrophobic effect.
J Am Chem Soc 129(15):4847–4852.

19. Xu L, Molinero V (2010) Liquid-vapor oscillations of water nanoconfined between
hydrophobic disks: Thermodynamics and kinetics. J Phys Chem B 114(21):7320–7328.

20. Sharma S, Debenedetti PG (2012) Evaporation rate of water in hydrophobic con-
finement. Proc Natl Acad Sci USA 109(12):4365–4370.

21. Sharma S, Debenedetti PG (2012) Free energy barriers to evaporation of water in
hydrophobic confinement. J Phys Chem B 116(44):13282–13289.

22. Mastropietro DJ, Ducker WA (2012) Forces between hydrophobic solids in concen-
trated aqueous salt solution. Phys Rev Lett 108(10):106101.

23. Yu N, Hagan MF (2012) Simulations of HIV capsid protein dimerization reveal the
effect of chemistry and topography on the mechanism of hydrophobic protein as-
sociation. Biophys J 103(6):1363–1369.

24. Discher BM, et al. (1999) Polymersomes: Tough vesicles made from diblock co-
polymers. Science 284(5417):1143–1146.

25. Percec V, et al. (2010) Self-assembly of Janus dendrimers into uniform dendrimersomes
and other complex architectures. Science 328(5981):1009–1014.

26. Chen Q, Yan J, Zhang J, Bae SC, Granick S (2012) Janus and multiblock colloidal
particles. Langmuir 28(38):13555–13561.

27. Vargo KB, Parthasarathy R, Hammer DA (2012) Self-assembly of tunable protein su-
prastructures from recombinant oleosin. Proc Natl Acad Sci USA 109(29):11657–11662.

28. Berne BJ, Weeks JD, Zhou R (2009) Dewetting and hydrophobic interaction in physical
and biological systems. Annu Rev Phys Chem 60:85–103.

29. Cerdeirina CA, Debenedetti PG, Rossky PJ, Giovambattista N (2011) Evaporation
length scales of confined water and some common organic liquids. J. Phys. Chem.
Lett. 2:1000–1003.

30. Giovambattista N, Rossky PJ, Debenedetti PG (2012) Computational studies of pres-
sure, temperature, and surface effects on the structure and thermodynamics of
confined water. Annu Rev Phys Chem 63:179–200.

31. Patel AJ, Varilly P, Chandler D (2010) Fluctuations of water near extended hydro-
phobic and hydrophilic surfaces. J Phys Chem B 114(4):1632–1637.

32. Patel AJ, Varilly P, Chandler D, Garde S (2011) Quantifying density fluctuations in
volumes of all shapes and sizes using indirect umbrella sampling. J Stat Phys 145(2):
265–275.

33. Hummer G, Garde S, García AE, Pohorille A, Pratt LR (1996) An information theory
model of hydrophobic interactions. Proc Natl Acad Sci USA 93(17):8951–8955.

34. Garde S, Hummer G, García AE, Paulaitis ME, Pratt LR (1996) Origin of entropy con-
vergence in hydrophobic hydration and protein folding. Phys Rev Lett 77(24):
4966–4968.

35. Rajamani S, Truskett TM, Garde S (2005) Hydrophobic hydration from small to large
lengthscales: Understanding and manipulating the crossover. Proc Natl Acad Sci USA
102(27):9475–9480.

36. Mittal J, Hummer G (2008) Static and dynamic correlations in water at hydrophobic
interfaces. Proc Natl Acad Sci USA 105(51):20130–20135.

37. Godawat R, Jamadagni SN, Garde S (2009) Characterizing hydrophobicity of in-
terfaces by using cavity formation, solute binding, and water correlations. Proc Natl
Acad Sci USA 106(36):15119–15124.

38. Varilly P, Patel AJ, Chandler D (2011) An improved coarse-grained model of solvation
and the hydrophobic effect. J Chem Phys 134(7):074109.

39. Jamadagni SN, Godawat R, Garde S (2011) Hydrophobicity of proteins and interfaces:
Insights from density fluctuations. Annu Rev Chem Biomol Eng 2:147–171.

40. Patel AJ, et al. (2011) Extended surfaces modulate hydrophobic interactions of
neighboring solutes. Proc Natl Acad Sci USA 108(43):17678–17683.

41. Rotenberg B, Patel AJ, Chandler D (2011) Molecular explanation for why talc surfaces
can be both hydrophilic and hydrophobic. J Am Chem Soc 133(50):20521–20527.

42. Patel AJ, et al. (2012) Sitting at the edge: How biomolecules use hydrophobicity to
tune their interactions and function. J Phys Chem B 116(8):2498–2503.

43. Remsing RC, Weeks JD (2013) Dissecting hydrophobic hydration and association.
J Phys Chem B 117(49):15479–15491.

44. Remsing RC, Patel AJ (2015) Water density fluctuations relevant to hydrophobic hy-
dration are unaltered by attractions. J Chem Phys 142(2):024502.

45. ten Wolde PR, Chandler D (2002) Drying-induced hydrophobic polymer collapse. Proc
Natl Acad Sci USA 99(10):6539–6543.

46. Maibaum L, Chandler D (2003) A coarse-grained model of water confined in a hy-
drophobic tube. J Phys Chem B 107:1189–1193.

47. Willard AP, Chandler D (2008) The role of solvent fluctuations in hydrophobic as-
sembly. J Phys Chem B 112(19):6187–6192.

48. Miller TF, 3rd, Vanden-Eijnden E, Chandler D (2007) Solvent coarse-graining and the
string method applied to the hydrophobic collapse of a hydrated chain. Proc Natl
Acad Sci USA 104(37):14559–14564.

49. Li J, Morrone JA, Berne BJ (2012) Are hydrodynamic interactions important in the
kinetics of hydrophobic collapse? J Phys Chem B 116(37):11537–11544.

50. Mondal J, Morrone JA, Berne BJ (2013) How hydrophobic drying forces impact the
kinetics of molecular recognition. Proc Natl Acad Sci USA 110(33):13277–13282.

51. Setny P, Baron R, Michael Kekenes-Huskey P, McCammon JA, Dzubiella J (2013) Sol-
vent fluctuations in hydrophobic cavity-ligand binding kinetics. Proc Natl Acad Sci
USA 110(4):1197–1202.

52. Ashbaugh HS (2013) Solvent cavitation under solvophobic confinement. J Chem Phys
139(6):064702.

53. Altabet YE, Debenedetti PG (2014) The role of material flexibility on the drying
transition of water between hydrophobic objects: A thermodynamic analysis. J Chem
Phys 141:18C531.

54. Guillemot L, Biben T, Galarneau A, Vigier G, Charlaix É (2012) Activated drying in
hydrophobic nanopores and the line tension of water. Proc Natl Acad Sci USA 109(48):
19557–19562.

55. Willard AP, Chandler D (2010) Instantaneous liquid interfaces. J Phys Chem B 114(5):
1954–1958.

56. Vega C, de Miguel E (2007) Surface tension of the most popular models of water by
using the test-area simulation method. J Chem Phys 126(15):154707.

57. Gebauer D, Völkel A, Cölfen H (2008) Stable prenucleation calcium carbonate clusters.
Science 322(5909):1819–1822.

58. Gebauer D, Cölfen H (2011) Prenucleation clusters and non-classical nucleation. Nano
Today 6:564–584.

59. Mittal J, Hummer G (2012) Pair diffusion, hydrodynamic interactions, and available
volume in dense fluids. J Chem Phys 137(3):034110.

60. Kumar V, Sridhar S, Errington JR (2011) Monte Carlo simulation strategies for com-
puting the wetting properties of fluids at geometrically rough surfaces. J Chem Phys
135(18):184702.

61. Berendsen HJC, Grigera JR, Straatsma TP (1987) The missing term in effective pair
potentials. J Phys Chem 91:6269–6271.

62. Varilly P, Chandler D (2013) Water evaporation: A transition path sampling study.
J Phys Chem B 117(5):1419–1428.

63. Bussi G, Donadio D, Parrinello M (2007) Canonical sampling through velocity rescal-
ing. J Chem Phys 126(1):014101.

64. Essmann U, et al. (1995) A smooth particle mesh ewald method. J Chem Phys 103:
8577–8593.

65. Ryckaert JP, Ciccotti G, Berendsen HJC (1977) Numerical integration of the cartesian
equations of motion of a system with constraints: Molecular dynamics of n-alkanes.
J Comput Phys 23:327–341.

8186 | www.pnas.org/cgi/doi/10.1073/pnas.1503302112 Remsing et al.

www.pnas.org/cgi/doi/10.1073/pnas.1503302112