Vapor nucleation paths in lyophobic nanopores

In recent years, technologies revolving around the use of lyophobic nanopores gained considerable attention in both fundamental and applied research. Owing to the enormous internal surface area, heterogeneous lyophobic systems (HLS), constituted by a nanoporous lyophobic material and a non-wetting liquid, are promising candidates for the efficient storage or dissipation of mechanical energy. These diverse applications both rely on the forced intrusion and extrusion of the non-wetting liquid inside the pores; the behavior of HLS for storage or dissipation depends on the hysteresis between these two processes, which, in turn, are determined by the microscopic details of the system. It is easy to understand that molecular simulations provide an unmatched tool for understanding phenomena at these scales. In this contribution we use advanced atomistic simulation techniques in order to study the nucleation of vapor bubbles inside lyophobic mesopores. The use of the string method in collective variables allows us to overcome the computational challenges associated with the activated nature of the phenomenon, rendering a detailed picture of nucleation in confinement. In particular, this rare event method efficiently searches for the most probable nucleation path(s) in otherwise intractable, high-dimensional free-energy landscapes. Results reveal the existence of several independent nucleation paths associated with different free-energy barriers. In particular, there is a family of asymmetric transition paths, in which a bubble forms at one of the walls; the other family involves the formation of axisymmetric bubbles with an annulus shape. The computed free-energy profiles reveal that the asymmetric path is significantly more probable than the symmetric one, while the exact position where the asymmetric bubble forms is less relevant for the free energetics of the process. A comparison of the atomistic results with continuum models is also presented, showing how, for simple liquids in mesoporous materials of characteristic size of ca. 4nm, the nanoscale effects reported for smaller pores have a minor role. The atomistic estimates for the nucleation free-energy barrier are in qualitative accord with those that can be obtained using a macroscopic, capillary-based nucleation theory.


Introduction
Nanoconfined cavitation is a topic of vivid interest for a number of reasons, mostly due to the promising technological applications associated with the subject [1] and to the exotic behavior of liquids in these extreme conditions [2]. For instance, understanding how heterogeneous nucleation takes place in nanoporous materials is of paramount importance for the design and characterization of mechanical energy storage or dissipation devices exploiting the intrusion of a non-wetting liquid in a lyophobic nanoporous material -the so-called heterogeneous lyophobic systems (HLS) [1,3]. This kind of applications is the object of a Contribution to the Topical Issue "Advances in Computational Methods for Soft Matter Systems" edited by Lorenzo Rovigatti, Flavio Romano, John Russo. a e-mail: atinti@is.mpg.de b e-mail: alberto.giacomello@uniroma1.it number of experimental studies, which showcased significant deviations from the behavior predicted using macroscopic arguments inspired to the classical nucleation theory (CNT) [2,4,5].
Molecular dynamics simulations are the tool of choice for the study of nanoconfined wetting and nucleation phenomenona as they are able to resolve the nanometric lengthscales characterizing the phenomena, allowing, for instance, for a microscopically-informed investigation of the wetting [6,7] or dewetting [8] of a single nanopore. The main computational challenge for these simulations is sampling the rare drying and filling processes of a hydrophobic nanochannel, which are thermally activated in nature. In other words, these phenomena have two distinct timescales: a fast one connected with the particle dynamics and a slow one characterizing the thermal fluctuations required to trigger cavitation. The long waiting times for a thermally activated (or rare) event can exceed by far the capabilities of the fastest available supercomputer. Therefore, advanced rare events sampling techniques are often needed in order to be able to access the wetting and dewetting processes inside a hydrophobic nanopore.
Among a wealth of rare events techniques [9], we adopt here the so-called string method in collective variables [10], which allows one to compute the most probable path associated with a rare event and the related free energy. This method is computationally efficient at sampling highdimensional free-energy landscapes, allowing considerable freedom in the choice of collective variables (CVs), i.e., the macroscopic descriptors of the process. The cost of usual rare event methods, such as umbrella sampling, instead, scales exponentially with the number of CVs [11] precluding their use in high-dimensional landscapes. For the particular case of wetting and nucleation phenomena in confinement, it has been shown that simulations with a single CV, e.g., the volume of the nucleating bubble, may be affected by significant artifacts [12,13]. In this context, the string method allows us to use the coarse-grained density field as CV; the string calculations with this CV yield the transition path(s) avoiding limitations connected with low-dimensional CVs [8,14].
In the present contribution we address the question whether one or more transition paths exist leading to the dewetting of a nanoscale cylindrical pore. We are able to identify via the string method at least three paths, with two distinct symmetries. Each path is then characterized in terms of its free-energy profile which, in turn, determines its probability and kinetics. A detailed comparison of these results with the confined water vapor nucleation computations of [8] and with macroscopic theories is also presented. This analysis underscores a clear distinction between water confined in pores with diameter of 2.6 nm and simple liquids in pores with diameter of ca. 12 σ (ca. 4 nm): nanoscale effects are significant in the first case while a macroscopic description is accurate in the latter.
The paper is organized as follows. In sect. 2 the simulation technique is briefly reviewed. Section 3 is devoted to the results of the string calculations, which are then discussed in sect. 4 and compared with macroscopic continuum models in sect. 5. Section 6 is left for conclusions.

Molecular dynamics simulations
In this contribution, we focus on the different paths for nucleation in lyophobic nanopores. Differently from previous works [8], here we employ Lennard-Jones (LJ) fluid and solid, with lyophobic interactions. Also the dimensions of the pore are larger than in the previous case.
The system under investigation consists of a single cylindrical pore of finite length excavated from a crystalline solid and surrounded by two relatively large liquid reservoirs. The pore has a diameter of 12.5 σ in LJ internal units, which would roughly correspond to a 4 nm pore in dimensional units (assuming the value of σ of the TIP4P/2005 model [15]). A snapshot of the simulated system is shown in fig. 1. The simulation comprises a total of  The main technical difference between these simulations and ref. [8] resides in the barostatting, as in the present LJ simulation we employ a Noosé-Hoover chain thermostat and barostat generating a dynamics compatible with the NPT ensemble [16] at a reduced temperature T = 0.8 /k B and at coexistence pressure (P ≈ 0). Given the small dimensions of the confined bubble as compared to the simulated system, this difference is not expected to have a significant impact on the results [17].
The liquid-liquid interactions are modeled as standard LJ interactions, with σ = = 1, while the liquid-solid interactions are modeled through the use of a modified LJ potential: where the parameter c tunes the attractive part of the potential in order to render hydrophilic or hydrophobic fluidliquid systems. The Young contact angle θ Y was measured by independent MD simulations of a cylindrical drop on a planar surface returning the value θ Y ≈ 113 • for the chosen value c = 0.6. The calculations were carried out using the MD engine LAMMPS [18] which was customized in order to implement the calculation of the collective variables. Similarly to the simulations discussed in [8], the CVs are related to the discretized density field inside the nanopore. More in detail, the CVs represent the number of fluid particles occupying each of the 160 parallelepipedic cells arranged in a 4 × 4 × 10 Cartesian grid enclosing the cylindrical pore.
Here the fundamental steps of the string method in CVs are briefly reviewed; for further details on the method we refer to the original papers and reviews on the subject [9,10,19,20]. The string method performs a local search for transition paths [19]; this implies that it is possible to find several paths for the same rare event which are characterized by a different probability. More in detail, the method is initialized with a tentative string, i.e., with a sequence of replicas (images) of the system spanning all the advancement of the rare event. In the present case, the images are replicas of the system shown in fig. 1 at different filling levels of the pore, from a liquid-filled nanopore to an empty one. The string is then evolved according to a first-order pseudo-dynamics along the perpendicular direction driven by the thermodynamic forces conjugate to the chosen CV. For each image and for each evolution step, the thermodynamic force is computed by an independent restrained molecular dynamics simulation [10,21]. At convergence, the pseudo-dynamics yields the local transition path. This scheme entails that, by choosing different initial conditions for the pseudo-dynamics, one can arrive at distinct nucleation paths. In addition, by integrating the thermodynamic forces, it is possible to reconstruct the free-energy profile connected with each path and to determine the most probable one by comparison.
Three complete string simulations were carried out using this procedure starting from different initial conditions. Each string consisted of 64 images that were evolved, following the optimization algorithm implied by the method [10,19], until convergence. This procedure typically requires 20 steps of evolution.

Results
Three different strategies were followed in order to initialize the string, as described below. a) A first string was initialized by generating a tentative path characterized by the formation of a cylindrical vapor bubble (annulus) spanning the whole middle section of the pore and growing symmetrically along the axis of the pore. We will refer to this simulation as the one initialized from a symmetric configuration. b) A second string calculation was initialized by extracting 64 configurations from an unbiased MD trajectory of a spontaneous dewetting event. This process, which at the thermodynamic conditions of the string simulations would be rare, was triggered by setting the barostat to negative pressures, so as to favor the observation bubble nucleation within a computationally accessible time. We refer to this calculation as the result of initializing the string using a spontaneous event. c) The last initial string was obtained similarly to the previous one, but with the difference that the dewetting event was triggered by adding a fictitious potential inducing the dewetting on one side of the middle section of the nanopore, in order to force the bubble nucleation from the geometrical center of the pore. These three strings described above were used in order to initialize independent evolutions of the pseudodynamics and converged to the three paths shown in fig. 2. The free-energy profiles associated with the three calculations are presented in fig. 3.

Discussion
A common feature of the three paths is that the maximum of the free energy (transition state) is connected to the breakup of the initial bubble, which gives birth to two symmetric menisci. These menisci are very similar in the three cases and have the typical shape prescribed by Young and Laplace laws: spherical caps meeting the pore walls with the prescribed contact angle. At the end of the nucleation process the menisci pin at the pore mouth, with the final state (free-energy minimum) characterized by a flat interface. The initial part of nucleation in which a bubble is formed, which determines the kinetics of cavitation, is instead significantly different for the three paths as described below. This difference is also reflected in the free energies ( fig. 3).  [22]). The profile connected to the formation of an annular bubble (green) is associated to a much larger nucleation barrier as compared to the other two nucleation paths involving the formation of an asymmetric bubble on one side of the pore (red and blue). Error bars represent the variance of the sample of FE profiles computed as explained in [23]: for clarity they are displayed only for the string initialized using a spontaneous nucleation event; similar values are expected for the other two profiles. We remark that, around the transition state, the errors are very small, comparable to the size of the symbols.
The strings initialized from spontaneous-like events (b) and c)) converged in few iterations of the pseudo-dynamics (∼ 20) to the nucleation path characterized by an asymmetric bubble with a shape resembling that of a potato chip similar to the one observed in the water simulations of ref. [8]. Nucleation originating from a non-axisymmetric bubble was also predicted by CNT arguments in [4], but there the CNT path was not continuous as in the present case. The main difference between path b) and c) resides in the position of the initial bubble: in the spontaneous case, the nucleation bubble forms relatively close to one of the two ends of the nanopore, while, in the other case, the bubble lies at the center due to the ad hoc triggering potential.
The two asymmetric nucleation paths result in analogous free-energy barriers of about 35 k B T ( fig. 3). This similarity may be ascribed to the azimuthal and translational symmetry of the non-axisymmetric bubble along the axis of the cylindrical pore which results in similar free energetics sufficiently far from the pore mouth. However, minor consequences of the different initial positions of the vapor bubble are still observable in the free-energy profiles: the transition state for the spontaneous bubble is anticipated, probably because the critical bubble is de-formed by pinning to the cavity mouth at one end of the pore. The free-energy profile in the initial phase of bubble formation seems also reminiscent of how the bubble is generated: more rounded if by imposing a small repulsive potential, steeper from a spontaneous fluctuation.
Differently from the previous cases, the string initialized using a cylindrical vapor annulus converged, after about 40 iterations, to a nucleation path that does not involve the loss of the axial symmetry. The nucleation bubble was formed in the shape of a cylindrical annulus growing larger and larger until the symmetric capillary neck breaks up giving rise to two independent menisci, as in the other two cases. The free-energy barrier connected with this path is much larger: The MD simulations presented in sect. 3 provide atomistic evidence of the existence of several nucleation paths inside a hydrophobic cylindrical pore. These paths can be classified in two different families: one in which nucleation starts from an axisymmetric bubble and a second from a non-axisymmetric one. The free-energy profiles of fig. 3 clearly show that the formation of an asymmetric bubble on one side of the pore is energetically more convenient as compared to the formation of an annular bubble. The main energetic benefit of this configuration comes from the fact that this shape reduces the costly liquid-vapor interface, while maximizing the favorable solid-vapor interface. It is important to remark that the time of the nucleation process scales exponentially with the free-energy barrier. If one assumes that the kinetic prefactors are similar for the two paths, a difference of 20k B T implies that the axisymmetric route is roughly 5·10 8 times slower than asymmetric one. In other words, in actual experiments one would never observe the axisymmetric case. This finding is in line with the break of symmetry predicted in other cases of confined wetting or nucleation, e.g., bubbles in rectangular or parallelepipedic grooves [14,[23][24][25].
The string in CVs is a quasi-equilibrium method, i.e., it assumes that the inertia of the process is negligible. In this limit, the nucleation paths coincide with those of the opposite phenomenon of liquid intrusion inside the nanopore, i.e., the process leading from the empty to the fully wet pore. At the thermodynamic conditions considered here, intrusion cannot happen, as it is from a low to a high free-energy state; however by increasing the pressure it is possible to induce it, while keeping the same ordering of the barrier related to the various paths [8,26]. As a consequence, the much higher free-energy barriers associated with the axisymmetric filling of the nanopore shown in fig. 3, provides a good justification of the asymmetric intrusion mechanism proposed in [8]. In sect. 5 we analyze with the aid of confined CNT predictions the energetic reasons of this break of symmetry (see also [14]).
The dynamics of asymmetric intrusion remains somewhat counterintuitive to picture when thinking in macroscopic terms of the most convenient way for the two facing liquid menisci to merge at the end of the intrusion process. At the nanoscale, however, the break of symmetry can be due to fluctuations in the interface shape. On larger scales, instead, inertial effects can become more significant and contribute to favor the (energetically unfavorable) sym-metric path; these subtle effects, however, need to be addressed with simulations techniques that go beyond the quasi-equilibrium approximation, e.g., with forward flux sampling [27].
The barriers associated with the vapor nucleation process are, in the most favorable case, represented by the nucleation of an asymmetric bubble, about 5-fold larger than the barriers reported in ref. [8] for the dewetting of a nanopore with diameter of 2.6 nm filled with water. In absolute terms, the LJ barriers are 35 times the thermal energy available to the system, k B T , which means that in the present system it is unlikely to observe dewetting at twophase coexistence; in particular, it is typically impossible within the timescales accessible to brute-force MD simulations. Due to the large nucleation and intrusion barriers, a significant intrusion/extrusion hysteresis is expected for a HLS with the present physical (diameter ca. 4 nm) and chemical (θ Y = 113 • ) characteristics: the intrusion and extrusion pressures are typically very different. In this case, application of the HLS for energy storage is not efficient, but it can be used as a damper [2] owing to the large dissipation involved in the intrusion/extrusion cycles. This difference between the present results and those of ref. [8] is mostly due to the slightly lower lyo/hydrophobicity of the solid walls (113 • vs. 119 • ) and to the size of the pore, which is about twice as large in the LJ case when compared to the characteristic correlation length of the fluids. The larger relative dimensions also account for the reduced deviations from continuum models which we discuss in the next section. Figure 4, for instance, shows how some nanoscale deviations remain in the magnitude of the free-energy barriers, but they appear less dramatic than in more confined cases.

Comparison with continuum theories
A very simple and informative approach to treat nucleation problems is represented by the classical nucleation theory (CNT), which assumes a macroscopic, sharpinterface description of the liquid-vapor-solid system. This CNT approach has been extended to investigate nucleation in complex confining geometries [4,28], with the final aim of determining the path and the free-energy profile for nucleation. In the simplest approach, a single CV is used -typically, the bubble volume. The approximate nucleation path is generated by the juxtaposition of bubble shapes minimizing the free energy for each value of the CV. A detailed analysis of confined CNT in a cylindrical geometry and of the limitations of such an approach is presented in refs. [4,14]. In fig. 4, we show the free energies related to the asymmetric "chip" solution, and to the two independent menisci. Even in this simplified treatment the energetic advantage of starting nucleation from an asymmetric bubble is apparent [4]. At least close to the initial bubble, the agreement with the string results is surprising.
For larger volumes, however, the confined CNT approach suffers from the limitations of using a single CV, e.g., it yields a discontinuous nucleation path [12,29]. In particular, close to the transition state, the transition from Fig. 4. Comparison between confined CNT predictions and atomistic data (red symbols) resulting from the string calculation initialized from the triggered nucleation event. The solid and dashed lines indicate the free-energy costs associated with the formation of a chip-like and of a cylindrical bubble, respectively, i.e., the bubble geometries considered in refs. [4,14]. The empty circles represent the values assumed by the confined CNT functional (for more details, see [8]), computed using the geometrical quantities extracted from the atomistic string calculation.
an asymmetric bubble to two independent menisci cannot be captured by this approach, hindering the evaluation of the free-energy barrier. In order to overcome these difficulties, it is possible to evaluate the free energies along the nucleation path computed via atomistic string calculations but using the sharp-interface free-energy functional. This approach allows one to identify which deviations are genuinely due to nanoscale effects, such as line contributions to the free energy [8]. The geometric parameters characterizing the interfaces are extracted from the average density fields of fig. 2, while the thermodynamic parameters (the pressure, surface tension, and contact angle) are computed via independent MD simulations. Figure 4 shows that CNT computed along the correct nucleation path is able to produce accurate predictions up to few nanometers, in fair agreement with the atomistic estimates: only a minor overestimation of the free-energy barrier is detectable. This result suggests that, for simple liquids at these scales, line and other nanoscale effects are negligibly small. In contrast, the results of ref. [8], obtained with the same procedure for water in more extreme confinement, showed that CNT models dramatically overestimate nucleation free-energy barriers and only line corrections and modified contact angles were able reconcile atomistic and CNT predictions.
The discrepancies between the water and LJ cases may have two distinct origins: the different levels of confinement and the different behavior of a simple liquid and of water. Quantitatively discriminating between these two effects would require a string simulation of the two liquids in equivalent nanopores, which is outside the scope of the present work. However, it is clear from fig. S6 of ref. [8] that the effect of confinement on the liquid structure within the nanopore is much more accentuated in the water case than in the LJ one: in the former case, layering effects are observed at all radial distances from the pore walls, while, for the latter, bulk properties are recovered at the center of the pore. It is difficult to identify other features stemming from the different nature of the two liquids other than the obvious difference in surface tensions. It is here useful to recall the case of water confined within hydrophobic cavities of characteristic size of 10 nm (but different geometries than the pore), for which the nucleation free-energy profiles were reported to be in fair agreement with macroscopic predictions [13,30]. These pieces of information together seem to suggest that, both for water and simple liquids, the size of confinement is the major contribution in determining deviations from macroscopic models.

Conclusions
In this contribution we presented extensive molecular dynamics simulations demonstrating the use of the string method in collective variables as a tool for the study of vapor nucleation in confinement. The simulation campaign focused on the formation of a LJ vapor bubble in a lyophobic mesopore, a subject of paramaount interest due to a number of technologically appealing applications relying on the peculiar physics of vapor nucleation at the nanoscale.
The inherent features of the string method in collective variables make it the ideal candidate for the study of the physical problem of our interest and allowed us to demonstrate the existence of several nucleation paths, without incurring in the artifacts associated to the use of a single collective variables and/or simplified capillary-based nucleation theories.
The primary merit of the calculations here presented is to clearly demonstrate the existence of several independent nucleation paths, associated with different freeenergy barriers ΔΩ † for nucleation and thus different nucleation rates, which scale as exp(−ΔΩ † /k B T ).
The free-energy data are critically compared to the studies available in the literature and to insights provided by simpler nucleation theories. In particular, comparing the present calculations with the TIP4P water simulations of ref. [8] (obtained using a similar simulation protocol on slightly smaller nanopores) it was possible to quantify how crucial the effect of pore size is on the nucleation dynamics: the system considered in the present work does not show significant deviations from the nucleation barriers predicted by continuum theories, resulting in larger freeenergy barriers and in a vanishing probability to observe spontaneous nucleation. The importance of these findings is readily explained by the fact that the use of HLS as energy storage devices or as dampers critically depends on their extrusion properties: the accelerated nucleation demonstrated in ref. [8] for extreme hydrophobic confinement is responsible for a highly efficient energy storage, while the present system is better suited for damping applications.
Open Access funding provided by Max Planck Society. The research leading to these results received funding from the European Research Council (ERC) under the European Union Seventh Framework Program (FP7/2007-2013)/ERC Grant Agreement 339446. We acknowledge PRACE for awarding us access to Fermi and Marconi at CINECA, Italy. We acknowledge the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support.

Author contribution statement
AG and CMC conceived the present simulations. AT performed the simulations. AT and AG analyzed the data. AT and AG wrote the first draft of the article. All authors discussed the results and contributed to the final manuscript.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.