Pathways to dewetting in hydrophobic confinement

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 extensive molecular simulations of water between two nanoscopic hydrophobic surfaces, in conjunction with advanced sampling techniques, here we show that for inter-surface separations that thermodynamically favor dewetting, the barrier to dewetting does not correspond to the formation of a (classical) critical vapor tube. Instead, it corresponds 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 fluctuations relative to Gaussian statistics - facilitates this non-classical behavior. By stabilizing isolated cavities relative to vapor tubes, enhanced water density fluctuations thus stabilize novel pathways, 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.

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 < d c [8,[10][11][12][13][14][15][16][18][19][20][21].These studies have highlighted that the bottleneck 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 thermodynamics, wherein the barrier arises primarily from a competition between the favorable solid-vapor and unfavorable liquid-vapor surface energies (Equation 1 and Figure 1).Thus, the classical mechanism for the dewetting transition prescribes that a vapor tube, which 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].
While it has been recognized that water density fluctuations [14,15] must play a crucial role in nucleating vapor tubes, 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 simulations 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][32][33][34][35][36][37][38][39][40][41][42][43].In particular, both simulations and theory have shown that the likelihood of observing low density fluctuations adjacent to extended hydrophobic surfaces is enhanced relative to Gaussian statistics [13,31,[36][37][38]42].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 confined region is decreased, the stability of isolated cavities relative to vapor tubes also decreases, and at a particular density, isolated cavities abruptly transition to vapor tubes.Surprisingly, for d d c , 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 critical vapor tubes predicted by classical theory.Because the newly formed vapor tube is super-critical, it grows spontaneously.Importantly, because the formation of this super-critical vapor tube involves a non-classical 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.1), suggesting that a vapor tube larger than a critical size must be nucleated before dewetting can proceed.

I. 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: ∆G th (r; d) = π[r 2 d∆P + 2rdγ + 2r 2 γ 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 conditions [29,44,45], whereas the line tension contribution can be important [20,46].The term containing cos θ is negative for hydrophobic surfaces and favors dewetting, whereas the term corresponding to formation of the vapor-liquid area is unfavorable.The functional form of ∆G th (r; d) given in Equation 1 is illustrated in Figure 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.

II. 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 between two, roughly square, hydrophobic surfaces of size L = 4 nm, separated by a distance, d, as shown in Figure 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,47].We characterize the statistics of water density fluctuations in the confined volume using Indirect Umbrella Sampling (IN-DUS) [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 Figure 2a for a range of separations, d, with the free energy of the liquid basin, N = N liq , 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 N kink .This discontinuity in slope is seen more clearly in the derivatives of the free energy, shown in Figure 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.

III. 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 N kink .We do so by building upon the instantaneous interface method of Willard and Chandler [1] to identify iso-surfaces that encompass the dewetted regions, that is, the regions from which water is absent.The details of the method are included in the SI.As illustrated in Figure 3a for d = 20 Å and N = N kink − 12, characteristic configurations with N N kink contain 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 region between the surfaces not shown in purple and are also present outside the confinement region.In contrast, for configurations with N N kink , isolated cavities are observed adjacent to as N is increased.The color scheme is the same as that in Figure 2. The value of N corresponding to h tube N = 0.5 (dashed line) is defined as N tube .(d) N tube is identical to the location of the kink in the free energy profiles, N kink , confirming that the kink demarcates conformations with and without vapor tubes.(e) Conditional average of the isolated cavity indicator function, hcav N , shows a sharp increase in the vicinity of N kink (the square symbols correspond to N tube ), followed by a gradual decrease at larger N -values, and eventually vanishing around N = N liq .(f) For the larger d-values, h tube N tube = hcav N tube = 0.5.However, for the smaller dvalues, hcav N tube < 0.5, suggesting the possibility of direct vapor tube nucleation without isolated cavities as intermediates.
one surface or the other; however, as seen in Figure 3b for N = N kink + 3, the cavities do not span the region between the surfaces to form vapor tubes.Movies corresponding to the configurations shown in Figures 3a  and 3b can be found in the SI.
The configurations shown in Figures 3a and 3b suggest that N = N kink marks the boundary between the vapor tube and isolated cavity ensembles.To put this notion on a quantitative footing, we define indicator functions, h tube and h cav , 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 h tube N , subject to the constraint that the number of water molecules in confinement is N , is shown in Figure 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 Figure 3c; details of the averaging procedure as well as the criteria employed in the definition of the indicator functions can be found in the SI.For a given separation, h tube N , 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 N tube to be the value of N at which h tube N undergoes a sharp transition; in particular, where h tube N crosses 0.5.As shown in Figure 3d, for the entire range of d-values studied here (11 Å ≤ d ≤ 25 Å), N tube is equal to N kink , formalizing the notion that the kink in ∆G(N ; d) demarcates configurations that display vapor tubes and those that do not.
Analogous to h tube N , the conditional average h cav N quantifies the fraction of configurations with N confined waters, which feature isolated cavities.For all separations, h cav N displays a sharp increase, followed by a gradual decrease; see Figure 3e.The sharp increase occurs in the vicinity of N = N kink 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 define N cav to be the value of N where h cav N crosses 0.5 (with a negative slope).Despite these common features in the functional form of h cav N , there are subtle differences in h cav N at small and large separations, which nevertheless have interesting consequences.For the largest separations, a well-defined plateau at h cav N = 1 separates the sharp increase in h cav N and its gradual decrease; this plateau demarcates the range of N -values, which reliably feature isolated cavities.In contrast, the plateau at h cav N = 1 is absent for the smaller separations.Instead, the sharp increase in h cav N and its gradual decrease overlap in their range of N -values.This overlap suggests that configurations corresponding to N = N kink are not limited to those with vapor tubes or isolated cavities, but could also be uniform.Given that h tube N = 0.5 at N = N tube by definition, a value of h cav N tube < 0.5 would correspond to a non-zero likelihood of observing uniform configurations with N = N tube = N kink .As shown in Figure 3f, that is indeed the case d ≤ 16 Å, suggesting that while 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.

IV. FREE ENERGETICS OF VAPOR TUBES AND ISOLATED CAVITIES
Given the abrupt change in the dewetted morphologies at N = N kink , we expect the functional form of the free energetics for N > N kink and N < N kink to be different.Figures 3c and 3d collectively show that configurations with N N kink feature a vapor tube spanning the confined region, consistent with classical arguments.While 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.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, The values of vapor tube radii thus obtained are consistent with the average radii of the vapor tubes observed in our simulations for N < N kink , as shown in the SI.Having a one-to-one relation 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 r kink < r < L/2; the corresponding free energies are shown as symbols in Figure 4a.The lines are fits to the ∆G(r)-data using the macroscopic expectation, ∆G fit (r) = ∆G th (r) + 2k B T ln(1 − 2r/L), where the logarithmic term corresponds to the translational entropy 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 k B T /nm 2 , comparable to the reported value of 14.5 k B T /nm 2 for the water model that we employ [49].Our fits yield −λ/γ in the 6.5 − 7.5 Å range, in accord with a recently reported experimental value of λ = −30 pN [46], which yields −λ/γ = 4.2 Å.These agreements are remarkable considering the simplicity of the model that we employ as well as the assumptions that we make (cylindrical vapor tube shape, constant surface tension independent of vapor tube curvature, etc.), suggesting that the energetics of the vapor tube are well-described by classical macroscopic theory.Further 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 G th (r) can be found in the SI.
To investigate the free energetics of the isolated cavity ensemble, in Figure 4b, we focus on ∆G(N > N kink ; d).In the liquid basin, that is, in the vicinity of N = N liq , the free energy (symbols) is parabolic (solid lines), indicating that the underlying density fluctuations are Gaussian.While ∆G remains harmonic for N > N liq , it crosses over to being roughly linear (dashed lines) for N < N liq .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 Figure 4b, the location of the crossover agrees well with the corresponding N cav -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 N kink + 20 ≤ N ≤ N liq − 50); the dashed lines shown in Figure 4b are linear fits with a slope of −0.396 k B T per water.Additional details of the fitting procedure and the values of the parameters obtained are provided in the SI.The difference between the values of N liq and the xintercepts 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, N liq − 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 inter-surface separation, d.As a result, the location of the kink, where isolated cavities become metastable with respect to vapor tubes, also depends on d.

V. NON-CLASSICAL DEWETTING MECHANISM REDUCES BARRIERS
Figure 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 Figure 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 = N kink , where the two free energy profiles intersect, 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 > N kink ) can result in smaller barriers to dewetting than anticipated by classical theory, as shown in Figure 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 = N kink ; 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, super-critical vapor tubes, which is clearly smaller than the classical barrier.In Figure 5c, we illustrate that the nascent vapor tubes are not super-critical for all separations; for d = 23 Å, the classical barrier appears in the stable segment of the free energy profile.Thus, while 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.the free energy to form a vapor tube of radius, r.The points were obtained from the simulated free energy profiles by employing the relation πr 2 /L 2 = (N liq − N )/N liq in the region r kink < r < L/2, and the lines are fits to macroscopic theory.(b) The portion of the free energy corresponding to the liquid basin (N ≥ N kink ) 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 hcav N is 0.5 with a negative slope.The linear regions have roughly the same slope for all separations.
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 Figure 5d ∆G at N max < N kink , suggesting that the newly formed vapor tubes are smaller than the corresponding critical vapor tubes.Interestingly, the separation at which the system transitions from non-classical to classical dewetting barriers is close to the separation at which there is coexistence between the liquid and vapor, d c .Thus, for the separations with a thermodynamically favorable driving force for dewetting, the mechanism for dewetting is manifestly non-classical, corresponding to the formation of super-critical vapor tubes from isolated cavities, and requiring a smaller free energetic barrier than anticipated by macroscopic theory.

VI. DISCUSSION AND OUTLOOK
Our results highlight that water density fluctuations play a central role in the pathways to dewetting in hydrophobic confinement.Enhanced water density fluctuations in the vicinity of hydrophobic surfaces stabilize isolated cavities relative to vapor tubes for N > N kink ; the system resides in the classical vapor tube ensemble only for N < N kink .While the free energetics of both the isolated cavity and vapor tube ensembles are welldescribed 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 landscape 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 .While beyond the scope of the present work, it will be interesting to uncover these additional order parameters that define the tran-sition state ensemble.It is conceivable that additional barriers may present themselves in the parameters that are orthogonal to N .
Dewetting in nanoscopic hydrophobic confinement plays an important 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 [50][51][52][53].Our results show that water density fluctuations stabilize non-classical pathways, which reduce the barriers along the N -coordinate, and should therefore enhance the kinetics of dewetting-mediated biophysical phenomena.
Dewetting in hydrophobic confinement is also important in a host of non-biological phenomena, ranging from heterogeneous nucleation of vapor bubbles and contact line pinning, to the Cassie-Wenzel transitions on textured surfaces [54].These phenomena involve intricate confinement geometries, which could result in complex pathways involving one or more transitions between various dewetted morphologies, such as the isolated cavities and vapor tubes that we observe here.Fluctuation-mediated pathways ought to similarly reduce dewetting barriers associated with these diverse phenomena.

VII. MATERIALS AND METHODS
We simulate SPC-E water in confinement between two square hydrophobic surfaces of size L = 4 nm, for a range of separations, d (See Figure 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 1008 atoms each, and are arranged on a hexagonal lattice with a spacing of 1.4 Å.The surface atoms interact with the water oxygens through the Lennard-Jones potential with the parameters, σ = 3.283 Å and = 0.121 kJ/mol; see refs.[20,21] for further details.These interactions result in a water droplet contact angle θ ≈ 120 • on the hydrophobic surface, as shown in the SI.We have chosen the SPC-E model of water [55] since 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,56].Our simulations contain roughly 10,000 to 15,000 water molecules and were performed in the NVT ensemble, thermostatted at T = 300 K using the canonical velocity rescaling thermostat [57].We employ a periodic simulation box with the hydrophobic 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 300K [31,47]; 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].Shortranged interactions were truncated at 1 nm, whereas long-ranged electrostatic interactions were computed using the particle mesh Ewald method [58].The bonds in water were constrained using SHAKE [59].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 indirect umbrella sampling (INDUS) method [31,32].Each biased simulation was run for 6 ns and the first 1 ns was discarded for equilibration.
To demonstrate the existence of kinks in ∆G(N ; d) more clearly, we plot smoothed derivatives of ∆G(N ; d) in Figure 2b of the main text.Directly computing derivatives of ∆G(N ) using finite difference, for example, using ∂∆G/∂N = ∆G(N + 1) − ∆G(N ), results in noisy estimates due to the numerical error in ∆G(N ) (Figure S1).To minimize such effects, we first smooth ∆G(N ) using a rectangular window function, such that is the smoothed free energy.Smoothed derivatives are then estimated from finite differences of such smoothed free energies, (S3) Derivatives of both the unsmoothed and the smoothed (using a width of ∆N = 10) free energies for select values of d are compared in Figure S1.
An analogous smoothing of the first derivative was also performed in order to obtain the second derivative of the free energy, which in turn enables us to robustly estimate N kink as the location of the minimum in ∂ 2 β∆G/∂N 2 .The error in N kink was estimated using block averaging; the trajectories were divided into five blocks of 1 ns each, and N kink was obtained for each of the five blocks.As an example, the second derivatives for five such blocks are shown in Figure S2 for d = 14 • A .* RCR and EX are joint first authors who contributed equally to this work † amish.patel@seas.upenn.edu

II. OBTAINING INSTANTANEOUS INTERFACES THAT ENVELOP DEWETTED REGIONS
Here, we closely follow and build upon the approach for estimating instantaneous interfaces, originally developed by Willard and Chandler [1].In addition to considering the coarse-grained water density, we also include a coarse-grained density arising from the plate atoms into an overall, normalized coarse-grained density at (x, y, z) as follows: ρtotal (x, y, z; R) = ρwater (x, y, z; R) (S4) Here, R represents the positions of all the heavy atoms in the system in a given configuration, the coarse-grained water density is normalized by the corresponding bulk density, and the coarse-grained plate density by the maximum coarse-grained plate density that occurs at the center of the plate.Defined this way, ρtotal is approximately equal to unity everywhere in the liquid state and near unity at the center of the plates and in the interfacial region.Configurations containing dewetted regions (cavities) will have significantly smaller values of ρtotal that approach zero in the vicinity of the cavity.Therefore, the ρtotal = 0.5 iso-density surface serves as a convenient definition of the instantaneous interface, allowing us to readily visualize the position, size, and shape of the cavity; we use the Marching Cube algorithm [2] to identify the instantaneous interface.The coarse-grained density fields of the individual species are estimated using ρα (x, y, z; R) where α represents either the water oxygen atoms or the plate atoms, N α is the number of atoms of type α, and (x i , y i , z i ) correspond to the coordinates of atom i.
For each configuration obtained from our simulations, we set up a three-dimensional grid to compute the coarsegrained density field with 0.1 nm spacing in each dimension.The coarse-graining function φ(x) is chosen to be a Gaussian with a width of 0.24 nm, which was truncated at 0.7 nm, shifted down to make it continuous, and normalized.The particular characteristics of the instantaneous interfaces thus computed, as well as those of the dewetted regions, such as their exact shapes and volumes, will depend on the choices made in Equations S4 and S5, as well as the parameters chosen.A discussion of how these choices affect the instantaneous interface calculation and which choices are judicious is beyond the scope of this work, and will be the subject of a separate publication.However, it is important to note that the qualitative insights that we obtain in this work are not sensitive to the particular choices that we make here.

III. MOVIES
The dynamical nature of the vapor tubes described in the main text can be observed in the movie shown in Figure S3.There, we show a top and a side view of a plate-spanning vapor tube from a simulation with a biasing potential, κ( Ñ − Ñ * ) 2 /2; κ = 0.12 kJ/mol and Ñ * = 650.The average value of Ñ in the presence of the biasing potential was Ñ = 653.The plate atoms are shown as cyan spheres, water molecules are not shown for clarity, and the purple mesh corresponds to the instantaneous interface enveloping the vapor region.In essence, the purple mesh defines regions devoid of water molecules, such that the space outside the mesh is high in water density.
We show similar movies for an isolated cavity at the surface(s) of the hydrophobic plate(s) in Figure S4.The coloring scheme is the same as in Figure S3, and the configurations are taken from a biased simulation with Ñ * = 660 and κ = 0.12 kJ/mol; the corresponding Ñ = 668.There, an isolated cavity can be observed fluctuating in shape and size, and even moving from one plate to the other.The mechanism for bubble migration from one plate to another is an open question relevant to dewetting transitions and deserves further investigation.

IV. DEFINITION OF THE TUBE AND CAVITY INDICATOR FUNCTIONS
The tube indicator function h tube is determined by examining the coarse-grained density field, ρtotal , between the two plates.We define the x-coordinate to be perpendicular to the plates, such that the two plates are located at x 1 and x 2 , respectively, with x 1 < x 2 .We additionally define a buffer region b = 0.4 nm from the center of each plate to avoid the region where the coarse-grained density originating from the plate atoms is larger than or close to 0.5.If at some location, (y * , z * ), the total coarse-grained density is below 0.5 at all x-values between x 1 + b and x 2 − b, we assign h tube = 1, i.e.
The cavity indicator function h cav is similarly determined from ρtotal .However, h cav is equal to unity if the coarse-grained density is less than a half at some location between the plates, but not for all x-values between the plates.In other words, : {ρ total (x * , y * , z * ) < 0.5} , h tube = 0 = 0 otherwise. (S7)

V. PROCEDURE FOR CALCULATING UNBIASED ENSEMBLE AVERAGES
All umbrella sampling simulations performed in this work employ a biasing potential.Therefore, when performing ensemble averages of any observable, care must be exercised in accounting for the bias introduced by the potential.This is done using the WHAM/MBAR formalism [3][4][5], such that the unbiased ensemble average of any observable, A(R), which can be expressed as a function of the configuration vector, R, is given by where Here, F k is the free energy of the kth biased ensemble, K is the number of biasing potentials or windows used, N k is the number of samples in window k, V k is the biasing potential for window k, and R j,n refers to configuration n in window j.Equation S8 is used here to calculate the ensemble average of h tube , conditioned on the number of waters in confinement being N , according to where δ N,N (Rj,n) is the Kronecker delta function.The ensemble average of h cav as a function of N is calculated in an analogous manner by replacing h tube (R) with h cav (R) in Equation S10.

VI. AVERAGE SHAPE OF THE VAPOR TUBE
The average shape of a vapor tube is identified by first taking the average of the coarse grained density, ρtot (x, y, z) ≡ ρtotal (x, y, z; R) , over 5000 frames.Fig-

ure S5 displays the averaged density of the d = 20
• A with an average of 569 waters (top) and 474 waters (bottom) between the plates, obtained from biased simulations with Ñ * = 570 and κ = 0.03 kJ/mol and Ñ * = 480 and κ = 0.03 kJ/mol, respectively.From left to right in each row presents a front view, side view and threedimensional view of the averaged shape of the vapor tube.The three-dimensional rendering of the vapor tube is obtained as the iso-density surface ρtot = 0.5, and this surface provides an accurate description of the average shape of the vapor tube.Additionally, the radii obtained directly from these isodensity surfaces, 1.03 nm and 1.30 nm, respectively for the top and bottom panels of Figure S5, are in reasonably good agreement with those obtained from the simple relation between r and N used in the main text (Equation S13), 1.22 nm and 1.45 nm, respectively.Note that the vapor tubes are cylindrical to a good approximation, in accord with macroscopic theory.However, further corrections to such theories could be obtained by taking into account finer details of the "hour-glass" shape of the vapor tubes.Because the precise shape of the vapor tubes depends on the details of the instantaneous interface calculation and the parameters employed, we do not attempt such corrections here.

VII. CONTACT ANGLE DETERMINATION
The contact angle of the surface was determined by performing a 2 ns simulation of a cylindrical droplet containing 4142 water molecules.This cylindrical droplet was divided into five slabs, each 1 nm in width, and the density as a function of the radius y and the height z of each slab was computed.The droplet profile is then defined as the point where the density of the droplet is equal to half that of the bulk density.
The contact angle was then determined by fitting this droplet profile to a circle for z > 0.7 nm, as shown in Figure S6, where the fit function is given by Using this functional form, the contact angle, θ, can be obtained from .
(S12) The contact angle was determined independently for each of the five slabs and averaged to yield a contact angle of θ = 120.2• ± 0.5 • , or cos θ = −0.503± 0.008.

VIII. DETAILS OF VAPOR BASIN FITS
A. Fitting used in the main text In order to fit the simulated free energies to macroscopic theory, we consider the formation of a vapor tube of radius r, which depends on the number of waters, N , between the hydrophobic plates as Note that Equation S13 is an approximate expression; however, this simple approximation captures the salient features of vapor tube formation, as detailed below and in the main text.The free energy as a function of r as predicted by macroscopic theory is (S14) where γ is the liquid-vapor surface tension, θ is the contact angle (determined as described above), and λ is the line tension.The first term in Equation S14 is the free energy of vapor tube formation, ∆G th , described in the main text, and the last term accounts for the translational entropy of the vapor tube.The effective distance between the plates, d eff , is obtained by subtracting a constant offset, ξ, from d, that is d eff = d−ξ.The x-intercept of a linear fit of the simulated d-dependence of N liq is equal to ξ, as shown in Figure S7, so that a plot of N liq vs d eff passes through the origin.Here we find a value of the offset to be ξ = 0.4964 nm.The fits shown in Figure 4a of the main text were obtained by fitting the simulated free energies for N < N kink and r < L/2 to the parameters, βγ and λ/γ, using Equation S14; the fit parameters obtained are listed in Table I.

B. Curvature Corrections
We also explored a number of other possible fits, the first of which includes curvature corrections to the surface tension using the Tolman length, δ, such that the curvature corrected surface tension is γ(r) = γ(1 − δ/r).The fit equation then becomes However, by referencing the simulated free energies for a given d-value to the free energy in the liquid basin, that is ∆G(N liq ) = 0, we force ∆G(r) to be 0 at r = 0. To be consistent with this convention, we neglect the constant From Equation S16, we see that the curvature correction δ acts in a manner analogous to the line tension λ, i.e. both are coefficients in the term linear in r.Therefore, we can simply relate the fit parameters in Equation S16 to those in Equation S14, according to where the primed quantities indicate those obtained from Equation S14.The value of λ/γ obtained from Equation S16 can then be readily predicted with knowledge of the Tolman length.Recent estimates of this length for SPC/E water at 300 K yield δ ≈ −0.1 nm [6,7].Therefore, the value of λ/γ changes by ten to sixteen percent through the inclusion of the Tolman length in our fitting procedure.

C. Omission of Line Tension
To ascertain the importance of line tension, we attempt to fit the simulated free energies without including the term containing λ in our expression for the macroscopic theory, and instead using δ as a fit parameter according to As one may anticipate from the above discussion, the data is fit equally well by Equation S18.However, this procedure yields unphysical estimates of the Tolman length, with δ ∼ −1 nm.Therefore, we can conclude that line tension is necessary to provide a physically accurate description of the vapor tube formation and growth that facilitates capillary evaporation between nanoscopic hydrophobic plates.

IX. LIQUID BASIN FITS
In order to characterize the free energetics of isolated cavities, we fit the simulated free energies in the range, N kink + 20 < N < N liq − 50, to straight lines.We first fit the d = 25 • A free energy, and use the corresponding slope for all d-values.Intercepts are then obtained by fitting the data for each d-value separately.For d ≤ 16 • A , data in the N kink + 20 < N < N liq − 50-range are insufficient to be fit reliably, so the fitting was not performed.The regions fit for each d-value are shown in Figure S8.In order to estimate errors, we divided the data into five blocks of equal length, and fit each block separately (similar block averaging analysis was performed for all quantities).Therefore, for each d, we show five data sets and five fits, although they are nearly indistinguishable.
From these fits, we can examine the behavior of N liq (d) − N int (d), where N int (d) is the x-intercept of the linear fit to the fat tails (for d > 16 • A ).Both N liq (d) and N int (d) are roughly linear in d, see Figures S7 and S8, respectively.However, as shown in Figure S9, the quantity N liq − N int , which represents the distance in N from the liquid basin one must move to observe non-Gaussian tails in ∆G(N ; d), or equivalently, to form isolated cavities, is largely insensitive to d.

X. BARRIER LOCATION AND HEIGHT: SIMULATIONS VS MACROSCOPIC THEORY
In this section, we focus on the location as well as the height of the barrier, as predicted by macroscopic theory, and compare them with the corresponding simulated values.Macroscopic theory predicts a linear dependence of the location of the barrier, r * , on d, where we have neglected the logarithmic term in Equation S14.This behavior is captured by the data shown in Figure S10, where the r * -values are obtained from the fits of the vapor tubes free energetics to macroscopic theory.The slope of the fitted line corresponds to cos θ = −0.547,which is slightly larger in magnitude than that obtained from direct simulations.The ratio λ/γ = −0.16nm obtained from the y-intercept of the linear fit lies in the same range as those predicted from fitting the free energies to Equation S14.Macroscopic theory also predicts ∆G(r * ; d) to be quadratic in d; this behavior is similarly captured as illustrated by fitting the data to a parabola (solid line).
In addition to the position and height of the maximum of the fits, we also include the position and the height of the free energy barrier obtained from the simulated free energy profiles, r max (d) and ∆G(r max ; d), respectively.For large d-values, we expect r max ≈ r * , however, this is not true at smaller plate separations because the critical vapor tube is in the metastable branch of the vapor tube free energy; for these d-values, the barrier corresponds to the kink in the free energy.Indeed, r * = r max at these separations, and the dependence of r max on d can not be described by a straight line over the entire range of d.Similar behavior is observed in ∆G(r max ; d), albeit to a lesser extent.

Figure 1 .
Figure 1.(a -c) Simulation snapshots of water (shown in red/white) in confinement 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 (Equation1), suggesting that a vapor tube larger than a critical size must be nucleated before dewetting can proceed.

Figure 2 .
Figure 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 N kink .

Figure 3 .
Figure 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 = N kink − 12, and (b) an isolated cavity for N = N kink + 3. Water molecules not shown for clarity.(c) Average of the binary vapor tube indicator function, h tube , conditioned on the number of waters in confinement being N , displays a sharp transition from 1 to 0as N is increased.The color scheme is the same as that in Figure2.The value of N corresponding to h tube N = 0.5 (dashed line) is defined as N tube .(d) N tube is identical to the location of the kink in the free energy profiles, N kink , confirming that the kink demarcates conformations with and without vapor tubes.(e) Conditional average of the isolated cavity indicator function, hcav N , shows a sharp increase in the vicinity of N kink (the square symbols correspond to N tube ), followed by a gradual decrease at larger N -values, and eventually vanishing around N = N liq .(f) For the larger d-values, h tube N tube = hcav N tube = 0.5.However, for the smaller dvalues, hcav N tube < 0.5, suggesting the possibility of direct vapor tube nucleation without isolated cavities as intermediates.

Figure 4 .
Figure 4. (a) β∆G(N ≤ N kink ; 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 employing the relation πr 2 /L 2 = (N liq − N )/N liq in the region r kink < r < L/2, and the lines are fits to macroscopic theory.(b) The portion of the free energy corresponding to the liquid basin (N ≥ N kink ) 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 hcav N is 0.5 with a negative slope.The linear regions have roughly the same slope for all separations.

Figure 5 .
Figure 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 Figure4.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 sub-critical, and has to grow larger for dewetting to proceed.(d) Comparison of the location of the kink, N kink , 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 ≈ N kink .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 N < N kink .The transition from non-classical to classical behavior occurs in the vicinity of the coexistence separation, dc.

Figure S3 .
Figure S3.Movies (click to view) illustrating a platespanning vapor tube in a biased simulation with Ñ * = 650 and κ = 0.12 kJ/mol, with Ñ = 653.Plate atoms are shown as spheres, while the purple mesh corresponds to the instantaneous interface enveloping the vapor region (water molecules have been omitted for clarity).The movie is shown from directions (top) perpendicular and (bottom) parallel to the plane of the plates.The duration of the movie corresponds to 250 ps of simulation time.Adobe Reader may be needed to view.Representative snapshots are shown here; movie files are available upon request and will appear in the final published version.

Figure S4 .
Figure S4.Movies (click to view) illustrating vapor bubbles formed at the surface of a hydrophobic plate in a biased simulation with Ñ * = 660 and κ = 0.12 kJ/mol, with Ñ = 668.Plate atoms are shown as spheres, while the purple mesh corresponds to the instantaneous interface enveloping the vapor region (water molecules have been omitted for clarity).The movie is shown from directions (top) perpendicular and (bottom) parallel to the plane of the plates.The duration of the movie corresponds to 250 ps of simulation time.Adobe Reader may be needed to view.Representative snapshots are shown here; movie files are available upon request and will appear in the final published version.

Figure S5 .Figure S6 .Figure S7 .Figure S8 .Figure S9 .
Figure S5.The average coarse-grained density, ρtot(x, y, z), for d = 20 Figure S10.(top) Barrier position r * and (bottom) barrier height β∆G(r * ) of the fitted vapor tube free energies, are wellfit by linear and parabolic functions of the inter-plate separation, d, respectively, in agreement with macroscopic theory.In contrast, the simulated barrier position rmax and barrier height β∆G(rmax) deviate from those functional forms.

Figure S14 .
Figure S14.Free energies in the vapor basin (N < N kink ) and the corresponding fits to Equation S14.Simulation data is shown as points and the fits are shown as solid lines.The arrow indicates the direction of increasing d, from d = 11 • A to d = 25 • A in increments of 1

Figure S15 .
Figure S15.Free energies in the liquid basin (N > N kink ) and the corresponding parabolic (solid) and linear (dashed) fits to the right side of the minimum and to the fat tail regions, respectively.

TABLE I .
Parameters obtained from fitting the vapor tube portion of the free energy to Equation S14