Free energy calculations along entropic pathways. III. Nucleation of capillary bridges and bubbles

The stability of liquid bridges between solid surfaces, as well as the mechanism by which they form, is a central phenomenon in interface science and adhesion.^1–7 The formation of these bridges, as well as the nucleation of bubbles, actually occurs in a wide range of systems that span many length scales, from colloidal systems to nanosized capillaries.^8–10 This has been shown, e.g., by atomic force microscopy and surface force apparatus experiments, which have revealed how liquid contacts can form and snap off, paving the way for applications in nanotribology and nanolithography.^11–20 In recent years, the use of molecular modeling and simulation to provide a molecular level understanding of this phenomenon has drawn considerable interest.^21–28 Density functional theory calculations,^29–32 as well as molecular dynamics^33–37 and Monte Carlo simulations,^38–50 have started to shed light on the formation of solid-like and liquid-like junctions in nanopores and on the dynamics of cavitation between solid surfaces. In particular, Neimark and co-workers^51–57 have examined the capillary condensation process in cylindrical nanopores and have shown that the scenario proposed by Everett and Haynes⁵⁸ for a macroscopic capillary could also extend to nanosized capillaries. This suggested a very complex and intriguing pathway for the capillary condensation process, involving a series of structural changes in the adsorbed fluid and the formation of a liquid bridge across the pore section as an intermediate towards the high density adsorbed phase.

Condensation and evaporation processes are challenging processes to model and simulate. This is due to the large free energy barrier of nucleation that the system has to overcome to complete the phase transition.^59–66 Simulating the pathways underlying condensation and evaporation in nanoscale capillaries requires suitable simulation techniques that allow for the sampling of these rare events by driving the formation of the new phase. In this work, we focus on elucidating both the capillary condensation and evaporation pathways and on shedding light on the relation between the structure of the nanoconfined fluid and its entropy. For this purpose, we extend the recently developed $μVT−S$ simulation method^67,68 to study these two phenomena. The $μVT−S$ approach is implemented within the grand-canonical ensemble, in which the chemical potential μ, the temperature T, and the volume V of the system are held fixed. As shown in the case of the nucleation of liquid droplets for single component systems⁶⁷ and binary mixtures,⁶⁸ this method provides a direct way to calculate the entropy of the system throughout the nucleation process and, therefore, to use it as a reaction coordinate for such activated processes. In this work, we carry out $μVT−S$ simulations with the aim of (i) simulating the adsorption and desorption of the fluid in the grand-canonical statistical ensemble, that is especially suited to study adsorption phenomena and (ii) using the entropy of the confined fluid as the reaction coordinate to elucidate the two pathways underlying nanoscale capillary condensation and evaporation. To achieve this, we impose a target value for the entropy of the system through an umbrella sampling bias potential.⁶⁹ Then, by varying the value of the entropy of the system, we are able to sample the entire pathway connecting the two nanoconfined phases, starting from the metastable confined vapor and ending with the nanoconfined liquid in the case of capillary condensation. Doing so, we are able to identify the different intermediate structures, specifically bridges and bubbles, that form in the nanopores during the capillary condensation and evaporation processes.

The paper is organized as follows. We start by discussing how we extend the $μVT−S$ approach to study the condensation and evaporation processes of argon in MCM-41 silica mesoporous molecular sieves. We also present the potential models used for the adsorbate and its interaction with the cylindrical nanopore and provide the technical details for the simulations. We then examine the capillary condensation process for argon in the nanopore and determine the free energy barrier associated with the process. We carry out a detailed structural analysis to identify the formation of the liquid bridge along this pathway. We also discuss the capillary evaporation process that occurs during the desorption of the confined fluid and apply the same analyses as for the reverse pathway. We finally draw our main conclusions in Sec. V.

We study the capillary condensation and evaporation of argon in a cylindrical pore of 10 atomic diameters. The geometry is defined as follows. The cylindrical pore is aligned with the z-axis, and we define the length, or lateral dimension, as L_z. Simulations are carried out at the boiling point for argon (T = 87.3 K) within a nanopore found in MCM-41 silica mesoporous molecular sieves.⁵⁴ We model argon with a Lennard-Jones potential and take the following parameters for the exclusion diameter $σ=3.4$ Å and for the potential well-depth $ϵ/kB=119.8$ K. The interactions between the adsorbed argon atoms and the nanopore are modeled with the following functional form:^52–57 

in which r is the radial coordinate of the Ar atom adsorbed in the pore, R is the pore radius (here $5σ$⁠), $ρs$ is the surface density of adsorption centers, and $F(α,β;γ;δ)$ is the hypergeometric series. The parameters for the solid-fluid interactions were taken as $ρsϵsf=2253$ K/nm² and $σsf=3.17$ Å. This functional form accurately models the interaction between the adsorbate and the structureless cylindrical layer of adsorption centers on the pore wall as shown by Ravikovitch et al.⁷⁰ Following Vishnyakov and Neimark,⁵⁴ we carry out simulations of capillary condensation and evaporation in nanopores with a long lateral dimension $Lz=30σ$ to allow for the sampling of symmetry breaking configurations containing bubbles and liquid bridges. Periodic boundary conditions are applied along this lateral direction z. We also calculate explicitly the interactions between argon atoms up to a distance of $5σ$ and neglect the fluid-fluid interactions beyond that cutoff distance. Finally, in the rest of this work, we use the conventional set of reduced units,⁷¹ with respect to the Lennard-Jones parameters of the fluid.

To span the entire pathway underlying the capillary condensation and evaporation processes, we carry out a series of $μVT−S$ simulations and gradually vary the entropy of the system to drive the phase transition within the confined fluid. We briefly outline here the principles of the $μVT−S$ method (more details may be found in previous work^67,68). $μVT−S$ simulations are carried out in the grand-canonical ensemble $(μ,V,T)$⁠, i.e., at constant chemical potential μ, temperature T, and volume V. This statistical ensemble is especially well suited to study the adsorption phenomena, since the number of atoms N adsorbed in the nanopore is allowed to vary in the $(μ,V,T)$ ensemble. Thus, this ensemble mimics what is observed experimentally during the adsorption/desorption processes. This ensemble has also the advantage of providing a direct way of estimating the entropy of the confined fluid. Taking the total Helmholtz free energy A = U − TS of the confined fluid to be equal to $Nμ$⁠, we obtain the following equation for the entropy of the adsorbed fluid S_loading:

In this equation, S_loading denotes the total entropy of the system and, as such, increases with the total number of atoms N adsorbed in the nanopore. We add that the pV term is omitted in Eq. (2) since it is negligible when compared to the other terms (pV actually accounts for less than 0.5% of the lowest values sampled for TS_loading). This approximation is similar to the one made by Waghe, Rasaiah, and Hummer in their calculation of the entropy of water adsorbed in carbon nanotubes.⁷² To sample the transition pathway, we use S_loading as the reaction coordinate and drive the nucleation event through an umbrella sampling approach.⁶⁹ For this purpose, we define the following bias energy:

in which S₀ is the target value for the total entropy. This bias energy is added to the total potential energy of the system, which appears in the Metropolis criteria used in the Monte Carlo simulations. In practice, for capillary condensation, we perform a series of $μVT−S$ simulations with increasing values for the target entropy S₀. This promotes the uptake of additional Ar atoms since the total entropy of the system is a function of N. Similarly, the pathway for the desorption process is sampled by carrying out successive $μVT−S$ simulations, with decreasing values for the target entropies S₀.

Previous work⁷³ has shown that the conditions for liquid-vapor equilibrium in the pore are obtained for $μ=−10.53ϵ$⁠. To observe the capillary condensation/evaporation processes, the chemical potential must be close enough to the chemical potential at equilibrium. More specifically, to simulate the capillary condensation process, we carry out $μVT−S$ simulations for a chemical potential of $μ=−10.48ϵ$⁠, while for the evaporation process, we use a chemical potential of $μ=−10.54ϵ$⁠. These two values for μ show that our findings are consistent with the results from Peterson and Gubbins.⁷³ As expected, these two values of μ bracket the estimate of Peterson and Gubbins for the chemical potential at the liquid-vapor equilibrium, with μ for the capillary condensation process being slightly above the chemical potential at coexistence and μ for the capillary evaporation process being slightly below. A total of 35 umbrella sampling windows are performed to connect the vapor and liquid phases for the confined fluid and to sample the entire condensation and evaporation pathways. For each window, we first carry out an equilibration run of $1×108$ MC steps to allow the system to relax towards the target value for the entropy. We then perform a production run of $2×108$ MC steps, during which averages are collected and structural analyses are run.

We start by discussing the results obtained for the capillary condensation process. To sample the entropic pathway for this process, we carry out a series of $μVT−S$ simulations at $μ=−10.48ϵ$⁠. We gradually increase the value for the target entropy, as we go from one umbrella sampling window to another, and, as a result, increase S_loading for the confined fluid. We show in Fig. 1 how the system responds to this increase in the entropy of the adsorbed fluid. The top panel of Fig. 1 consists of a plot of the number of argon atoms N adsorbed in the nanopore as a function of the reduced entropy of the adsorbed fluid (noted as $Sloading*$⁠). We observe the following trends. First, we find that N increases with $Sloading*$⁠, starting from about 900 Ar atoms adsorbed for $Sloading*=3300$⁠, and reaches a value of about 1600 Ar atoms for $Sloading*=6800$⁠. This plot shows that the successive $μVT−S$ simulations, with increasing values for the total entropy of the system, allow to simulate the phase transition process from a metastable phase of low density to the stable phase of high density. The bottom panel of Fig. 1 shows the corresponding variation of the interaction energy within the confined fluid during the $μVT−S$ simulations. Overall, we observe a decrease in the interaction energy as $Sloading*$ increases. This behavior is consistent with the increase in the amount of Ar adsorbed in the cylindrical nanopore observed when $Sloading*$ increases.

We plot in Fig. 2 the free energy profile obtained for the capillary condensation process. The starting point for the $μVT−S$ simulations is a metastable phase of low density adsorbed in the nanopore, with a reduced entropy of $Sloading*=3300$⁠. As shown on the left of Fig. 2, the free energy profile exhibits a local minimum for this metastable phase. For convenience, we have chosen to assign the origin for the free energy to this local minimum. Several regimes can be identified in this plot of the free energy as a function of the entropy of the adsorbed phase. First, the free energy profile exhibits a steep increase, over entropies ranging from $Sloading*=3300$ to $Sloading*=4800$⁠. This increase in the free energy corresponds to the steady increase in N, and accordingly to the steady decrease of $Uinteraction*$⁠, over the same range of entropies. Then, for $Sloading*=4800$ to $Sloading*=6000$⁠, the free energy profile exhibits a second regime which can be characterized as an almost flat top that very slowly decreases with $Sloading*$⁠. This second regime is associated with changes in both the slopes for the variations of N and $Uinteraction*$ in Fig. 1. Then, for $Sloading*=6000$ to $Sloading*=6850$⁠, we observe a third regime for the free energy profile and a steady decrease in free energy. During this stage, the system completes the phase transition towards the stable phase of high density. As expected, we obtain a free energy for the stable phase that is below, by 11 k_BT, the free energy of the metastable phase of low density. This establishes that the end point for the $μVT−S$ simulations is indeed the stable phase. As for the previous two regimes in the free energy profile, we find that the decrease for the free energy, observed during the third regime, is connected to sharp changes in the variations of both N and $Uinteraction*$ for the adsorbed phase. This indicates that significant structural changes take place within the fluid confined in the nanopore.

To elucidate the mechanism underlying capillary condensation, we carry out a series of structural analyses. We begin by determining the density profiles for the adsorbed fluid across the nanopore. The results obtained during several of the $μVT−S$ simulations are shown in Fig. 3. At the start of the capillary condensation process, i.e., for the metastable phase of low density, the adsorbed fluid mostly consists of two layers close to the wall, as shown by the two peaks on the density profile for radii of about 4 σ (fluid layer closest to the wall) and 3 σ (second layer). As $Sloading*$ increases, Ar atoms start to populate the center of the nanopore, and a third peak for a radius around 2 σ develops (see the density profile in red for $Sloading*=4290$⁠). We detect the formation of a fourth and then of a fifth peak when the system reaches the top of the free energy barrier (see density profile in green). Then, as $Sloading*$ further increases, the height of the peaks associated with the innermost fluid layers continue to increase, until the nanopore has completely filled and the system reaches the stable phase of high density.

We now more closely look at the radial density profiles of Fig. 3 to determine why, for intermediate values of $Sloading*$⁠, we find lower peaks for the innermost fluid layers. This prompts the following question. Are these lower peaks associated with a uniform distribution of Ar atoms along the side of the tube (along the z direction)? Alternatively, do these lower peaks signal the formation of a liquid contact within the nanotube? To answer this question, we need to determine how the distribution of atoms evolves along the side of the cylindrical nanopore, i.e., along the z axis. For this purpose, we refine our structural analysis by computing the density profiles against both the radius and the z coordinate and show how the radial density varies along the nanopore. We show in Fig. 4(a) the results obtained for the metastable phase of low density. The plot for the density profile shows that the atoms adsorbed are very predominantly located close to the wall. More specifically, the fluid is organized into two layers close to the wall as shown by the two bright lines (i.e., regions of high density) on the density profile of Fig. 4(a). This is consistent with the profiles of Fig. 3 that showed the formation of two layers close to the wall. The additional information brought by Fig. 4(a) is the fact that Ar atoms in these two layers are distributed uniformly along z. We also show in Fig. 4(a) a snapshot of a configuration of the system for the metastable phase of low density and highlight in cyan the region where the void dominates (we define void as any region of space for which the reduced density is below a threshold value). We carry out the same analysis along the rest of the capillary condensation pathway. At the top of the free energy barrier, we obtain the density profile and snapshot shown in Fig. 4(b). Both show that, at this stage, the Ar atoms located in the inner part of the nanopore are not uniformly distributed along z and that a liquid-like region has developed for a value of z around 20 σ on the density profile. This can best be seen on the density profile with the onset, for z = 20, of 3 bright spots for r around 0, 1, and 2, that indicate the formation of 3 partial liquid layers. The snapshot of a configuration of the system at the top of the free energy barrier (see Fig. 4(b)) confirms that a liquid bridge has developed across the nanopore for this value of z = 20. As the capillary condensation further proceeds, the liquid bridge becomes wider and wider along the nanopore, resulting in a recess of the void region. Finally, as the system reaches the stable high density phase, we see that the structure of the adsorbed fluid exhibits a uniform structure with multiple layers (see Fig. 4(c)). The mechanism observed for capillary condensation is thus consistent with the scenario proposed by Everett and Haynes⁵⁸ for a macroscopic capillary and the results from gauge cell Monte Carlo simulations from Neimark and Vishnyakov,⁵⁴ where both pointed to a liquid bridge-mediated capillary condensation process. We add that the free energy barriers are of the same order as those previously obtained for the capillary condensation process by Vishnyakov and Neimark⁵⁴ using the gauge cell method.

We now turn to the results obtained for the capillary evaporation process (⁠$μ=−10.54ϵ$⁠). We show in Fig. 5 the free energy profile obtained as the system moves from a metastable phase of high density to the stable phase of low density. In line with our results for the capillary condensation process, we set the origin for the free energy axis to the free energy of the metastable phase. The main difference with the case of the condensation process, studied in the first part of the paper, is that we now carry out a series of $μVT−S$ simulations with decreasing values for the total entropy of the system. Therefore, $Sloading*$ now varies in the opposite sense to the one followed for the condensation process. Looking at the free energy profile as a function of $Sloading*$⁠, we find three successive regimes for the variations of the free energy of the system. First, for $Sloading*$ varying from 7000 (metastable phase of high density) to 6500 (top of the free energy barrier), we observe a steep increase in free energy. Then, for $Sloading*$ between 6500 and 4800, the free energy profile exhibits a much more gradual behavior, as the free energy slowly decreases with the entropy of the system, in an almost linear fashion. The third regime is observed for lower values of $Sloading*$ (between 4800 and 2900). For this range of entropies, the free energy decreases more rapidly with entropy, until the system reaches the stable phase of low density. At this point, the free energy of the system is 55 k_BT lower than the metastable state from which the evaporation pathway started. This confirms that the system has reached a thermodynamically stable phase at the end of the $μVT−S$ simulations.

We now analyze the structure of the adsorbed fluid along the evaporation pathway. Fig. 6 shows the density profiles across the cylindrical nanopore for different values of the entropy. The starting point for the $μVT−S$ simulations displays the expected multi-peaked density profile. This density profile is consistent with the organization of the metastable phase of high density in multiple fluid layers within the nanopore (⁠$Sloading*=6980$⁠). Close to the top of the free energy barrier (⁠$Sloading*=4761$⁠), we observe a small, but noticeable, increase in the number of Ar atoms adsorbed in the nanopore. This indicates that the capillary evaporation process starts with the adsorption of a few extra Ar atoms, amounting to about 1%–2% of the overall number of Ar atoms within the system. These extra Ar atoms destabilize the organization within the nanopore and trigger the desorption process, as shown by the decrease in the amplitude of the peaks in the density profile that follows (see, e.g., the density profile obtained for $Sloading*=5750$⁠). This leads to a decrease in the free energy as $Sloading*$ further decreases. Then, at the end of the second regime found in the free energy plot, we observe that the density profile exhibits two main peaks, the third inner most layer showing a reduced fluid density (⁠$Sloading*=4761$⁠). Finally, we recover the structure of the metastable phase of low density (see the density profile for at $Sloading*=2887$⁠), for which the Ar atoms are very predominantly adsorbed close to the wall of the cylindrical nanopore.

We finally comment on the local structure within the fluid along the side of the nanopore. We observe that once the organization of the high density phase has been destabilized, capillary evaporation proceeds through the formation of a bubble within the nanopore. This can best be seen in Fig. 7 for $Sloading*=5750$⁠. The density profile of Fig. 7 shows that a dark spot develops for values of z around 15 and that the innermost layers of the adsorbed fluid start getting depleted. Furthermore, this plot establishes that a non-uniform distribution of Ar atoms has formed along the side of the nanopore. In this case, the dark spots on the density profile on the left of Fig. 7, as well as the reduction in the height of the peaks for the innermost layers in Fig. 6, are due to the formation of a bubble, as shown through the low density region indicated on the snapshot on the right of Fig. 7.

In this work, we carry out $μVT−S$ simulations to elucidate the capillary condensation and evaporation processes. Using the total entropy for the adsorbate, we are able to drive these processes along the entropic pathways underlying these processes. This allows us to shed light on the structural changes that occur within the confined fluid and give rise to the phase transition. Considering the example of argon adsorbed in a smooth cylindrical nanopore typical of the MCM-41 silica adsorbent, we start by focusing on capillary condensation. Our simulations allow us to identify a complex free energy profile, corresponding to the successive stages in the condensation process. Our results show that capillary condensation from a metastable phase of low density starts with the nucleation of a liquid bridge within the nanopore. This liquid bridge then expands, becomes wider and wider, as the total entropy of the system increases, and finally takes over the system to yield the stable phase of high density. Applying this approach to the phenomenon of capillary evaporation similarly uncovers a multi-stage process. In this case, the layered structure of the metastable phase of high density is first destabilized by the adsorption of a few extra atoms. This is then followed by the nucleation of a bubble within the nanopore. Overall, the mechanism identified for capillary condensation is consistent with the scenario of Everett and Haynes⁵⁸ for a macroscopic system and with gauge cell Monte Carlo simulations by Neimark and Vishnyakov.⁵⁴ The results obtained here further establish the key role played by the nucleation of liquid bridges in the capillary condensation process, and of the nucleation of bubbles during capillary evaporation, and shed light on the interplay between the entropy and the structure of the adsorbed fluid throughout these processes.

Partial funding for this research was provided by NSF through CAREER Award No. DMR-1052808. Acknowledgement is made to the Donors of the American Chemical Society Petroleum Research Fund for partial support of this research through Grant No. 548002-ND10.