Classical nucleation of vapor between hydrophobic plates

Nucleation of bubbles in nanoconfined liquids has a major importance in diverse fields such as physics, chemistry, engineering, and biology: the proximity of two hydrophobic surfaces promotes the evaporation of the intervening liquid,¹ which, for instance, can be exploited in self-recovery superhydrophobic surfaces^2–6 and is crucial to achieve extrusion of liquids entrapped within nanoporous materials^7–15 and in nanofluidic circuits.¹⁶ In biology, the formation of localized vapor pockets gives rise to attractive forces between hydrophobic protein residues,^17–19 playing a major role in protein folding^20,21 and gating of ion channels.^22–26 

In this work, we devise an exact classical nucleation theory (CNT) for the formation of vapor bubbles confined between two infinite parallel plates, showing that different bubble morphologies are possible and dominate at different pressures. This theory is then used to interpret the existing simulation results of water cavitation between plates at nanometer separation, discriminating purely “classical,” i.e., that can be explained in terms of CNT, from “non-classical” effects. Unexpectedly, we found quantitative agreement of CNT down to plate separations of 2 nm, suggesting that a 150-year-old theory is still fresh and insightful and, in perspective, can be used to model and understand more complex phenomena such as nucleation within flexible boundaries such as the hydrophobic regions of membrane proteins.

The objective of the present section is to offer a brief overview of the heterogeneous CNT that we make use of and that we are going to identify as “confinement CNT.” In the classical capillarity framework, the system is supposed to be composed of three phases—solid, liquid, and vapor. It is assumed that these phases are divided by geometrical surfaces (which we denote, using subscripts, as Σ_lv, Σ_sl, Σ_sv) and that phases display bulk properties up to the interface (sharp interface approximation). The convenient free energy (i.e., the grand canonical potential²⁷) of the system can now be expressed as a sum of the bulk and interface terms of the various phases. Since the nanoscale confinement line tension has been shown to be relevant,^9,11 we also include free-energy terms scaling with the length of the triple line (i.e., the line where the three phases meet).

In our treatment, we are going to consider a rigid confinement such that the solid surface (Σ_s = Σ_sl ∪ Σ_sv) does not change during nucleation. Due to this assumption, the bulk term in the free energy associated with the bulk solid phase is constant and can be set to zero without the loss of generality. For the same reason, the sum of the areas of the solid–vapor and solid–liquid interfaces is constant A_sv + A_sl = Area(Σ_sl ∪ Σ_sv) = const.

In the confined case, Laplace’s equation has to be supplemented with a boundary condition prescribing the angle at which liquid–vapor surfaces meet with the solid walls. We will return to this point in Subsection II C when discussing the general case in the presence of a line tension.

In this article, we focus on axisymmetric vapor bubbles,^34,35 confined between planar surfaces.³⁶ 

The geometry that will be considered is one composed of two infinite parallel plates. This geometry is not exclusive to vapor nucleation studies, but it is also common to investigations of the (thermodynamic) Casimir forces mediated by thin confined films close to the critical point (see, e.g. Ref. 37). Within the present CNT framework, in order to comply with force balance at the boundary, CMC bubble surfaces have to satisfy the boundary conditions dictating the specific shape of the various bubbles. In the absence of a line tension, this boundary condition is represented by the familiar Young equation prescribing the contact angle at which vapor bubbles approach the solid surface. In addition, this latter equation can be obtained variationally.³¹ 

When considering the presence of a line tension, it can be shown, using variational arguments that are not dissimilar from those previously discussed, how the addition of a triple line tension results in an additional term in the Young equation,³⁸ which is proportional to the geodesic curvature of the triple line, when thought as a curve on the surface of the solid boundary with a proportionality factor that amounts to the ratio τ/γ_lv.

From a geometrical perspective, the addition of a line tension to the model is not directly able to change the general CMC shape of droplets, whose shape, within CNT, is prescribed by Eq. (4). Yet, the presence of line tension results in a complicated boundary condition prescribing a size-dependent contact angle [see Fig. 5(d)], which is, in general, different from the Young one. It should be noticed that, by introducing a line tension, one also introduces an intrinsic lengthscale, defined by the ratio τ/γ_lv, disrupting the general property of capillary shapes to be scale invariant.

Restricting our attention to axisymmetric surfaces makes it possible to considerably simplify the problem of determining the bubble shapes, making it tractable analytically or semi-analytically. This is due to the fact, proven almost two centuries ago by Delaunay,⁴⁰ that there exist only a limited number of families of axisymmetric CMC surfaces.

Such surfaces have in common the fact that they are generated by the revolution along the x-axis of curves that are rouloids of conic sections. These curves, known as the undulary, the nodary, and the catenary, represent the locus of the points spanned by the foci of different conic sections when rolling without slipping along the x axis. More specifically, the undulary is the locus of points occupied by the focus of an ellipse while rolling along the x axis, the catenary is similarly generated by the focus of a parabola while the nodary is generated by the focus of a hyperbola. The rouloid curves, together with their surfaces of revolution, are shown in Fig. 1. The more familiar CMC surfaces represented by the cylinder and the sphere can also be rationalized in the same framework by considering these shapes as the result of the revolution, respectively, of a straight line of a semi-circle (i.e., the rouloids for a circle and a degenerate ellipse). Adding to the surfaces mentioned previously, the plane is finally the last member of the Delaunay family of CMC surfaces of revolution.

The six families: cylinders, spheres, catenoids, unduloids, nodoids, and planes correspond to the totality of axisymmetric CMC surfaces. We refer to Refs. 41 and 68 for an enlightening review of this subject.

Although the classification of axisymmetric CMC surfaces dates back to Delaunay, it was not until 50 years ago that Kenmotsu derived convenient closed-form expressions for such surfaces in terms of elliptic integrals.⁴² Variations of such expressions are reported in Subsection 2 of the  Appendix.

It is now clear how the candidate solutions to our variational problem for the minimal bubble shapes in confinement, associated with the given value of the vapor volume, should be singled out by selecting suitable portions of a CMC surface Σ_lv, compatible with the prescribed curvatures, vapor volumes, and boundary conditions [i.e., the contact angle prescribed by the generalized Young law, Eq. (5)]. This procedure amounts to determining the coefficients (e.g., B, or a and c for the parametrizations given in Subsection 2 of the  Appendix) in the expressions for the constant mean curvature surfaces. This procedure can be performed analytically for the simplest cases (e.g., in the absence of line tension) or numerically in the other cases.

From the previous discussion, it is clear how finding a section of a mean curvature surface satisfying prescribed boundary conditions corresponds to the solution of a boundary value differential problem. As such uniqueness of the result is not guaranteed and, in general, several shapes exist satisfying the CMC and boundary conditions. In particular, unduloids and nodoids are periodic functions and, typically, more than one section of these surfaces is compatible with a given bubble volume and contact angle (see Sec. I C). In this case, CNT suggests that the most probable bubble configuration is the least energetic one. We will, therefore, screen candidate (constrained) minima via Eq. (2).

In the following discussion, we focus on nucleation occurring between two hydrophobic parallel plates (see, for instance, the right panels of Fig. 2). Two main families of bubble geometries are relevant in such confines: a spherical-cap bubble in contact with one wall and an hourglass-shaped bubble bridging the two walls [see Figs. 5(b) and 5(c), respectively]. As the hourglass-shaped bubble grows along its nucleation path and its volume changes, the mean curvature H varies proportionally to the Lagrange multiplier λ and bubble shapes change with continuity. In this discussion, we have neglected the possibility of spherical bubbles nucleating in the bulk, and of multiple surface-attached spherical caps, which can be shown to be non-minimal in terms of ΔΩ.

Within the present framework, it is, particularly, instructive to single out zero mean curvature bubbles, which at (bulk) liquid–vapor coexistence (ΔP = 0), correspond to unconstrained extrema (i.e., maxima in this case) of the grand potential. These shapes correspond naturally to the most energetic bubbles along the nucleation path, that is, the transition state.

An analytic solution to the problem of finding an axisymmetric surface with zero curvature dates back to Meusnier,⁴³ who in the late XVIII century proved that the only minimal surfaces of revolution are the plane and the catenoid, i.e., the surface obtained through a rotation around the x axis of the graph of a catenary.

Summarizing, there are six axisymmetric morphologies of constant curvature that can describe bubble shapes compatible with CNT, as specified by Eqs. (4) and (5): the sphere/spherical cap, the plane, the cylinder, the catenoid, the unduloid, and the nodoid. In particular, the nucleation path is constituted by the succession of bubbles having the minimum grand potential at fixed V_v; these bubble shapes are CMC surfaces that are sought for among the aforementioned bubble morphologies. The values of the parameters describing the surfaces are eventually determined by the volume and contact angle of the bubble. For several cases (e.g., τ = 0 and ΔP = 0), the values of the parameters can be determined analytically in a simple way. In the other cases, the parameters can be obtained numerically by minimizing the grand potential in which all geometrical quantities are expressed through Eq. (A3) or simplified versions of it.

In this section, we exploit the CNT summarized in Sec. II C to interpret the nucleation free-energy barriers obtained by means of accurate simulations in a sequence of papers^45–47 that focused on the confined nucleation of vapor between hydrophobic plates. This procedure allowed us to rationalize the scaling of free energy barriers with the plate distance and to narrow down a range of values for the line tension τ.

Since our CNT is based on the sharp interface assumption, it is clear that some convention on the positioning of the interfaces has to be assumed when interpreting atomistic simulation data; this is especially relevant for the positioning of the plates. In the following discussion, we assume the solid interfaces to be shifted by σ = σ_O–O away from the nominal positions of the wall atoms, accounting for the excluded volume at the plates (see right panels of Fig. 2). We are going to refer to this measure for the plate distance as D_eff, whereas we indicate the nominal plate distance, which is typically addressed in the simulation work, with the symbol D_nom. The same convention was adopted also when independently measuring the contact angle (see Fig. 9).

Figure 2 directly compares the CNT results using our convention for the plate distance D and the atomistic free-energy barriers from Refs. 45–47. Thermodynamic and material properties of the atomistic system are taken from the literature or computed via independent molecular dynamics (MD) simulations: ΔP = 0, γ = 0.0636 N/m, and θ_Y = 135° (for details, see Subsection 1 of the  Appendix). We remark that there are no fitting parameters, besides the uneludible need of choosing a sensible convention to express the atomistic plate distances. With these positions, it is seen that τ = 0 pN yields a quantitative accord with the atomistic data. The bars further indicate that all simulation results seem to fall within the narrow range τ = 0 ± 5 pN, a small value that is smaller than previous reports, but which is found to be consistent with recent simulations performed by Bey et al.⁴⁸ In this recent paper, the authors performed a direct measurement (based on the Navascués and Tarazona theory⁴⁹) of the line tension of SPC/E water on similar model hydrophobic substrates obtaining values of the line tension in the same range.

Within our new estimate of τ, it is clear that the dominant contribution to the free-energy barriers, and thus to the nucleation kinetics, comes from the surface terms in Eq. (2). In other words, even at two-phase coexistence (ΔP = 0), the competition between the cost of forming a liquid–vapor interface and the free-energy gain obtained by dewetting a portion of the hydrophobic walls is sufficient to determine the thermally activated drying of the two plates even without invoking major line tension effects. Accordingly, the CNT expression for the free-energy barriers is ΔΩ^† = α(θ_Y)γD², in agreement with the scaling predicted in Refs. 50 and 51 for which we have provided the full the analytical expression, Eq. (7).

CNT in the absence of line tension offers an additional interpretation on the ΔΩ^† ∼ D² scaling: if the free energy is dominated only by surface terms, all the transition-state bubbles are the same up to a simple geometrical rescaling. When τ ≠ 0, the intrinsic lengthscale imposed by the τ/γ term in the generalized Young equation results in a more complicated scenario, where these simple scaling arguments do not apply anymore, due to the size-dependent shape of the transition-state bubbles which have different contact angles. Figure 3 clearly shows the different scalings of ΔΩ^†(D) with and without the line tension: in a log–log scale $ΔΩ†(Deff)τ=0$ is represented by a straight line, indicative of the expected power law behavior. When, instead, τ ≠ 0, deviations are observed which become more prominent at small plate distances (D_nom ∼ 2 nm). Interestingly, the largest deviations are observed for negative line tensions, which account for effectively more hydrophobic plate behavior.

As in the case of Fig. 2, the analysis of nucleation rates underscores how smaller to vanishing line tension values τ = 0 ± 5 pN yield good agreement between the atomistic data and CNT.

It is somewhat surprising that the simple and generic sharp-interface CNT is in quantitative agreement with simulation data reported by different authors and obtained with different simulation techniques for the nucleation of water between simple hydrophobic planar walls.

At the same time, it should be noted that our result does not rule out line tension in the general case, for its value might depend on the confining geometry.⁵⁶ For example, experimental data on water in hydrophobized MCM-41 nanopores reported τ = −30 pN⁹ and a similar figure are obtained in MD estimates.¹¹ In this latter system, the cylindrical geometry of the pores imposes bubble shapes in which the triple line terms are magnified, especially close to the transition state. Altogether, our study brings to the conclusion that accurate sharp-interface models are necessary to interpret experimental and computational data to extract the elusive line tension.

Summarizing, the exact CNT presented here pushes the value of line corrections extracted from the data of water cavitation between plates down to values smaller than 5 pN in the absolute value. Importantly, this result suggests that the classical budget of the volume and surface terms in (2) yields accurate estimates of the free-energy barriers for nucleation in confinement between plates. This tool is conceptually simple and computationally inexpensive and thus shows promise for investigating wetting and drying phenomena at the nanoscale.

In order to further test the predictions of the CNT hereby presented, it is instructive to analyze the full nucleation pathway. Earlier simulation reports, e.g., Ref. 47 have established that, at thermodynamic conditions close to liquid–vapor coexistence, the typical nucleation pathway is characterized by initial density fluctuations that eventually result in the formation of a vapor bubble attached to one of the walls. This spherical-cap bubble grows until an “hourglass-like” bubble becomes more energetically favorable. The maximum of the free energy (and thus the transition state) is found in correspondence with this hourglass morphology. Subsequent growth of the bubble results in the spontaneous growth of the hourglass-like shape, until the finite size of the simulated plates induces the pinning of the triple line not considered in the present theory. Figure 5 shows that the current CNT predicts a similar path formed by the juxtaposition of two free-energy branches corresponding to the spherical-cap and hourglass morphologies which dominate at small and large bubble volumes V_v, respectively. Bubble morphologies can be expected to depart from the atomistic ones for finite plates at low V_v, because water density fluctuations are not captured by CNT, and at large V_v, because the infinite plates of our CNT do not induce pinning.

The qualitative agreement between the atomistic and CNT nucleation mechanisms further motivated us to quantitatively test the CNT predictions for the entire free-energy profile against available simulation data.⁴⁷ The detailed comparison is presented in Fig. 6. In general, the agreement of the free-energy profiles appears fair, with major deviations of CNT from atomistic data observed only for the system with the smaller plate distance D_nom = 14 Å; in this case, the effective distance between the plates D_eff < 1 nm. For the other systems (D_nom = 17, 20, and 23 Å), the central part of the free-energy profiles, including the transition state, is in quantitative accord with simulations, while deviations appear at small and large V_v for different reasons, explained below.

At small V_v, atomistic effects are observed to dominate the simulation free energies, accounting for the deviation between simulations and CNT. In particular, Gaussian density fluctuations in the confined liquid without the formation of an actual interface are expected at the molecular scale,^17,47,57 while CNT is a mean-field theory according to which the smallest bubbles possess a well-defined and sharp interface. These behaviors are known to give rise to opposite curvatures⁵⁸ in the free energy profile (see the inset in Fig. 6). Since the volume of the vapor bubble in the fluctuation-dominated regime is actually zero, the comparison in Fig. 6 is made by adding a small horizontal shift to the CNT profiles, which allows the superimposition of the free-energy branch corresponding to the spherical cap, not dissimilar from what can be observed in Ref. 47, Fig. 5. The value ΔV = 2200 ų, which may look arbitrary, can be rationalized by measuring, in independent simulations, the actual volume associated with the liquid phase intruding the plates. The small ΔV emerges as the difference between the inter-plate volume measured according to the effective plate distance convention and the actual volume occupied by water molecules. Using a convex hull to estimate by excess this latter quantity, we were able to measure, in independent simulations, a ΔV of roughly 1900 ų. It should also be noticed how a certain degree of approximation in the mapping between CNT and atomistic simulation is inherent in the fact that atomistic simulation free energies are obtained using the number of atoms N within a control volume as the order parameter, while CNT uses the bubble volume V_v. In the present comparison, a constant density was used to convert N into V_v; this approximation is good for larger volumes, yet it deteriorates for small values V_v. The agreement, with the volume shift explained previously, is excellent for the spherical cap regime and in the vicinity of the transition state, especially for the cases D = 23 and 20 Å in which a proper bubble is formed before the transition state (see below). It is also useful to notice how the slope for the spherical bubble branch resulting from CNT exactly matches the slope obtained in simulations without any fit, which reassures on the horizontal shift procedure discussed previously.

For large V_v, the atomistic and CNT free-energy profiles would tend to diverge for different reasons, related to the finite size of the plates used in simulations. As a matter of fact, in simulations, the hourglass bubble becomes pinned at the plate boundary leading to the emergence of a second minimum in the free energy associated with the metastably dry plates. For this very reason shown in Fig. 6, we opted to report the CNT curves associated with bubbles that are at least σ_O–O away from the plate boundaries, to reflect the fact that pinning effects have not been considered in our CNT model, based on the assumption of infinite plates. Such finite-size effects not only induce bubble pinning at the edges but also influence indirectly thermodynamic parameters, such as θ_Y, when the liquid–vapor interface is close to the plate boundaries.

As already observed by Remsing et al.,⁴⁷ there is a particular bubble volume at which the free energy associated with the formation of the spherical-cap and the hourglass bubbles coincide, resulting in a “kink” in the free energy profile ΔΩ(V_v), associated with the change in the different slopes of the two branches. It is known that kink-like features are a signature that the order parameters used in CNT and in simulation^11,59 are not able to locally describe a morphological change in the bubble and discontinuously switch between two minimal bubble morphologies. More in detail, ordinary CNT assumes that the nucleation pathway is a sequence of quasi-equilibrium, volume-constrained vapor bubbles; similar hypotheses are also used in simulations that use the number of water molecules N or V_v as collective variables. Since these approximations are equivalent for CNT and the reported simulations, both feature a kink in the free-energy profile. At pressure conditions far from liquid–vapor coexistence, it is possible that the maximum of the free energy profile (transition state) coincides with the kink, as shown in Fig. 5. While the detailed exploration of this regime is outside the scope of the present work, it would be desirable in the future to deploy dedicated methods to overcome the subtleties associated with a single order parameter/collective variable.^11,12,59,60

The agreement with simulation data suggests that the CNT model can be used as a flexible and computationally convenient tool to predict, at least approximately, nucleation free-energy barriers, rates, and mechanism. For instance, within the CNT model, it is straightforward to continuously vary geometrical (D) and thermodynamic parameters (ΔP, θ_Y, τ, etc.) and recompute, for each set, the nucleation path and free-energy profile. This parameter scan allows one to identify different nucleation regimes, which, in turn, determine non-trivially the nucleation kinetics via Eq. (8). Here by “regimes” we refer to the bubble configurations connected to the free-energy maximum, which change depending on confinement and thermodynamic conditions.

The CNT predictions for the nucleation regimes are reported in Fig. 7(a). Recalling Fig. 5, it is seen that decreasing ΔP shifts the maximum to lower V_v, first to the kink and then to the spherical-cap configuration. Decreasing the distance D, instead, causes the free-energy profile connected with the hourglass bubble to become steeper, while shifting the maximum toward the point where the two free-energy branches meet (Fig. 6). Intuitively, as D increases, one moves toward the single-plate limit, where only the spherical-cap regime exists. Finally, Fig. 7(b) clearly shows that line tension τ changes the height of the barriers and the critical volume but has a minor influence on the nucleation regime, slightly altering the stability of hourglass-shaped transition states, in which two contact lines exist.

A partial comparison of the nucleation regimes is possible with the mesoscale predictions of Ref. 44, remembering that the main case considered there is that of condensation. Both theories predict that, close to coexistence (ΔP ≈ 0 or unit supersaturation in Ref. 44), the hourglass regime is favored over the spherical-cap one, while, away from this condition, the opposite happens; similarly, both the macroscopic and the mesoscopic theories predict that increasing D facilitates the spherical-cap regime. A conceptual novelty introduced by the present CNT approach is that one should consider the entire nucleation path and not only the transition states: in this case, the construction in Fig. 5 clearly shows that the nucleation process must always start with a spherical cap, simply because for vanishingly small V_v, a gap-spanning vapor tube meeting the plates with the prescribed contact angle cannot exist (in geometrical terms, no hourglass can be constructed for V_v → 0). The consequence of this approach is that an intermediate regime between the spherical cap and the hourglass exists, the kink regime, which actually occupies a substantial domain in the ΔP vs D plane (Fig. 7). We would like to underline how this prediction is consistent with the results of atomistic simulations.⁴⁷ 

Summarizing, the present CNT is capable of explaining the different nucleation regimes between two parallel plates for a variety of surface chemistries and fluids and in a wide range of thermodynamic conditions. A favorable comparison is found between CNT predictions and a number of simulation data obtained by various authors with different techniques. Surprisingly, the agreement between the free-energy barriers computed by MD and by CNT holds down to very small plate separations of the order of roughly 2 nm. A significant non-classical effect due to equilibrium density fluctuations was only observed in the confined liquid close to the fully wet state. These phenomena are well outside the scope of the present model inspired by macroscopic capillarity, yet in this regard, the exact CNT solution may help interpreting deviations of the free energy from this fluctuation-dominated basin as the emergence of a fully formed spherical bubble.⁵⁹ 

In this work, we have shown the power of an exact classical nucleation theory at interpreting experiments and simulations of nanoconfined water. Using this theory as a reference, we have identified from the results of recent molecular dynamics simulation^{46,47,52,61,62} non-classical effects, that is, those which deviate from the hypotheses of classical nucleation theory, for example, due to nanoscale confinement. Within the CNT realm, line tension is found to be less important than previously assumed in the plate geometry, resulting in a reasonable collapse of simulation data for τ = 0, even for distances between plates close to 2 nm. Additionally, CNT was able to rationalize the different nucleation regimes reported in simulations⁴⁷ and mean-field calculations,⁴⁴ with the transition state assuming different shapes. Non-classical effects related to fluctuations of the confined liquid are important in the free-energy profiles only close to the fully wet state. However, these deviations from CNT do not seem to significantly affect the estimates of the free-energy barriers and of the nucleation rates, demonstrating the capability of simple CNT to predict the energetics of nucleation in extreme hydrophobic confinement. In more detail, the computational strategy suggested by the present results is to use simple microscopic calculations (e.g., molecular dynamics) only to measure the relevant thermodynamic parameters (ΔP, γ, and θ_Y, as explained in Subsection 1 of the  Appendix) and computationally inexpensive CNT calculations to treat nucleation phenomena, in the spirit of morphometric thermodynamics.⁶³ This powerful and flexible approach may prove an invaluable tool to investigate phenomena of interest in nanotechnology, surface science, and biology far beyond what can be accessed today by molecular dynamics alone. Future investigations could be devoted to the investigation of the role of thermal fluctuations for the confined menisci of present interest via minimization of a non-local energy functional.⁶⁴ 

In conclusion, we observed that in the highly sophisticated simulation data available, only relatively minor non-classical effects can be detected as to vapor nucleation between extended planar surfaces even at nanoscale distance.

The authors acknowledge thoughtful discussion with Matteo Amabili. This project received funding from the European Research Council (ERC) under the European Union’s Seventh Framework Programme (Grant No. FP7/2007-2013)/ERC Grant Agreement No. 339446 and Horizon 2020 research and innovation program (Grant Agreement No. 803213).

S.M. acknowledges financial support from the Sapienza University of Rome for financial support via Grant No. RG11715C81D4F43C and from the University of Ferrara under the program Finanziamento per l’Incentivo alla Ricerca (FIR) 2020.

The authors have no conflicts to disclose.

Antonio Tinti: Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Alberto Giacomello: Methodology (equal); Supervision (equal); Writing – review & editing (equal). Simone Meloni: Supervision(equal); Writing – review & editing (equal). Carlo Massimo Casciola: Supervision (equal); Writing – review & editing (equal).

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

In order to identify non-classical effects, the thermodynamic parameters appearing in the CNT should be accurately matched to the atomistic ones. For the liquid–vapor surface tension, we use the result γ = 0.0636 N/m from Ref. 65; the pressure, instead, is the coexistence one ΔP = 0 for the considered simulations.^47,52 In the CNT augmented with line tension, τ enters only through the boundary condition, modifying the contact angle θ as compared to the Young law, Eq. (5). Therefore, if one wants to compute τ by comparing CNT with simulations, as in the present case, it is crucial to obtain an accurate estimate of θ_Y avoiding small-size artifacts. With this aim, we performed 6ns-long MD simulations (∼250 000 Atoms) of a large cylindrical drop with a radius of ∼100 Å (see Fig. 8), with the same inter-atomic potentials of Refs. 47 and 52 obtaining θ_Y ≈ 135°, a value remarkably close to the one in Ref. 45. The cylindrical droplet allowed us to avoid line tension contributions to the contact angle (see for instance Ref. 66). A bulk-like layer of water is also present below the wall in order to mimic the actual conditions presented in the literature.

The contact angle was then obtained by calculating the slope of the circle fitting the density isoline ρ/ρ_bulk = 0.5 at a distance of 1 σ from the center of the solid atoms (Fig. 9). The distance D between the plates was estimated by subtracting from the distance between the centers of the wall atoms twice the depleted layer of ∼1  σ_O–O. It is important that the same convention is used for the measurement of the contact angle, for the distance D, and for all geometrical quantities appearing in the grand potential.

The parameters appearing in the aforementioned equations, μ = 2/(a + c) = 2H, k² = (c² − a²)/c², m = (c² − a²)/2, n = (c² + a²)/2, and finally $BH≠0=1−2aa+c$⁠, solely depend on the real parameters a and c, which, in turn, control the shape of the surface; the inequality c > a is assumed for non-degenerate unduloids.⁶⁷ More specifically, from the functional form of y(u), it is possible to quickly establish that the distance from the symmetry axis is bounded between a and c. The surfaces of revolution are found to be periodic along the symmetry axis, with a distance π(a + c)/4 occurring between the points of minimum distance from the axis. The mean curvature of the surface obtained via revolution of the undulary about the x axis is now conveniently expressed by the constant H = 1/(a + c).

In conclusion, by selecting appropriate values of the parameters a and c, it is possible to obtain curves whose revolution generates several families of CMC surfaces, including spheres, cylinders, unduloids, and nodoids. Nodoids are, in general, self-intersecting surfaces, and special care must be taken when branches of nodoids are meant to represent liquid–vapor dividing surfaces.

Due to the nature of the energy minimization process, which can be recast as a differential boundary value problem, the uniqueness of the solution is not guaranteed. This reasoning can be further specialized by considering that unduloids and nodoids are, in general, periodic shapes and as such a potentially infinite number of unduloid or nodoid branches can be found, with the only constraint that they comply with the generalized Young boundary conditions. Without recurring to sophisticated stability arguments, we resolved this difficulty by noticing that the free energy Ω of these shapes, which still are extremants for our variational problem, is larger than the one associated with the uncorrugated shapes. These latter shapes can be numerically extracted by iteratively perturbing the analytic catenoid solution at zero mean curvature to yield the full nucleation path. An example of the multiple shapes with the same mean curvature that can be computed is provided in Fig. 10. It should also be remarked how, in this work, out of simplicity, only solutions that are symmetric with respect to the plane equidistant from the plates have been investigated. This symmetry is implicit in the case of identical plates and no line tension, and it was retained for small perturbations of the τ = 0 solutions that are presently reported.